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