barrel_ecal/scripts/emcal_barrel_pi0_analysis.cxx

0001 ////////////////////////////////////////
0002 // Read reconstruction ROOT output file
0003 // Plot variables
0004 ////////////////////////////////////////
0005
0006 #include "ROOT/RDataFrame.hxx"
0007 #include <iostream>
0008
0009 #include "edm4hep/MCParticleCollection.h"
0010 #include "edm4hep/SimCalorimeterHitCollection.h"
0011
0012 #include "common_bench/benchmark.h"
0013 #include "common_bench/mt.h"
0014 #include "common_bench/util.h"
0015
0016 R__LOAD_LIBRARY(libfmt.so)
0017 #include "fmt/core.h"
0018
0019 #include "TCanvas.h"
0020 #include "TStyle.h"
0021 #include "TMath.h"
0022 #include "TH1.h"
0023 #include "TF1.h"
0024 #include "TH1D.h"
0025 #include "TFitResult.h"
0026 #include <nlohmann/json.hpp>
0027
0028 using ROOT::RDataFrame;
0029 using namespace ROOT::VecOps;
0030 using json = nlohmann::json;
0031
0032 void emcal_barrel_pi0_analysis(
0033                                 const char* input_fname = "sim_output/sim_emcal_barrel_pi0.edm4hep.root"
0034                                 //const char* input_fname = "../sim_output/sim_emcal_barrel_uniform_pi0.edm4hep.root"
0035                                 )
0036 {
0037   // Setting for graphs
0038   gROOT->SetStyle("Plain");
0039   gStyle->SetOptFit(1);
0040   gStyle->SetLineWidth(2);
0041   gStyle->SetPadTickX(1);
0042   gStyle->SetPadTickY(1);
0043   gStyle->SetPadGridX(1);
0044   gStyle->SetPadGridY(1);
0045   gStyle->SetPadLeftMargin(0.14);
0046   gStyle->SetPadRightMargin(0.14);
0047
0048   ROOT::EnableImplicitMT();
0049   ROOT::RDataFrame d0("events", input_fname);
0050
0051   // Sampling Fraction grabbed from json file
0052   json j;
0053   std::ifstream prev_steps_ifstream("results/emcal_barrel_electron_calibration.json");
0054   prev_steps_ifstream >> j;
0055   double samp_frac = j["electron"]["sampling_fraction"];
0056
0057   // Thrown Energy [GeV]
0058   auto Ethr = [](std::vector<edm4hep::MCParticleData> const& input) {
0059     return TMath::Sqrt(input[2].momentum.x*input[2].momentum.x + input[2].momentum.y*input[2].momentum.y + input[2].momentum.z*input[2].momentum.z + input[2].mass*input[2].mass);
0060   };
0061
0062   // Number of hits
0063   auto nhits = [] (const std::vector<edm4hep::SimCalorimeterHitData>& evt) {return (int) evt.size(); };
0064
0065   // Energy deposition [GeV]
0066   auto Esim = [](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0067     auto total_edep = 0.0;
0068     for (const auto& i: evt){
0069       total_edep += i.energy;
0070     }
0071     return total_edep;
0072   };
0073
0074   // Sampling fraction = Esampling / Ethrown
0075   auto fsam = [](const double& sampled, const double& thrown) {
0076     return sampled / thrown;
0077   };
0078
0079   // Energy Resolution = Esampling/Sampling_fraction - Ethrown
0080   auto eResol = [&](const double& sampled, const double& thrown){
0081     return sampled / samp_frac - thrown;
0082   };
0083
0084   // Relative Energy Resolution = (Esampling/Sampling fraction - Ethrown)/Ethrown
0085   auto eResol_rel = [&](const double& sampled, const double& thrown){
0086     return (sampled / samp_frac - thrown) / thrown;
0087   };
0088
0089   // Returns the pdgID of the particle
0090   auto getpid = [](std::vector<edm4hep::MCParticleData> const& input) {
0091     return input[2].PDG;
0092   };
0093
0094   // Returns number of particle daughters
0095   auto getdau = [](std::vector<edm4hep::MCParticleData> const& input) {
0096     return input[2].daughters_begin;
0097   };
0098
0099   // Define variables
0100   auto d1 = ROOT::RDF::RNode(
0101     d0.Define("Ethr", Ethr, {"MCParticles"})
0102       .Define("pid", getpid, {"MCParticles"})
0103       .Define("dau", getdau, {"MCParticles"})
0104   );
0105
0106   auto Ethr_max = 7.5;
0107   auto fsam_est = 1.0;
0108   if (d1.HasColumn("EcalBarrelScFiHits")) {
0109     d1 = d1.Define("nhits", nhits, {"EcalBarrelImagingHits"})
0110            .Define("EsimImg", Esim, {"EcalBarrelImagingHits"})
0111            .Define("EsimScFi", Esim, {"EcalBarrelScFiHits"})
0112            .Define("Esim", "EsimImg+EsimScFi")
0113            .Define("fsamImg", fsam, {"EsimImg", "Ethr"})
0114            .Define("fsamScFi", fsam, {"EsimScFi", "Ethr"})
0115            .Define("fsam", fsam, {"Esim", "Ethr"});
0116     fsam_est = 1.2*samp_frac;
0117   } else {
0118     d1 = d1.Define("nhits", nhits, {"EcalBarrelSciGlassHits"})
0119            .Define("Esim", Esim, {"EcalBarrelSciGlassHits"})
0120            .Define("fsam", fsam, {"Esim", "Ethr"});
0121     fsam_est = 1.0;
0122   }
0123   d1 = d1.Define("dE",         eResol,        {"Esim","Ethr"})
0124          .Define("dE_rel",     eResol_rel,    {"Esim","Ethr"});
0125
0126   // Define Histograms
0127   std::vector <std::string> titleStr = {
0128                                         "Thrown Energy; Thrown Energy [GeV]; Events",
0129                                         "Number of hits per events; Number of hits; Events",
0130                                         "Energy Deposit; Energy Deposit [GeV]; Events",
0131                                         "dE Relative; dE Relative; Events"
0132   };
0133
0134   std::vector<std::vector<double>> range = {{0, Ethr_max}, {0, 2000}, {0, fsam_est * Ethr_max}, {-3, 3}};
0135   std::vector<std::string> col           = {"Ethr",   "nhits",   "Esim", "dE_rel"};
0136
0137   double meanE  = 5;
0138   int nCol = range.size();
0139   for (int i = 0; i < nCol; i++){
0140     int binNum = 100;
0141     auto h = d1.Histo1D({"hist", titleStr[i].c_str(), binNum, range[i][0], range[i][1]}, col[i].c_str());
0142     if (col[i] == "Ethr"){
0143       meanE = h->GetMean();
0144     }
0145     auto *c = new TCanvas("c", "c", 700, 500);
0146     c->SetLogy(1);
0147     auto h1 = h->DrawCopy();
0148     h1->GetYaxis()->SetTitleOffset(1.4);
0149     h1->SetLineWidth(2);
0150     h1->SetLineColor(kBlue);
0151     c->SaveAs((fmt::format("results/emcal_barrel_pi0_{}.png", col[i])).c_str());
0152     c->SaveAs((fmt::format("results/emcal_barrel_pi0_{}.pdf", col[i])).c_str());
0153     //std::printf("Generic %d\n", i);
0154   }
0155
0156   // Resolution Plots
0157   titleStr = {
0158               "Sampling Fraction; Sampling Fraction; Events",
0159               "dE; dE[GeV]; Events"
0160   };
0161   range = {{0,fsam_est}, {-3, 3}};
0162   col   = {"fsam",   "dE"};
0163   nCol  = range.size();
0164   std::printf("Here %d\n", 10);
0165   std::vector<std::vector<double>> fitRange = {{0.005, fsam_est}, {-3, 3}};
0166   double sigmaOverE = 0;
0167
0168   auto hr = d1.Histo1D({"histr", titleStr[0].c_str(), 150, range[0][0], range[0][1]}, col[0].c_str());
0169   auto *c = new TCanvas("c", "c", 700, 500);
0170   c->SetLogy(1);
0171   auto h2 = hr->DrawCopy();
0172   h2->GetYaxis()->SetTitleOffset(1.4);
0173   h2->SetLineWidth(2);
0174   h2->SetLineColor(kBlue);
0175   h2->Fit("gaus","","", fitRange[0][0], fitRange[0][1]);
0176   h2->GetFunction("gaus")->SetLineWidth(2);
0177   h2->GetFunction("gaus")->SetLineColor(kRed);
0178
0179   c->SaveAs((fmt::format("results/emcal_barrel_pi0_{}.png", col[0])).c_str());
0180   c->SaveAs((fmt::format("results/emcal_barrel_pi0_{}.pdf", col[0])).c_str());
0181   std::printf("Resolution %d\n", 0);
0182
0183
0184   auto hs = d1.Histo1D({"hists", titleStr[1].c_str(), 100, range[1][0], range[1][1]}, col[1].c_str());
0185   auto *c1 = new TCanvas("c1", "c1", 700, 500);
0186   c1->SetLogy(1);
0187   auto h3 = hs->DrawCopy();
0188   h3->GetYaxis()->SetTitleOffset(1.4);
0189   h3->SetLineWidth(2);
0190   h3->SetLineColor(kBlue);
0191   auto fit = h3->Fit("gaus","","", fitRange[1][0], fitRange[1][1]);
0192   if (fit == 0) {
0193     double* res = h3->GetFunction("gaus")->GetParameters();
0194     sigmaOverE  = res[2] / meanE;
0195   } else {
0196     std::printf("Fit failed\n");
0197     sigmaOverE  = h3->GetStdDev() / h3->GetMean();
0198   }
0199   c1->SaveAs((fmt::format("results/emcal_barrel_pi0_{}.png", col[1])).c_str());
0200   c1->SaveAs((fmt::format("results/emcal_barrel_pi0_{}.pdf", col[1])).c_str());
0201   std::printf("Resolution %d\n", 1);
0202
0203   // Energy Resolution Calculation
0204   std::string test_tag = "Barrel_emcal_pi0";// TODO: Change test_tag to something else
0205   std::string detEle   = "Barrel_emcal";
0206
0207   // Energy resolution in the barrel region (-1 < eta < 1)
0208   // Taken from : Initial considerations for EMCal of the EIC detector by A. Bazilevsky
0209   // sigma_E / E = 12% / E^0.5 convoluted with 2%
0210   // sigma_E / E = [ (0.12/E^0.5)^2 + 0.02^2]^0.5, with E in [GeV]
0211   double resolutionTarget = TMath::Sqrt(0.12 * 0.12 / meanE + 0.02 * 0.02);
0212
0213   common_bench::Test pi0_energy_resolution{
0214    {{"name", fmt::format("{}_energy_resolution", test_tag)},
0215    {"title", "Pi0 Energy resolution"},
0216    {"description",
0217     fmt::format("Pi0 energy resolution for {}, estimated using a Gaussian fit.", detEle)},
0218    {"quantity", "resolution (in %)"},
0219    {"target", std::to_string(resolutionTarget)}}
0220   };
0221
0222   //// Pass/Fail
0223   sigmaOverE <= resolutionTarget ? pi0_energy_resolution.pass(sigmaOverE) : pi0_energy_resolution.fail(sigmaOverE);
0224   common_bench::write_test({pi0_energy_resolution}, fmt::format("results/{}_pi0.json", detEle));
0225 }