From a38676b568ae8a0a297faaa76432587cfbfea017 Mon Sep 17 00:00:00 2001
From: Maria Zurek <zurek@anl.gov>
Date: Wed, 14 Jul 2021 18:26:34 +0000
Subject: [PATCH] Resolve "Add information from SciFi layers to benchmarks"
---
.../scripts/emcal_barrel_electrons.cxx | 1 +
.../emcal_barrel_electrons_analysis.cxx | 153 ------------------
.../scripts/emcal_barrel_electrons_reader.cxx | 1 +
.../emcal_barrel_energy_scan_analysis.cxx | 12 +-
.../emcal_barrel_particles_analysis.cxx | 4 +-
.../scripts/emcal_barrel_pi0_analysis.cxx | 4 +-
.../scripts/emcal_barrel_pions_analysis.cxx | 53 ++----
.../emcal_barrel_pions_electrons_analysis.cxx | 4 +-
8 files changed, 36 insertions(+), 196 deletions(-)
delete mode 100644 benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_analysis.cxx
diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons.cxx
index 12a78da3..217b2f55 100644
--- a/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons.cxx
+++ b/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons.cxx
@@ -78,3 +78,4 @@ void emcal_barrel_electrons(int n_events = 1e6, double e_start = 0.0, double e_e
hepmc_output.close();
std::cout << "Events parsed and written: " << events_parsed << std::endl;
}
+
diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_analysis.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_analysis.cxx
deleted file mode 100644
index 22cb118c..00000000
--- a/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_analysis.cxx
+++ /dev/null
@@ -1,153 +0,0 @@
-////////////////////////////////////////
-// 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"
-
-#include "emcal_barrel_common_functions.h"
-#include <cctype>
-
-using ROOT::RDataFrame;
-using namespace ROOT::VecOps;
-
-void emcal_barrel_electrons_analysis(const char* input_fname = "sim_output/sim_emcal_barrel_uniform_electrons.root")
-{
- input_fname = "temp_pions_electrons.root";
- // 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);
-
- ROOT::EnableImplicitMT();
- ROOT::RDataFrame d0("events", input_fname);
-
- // Environment Variables
- std::string detector_path = "";
- std::string detector_name = "topside";
- if(std::getenv("DETECTOR_PATH")) {
- detector_path = std::getenv("DETECTOR_PATH");
- }
- if(std::getenv("JUGGLER_DETECTOR")) {
- detector_name = std::getenv("JUGGLER_DETECTOR");
- }
-
- // 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");
-
- addDetectorName(detector_name, hEthr.GetPtr());
- addDetectorName(detector_name, hEsim.GetPtr());
- addDetectorName(detector_name, hfsam.GetPtr());
-
- // 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);
- c1->SaveAs("results/emcal_barrel_electrons_Ethr.png");
- c1->SaveAs("results/emcal_barrel_electrons_Ethr.pdf");
- }
- 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);
- c2->SaveAs("results/emcal_barrel_electrons_nhits.png");
- c2->SaveAs("results/emcal_barrel_electrons_nhits.pdf");
- }
-
- {
- 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);
- c3->SaveAs("results/emcal_barrel_electrons_Esim.png");
- c3->SaveAs("results/emcal_barrel_electrons_Esim.pdf");
- }
-
- std::cout << "derp4\n";
- {
- 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);
- //h->GetFunction("gaus")->SetLineWidth(2);
- //h->GetFunction("gaus")->SetLineColor(kRed);
- c4->SaveAs("results/emcal_barrel_electrons_fsam.png");
- c4->SaveAs("results/emcal_barrel_electrons_fsam.pdf");
- }
-}
diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_reader.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_reader.cxx
index 75c2dae9..90622139 100644
--- a/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_reader.cxx
+++ b/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_reader.cxx
@@ -123,3 +123,4 @@ void emcal_barrel_electrons_reader(double e_start = 0.0, double e_end = 30.0, co
c6->SaveAs("results/input_emcal_barrel_electrons_p.png");
c6->SaveAs("results/input_emcal_barrel_electrons_p.pdf");
}
+
diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_energy_scan_analysis.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_energy_scan_analysis.cxx
index 04fc024f..61041363 100644
--- a/benchmarks/barrel_ecal/scripts/emcal_barrel_energy_scan_analysis.cxx
+++ b/benchmarks/barrel_ecal/scripts/emcal_barrel_energy_scan_analysis.cxx
@@ -86,7 +86,9 @@ std::tuple <double, double, double, double> extract_sampling_fraction_parameters
// Define variables
auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
.Define("nhits", nhits, {"EcalBarrelHits"})
- .Define("Esim", Esim, {"EcalBarrelHits"})
+ .Define("EsimImg", Esim, {"EcalBarrelHits"})
+ .Define("EsimScFi", Esim, {"EcalBarrelScFiHits"})
+ .Define("Esim", "EsimImg+EsimScFi")
.Define("fsam", fsam, {"Esim", "Ethr"});
// Define Histograms
@@ -105,7 +107,7 @@ std::tuple <double, double, double, double> extract_sampling_fraction_parameters
"fsam");
// Number of layers to read the edep from
- int nlayers = 20;
+ int nlayers = 7;
TGraphErrors gr_no_edep(nlayers);
TGraphErrors gr_edep_mean(nlayers);
@@ -125,10 +127,12 @@ std::tuple <double, double, double, double> extract_sampling_fraction_parameters
return layer_edep;
};
- auto d2 = d0.Define(fmt::format("Esim_layer_{}",layer).c_str(), Esim_layer, {"EcalBarrelHits"});
+ auto d2 = d0.Define(fmt::format("EsimImg_layer_{}",layer).c_str(), Esim_layer, {"EcalBarrelHits"})
+ .Define(fmt::format("EsimScFi_layer_{}",layer).c_str(), Esim_layer, {"EcalBarrelScFiHits"})
+ .Define(fmt::format("Esim_layer_{}",layer).c_str(), fmt::format("EsimImg_layer_{}+EsimScFi_layer_{}",layer,layer).c_str());
std::cout << "Layer to process: " << layer << std::endl;
- 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(), 800, 0.0, 0.04}, fmt::format("Esim_layer_{}",layer));
+ 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));
clayer->cd(layer);
gPad->SetLogy();
diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_particles_analysis.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_particles_analysis.cxx
index 4d52a6e5..0e35ac33 100644
--- a/benchmarks/barrel_ecal/scripts/emcal_barrel_particles_analysis.cxx
+++ b/benchmarks/barrel_ecal/scripts/emcal_barrel_particles_analysis.cxx
@@ -94,7 +94,9 @@ void emcal_barrel_particles_analysis(std::string particle_name = "electron", boo
// Define variables
auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
.Define("nhits", nhits, {"EcalBarrelHits"})
- .Define("Esim", Esim, {"EcalBarrelHits"})
+ .Define("EsimImg", Esim, {"EcalBarrelHits"})
+ .Define("EsimScFi", Esim, {"EcalBarrelScFiHits"})
+ .Define("Esim", "EsimImg+EsimScFi")
.Define("fsam", fsam, {"Esim", "Ethr"});
// Define Histograms
diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_analysis.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_analysis.cxx
index 5034c1b2..b60a0eed 100644
--- a/benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_analysis.cxx
+++ b/benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_analysis.cxx
@@ -99,7 +99,9 @@ void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_b
// Define variables
auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
.Define("nhits", nhits, {"EcalBarrelHits"})
- .Define("Esim", Esim, {"EcalBarrelHits"})
+ .Define("EsimImg", Esim, {"EcalBarrelHits"})
+ .Define("EsimScFi", Esim, {"EcalBarrelScFiHits"})
+ .Define("Esim", "EsimImg+EsimScFi")
.Define("fsam", fsam, {"Esim","Ethr"})
.Define("pid", getpid, {"mcparticles"})
.Define("dau", getdau, {"mcparticles"})
diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_analysis.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_analysis.cxx
index a171175e..343eb26c 100644
--- a/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_analysis.cxx
+++ b/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_analysis.cxx
@@ -45,9 +45,9 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
// Thrown Energy [GeV]
auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
- std::vector<double> result;
- result.push_back(TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass));
- return result;
+ 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
@@ -55,45 +55,25 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
// Energy deposition [GeV]
auto Esim = [](const std::vector<dd4pod::CalorimeterHitData>& evt) {
- std::vector<double> result;
auto total_edep = 0.0;
for (const auto& i: evt)
total_edep += i.energyDeposit;
- result.push_back(total_edep);
- return result;
+ return total_edep;
};
// Sampling fraction = Esampling / Ethrown
- auto fsam = [](const std::vector<double>& sampled, const std::vector<double>& thrown) {
- std::vector<double> result;
- auto it_sam = sampled.cbegin();
- auto it_thr = thrown.cbegin();
- for (; it_sam != sampled.end() && it_thr != thrown.end(); ++it_sam, ++it_thr) {
- result.push_back(*it_sam / *it_thr);
- }
- return result;
+ auto fsam = [](const double sampled, const double thrown) {
+ return sampled / thrown;
};
// Energy Resolution = Esampling/Sampling_fraction - Ethrown
- auto eResol = [samp_frac](const std::vector<double>& sampled, const std::vector<double>& thrown) {
- std::vector<double> result;
- auto it_sam = sampled.cbegin();
- auto it_thr = thrown.cbegin();
- for (; it_sam != sampled.end() && it_thr != thrown.end(); ++it_sam, ++it_thr) {
- result.push_back(*it_sam / samp_frac - *it_thr);
- }
- return result;
+ auto eResol = [samp_frac](double sampled, double thrown) {
+ return sampled / samp_frac - thrown;
};
// Relative Energy Resolution = (Esampling/Sampling fraction - Ethrown)/Ethrown
- auto eResol_rel = [samp_frac](const std::vector<double>& sampled, const std::vector<double>& thrown) {
- std::vector<double> result;
- auto it_sam = sampled.cbegin();
- auto it_thr = thrown.cbegin();
- for (; it_sam != sampled.end() && it_thr != thrown.end(); ++it_sam, ++it_thr) {
- result.push_back((*it_sam / samp_frac - *it_thr) / *it_thr);
- }
- return result;
+ auto eResol_rel = [samp_frac](double sampled, double thrown) {
+ return (sampled / samp_frac - thrown) / thrown;
};
// Returns the pdgID of the particle
@@ -107,12 +87,13 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
};
// Define variables
- auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
- .Define("nhits", nhits, {"EcalBarrelHits"})
- .Define("Esim", Esim, {"EcalBarrelHits"})
- .Define("fsam", fsam, {"Esim","Ethr"})
- .Define("pid", getpid, {"mcparticles"})
- ;
+ auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
+ .Define("nhits", nhits, {"EcalBarrelHits"})
+ .Define("EsimImg", Esim, {"EcalBarrelHits"})
+ .Define("EsimScFi", Esim, {"EcalBarrelScFiHits"})
+ .Define("Esim", "EsimImg+EsimScFi")
+ .Define("fsam", fsam, {"Esim", "Ethr"})
+ .Define("pid", getpid, {"mcparticles"});
// Define Histograms
auto hEthr = d1.Histo1D({"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 7.5}, "Ethr");
diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_analysis.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_analysis.cxx
index 3bf036cd..e38a5c20 100644
--- a/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_analysis.cxx
+++ b/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_analysis.cxx
@@ -140,7 +140,9 @@ void emcal_barrel_pions_electrons_analysis(const char* input_fname1 = "sim_outpu
auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
.Define("Pthr", Pthr, {"mcparticles"})
.Define("nhits", nhits, {"EcalBarrelHits"})
- .Define("Esim", Esim, {"EcalBarrelHits"})
+ .Define("EsimImg", Esim, {"EcalBarrelHits"})
+ .Define("EsimScFi", Esim, {"EcalBarrelScFiHits"})
+ .Define("Esim", "EsimImg+EsimScFi")
.Define("fsam", fsam, {"Esim","Ethr"})
.Define("pid", getpid, {"mcparticles"})
.Define("dau", getdau, {"mcparticles"})
--
GitLab