barrel_ecal/scripts/emcal_barrel_energy_scan_analysis.cxx

0001 ////////////////////////////////////////
0002 // Read reconstruction ROOT output file
0003 // Plot variables
0004 // M.Żurek
0005 ////////////////////////////////////////
0006
0007 #include "ROOT/RDataFrame.hxx"
0008 #include <iostream>
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 #include "DD4hep/Detector.h"
0023 #include "DDG4/Geant4Data.h"
0024 #include "DDRec/CellIDPositionConverter.h"
0025
0026 using ROOT::RDataFrame;
0027 using namespace ROOT::VecOps;
0028
0029
0030 void save_canvas(TCanvas* c, std::string label)
0031 {
0032   c->SaveAs(fmt::format("results/energy_scan/{}.png",label).c_str());
0033 }
0034
0035 void save_canvas(TCanvas* c, std::string label, double E)
0036 {
0037   std::string label_with_E = fmt::format("{}/{}", E, label);
0038   save_canvas(c, label_with_E);
0039 }
0040 void save_canvas(TCanvas* c, std::string label, std::string E_label)
0041 {
0042   std::string label_with_E = fmt::format("{}/{}", E_label, label);
0043   save_canvas(c, label_with_E);
0044 }
0045 void save_canvas(TCanvas* c, std::string var_label, std::string E_label, std::string particle_label)
0046 {
0047   std::string label_with_E = fmt::format("{}/emcal_barrel_{}_{}", E_label, particle_label, var_label);
0048   save_canvas(c, label_with_E);
0049 }
0050
0051 std::tuple <double, double, double, double> extract_sampling_fraction_parameters(std::string particle_label, std::string E_label, dd4hep::Detector& detector)
0052 {
0053   std::string input_fname = fmt::format("sim_output/energy_scan/{}/sim_emcal_barrel_{}.edm4hep.root", E_label, particle_label);
0054   ROOT::EnableImplicitMT();
0055   ROOT::RDataFrame d0("events", input_fname);
0056
0057   // Thrown Energy [GeV]
0058   auto Ethr = [](std::vector<edm4hep::MCParticleData> const& input) {
0059     auto p = input[2];
0060     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);
0061     return energy;
0062   };
0063
0064   // Number of hits
0065   auto nhits = [] (const std::vector<edm4hep::SimCalorimeterHitData>& evt) {return (int) evt.size(); };
0066
0067   // Energy deposition [GeV]
0068   auto Esim = [](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0069     auto total_edep = 0.0;
0070     for (const auto& i: evt)
0071       total_edep += i.energy;
0072     return total_edep;
0073   };
0074
0075   // Sampling fraction = Esampling / Ethrown
0076   auto fsam = [](const double sampled, const double thrown) {
0077     return sampled / thrown;
0078   };
0079
0080   // Energy deposited in layers
0081   auto decoder = detector.readout("EcalBarrelImagingHits").idSpec().decoder();
0082   fmt::print("{}\n", decoder->fieldDescription());
0083   auto layer_index = decoder->index("layer");
0084   fmt::print(" layer index is {}.\n", layer_index);
0085
0086   // Define variables
0087   auto d1 = d0.Define("Ethr", Ethr, {"MCParticles"})
0088                 .Define("nhits", nhits, {"EcalBarrelImagingHits"})
0089                 .Define("EsimImg", Esim, {"EcalBarrelImagingHits"})
0090                 .Define("EsimScFi", Esim, {"EcalBarrelScFiHits"})
0091                 .Define("Esim", "EsimImg+EsimScFi")
0092                 .Define("fsam", fsam, {"Esim", "Ethr"});
0093
0094   // Define Histograms
0095   auto hEthr = d1.Histo1D(
0096       {"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 25.0},
0097       "Ethr");
0098   auto hNhits =
0099       d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events",
0100                   100, 0.0, 2000.0},
0101                  "nhits");
0102   auto hEsim = d1.Histo1D(
0103       {"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 500, 0.0, 0.5},
0104       "Esim");
0105   auto hfsam = d1.Histo1D(
0106       {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 800, 0.0, 0.2},
0107       "fsam");
0108
0109   // Number of layers to read the edep from
0110   int nlayers = 7;
0111
0112   TGraphErrors gr_no_edep(nlayers);
0113   TGraphErrors gr_edep_mean(nlayers);
0114
0115   // Draw energy per layer
0116   TCanvas* clayer = new TCanvas("clayer", "clayer", 1250, 1000);
0117   clayer->Divide(4,5);
0118
0119   for(int layer=1; layer<nlayers+1; layer++) {
0120     auto Esim_layer = [=](const std::vector<edm4hep::SimCalorimeterHitData>& evt) {
0121       auto layer_edep = 0.0;
0122       for (const auto& i: evt) {
0123         if (decoder->get(i.cellID, layer_index) == layer) {
0124           layer_edep += i.energy;
0125         }
0126       }
0127       return layer_edep;
0128     };
0129
0130     auto d2 = d0.Define(fmt::format("EsimImg_layer_{}",layer).c_str(), Esim_layer, {"EcalBarrelImagingHits"})
0131                   .Define(fmt::format("EsimScFi_layer_{}",layer).c_str(), Esim_layer, {"EcalBarrelScFiHits"})
0132                   .Define(fmt::format("Esim_layer_{}",layer).c_str(), fmt::format("EsimImg_layer_{}+EsimScFi_layer_{}",layer,layer).c_str());
0133     std::cout << "Layer to process: " << layer << std::endl;
0134
0135     auto hEsim_layer = d2.Histo1D({fmt::format("hEsim_layer_{}",layer).c_str(),fmt::format("Energy Deposit in layer {}; Energy Deposit [Gev]; Events",layer).c_str(), 200, 0.0, 0.04}, fmt::format("Esim_layer_{}",layer));
0136
0137     clayer->cd(layer);
0138     gPad->SetLogy();
0139     auto h = hEsim_layer->DrawCopy();
0140
0141     auto up_range = h->GetMean() + 3*h->GetStdDev();
0142     auto down_range = h->GetMean() - 3*h->GetStdDev();
0143     h->SetLineWidth(2);
0144     h->SetLineColor(kBlue);
0145
0146     h->GetXaxis()->SetRange(h->GetXaxis()->GetBinUpEdge(1), up_range); // skip 0th bin
0147     auto mean_layer = h->GetMean();
0148     auto rms_layer = h->GetStdDev();
0149     h->GetXaxis()->SetRange(); // reset the range
0150     h->GetXaxis()->SetRangeUser(0.,up_range);
0151
0152     auto no_edep = (h->GetBinContent(1)/h->GetEntries())*100;
0153     gr_no_edep.SetPoint(gr_no_edep.GetN(),layer,no_edep);
0154     gr_edep_mean.SetPoint(gr_edep_mean.GetN(),layer, mean_layer);
0155     gr_edep_mean.SetPointError(gr_edep_mean.GetN()-1,0, rms_layer);
0156   }
0157   save_canvas(clayer, "Esim_layer", E_label, particle_label);
0158
0159   // Event Counts
0160   auto nevents_thrown = d1.Count();
0161   std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
0162
0163   // Draw Histograms and Graphs for every energy
0164   TCanvas* cnoedep = new TCanvas("c_noedep", "c_noedep", 700, 500);
0165   cnoedep->cd();
0166   gr_no_edep.SetTitle("% of events with no dE;Layer;Events with no dE [%]");
0167   gr_no_edep.SetMarkerStyle(20);
0168   gr_no_edep.Draw("AP");
0169   save_canvas(cnoedep, "Layer_nodep", E_label, particle_label);
0170
0171   TCanvas* cmean = new TCanvas("c_mean", "c_mean", 700, 500);
0172   cmean->cd();
0173   gr_edep_mean.SetTitle("Mean and RMS of energy deposit;Layer;Mean dE [GeV]");
0174   gr_edep_mean.GetYaxis()->SetTitleOffset(1.4);
0175   gr_edep_mean.SetFillColor(4);
0176   gr_edep_mean.SetFillStyle(3010);
0177   gr_edep_mean.Draw("a3P");
0178   save_canvas(cmean, "Layer_Esim_mean", E_label, particle_label);
0179
0180   {
0181     TCanvas* c1 = new TCanvas("c1", "c1", 700, 500);
0182     c1->SetLogy(1);
0183     auto h = hEthr->DrawCopy();
0184     h->SetLineWidth(2);
0185     h->SetLineColor(kBlue);
0186     save_canvas(c1, "Ethr", E_label, particle_label);
0187   }
0188
0189   {
0190     TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
0191     c2->SetLogy(1);
0192     auto h = hNhits->DrawCopy();
0193     h->SetLineWidth(2);
0194     h->SetLineColor(kBlue);
0195     save_canvas(c2, "nhits", E_label, particle_label);
0196   }
0197
0198   {
0199     TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
0200     c3->SetLogy(1);
0201     auto h = hEsim->DrawCopy();
0202     h->SetLineWidth(2);
0203     h->SetLineColor(kBlue);
0204     double up_fit = h->GetMean() + 5*h->GetStdDev();
0205     double down_fit = h->GetMean() - 5*h->GetStdDev();
0206     h->GetXaxis()->SetRangeUser(0.,up_fit);
0207     save_canvas(c3, "Esim", E_label, particle_label);
0208   }
0209
0210   {
0211     TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
0212     auto h = hfsam->DrawCopy();
0213     h->SetLineWidth(2);
0214     h->SetLineColor(kBlue);
0215     double up_fit = h->GetMean() + 5*h->GetStdDev();
0216     double down_fit = h->GetMean() - 5*h->GetStdDev();
0217     if(down_fit <=0 ) down_fit = h->GetXaxis()->GetBinUpEdge(1);
0218     h->Fit("gaus", "", "", down_fit, up_fit);
0219     h->GetXaxis()->SetRangeUser(0.,up_fit);
0220     TF1 *gaus = h->GetFunction("gaus");
0221     gaus->SetLineWidth(2);
0222     gaus->SetLineColor(kRed);
0223     double mean = gaus->GetParameter(1);
0224     double sigma = gaus->GetParameter(2);
0225     double mean_err = gaus->GetParError(1);
0226     double sigma_err = gaus->GetParError(2);
0227     save_canvas(c4, "fsam", E_label, particle_label);
0228     return std::make_tuple(mean, sigma, mean_err, sigma_err);
0229   }
0230 }
0231
0232 std::vector<std::string> read_scanned_energies(std::string input_energies_fname)
0233 {
0234     std::vector<std::string> scanned_energies;
0235     std::string E_label;
0236     ifstream E_file (input_energies_fname);
0237     if (E_file.is_open())
0238     {
0239         while (E_file >> E_label)
0240         {
0241             scanned_energies.push_back(E_label);
0242         }
0243         E_file.close();
0244         return scanned_energies;
0245     }
0246     else
0247     {
0248         std::cout << "Unable to open file " << input_energies_fname << std::endl;
0249         abort();
0250     }
0251 }
0252
0253 void emcal_barrel_energy_scan_analysis(std::string particle_label = "electron")
0254 {
0255
0256   // Setting for graphs
0257   gROOT->SetStyle("Plain");
0258   gStyle->SetOptFit(1);
0259   gStyle->SetLineWidth(2);
0260   gStyle->SetPadTickX(1);
0261   gStyle->SetPadTickY(1);
0262   gStyle->SetPadGridX(1);
0263   gStyle->SetPadGridY(1);
0264   gStyle->SetPadLeftMargin(0.14);
0265   gStyle->SetPadRightMargin(0.14);
0266
0267   auto scanned_energies = read_scanned_energies(fmt::format("sim_output/emcal_barrel_energy_scan_points_{}.txt", particle_label));
0268
0269   //Take detector layers
0270   std::string detector_path = "";
0271   std::string detector_name = "athena";
0272   if(std::getenv("DETECTOR_PATH")) {
0273     detector_path = std::getenv("DETECTOR_PATH");
0274   }
0275   if(std::getenv("DETECTOR_CONFIG")) {
0276     detector_name = std::getenv("DETECTOR_CONFIG");
0277   }
0278
0279   dd4hep::Detector& detector = dd4hep::Detector::getInstance();
0280   detector.fromCompact(fmt::format("{}/{}.xml", detector_path, detector_name));
0281
0282   TGraphErrors gr_fsam(scanned_energies.size()-1);
0283   TGraphErrors gr_fsam_res(scanned_energies.size()-1);
0284
0285   for (const auto& E_label : scanned_energies) {
0286       auto [fsam, fsam_res, fsam_err, fsam_res_err] = extract_sampling_fraction_parameters(particle_label, E_label, detector);
0287       auto E = std::stod(E_label);
0288
0289       gr_fsam.SetPoint(gr_fsam.GetN(),E,100*fsam);
0290       gr_fsam.SetPointError(gr_fsam.GetN()-1,0., 100*fsam_err);
0291       gr_fsam_res.SetPoint(gr_fsam_res.GetN(),E,100.0*(fsam_res/fsam));
0292       auto fsam_res_rel_err = 100.0*(sqrt(pow((fsam_res_err/fsam),2)+pow((fsam_err*fsam_res)/(fsam*fsam),2)));
0293       gr_fsam_res.SetPointError(gr_fsam_res.GetN()-1,0.,fsam_res_rel_err);
0294   }
0295
0296   TCanvas* c5 = new TCanvas("c5", "c5", 700, 500);
0297   c5->cd();
0298   gr_fsam.SetTitle("Sampling Fraction Scan;True Energy [GeV];Sampling Fraction [%]");
0299   gr_fsam.SetMarkerStyle(20);
0300   gr_fsam.Fit("pol0", "", "", 2., 20.);
0301   gr_fsam.Draw("APE");
0302   save_canvas(c5, fmt::format("emcal_barrel_{}_fsam_scan", particle_label));
0303
0304   TCanvas* c6 = new TCanvas("c6", "c6", 700, 500);
0305   c6->cd();
0306   TF1* func_res = new TF1("func_res", "[0]/sqrt(x) + [1]", 0.25, 20.);
0307   func_res->SetLineWidth(2);
0308   func_res->SetLineColor(kRed);
0309   gr_fsam_res.SetTitle("Energy Resolution;True Energy [GeV];#Delta E/E [%]");
0310   gr_fsam_res.SetMarkerStyle(20);
0311   gr_fsam_res.Fit(func_res,"R");
0312   gr_fsam_res.Draw("APE");
0313   save_canvas(c6,fmt::format("emcal_barrel_{}_fsam_scan_res", particle_label));
0314
0315 }