barrel_ecal/scripts/emcal_barrel_pions_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 #include "TCanvas.h"
0017 #include "TStyle.h"
0018 #include "TMath.h"
0019 #include "TH1.h"
0020 #include "TF1.h"
0021 #include "TH1D.h"
0022 #include "TFitResult.h"
0023
0024 using ROOT::RDataFrame;
0025 using namespace ROOT::VecOps;
0026
0027 void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal_barrel_piplus.edm4hep.root")
0028 {
0029   // Setting for graphs
0030   gROOT->SetStyle("Plain");
0031   gStyle->SetOptFit(1);
0032   gStyle->SetLineWidth(2);
0033   gStyle->SetPadTickX(1);
0034   gStyle->SetPadTickY(1);
0035   gStyle->SetPadGridX(1);
0036   gStyle->SetPadGridY(1);
0037   gStyle->SetPadLeftMargin(0.14);
0038   gStyle->SetPadRightMargin(0.14);
0039
0040   ROOT::EnableImplicitMT();
0041   ROOT::RDataFrame d0("events", input_fname);
0042
0043   // Sampling Fraction
0044   double samp_frac = 0.0136;
0045
0046   // Thrown Energy [GeV]
0047   auto Ethr = [](std::vector<edm4hep::MCParticleData> const& input) {
0048     auto p = input[2];
0049     auto energy = TMath::Sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y + p.momentum.z * p.momentum.z + p.mass * p.mass);
0050     return energy;
0051   };
0052
0053   // Number of hits
0054   auto nhits = [] (const std::vector<edm4hep::SimCalorimeterHitData>& evt) {return (int) evt.size(); };
0055
0056   // Energy deposition [GeV]
0057   auto Esim = [](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0058     auto total_edep = 0.0;
0059     for (const auto& i: evt)
0060       total_edep += i.energy;
0061     return total_edep;
0062   };
0063
0064   // Sampling fraction = Esampling / Ethrown
0065   auto fsam = [](const double sampled, const double thrown) {
0066     return sampled / thrown;
0067   };
0068
0069   // Energy Resolution = Esampling/Sampling_fraction - Ethrown
0070   auto eResol = [samp_frac](double sampled, double thrown) {
0071     return sampled / samp_frac - thrown;
0072   };
0073
0074   // Relative Energy Resolution = (Esampling/Sampling fraction - Ethrown)/Ethrown
0075   auto eResol_rel = [samp_frac](double sampled, double thrown) {
0076       return (sampled / samp_frac - thrown) / thrown;
0077   };
0078
0079   // Returns the pdgID of the particle
0080   auto getpid = [](std::vector<edm4hep::MCParticleData> const& input) {
0081     return input[2].PDG;
0082   };
0083
0084   // Returns number of particle daughters
0085   auto getdau = [](std::vector<edm4hep::MCParticleData> const& input) {
0086     return input[2].daughters_begin;
0087   };
0088
0089   // Define variables
0090   auto d1 = ROOT::RDF::RNode(
0091     d0.Define("Ethr", Ethr, {"MCParticles"})
0092       .Define("pid",    getpid,     {"MCParticles"})
0093   );
0094
0095   auto Ethr_max = 7.5;
0096   auto fsam_est = 1.0;
0097   if (d1.HasColumn("EcalBarrelScFiHits")) {
0098     d1 = d1.Define("nhits", nhits, {"EcalBarrelImagingHits"})
0099            .Define("EsimImg", Esim, {"EcalBarrelImagingHits"})
0100            .Define("EsimScFi", Esim, {"EcalBarrelScFiHits"})
0101            .Define("Esim", "EsimImg+EsimScFi")
0102            .Define("fsamImg", fsam, {"EsimImg", "Ethr"})
0103            .Define("fsamScFi", fsam, {"EsimScFi", "Ethr"})
0104            .Define("fsam", fsam, {"Esim", "Ethr"});
0105     fsam_est = 0.1;
0106   } else {
0107     d1 = d1.Define("nhits", nhits, {"EcalBarrelSciGlassHits"})
0108            .Define("Esim", Esim, {"EcalBarrelSciGlassHits"})
0109            .Define("fsam", fsam, {"Esim", "Ethr"});
0110     fsam_est = 1.0;
0111   }
0112
0113   // Define Histograms
0114   auto hEthr  = d1.Histo1D({"hEthr",  "Thrown Energy; Thrown Energy [GeV]; Events",        100,  0.0, Ethr_max}, "Ethr");
0115   auto hNhits = d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events", 100,  0.0,   2000.0}, "nhits");
0116   auto hEsim  = d1.Histo1D({"hEsim",  "Energy Deposit; Energy Deposit [GeV]; Events",      100,  0.0,      1.0}, "Esim");
0117   auto hfsam  = d1.Histo1D({"hfsam",  "Sampling Fraction; Sampling Fraction; Events",      100,  0.0, fsam_est}, "fsam");
0118   auto hpid   = d1.Histo1D({"hpid",   "PID; PID; Count",                                   100,  -220,     220}, "pid");
0119
0120   // Event Counts
0121   auto nevents_thrown      = d1.Count();
0122   std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
0123
0124   // Draw Histograms
0125   TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
0126   c1->SetLogy(1);
0127   hEthr->GetYaxis()->SetTitleOffset(1.4);
0128   hEthr->SetLineWidth(2);
0129   hEthr->SetLineColor(kBlue);
0130   hEthr->DrawClone();
0131   c1->SaveAs("results/emcal_barrel_pions_Ethr.png");
0132   c1->SaveAs("results/emcal_barrel_pions_Ethr.pdf");
0133
0134   TCanvas *c2 = new TCanvas("c2", "c2", 700, 500);
0135   c2->SetLogy(1);
0136   hNhits->GetYaxis()->SetTitleOffset(1.4);
0137   hNhits->SetLineWidth(2);
0138   hNhits->SetLineColor(kBlue);
0139   hNhits->DrawClone();
0140   c2->SaveAs("results/emcal_barrel_pions_nhits.png");
0141   c2->SaveAs("results/emcal_barrel_pions_nhits.pdf");
0142
0143   TCanvas *c3 = new TCanvas("c3", "c3", 700, 500);
0144   c3->SetLogy(1);
0145   hEsim->GetYaxis()->SetTitleOffset(1.4);
0146   hEsim->SetLineWidth(2);
0147   hEsim->SetLineColor(kBlue);
0148   hEsim->DrawClone();
0149   c3->SaveAs("results/emcal_barrel_pions_Esim.png");
0150   c3->SaveAs("results/emcal_barrel_pions_Esim.pdf");
0151
0152   TCanvas *c4 = new TCanvas("c4", "c4", 700, 500);
0153   c4->SetLogy(1);
0154   hfsam->GetYaxis()->SetTitleOffset(1.4);
0155   hfsam->SetLineWidth(2);
0156   hfsam->SetLineColor(kBlue);
0157   hfsam->DrawClone();
0158   c4->SaveAs("results/emcal_barrel_pions_fsam.png");
0159   c4->SaveAs("results/emcal_barrel_pions_fsam.pdf");
0160
0161   TCanvas *c5 = new TCanvas("c5", "c5", 700, 500);
0162   c5->SetLogy(1);
0163   hpid->GetYaxis()->SetTitleOffset(1.4);
0164   hpid->SetLineWidth(2);
0165   hpid->SetLineColor(kBlue);
0166   hpid->DrawClone();
0167   c5->SaveAs("results/emcal_barrel_pions_pid.png");
0168   c5->SaveAs("results/emcal_barrel_pions_pid.pdf");
0169
0170 }