barrel_hcal/scripts/hcal_barrel_energy_scan_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 #include "TGraphErrors.h"
0021
0022 using ROOT::RDataFrame;
0023 using namespace ROOT::VecOps;
0024
0025
0026 void save_canvas(TCanvas* c, std::string label)
0027 {
0028   c->SaveAs(fmt::format("results/energy_scan/{}.png",label).c_str());
0029   c->SaveAs(fmt::format("results/energy_scan/{}.pdf",label).c_str());
0030 }
0031
0032 void save_canvas(TCanvas* c, std::string label, double E)
0033 {
0034   std::string label_with_E = fmt::format("{}/{}", E, label);
0035   save_canvas(c, label_with_E);
0036 }
0037 void save_canvas(TCanvas* c, std::string label, std::string E_label)
0038 {
0039   std::string label_with_E = fmt::format("{}/{}", E_label, label);
0040   save_canvas(c, label_with_E);
0041 }
0042 void save_canvas(TCanvas* c, std::string var_label, std::string E_label, std::string particle_label)
0043 {
0044   std::string label_with_E = fmt::format("{}/hcal_barrel_{}_{}", E_label, particle_label, var_label);
0045   save_canvas(c, label_with_E);
0046 }
0047
0048 std::tuple <double, double, double, double> extract_sampling_fraction_parameters(std::string particle_label, std::string E_label)
0049 {
0050   std::string input_fname = fmt::format("sim_output/energy_scan/{}/sim_hcal_barrel_{}.edm4hep.root", E_label, particle_label);
0051   ROOT::EnableImplicitMT();
0052   ROOT::RDataFrame d0("events", input_fname);
0053
0054   // Thrown Energy [GeV]
0055   auto Ethr = [](std::vector<edm4hep::MCParticleData> const& input) {
0056     auto p = input[2];
0057     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);
0058     return energy;
0059   };
0060
0061   // Number of hits
0062   auto nhits = [] (const std::vector<edm4hep::SimCalorimeterHitData>& evt) {return (int) evt.size(); };
0063
0064   // Energy deposition [GeV]
0065   auto Esim = [](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0066     auto total_edep = 0.0;
0067     for (const auto& i: evt)
0068       total_edep += i.energy;
0069     return total_edep;
0070   };
0071
0072   // Sampling fraction = Esampling / Ethrown
0073   auto fsam = [](const double sampled, const double thrown) {
0074     return sampled / thrown;
0075   };
0076
0077   // Define variables
0078   auto d1 = d0.Define("Ethr", Ethr, {"MCParticles"});
0079
0080   // Define assumptions
0081   auto Ethr_max = 25.0;
0082   auto fsam_estimate = 1.0;
0083   if (d1.HasColumn("EcalBarrelScFiHits")) {
0084     d1 = d1.Define("nhits", nhits, {"EcalBarrelImagingHits"})
0085            .Define("Esim", Esim, {"EcalBarrelImagingHits"})
0086            .Define("fsam", fsam, {"Esim", "Ethr"});
0087     fsam_estimate = 0.1;
0088   } else {
0089     d1 = d1.Define("nhits", nhits, {"EcalBarrelSciGlassHits"})
0090            .Define("Esim", Esim, {"EcalBarrelSciGlassHits"})
0091            .Define("fsam", fsam, {"Esim", "Ethr"});
0092     fsam_estimate = 1.0;
0093   }
0094
0095   // Define Histograms
0096   auto hEthr = d1.Histo1D(
0097       {"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, Ethr_max},
0098       "Ethr");
0099   auto hNhits =
0100       d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events",
0101                   100, 0.0, 2000.0},
0102                  "nhits");
0103   auto hEsim = d1.Histo1D(
0104       {"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 500, 0.0, fsam_estimate * Ethr_max},
0105       "Esim");
0106   auto hfsam = d1.Histo1D(
0107       {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 200, 0.0, 2.0 * fsam_estimate},
0108       "fsam");
0109
0110   // Event Counts
0111   auto nevents_thrown = d1.Count();
0112   std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
0113
0114   // Draw Histograms
0115   {
0116     TCanvas* c1 = new TCanvas("c1", "c1", 700, 500);
0117     c1->SetLogy(1);
0118     auto h = hEthr->DrawCopy();
0119     //h->GetYaxis()->SetTitleOffset(1.4);
0120     h->SetLineWidth(2);
0121     h->SetLineColor(kBlue);
0122     save_canvas(c1, "Ethr", E_label, particle_label);
0123   }
0124
0125   {
0126     TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
0127     c2->SetLogy(1);
0128     auto h = hNhits->DrawCopy();
0129     //h->GetYaxis()->SetTitleOffset(1.4);
0130     h->SetLineWidth(2);
0131     h->SetLineColor(kBlue);
0132     save_canvas(c2, "nhits", E_label, particle_label);
0133   }
0134
0135   {
0136     TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
0137     c3->SetLogy(1);
0138     auto h = hEsim->DrawCopy();
0139     //h->GetYaxis()->SetTitleOffset(1.4);
0140     h->SetLineWidth(2);
0141     h->SetLineColor(kBlue);
0142     double up_fit = h->GetMean() + 5*h->GetStdDev();
0143     double down_fit = h->GetMean() - 5*h->GetStdDev();
0144     h->GetXaxis()->SetRangeUser(0.,up_fit);
0145     save_canvas(c3, "Esim", E_label, particle_label);
0146   }
0147
0148   {
0149     TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
0150     //c4->SetLogy(1);
0151     auto h = hfsam->DrawCopy();
0152     //h->GetYaxis()->SetTitleOffset(1.4);
0153     h->SetLineWidth(2);
0154     h->SetLineColor(kBlue);
0155     double up_fit = h->GetMean() + 5*h->GetStdDev();
0156     double down_fit = h->GetMean() - 5*h->GetStdDev();
0157     h->Fit("gaus", "", "", down_fit, up_fit);
0158     h->GetXaxis()->SetRangeUser(0.,up_fit);
0159     TF1 *gaus = h->GetFunction("gaus");
0160     gaus->SetLineWidth(2);
0161     gaus->SetLineColor(kRed);
0162     double mean = gaus->GetParameter(1);
0163     double sigma = gaus->GetParameter(2);
0164     double mean_err = gaus->GetParError(1);
0165     double sigma_err = gaus->GetParError(2);
0166     save_canvas(c4, "fsam", E_label, particle_label);
0167     return std::make_tuple(mean, sigma, mean_err, sigma_err);
0168   }
0169 }
0170
0171 std::vector<std::string> read_scanned_energies(std::string input_energies_fname)
0172 {
0173     std::vector<std::string> scanned_energies;
0174     std::string E_label;
0175     ifstream E_file (input_energies_fname);
0176     if (E_file.is_open())
0177     {
0178         while (E_file >> E_label)
0179         {
0180             scanned_energies.push_back(E_label);
0181         }
0182         E_file.close();
0183         return scanned_energies;
0184     }
0185     else
0186     {
0187         std::cout << "Unable to open file " << input_energies_fname << std::endl;
0188         abort();
0189     }
0190 }
0191
0192 void hcal_barrel_energy_scan_analysis(std::string particle_label = "electron")
0193 {
0194   // Setting for graphs
0195   gROOT->SetStyle("Plain");
0196   gStyle->SetOptFit(1);
0197   gStyle->SetLineWidth(2);
0198   gStyle->SetPadTickX(1);
0199   gStyle->SetPadTickY(1);
0200   gStyle->SetPadGridX(1);
0201   gStyle->SetPadGridY(1);
0202   gStyle->SetPadLeftMargin(0.14);
0203   gStyle->SetPadRightMargin(0.14);
0204
0205     auto scanned_energies = read_scanned_energies(fmt::format("sim_output/hcal_barrel_energy_scan_points_{}.txt", particle_label));
0206
0207     TGraphErrors gr_fsam(scanned_energies.size()-1);
0208     TGraphErrors gr_fsam_res(scanned_energies.size()-1);
0209
0210     for (const auto& E_label : scanned_energies) {
0211         auto [fsam, fsam_res, fsam_err, fsam_res_err] = extract_sampling_fraction_parameters(particle_label, E_label);
0212         auto E = std::stod(E_label);
0213
0214         gr_fsam.SetPoint(gr_fsam.GetN(),E,100*fsam);
0215         gr_fsam.SetPointError(gr_fsam.GetN()-1,0., 100*fsam_err);
0216         gr_fsam_res.SetPoint(gr_fsam_res.GetN(),E,100.0*(fsam_res/fsam));
0217         auto fsam_res_rel_err = 100.0*(sqrt(pow((fsam_res_err/fsam),2)+pow((fsam_err*fsam_res)/(fsam*fsam),2)));
0218         gr_fsam_res.SetPointError(gr_fsam_res.GetN()-1,0.,fsam_res_rel_err);
0219     }
0220
0221     TCanvas* c5 = new TCanvas("c5", "c5", 700, 500);
0222     c5->cd();
0223     gr_fsam.SetTitle("Sampling Fraction Scan;True Energy [GeV];Sampling Fraction [%]");
0224     gr_fsam.SetMarkerStyle(20);
0225     gr_fsam.Fit("pol0", "", "", 2., 20.);
0226     gr_fsam.Draw("APE");
0227     save_canvas(c5, fmt::format("hcal_barrel_{}_fsam_scan", particle_label));
0228
0229     TCanvas* c6 = new TCanvas("c6", "c6", 700, 500);
0230     c6->cd();
0231     TF1* func_res = new TF1("func_res", "[0]/sqrt(x) + [1]", 0.25, 20.);
0232     func_res->SetLineWidth(2);
0233     func_res->SetLineColor(kRed);
0234     gr_fsam_res.SetTitle("Energy Resolution;True Energy [GeV];#Delta E/E [%]");
0235     gr_fsam_res.SetMarkerStyle(20);
0236     gr_fsam_res.Fit(func_res,"R");
0237     gr_fsam_res.Draw("APE");
0238     save_canvas(c6,fmt::format("hcal_barrel_{}_fsam_scan_res", particle_label));
0239 }