diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_analysis.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_analysis.cxx index 59c2d7d7d337e16ff0809eb9dcb4e422a9db8aca..7436c67b3846823b98e4d58e4e2519368805a833 100644 --- a/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_analysis.cxx +++ b/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_analysis.cxx @@ -83,7 +83,7 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal return result; }; - // Layer Sampling fraction = Esampling / Ethrown + // Sampling fraction = Esampling / Ethrown auto fsam = [](const std::vector<double>& sampled, const std::vector<double>& thrown) { std::vector<double> result; auto it_sam = sampled.cbegin(); @@ -116,6 +116,11 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal return result; }; + // Relative Energy Resolution = (Esampling - Ethrown)/Ethrown + auto getpid = [](std::vector<dd4pod::Geant4ParticleData> const& input) { + return input[2].pdgID(); + }; + // Define variables auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"}) .Define("nhits", nhits, {"EcalBarrelHits"}) @@ -123,13 +128,15 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal .Define("fsam", fsam, {"Esim","Ethr"}) .Define("dE", eResol, {"Esim","Ethr"}) .Define("dE_rel", eResol_rel, {"Esim","Ethr"}) + .Define("pid", getpid, {"mcparticles"}) ; // Define Histograms auto hEthr = d1.Histo1D({"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 7.5}, "Ethr"); auto hNhits = d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events", 100, 0.0, 2000.0}, "nhits"); auto hEsim = d1.Histo1D({"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 100, 0.0, 1.0}, "Esim"); - auto hfsam = d1.Histo1D({"hfsam", "Layer Sampling Fraction; Layer Sampling Fraction; Events", 100, 0.0, 0.1}, "fsam"); + auto hfsam = d1.Histo1D({"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1}, "fsam"); + auto hpid = d1.Histo1D({"hpid", "PID; PID; Count", 100, -220, 220}, "pid"); // Event Counts auto nevents_thrown = d1.Count(); @@ -175,6 +182,15 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal c4->SaveAs("results/emcal_barrel_pions_fsam.png"); c4->SaveAs("results/emcal_barrel_pions_fsam.pdf"); + TCanvas *c5 = new TCanvas("c5", "c5", 700, 500); + c5->SetLogy(1); + hpid->GetYaxis()->SetTitleOffset(1.4); + hpid->SetLineWidth(2); + hpid->SetLineColor(kBlue); + hpid->DrawClone(); + c5->SaveAs("results/emcal_barrel_pions_pid.png"); + c5->SaveAs("results/emcal_barrel_pions_pid.pdf"); + //Energy Resolution Work auto hdE = d1.Histo1D({"hdE", "dE; dE[GeV]; Events", 20, -7.5, 7.5}, "dE");//changed from 100 auto hdE_rel = d1.Histo1D({"hdE_rel", "dE Relative; dE Relative; Events", 20, -2, 2}, "dE_rel");//changed from 100