diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons.cxx
index 12a78da363d28a2725f022685a8306aa8532ec2f..217b2f55af9e452112d8b720d6c47993c4dc89cf 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 22cb118ce62847bd6ebef8e313537edbd04e3f80..0000000000000000000000000000000000000000
--- 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 75c2dae9c9218f036bbca7dd65f65e85bd67eca2..90622139975845b45f359fada3d2debf79145254 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 04fc024fdb8eb4652c8b8011d9c7d41871f3287b..61041363041ba8f110bb09f22eef3cf987db00ec 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 4d52a6e5e8ddd47e0c24747591d529bc50070b3c..0e35ac33bc6c380d0ebc69df2dd6c113e0f337a5 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 5034c1b2791a19f3b759a0974ff9e59b19a01979..b60a0eed282d40c74131268ce2ae1ea6d5e4983c 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 a171175e5a7dea65f792ba6362401a0ae3075d37..343eb26cc918d9c178a9be9f807ab6892fc3821a 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 3bf036cde29b3900cee9465c9e51b2cd178046ae..e38a5c207415ea7f1e631c06616158016bb2c21c 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"})