diff --git a/benchmarks/imaging_ecal/scripts/emcal_barrel_energy_scan_analysis_pi_e.cxx b/benchmarks/imaging_ecal/scripts/emcal_barrel_energy_scan_analysis_pi_e.cxx
index 8dc61a822de2c3e1fae092972d7181a2ce75bfd8..18f58681bf2db704d7400b7ad1b27b9ba80633cb 100644
--- a/benchmarks/imaging_ecal/scripts/emcal_barrel_energy_scan_analysis_pi_e.cxx
+++ b/benchmarks/imaging_ecal/scripts/emcal_barrel_energy_scan_analysis_pi_e.cxx
@@ -33,6 +33,25 @@
using ROOT::RDataFrame;
using namespace ROOT::VecOps;
+struct infoHist1D {
+ string name; // histo name
+ string column; // column name
+ string titles; // histo and axis titles
+ int nbins; // Nb of bins
+ double x0; // Lower hist limit
+ double x1; // Upper hist limit
+};
+
+struct infoHist2D {
+ string name; // Plot variable name
+ int nbinsx; // Nb of bins in x
+ double x0; // Lower hist limit in x
+ double x1; // Upper hist limit in x
+ int nbinsy; // Nb of bins in y
+ double y0; // Lower hist limit in y
+ double y1; // Upper hist limit in y
+};
+
void save_canvas(TCanvas* c, std::string label)
{
c->SaveAs(fmt::format("results/{}.png",label).c_str());
@@ -64,12 +83,19 @@ void set_histo_range(TH1D* h)
std::tuple <double, double> extract_sampling_fraction_parameters(std::string E_label)
{
std::string input_fname_electron = fmt::format("sim_output/rec_emcal_barrel_{}_{}.root", "electron", E_label);
- ROOT::EnableImplicitMT();
- ROOT::RDataFrame d0("events", input_fname_electron);
-
std::string input_fname_pionm = fmt::format("sim_output/rec_emcal_barrel_{}_{}.root", "pion-", E_label);
+
ROOT::EnableImplicitMT();
- ROOT::RDataFrame d2("events", input_fname_pionm);
+ ROOT::RDataFrame d0("events", {input_fname_electron, input_fname_pionm});
+
+
+ auto isElectron = [](std::vector<dd4pod::Geant4ParticleData> const& input){
+ return (input[2].pdgID == 11 ? true : false);
+ };
+
+ auto isPion = [](std::vector<dd4pod::Geant4ParticleData> const& input){
+ return (input[2].pdgID == -211 ? true : false);
+ };
// Thrown Energy [GeV]
auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
@@ -118,126 +144,149 @@ std::tuple <double, double> extract_sampling_fraction_parameters(std::string E_l
};
// Define variables
- auto delectron = d0.Define("Ethr", Ethr, {"mcparticles2"})
- .Define("mom", momentum, {"mcparticles2"})
- .Define("mom_smeared", momentum_smeared, {"mom"})
+ auto d1 = d0.Define("Ethr", Ethr, {"mcparticles2"})
+ .Define("Mom", momentum, {"mcparticles2"})
+ .Define("MomSmeared", momentum_smeared, {"mom"})
.Define("ErecImg", Erec, {"RecoEcalBarrelImagingHits"})
.Define("ErecScFi", Erec, {"EcalBarrelScFiHitsReco"})
.Define("NClusterImg", ncluster, {"EcalBarrelImagingClusters"})
.Define("NClusterScFi", ncluster, {"EcalBarrelScFiClusters"})
.Define("EClusterImg", Ecluster, {"EcalBarrelImagingClusters"})
.Define("EClusterScFi", Ecluster, {"EcalBarrelScFiClusters"})
- .Define("E_over_p_ClusterScFi", E_over_mom, {"EClusterScFi", "mom_smeared"})
- .Define("E_over_p_RecoScFi", E_over_mom, {"ErecScFi", "mom_smeared"});
-
- auto dpion = d2.Define("Ethr", Ethr, {"mcparticles2"})
- .Define("mom", momentum, {"mcparticles2"})
- .Define("mom_smeared", momentum_smeared, {"mom"})
- .Define("ErecImg", Erec, {"RecoEcalBarrelImagingHits"})
- .Define("ErecScFi", Erec, {"EcalBarrelScFiHitsReco"})
- .Define("NClusterImg", ncluster, {"EcalBarrelImagingClusters"})
- .Define("NClusterScFi", ncluster, {"EcalBarrelScFiClusters"})
- .Define("EClusterImg", Ecluster, {"EcalBarrelImagingClusters"})
- .Define("EClusterScFi", Ecluster, {"EcalBarrelScFiClusters"})
- .Define("E_over_p_ClusterScFi", E_over_mom, {"EClusterScFi", "mom_smeared"})
- .Define("E_over_p_RecoScFi", E_over_mom, {"ErecScFi", "mom_smeared"});
+ .Define("EoverPClusterScFi", E_over_mom, {"EClusterScFi", "MomSmeared"})
+ .Define("EoverPRecoScFi", E_over_mom, {"ErecScFi", "MomSmeared"});
+
+ auto dpions = d1.Filter(isPion, {"mcparticles2"});
+ auto delectrons = d1.Filter(isElectron, {"mcparticles2"});
+
+ enum typeColumns {Ethr, Mom, MomSmeared, EClusterScFi, NClusterScFi, EoverPScFi, EoverPRecScFi, nColumnsToPlot};
+
+ infoHist1D infoHists[nColumnsToPlot];
+ infoHists[Ethr] = {"hEthr", "Ethr", "Thrown Energy; Thrown Energy [GeV]; Events", 500, 0.0, 11.0};
+ infoHists[Mom] = {"hMom", "Mom", "Thrown Momentum; Thrown Momentum [GeV/c]; Events", 500, 0.0, 11.0};
+ infoHists[MomSmeared] = {"hMomSmeared", "MomSmeared", "Thrown Momentum Smeared; Thrown Momentum Smeared [GeV/c]; Events", 500, 0.0, 11.0};
+ infoHists[EClusterScFi] = {"hEClusterScFi", "EClusterScFi", "Cluster Energy; Cluster Energy [GeV]; Events", 500, 0.0, 11.0};
+ infoHists[NClusterScFi] = {"hNClusterScFi", "NClusterScFi", "Number of Clusters; # of Clusters; Events", 10, 0, 10};
+ infoHists[EoverPScFi] = {"hEoverPScFi", "EoverPClusterScFi", "Cluster Energy/Momentum; Cluster Energy/Momentum; Events", 480, 0., 1.2};
+ infoHists[EoverPRecScFi] = {"hEoverPScFi", "EoverPRecoScFi", "Reco Hits Energy/Momentum; Reco Hits Energy/Momentum; Events", 500, 0.0, 0.25};
+
+ // Define and draw Histograms
+ for(int col=0; col<nColumnsToPlot; col++){
+ auto infoHist = infoHists[col];
+ auto helectron = delectrons.Histo1D({fmt::format("{}_electron",infoHist.name).c_str(), infoHist.titles.c_str(), infoHist.nbins, infoHist.x0, infoHist.x1},infoHist.column.c_str());
+ auto hpion = dpions.Histo1D({fmt::format("{}_pion",infoHist.name).c_str(), infoHist.titles.c_str(), infoHist.nbins, infoHist.x0, infoHist.x1},infoHist.column.c_str());
+
+ auto c = new TCanvas(fmt::format("c{}",infoHist.comlumn).c_str(), fmt::format("c{}",infoHist.comlumn).c_str(), 700, 500);
+ helectron->GetYaxis()->SetTitleOffset(1.4);
+ helectron->SetLineWidth(2);
+ helectron->SetLineColor(kBlue);
+ helectron->DrawClone();
+ hpion->SetLineColor(kRed);
+ hpion->DrawClone("same");
+ save_canvas(c, infoHist.column.c_str(), E_label, "electron_pion");
+ }
- // Define Histograms
-
- auto hEthr_el = delectron.Histo1D({"hEthr_el", "Thrown Energy; Thrown Energy [GeV]; Events", 1000, 0.0, 11.0},"Ethr");
- auto hMom_el = delectron.Histo1D({"hMom_el", "Thrown Momentum; Thrown Momentum [GeV/c]; Events", 1000, 0.0, 11.0},"mom");
- auto hMomSmeared_el = delectron.Histo1D({"hMom_el", "Thrown Momentum Smeared; Thrown Momentum Smeared [GeV/c]; Events", 1000, 0.0, 11.0},"mom_smeared");
-
- // auto hEClusterImg_el = delectron.Histo1D({"hEClusterImg", "Cluster Energy; Cluster Energy [GeV]; Events", 200, 0.0, 25.0},"EClusterImg");
- // auto hNClusterImg_el = delectron.Histo1D({"hNClusterImg", "Number of Clusters; # of Clusters; Events", 100, 0.0, 100.0}, "NClusterImg");
- // auto hEoverPImg_el = delectron.Histo1D({"hfsamImg", "Cluster Energy/E true; Cluster Energy/E true; Events", 100, 0.0, 1.5},"fsamClusterImg");
- // auto hEoverPRecImg_el = delectron.Histo1D({"hfsamRecImg", "Reco Hits Energy/E true; Reco Hits Energy/E true; Events", 50, 0.0, 0.1},"fsamRecImg");
-
- auto hEClusterScFi_el = delectron.Histo1D({"hEClusterScFi_el", "Cluster Energy; Cluster Energy [GeV]; Events", 500, 0.0, 11.0},"EClusterScFi");
- auto hNClusterScFi_el = delectron.Histo1D({"hNClusterScFi_el", "Number of Clusters; # of Clusters; Events", 10, 0, 10}, "NClusterScFi");
- auto hEoverPScFi_el = delectron.Histo1D({"hEoverPScFi_el", "Cluster Energy/E true; Cluster Energy/E true; Events", 480, 0., 1.2},"E_over_p_ClusterScFi");
- auto hEoverPRecScFi_el = delectron.Histo1D({"hEoverPRecScFi_el", "Reco Hits Energy/E true; Reco Hits Energy/E true; Events", 100, 0.0, 0.25},"E_over_p_RecoScFi");
-
- auto hEthr_pion = dpion.Histo1D({"hEthr_pion", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 25.0},"Ethr");
- auto hMom_pion = dpion.Histo1D({"hMom_pion", "Thrown Momentum; Thrown Momentum [GeV/c]; Events", 100, 0.0, 25.0},"mom");
- auto hMomSmeared_pion = dpion.Histo1D({"hMomSmeared_pion", "Thrown Momentum Smeared; Thrown Momentum Smeared [GeV/c]; Events", 100, 0.0, 25.0},"mom_smeared");
-
- // auto hEClusterImg_pion = dpion.Histo1D({"hEClusterImg", "Cluster Energy; Cluster Energy [GeV]; Events", 200, 0.0, 25.0},"EClusterImg");
- // auto hNClusterImg_pion = dpion.Histo1D({"hNClusterImg", "Number of Clusters; # of Clusters; Events", 100, 0.0, 100.0}, "NClusterImg");
- // auto hEoverPImg_pion = dpion.Histo1D({"hfsamImg", "Cluster Energy/E true; Cluster Energy/E true; Events", 100, 0.0, 1.5},"fsamClusterImg");
- // auto hEoverPRecImg_pion = dpion.Histo1D({"hfsamRecImg", "Reco Hits Energy/E true; Reco Hits Energy/E true; Events", 50, 0.0, 0.1},"fsamRecImg");
+ auto nlayers = 20;
+ enum typeLayerColumns {ErecScFi, ErecImg, nhitsScFi, nhitsImg, nLayerColumnsToPlot};
+ std::string cname[nLayerColumnsToPlot] = {"ErecScFi_layer", "ErecImg_layer", "nhitsScFi_layer", "nhitsImg_layer"};
+ TCanvas *clayer[nLayerColumnsToPlot];
+ for(int col=0; col<nLayerColumnsToPlot; col++){
+ clayer[col] = new TCanvas(fmt::format("clayer{}",col).c_str(),fmt::format("clayer{}",col).c_str(),1250, 1000);
+ clayer[col]->Divide(4,5);
+ }
- auto hEClusterScFi_pion = dpion.Histo1D({"hEClusterScFi_pion", "Cluster Energy; Cluster Energy [GeV]; Events", 500, 0.0, 11.0},"EClusterScFi");
- auto hNClusterScFi_pion = dpion.Histo1D({"hNClusterScFi_pion", "Number of Clusters; # of Clusters; Events", 10, 0, 10}, "NClusterScFi");
- auto hEoverPScFi_pion = dpion.Histo1D({"hEoverPScFi_pion", "Cluster Energy/E true; Cluster Energy/E true; Events", 480, 0., 1.2},"E_over_p_ClusterScFi");
- auto hEoverPRecScFi_pion = dpion.Histo1D({"hEoverPRecScFi_pion", "Reco Hits Energy/E true; Reco Hits Energy/E true; Events", 50, 0.0, 0.25},"E_over_p_RecoScFi");
+ for(int layer = 1; layer < nlayers + 1; layer++) {
+ std::cout << "Layer to process: " << layer << std::endl;
+ // Energy reconstructed per layer
+ auto ERecInLayer = [=](const std::vector<dd4pod::CalorimeterHitData>& evt) {
+ auto layer_edep = 0.0;
+ for (const auto& i: evt) {
+ if (i.layer == layer) {
+ layer_edep += i.energyDeposit;
+ }
+ }
+ return layer_edep;
+ };
+
+ // Hits reconstructed per layer
+ auto nhitsInLayer = [=](const std::vector<dd4pod::CalorimeterHitData>& evt) {
+ auto nhits = 0;
+ for (const auto& i: evt) {
+ if (i.layer == layer) {
+ nhits++;
+ }
+ }
+ return nhits;
+ };
+
+ // dE up to the layer
+ auto ERecUpToLayer = [=](const std::vector<eic::CalorimeterHitData>& evt) {
+ double edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer < layer+1){edep += i.energy;}
+ }
+ return edep;
+ };
+
+ auto d2 = d0.Define(fmt::format("ErecImg_layer_{}",layer).c_str(), ERecInLayer, {"RecoEcalBarrelImagingHits"})
+ .Define(fmt::format("ErecScFi_layer_{}",layer).c_str(), ERecInLayer, {"EcalBarrelScFiHitsReco"})
+ .Define(fmt::format("Erec_layer_{}",layer).c_str(), fmt::format("EsimImg_layer_{}+EsimScFi_layer_{}",layer,layer).c_str())
+ .Define(fmt::format("nhitsImg_layer_{}",layer).c_str(), nhitsInLayer, {"RecoEcalBarrelImagingHits"})
+ .Define(fmt::format("nhitsScFi_layer_{}",layer).c_str(), nhitsInLayer, {"EcalBarrelScFiHitsReco"})
+ .Define(fmt::format("ErecSumImg_layer_{}",layer).c_str(), ERecUpToLayer, {"RecoEcalBarrelImagingHits"})
+ .Define(fmt::format("ErecSumScFi_layer_{}",layer).c_str(), ERecUpToLayer, {"EcalBarrelScFiHitsReco"});
+
+ auto dpionslayer = d2.Filter(isPion,{"mcparticles2"});
+ auto delectronslayer = d2.Filter(isElectron,{"mcparticles2"});
+
+ infoHist1D infoLayerHists[nLayerColumnsToPlot];
+
+ infoLayerHists[ErecImg] = {fmt::format("hErecImg_layer_{}",layer), fmt::format("ErecImg_layer_{}",layer), fmt::format("Energy Deposit in layer {}; Energy Deposit [Gev]; Events",layer), 200, 0.0, 0.04};
+ infoLayerHists[ErecScFi] = {fmt::format("hErecScFi_layer_{}",layer), fmt::format("ErecScFi_layer_{}",layer), fmt::format("Energy Deposit in layer {}; Energy Deposit [Gev]; Events",layer), 200, 0.0, 0.1};
+ infoLayerHists[nhitsImg] = {fmt::format("hnhitsImg_layer_{}_el",layer), fmt::format("nhitsImg_layer_{}",layer), fmt::format("Number of Reconstructed Hits in layer (); Number of hits; Events",layer), 500, 0, 500};
+ infoLayerHists[nhitsScFi] = {fmt::format("hnhitsScFi_layer_{}_el",layer), fmt::format("nhitsScFi_layer_{}",layer), fmt::format("Number of Reconstructed Hits in layer (); Number of hits; Events",layer), 500, 0, 500};
+
+ for(int col=0; col<nLayerColumnsToPlot; col++){
+ auto infoHist = infoLayerHists[col];
+ auto helectron = delectronslayer.Histo1D({fmt::format("{}_electron",infoHist.name).c_str(), infoHist.titles.c_str(), infoHist.nbins, infoHist.x0, infoHist.x1},infoHist.column.c_str());
+ auto hpion = dpionslayer.Histo1D({fmt::format("{}_pion",infoHist.name).c_str(), infoHist.titles.c_str(), infoHist.nbins, infoHist.x0, infoHist.x1},infoHist.column.c_str());
+
+ clayer[col]->cd(layer);
+ gPad->SetLogy();
+ helectron->GetYaxis()->SetTitleOffset(1.4);
+ helectron->SetLineWidth(2);
+ helectron->SetLineColor(kBlue);
+ helectron->DrawClone();
+ hpion->SetLineColor(kRed);
+ hpion->DrawClone("same");
+
+ // auto h = hEsim_layer->DrawCopy();
+
+ // auto up_range = h->GetMean() + 3*h->GetStdDev();
+ // auto down_range = h->GetMean() - 3*h->GetStdDev();
+ // h->SetLineWidth(2);
+ // h->SetLineColor(kBlue);
+
+ // h->GetXaxis()->SetRange(h->GetXaxis()->GetBinUpEdge(1), up_range); // skip 0th bin
+ // auto mean_layer = h->GetMean();
+ // auto rms_layer = h->GetStdDev();
+ // h->GetXaxis()->SetRange(); // reset the range
+ // h->GetXaxis()->SetRangeUser(0.,up_range);
+
+ // auto no_edep = (h->GetBinContent(1)/h->GetEntries())*100;
+ // gr_no_edep.SetPoint(gr_no_edep.GetN(),layer,no_edep);
+ // gr_edep_mean.SetPoint(gr_edep_mean.GetN(),layer, mean_layer);
+ // gr_edep_mean.SetPointError(gr_edep_mean.GetN()-1,0, rms_layer);
+ }
+ for(int col=0; col<nLayerColumnsToPlot; col++) save_canvas(clayer[col], cname[col] , E_label, "electron_pion");
+
+
// Event Counts
auto nevents_thrown = delectron.Count();
std::cout << "Number of Thrown Electron Events: " << (*nevents_thrown) << "\n";
- // Draw Histograms
- {
- TCanvas* c1 = new TCanvas("c1", "c1", 700, 500);
- c1->Divide(2,2);
- c1->SetLogy(1);
- c1->cd(1);
- auto h = hEthr_el->DrawCopy();
- h->GetYaxis()->SetTitleOffset(1.4);
- h->SetLineWidth(2);
- h->SetLineColor(kBlue);
- auto h_pion = hEthr_pion->DrawCopy("same");
- h_pion->SetLineColor(kRed);
- c1->cd(2);
- hMom_el->GetYaxis()->SetTitleOffset(1.4);
- hMom_el->SetLineWidth(2);
- hMom_el->SetLineColor(kBlue);
- hMom_el->DrawClone();
- auto up_fit = hMom_el->GetMean() + 5*hMom_el->GetStdDev();
- auto down_fit = hMom_el->GetMean() - 5*hMom_el->GetStdDev();
- hMom_el->GetXaxis()->SetRangeUser(down_fit,up_fit);
- hMom_pion->SetLineColor(kRed);
- hMom_pion->DrawClone("same");
- c1->cd(3);
- hMomSmeared_el->GetYaxis()->SetTitleOffset(1.4);
- hMomSmeared_el->SetLineWidth(2);
- hMomSmeared_el->SetLineColor(kBlue);
- hMomSmeared_el->DrawClone();
- hMomSmeared_el->GetXaxis()->SetRangeUser(down_fit,up_fit);
- hMomSmeared_pion->SetLineColor(kRed);
- hMomSmeared_pion->DrawClone("same");
- save_canvas(c1, "gen", E_label, "electron_pion");
- }
-
- TCanvas* c2 = new TCanvas("c2", "c2", 1400, 1000);
- hNClusterScFi_el->GetYaxis()->SetTitleOffset(1.4);
- hNClusterScFi_el->SetLineWidth(2);
- hNClusterScFi_el->SetLineColor(kBlue);
- hNClusterScFi_el->DrawClone();
- hNClusterScFi_pion->SetLineColor(kRed);
- hNClusterScFi_pion->DrawClone("same");
- save_canvas(c2, "NCluster", E_label, "electron_pion");
-
- TCanvas* c3 = new TCanvas("c3", "c3", 1400, 1000);
- hEoverPScFi_el->GetYaxis()->SetTitleOffset(1.4);
- hEoverPScFi_el->SetLineWidth(2);
- hEoverPScFi_el->SetLineColor(kBlue);
- hEoverPScFi_el->DrawClone();
- hEoverPScFi_pion->SetLineColor(kRed);
- hEoverPScFi_pion->DrawClone("same");
- save_canvas(c3, "EoverPClusterScFi", E_label, "electron_pion");
-
- TCanvas* c4 = new TCanvas("c4", "c4", 1400, 1000);
- hEoverPRecScFi_el->GetYaxis()->SetTitleOffset(1.4);
- hEoverPRecScFi_el->SetLineWidth(2);
- hEoverPRecScFi_el->SetLineColor(kBlue);
- hEoverPRecScFi_el->DrawClone();
- hEoverPRecScFi_pion->SetLineColor(kRed);
- hEoverPRecScFi_pion->DrawClone("same");
- save_canvas(c4, "EoverPRecScFi", E_label, "electron_pion");
-
-
return std::make_tuple(0., 0.);
}