Implement detector benchmark for the barrel calorimeter for energy resolution

benchmarks/barrel_ecal/scripts/emcal_barrel_energy_scan_analysis.cxx

+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/CalorimeterHitCollection.h"
+
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TMath.h"
+#include "TH1.h"
+#include "TF1.h"
+#include "TH1D.h"
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+  // Setting for graphs
+  gROOT->SetStyle("Plain");
+  gStyle->SetOptFit(1);
+  gStyle->SetLineWidth(2);
+  gStyle->SetPadTickX(1);
+  gStyle->SetPadTickY(1);
+  gStyle->SetPadGridX(1);
+  gStyle->SetPadGridY(1);
+  gStyle->SetPadLeftMargin(0.14);
+  gStyle->SetPadRightMargin(0.14);
+
+
+void save_canvas(TCanvas* c, std::string label, double E)
+{
+    c->SaveAs(std::format("results/{}_{}.png",label, E));
+    c->SaveAs(std::format("results/{}_{}.pdf",label, E));
+}
+
+void save_canvas(TCanvas* c, std::string label)
+{
+    c->SaveAs(std::format("results/{}.png",label));
+    c->SaveAs(std::format("results/{}.pdf",label));
+}
+
+std::tuple <double, double> extract_sampling_fraction_parameters(double E)
+{
+  auto input_fname = std::format("sim_output/energy_scan/{}/emcal_barrel_uniform_electrons.root", E);
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame d0("events", input_fname);
+
+  // Thrown Energy [GeV]
+  auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+    auto p = input[2];
+    auto energy = TMath::Sqrt(p.psx * p.psx + p.psy * p.psy + p.psz * p.psz + p.mass * p.mass);
+    return energy;
+  };
+
+  // Number of hits
+  auto nhits = [] (const std::vector<dd4pod::CalorimeterHitData>& evt) {return (int) evt.size(); };
+
+  // Energy deposition [GeV]
+  auto Esim = [](const std::vector<dd4pod::CalorimeterHitData>& evt) {
+    auto total_edep = 0.0;
+    for (const auto& i: evt)
+      total_edep += i.energyDeposit;
+    return total_edep;
+  };
+
+  // Sampling fraction = Esampling / Ethrown
+  auto fsam = [](const double sampled, const double thrown) {
+    return sampled / thrown;
+  };
+
+  // Define variables
+  auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
+                .Define("nhits", nhits, {"EcalBarrelHits"})
+                .Define("Esim", Esim, {"EcalBarrelHits"})
+                .Define("fsam", fsam, {"Esim", "Ethr"});
+
+  // 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", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1},
+      "fsam");
+
+  // Event Counts
+  auto nevents_thrown = d1.Count();
+  std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
+
+  // Draw Histograms
+  {
+    TCanvas* c1 = new TCanvas("c1", "c1", 700, 500);
+    c1->SetLogy(1);
+    auto h = hEthr->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    save_canvas(c1, "emcal_barrel_electrons_Ethr", E);  
+  }
+  std::cout << "derp1\n";
+
+  {
+    TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
+    c2->SetLogy(1);
+    auto h = hNhits->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    save_canvas(c2, "emcal_barrel_electrons_nhits", E);
+  }
+
+  {
+    TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
+    c3->SetLogy(1);
+    auto h = hEsim->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    save_canvas(c3, "emcal_barrel_electrons_Esim", E);
+  }
+
+  {
+    TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
+    c4->SetLogy(1);
+    auto h = hfsam->DrawCopy();
+    //h->GetYaxis()->SetTitleOffset(1.4);
+    h->SetLineWidth(2);
+    h->SetLineColor(kBlue);
+    h->Fit("gaus", "", "", 0.01, 0.1);
+    TF1 *gaus = h->GetFunction("gaus");
+    gaus->SetLineWidth(2);
+    gaus->SetLineColor(kRed);    
+    double mean = gaus->GetParameter(2);
+    double sigma = gaus->GetParameter(3);
+    save_canvas(c4, "emcal_barrel_electrons_fsam", E);
+    return std::make_tuple(mean, sigma);
+  }
+}
+
+std::vector<double> read_scanned_energies(std::string input_energies_fname)
+{
+    std::vector<double> scanned_energies;
+    double E;
+    ifstream E_file (input_energies_fname);
+    if (E_file.is_open())
+    {
+        while (E_file >> E)
+        {
+            scanned_energies.push_back(E);
+        }
+        myfile.close();
+        return scanned_energies;
+    }
+    else
+    {
+        std::cout << std::format("Unable to open file {}", input_energies_fname);
+        abort();
+    }
+};
+
+void emcal_barrel_electrons_energy_scan_analysis()
+{
+    vector<double> scanned_energies = read_scanned_energies("emcal_barrel_energy_scan_points.txt");
+    TGraph gr_fsam(scanned_energies.size());
+    TGraph gr_fsam_res(scanned_energies.size());
+    for (const auto& E : scanned_energies) {
+        auto [fsam, fsam_res] = extract_sampling_fraction_parameters(E);
+        gr_fsam.AddPoint(E,fsam);
+        gr_fsam_res.AddPoint(E,fsam_res);
+    }
+    TCanvas* c5 = new TCanvas("c5", "c5", 700, 500);
+    c5->cd();
+    gr_fsam.Draw("AP");
+    save_canvas(c5,"emcal_barrel_electrons_fsam_scan");
+
+    TCanvas* c6 = new TCanvas("c6", "c6", 700, 500);
+    c6->cd();
+    gr_fsam_res.Draw("AP");
+    save_canvas(c6,"emcal_barrel_electrons_fsam_scan_res");
+}