barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx

0001 ////////////////////////////////////////
0002 // Read reconstruction ROOT output file
0003 // Plot variables
0004 ////////////////////////////////////////
0005
0006 #include "ROOT/RDataFrame.hxx"
0007 #include <iostream>
0008 #include <fstream>
0009 #include <fmt/core.h>
0010
0011 #include "edm4hep/MCParticleCollection.h"
0012 #include "edm4hep/SimCalorimeterHitCollection.h"
0013
0014 #include "TCanvas.h"
0015 #include "TStyle.h"
0016 #include "TMath.h"
0017 #include "TH1.h"
0018 #include "TF1.h"
0019 #include "TH1D.h"
0020
0021 #include "hcal_barrel_common_functions.h"
0022
0023 using ROOT::RDataFrame;
0024 using namespace ROOT::VecOps;
0025
0026 void save_canvas(TCanvas* c, std::string label)
0027 {
0028   c->SaveAs(fmt::format("results/{}.png",label).c_str());
0029   c->SaveAs(fmt::format("results/{}.pdf",label).c_str());
0030 }
0031
0032 void save_canvas(TCanvas* c, std::string label, std::string particle_label)
0033 {
0034   std::string label_with_E = fmt::format("hcal_barrel_{}_{}", particle_label, label);
0035   save_canvas(c, label_with_E);
0036 }
0037
0038 void hcal_barrel_particles_analysis(std::string particle_name = "electron")
0039 {
0040   // Setting for graphs
0041   gROOT->SetStyle("Plain");
0042   gStyle->SetOptFit(1);
0043   gStyle->SetLineWidth(2);
0044   gStyle->SetPadTickX(1);
0045   gStyle->SetPadTickY(1);
0046   gStyle->SetPadGridX(1);
0047   gStyle->SetPadGridY(1);
0048   gStyle->SetPadLeftMargin(0.14);
0049   gStyle->SetPadRightMargin(0.14);
0050
0051   ROOT::EnableImplicitMT();
0052   std::string input_fname = fmt::format("sim_output/sim_hcal_barrel_{}.edm4hep.root", particle_name);
0053   ROOT::RDataFrame d0("events", input_fname);
0054
0055   // Thrown Energy [GeV]
0056   auto Ethr = [](std::vector<edm4hep::MCParticleData> const& input) {
0057     auto p = input[2];
0058     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);
0059     return energy;
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     return total_edep;
0071   };
0072
0073   // Sampling fraction = Esampling / Ethrown
0074   auto fsam = [](const double sampled, const double thrown) {
0075     return sampled / thrown;
0076   };
0077
0078   // Define variables
0079   auto d1 = d0.Define("Ethr", Ethr, {"MCParticles"})
0080                 .Define("nhits", nhits, {"HcalBarrelHits"})
0081                 .Define("Esim", Esim, {"HcalBarrelHits"})
0082                 .Define("fsam", fsam, {"Esim", "Ethr"});
0083
0084   // Define Histograms
0085   auto hEthr = d1.Histo1D(
0086       {"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 25.0},
0087       "Ethr");
0088   auto hNhits =
0089       d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events", 100, 0.0, 2000.0},
0090       "nhits");
0091   auto hEsim = d1.Histo1D(
0092       {"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 500, 0.0, 0.5},
0093       "Esim");
0094   auto hfsam = d1.Histo1D(
0095       {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 200, 0.0, 0.05},
0096       "fsam");
0097
0098   // Event Counts
0099   auto nevents_thrown = d1.Count();
0100   std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
0101
0102   // Draw Histograms
0103   {
0104     TCanvas* c1 = new TCanvas("c1", "c1", 700, 500);
0105     c1->SetLogy(1);
0106     auto h = hEthr->DrawCopy();
0107     h->SetLineWidth(2);
0108     h->SetLineColor(kBlue);
0109     save_canvas(c1,"Ethr",particle_name);
0110   }
0111
0112   {
0113     TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
0114     c2->SetLogy(1);
0115     auto h = hNhits->DrawCopy();
0116     h->SetLineWidth(2);
0117     h->SetLineColor(kBlue);
0118     save_canvas(c2,"nhits",particle_name);
0119   }
0120
0121   {
0122     TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
0123     c3->SetLogy(1);
0124     auto h = hEsim->DrawCopy();
0125     h->SetLineWidth(2);
0126     h->SetLineColor(kBlue);
0127     double up_fit = h->GetMean() + 5*h->GetStdDev();
0128     double down_fit = h->GetMean() - 5*h->GetStdDev();
0129     h->GetXaxis()->SetRangeUser(0.,up_fit);
0130     save_canvas(c3,"Esim",particle_name);
0131   }
0132
0133   {
0134     TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
0135     c4->SetLogy(1);
0136     auto h = hfsam->DrawCopy();
0137     h->SetLineWidth(2);
0138     h->SetLineColor(kBlue);
0139     double up_fit = h->GetMean() + 5*h->GetStdDev();
0140     double down_fit = h->GetMean() - 5*h->GetStdDev();
0141     h->Fit("gaus", "", "", down_fit, up_fit);
0142     h->GetXaxis()->SetRangeUser(0.,up_fit);
0143     TF1 *gaus = h->GetFunction("gaus");
0144     gaus->SetLineWidth(2);
0145     gaus->SetLineColor(kRed);
0146     save_canvas(c4,"fsam",particle_name);
0147   }
0148 }
0149