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