Skip to content
Snippets Groups Projects

Resolve "pi0 resolution for ECal barrel"

Merged Resolve "pi0 resolution for ECal barrel"
19-pi0-resolution-for-ecal-barrel into master
Merged Marshall Scott requested to merge 19-pi0-resolution-for-ecal-barrel into master
Compare
and
Show latest version
3 files
+ 54
− 48
Compare changes
  • Side-by-side
  • Inline
Files
3
benchmarks/barrel_ecal/scripts/emcal_barrel_pions_analysis.cxx
+ 33
− 48
@@ -9,9 +9,9 @@
#include "dd4pod/Geant4ParticleCollection.h"
#include "dd4pod/CalorimeterHitCollection.h"
//#include <detector_benchmarks/include/benchmark.h>
//#include "mt.h"
//#include "util.h"
#include "benchmark.h"
#include "mt.h"
#include "util.h"
#include "TCanvas.h"
#include "TStyle.h"
@@ -19,9 +19,7 @@
#include "TH1.h"
#include "TF1.h"
#include "TH1D.h"
#include <TFitResult.h>
#include "TSystem.h"
#include <fstream>
#include "TFitResult.h"
using ROOT::RDataFrame;
using namespace ROOT::VecOps;
@@ -44,15 +42,15 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
//TODO: Change test_tag to something else
std:string detector = "Barrel_emcal";
double resolutionTarget = 0.1;
/*
eic::util::Test pion0_Energy_resolution{
eic::util::Test pi0_energy_resolution{
{{"name", fmt::format("{}_energy_resolution", test_tag)},
{"title", "Pion0 Energy resolution"},
{"description",
fmt::format("Pion0 energy resolution with {}, estimated using a Gaussian fit.", detector)},
{"quantity", "resolution (in %)"},
{"target", std::to_string(resolutionTarget)}}};
*/
ROOT::EnableImplicitMT();
ROOT::RDataFrame d0("events", input_fname);
@@ -87,16 +85,6 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
// Sampling fraction = Esampling / Ethrown
auto fsam = [](const std::vector<double>& sampled, const std::vector<double>& thrown) {
std::vector<double> result;
for (const auto& E1 : thrown) {
for (const auto& E2 : sampled)
result.push_back(E2/E1);
}
return result;
};
// Sampling fraction = Esampling / Ethrown, revised
auto fsam2 = [](const std::vector<double>& sampled, const std::vector<double>& thrown) {
std::vector<double> result;
auto it_sam = sampled.cbegin();
auto it_thr = thrown.cbegin();
@@ -117,33 +105,30 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
return result;
};
// Relative Energy Resolution = (Esampling - Ethrown)/Esampling
// Relative Energy Resolution = (Esampling - Ethrown)/Ethrown
auto eResol_rel = [](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) / *it_sam);
result.push_back((*it_sam - *it_thr) / *it_thr);
}
return result;
};
// Thrown Energy [GeV]
auto charge = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
std::vector<double> result;
result.push_back(input[2].charge);
return result;
// Relative Energy Resolution = (Esampling - Ethrown)/Ethrown
auto getpid = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
return input[2].pdgID;
};
// Define variables
auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
.Define("nhits", nhits, {"EcalBarrelHits"})
.Define("Esim", Esim, {"EcalBarrelHits"})
//.Define("fsam", fsam, {"Esim","Ethr"})
.Define("fsam", fsam2, {"Esim","Ethr"})
.Define("dE", eResol, {"Esim", "Ethr"})
.Define("dE_rel", eResol_rel, {"Esim", "Ethr"})
.Define("charge", charge, {"mcparticles"})
.Define("fsam", fsam, {"Esim","Ethr"})
.Define("dE", eResol, {"Esim","Ethr"})
.Define("dE_rel", eResol_rel, {"Esim","Ethr"})
.Define("pid", getpid, {"mcparticles"})
;
// Define Histograms
@@ -151,7 +136,7 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
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");
auto hcharge = d1.Histo1D({"hcharge", "charge; charge; Count", 10, -2, 2}, "charge");
auto hpid = d1.Histo1D({"hpid", "PID; PID; Count", 100, -220, 220}, "pid");
// Event Counts
auto nevents_thrown = d1.Count();
@@ -198,27 +183,27 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
c4->SaveAs("results/emcal_barrel_pions_fsam.pdf");
TCanvas *c5 = new TCanvas("c5", "c5", 700, 500);
hcharge->GetYaxis()->SetTitleOffset(1.4);
hcharge->SetLineWidth(2);
hcharge->SetLineColor(kBlue);
hcharge->DrawClone();
c5->SaveAs("results/emcal_barrel_pions_charge.png");
c5->SaveAs("results/emcal_barrel_pions_chage.pdf");
c5->SetLogy(1);
hpid->GetYaxis()->SetTitleOffset(1.4);
hpid->SetLineWidth(2);
hpid->SetLineColor(kBlue);
hpid->DrawClone();
c5->SaveAs("results/emcal_barrel_pions_pid.png");
c5->SaveAs("results/emcal_barrel_pions_pid.pdf");
//Energy Resolution Work
auto hdE = d1.Histo1D({"hdE", "dE; dE[GeV]; Events", 20, -7.5, 7.5}, "dE");//changed from 100
auto hdE_rel = d1.Histo1D({"hdE_rel", "dE Relative; dE Relative; Events", 20, -7.5, 7.5}, "dE_rel");//changed from 100
//hdE_rel->Fit("gaus", "", "", -7.5, 7.5);
//double* res = hdE_rel->GetFunction("gaus")->GetParameters();
auto hdE_rel = d1.Histo1D({"hdE_rel", "dE Relative; dE Relative; Events", 20, -2, 2}, "dE_rel");//changed from 100
hdE_rel->Fit("gaus", "", "", -2, 2);
double* res = hdE_rel->GetFunction("gaus")->GetParameters();
//Pass/Fail
/*
if (res[2] <= resolutionTarget) {
pion0_energy_resolution.pass(res[2]);
pi0_energy_resolution.pass(res[2]);
} else {
pion0_energy_resolution.fail(res[2]);
pi0_energy_resolution.fail(res[2]);
}
*/
//std::printf("Energy Resolution is %f\n", res[2]);
//Energy Resolution Histogram Plotting
@@ -235,13 +220,13 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
hdE_rel->GetYaxis()->SetTitleOffset(1.4);
hdE_rel->SetLineWidth(2);
hdE_rel->SetLineColor(kBlue);
//hdE_rel->GetFunction("gaus")->SetLineWidth(2);
//hdE_rel->GetFunction("gaus")->SetLineColor(kRed);
hdE_rel->GetFunction("gaus")->SetLineWidth(2);
hdE_rel->GetFunction("gaus")->SetLineColor(kRed);
hdE_rel->DrawClone();
cdE_rel->SaveAs("results/emcal_barrel_pi0_dE_rel.png");
cdE_rel->SaveAs("results/emcal_barrel_pi0_dE_rel.pdf");
//eic::util::write_test({pion0_energy_resolution}, fmt::format("{}_pions.json", detector));
eic::util::write_test({pi0_energy_resolution}, fmt::format("{}_pions.json", detector));
}