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
1 file
+ 43
− 34
Compare changes
  • Side-by-side
  • Inline
benchmarks/barrel_ecal/scripts/emcal_barrel_pions_analysis.cxx
+ 43
− 34
@@ -41,7 +41,11 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
std::string test_tag = "Barrel_emcal_pi0";
//TODO: Change test_tag to something else
std:string detector = "Barrel_emcal";
double resolutionTarget = 0.1;
//Energy resolution in teh barrel region(-1 < eta < 1)
// sigma_E / E = 12% / E^0.5 convoluted with 2%
// sigma_E / E = [ (0.12/E^0.5)^2 + 0.02^2]^0.5, with E in [GeV]
double thrown_energy = 5; // Current thrown energy, will need to grab from json file
double resolutionTarget = TMath::Sqrt(0.12 * 0.12 / thrown_energy + 0.02 * 0.02);
eic::util::Test pi0_energy_resolution{
{{"name", fmt::format("{}_energy_resolution", test_tag)},
@@ -54,14 +58,9 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
ROOT::EnableImplicitMT();
ROOT::RDataFrame d0("events", input_fname);
/*
auto colNames = d0.GetColumnNames();
ofstream out;
out.open("results/column_list.txt");
for (auto &&col : colNames){out << col << endl; cout << col << endl;}
out.close();
*/
// Sampling Fraction
double samp_frac = 0.0136;
// Thrown Energy [GeV]
auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
@@ -94,49 +93,58 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
return result;
};
// Energy Resolution = Esampling - Ethrown
auto eResol = [](const std::vector<double>& sampled, const std::vector<double>& 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 - *it_thr);
result.push_back(*it_sam / samp_frac - *it_thr);
}
return result;
};
// Relative Energy Resolution = (Esampling - Ethrown)/Ethrown
auto eResol_rel = [](const std::vector<double>& sampled, const std::vector<double>& 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 - *it_thr) / *it_thr);
result.push_back((*it_sam / samp_frac - *it_thr) / *it_thr);
}
return result;
};
// Relative Energy Resolution = (Esampling - Ethrown)/Ethrown
// Returns the pdgID of the particle
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("dE", eResol, {"Esim","Ethr"})
.Define("dE_rel", eResol_rel, {"Esim","Ethr"})
.Define("pid", getpid, {"mcparticles"})
;
auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
.Define("nhits", nhits, {"EcalBarrelHits"})
.Define("Esim", Esim, {"EcalBarrelHits"})
.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");
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 hpid = d1.Histo1D({"hpid", "PID; PID; Count", 100, -220, 220}, "pid");
auto hpid = d1.Histo1D({"hpid", "PID; PID; Count", 100, -220, 220}, "pid");
// Set sampling Fraction, ideally this will be taken from a json file
cout << "The sampling fraction mean is " << samp_frac << endl;
samp_frac = hfsam -> GetMean();
cout << "The sampling fraction mean is " << samp_frac << endl;
auto d2 = d1.Define("dE", eResol, {"Esim","Ethr"})
.Define("dE_rel", eResol_rel, {"Esim","Ethr"})
;
// Event Counts
auto nevents_thrown = d1.Count();
@@ -191,18 +199,19 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
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, -2, 2}, "dE_rel");//changed from 100
hdE_rel->Fit("gaus", "", "", -2, 2);
double* res = hdE_rel->GetFunction("gaus")->GetParameters();
//Energy Resolution Calculation
auto hdE = d2.Histo1D({"hdE", "dE; dE[GeV]; Events", 20, -2.5, 2.5}, "dE");//changed from 100
auto hdE_rel = d2.Histo1D({"hdE_rel", "dE Relative; dE Relative; Events", 20, -2.5, 2.5}, "dE_rel");//changed from 100
hdE->Fit("gaus", "", "", -2.5, 2.5);
double* res = hdE->GetFunction("gaus")->GetParameters();
double sigmaOverE = res[2] / thrown_energy;
//Pass/Fail
if (res[2] <= resolutionTarget) {
pi0_energy_resolution.pass(res[2]);
if (sigmaOverE <= resolutionTarget) {
pi0_energy_resolution.pass(sigmaOverE);
} else {
pi0_energy_resolution.fail(res[2]);
pi0_energy_resolution.fail(sigmaOverE);
}
//std::printf("Energy Resolution is %f\n", res[2]);
@@ -212,6 +221,8 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
hdE->GetYaxis()->SetTitleOffset(1.4);
hdE->SetLineWidth(2);
hdE->SetLineColor(kBlue);
hdE->GetFunction("gaus")->SetLineWidth(2);
hdE->GetFunction("gaus")->SetLineColor(kRed);
hdE->DrawClone();
cdE->SaveAs("results/emcal_barrel_pi0_dE.png");
cdE->SaveAs("results/emcal_barrel_pi0_dE.pdf");
@@ -220,8 +231,6 @@ 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->DrawClone();
cdE_rel->SaveAs("results/emcal_barrel_pi0_dE_rel.png");
cdE_rel->SaveAs("results/emcal_barrel_pi0_dE_rel.pdf");