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
+ 49
− 21
Compare changes
  • Side-by-side
  • Inline
benchmarks/barrel_ecal/scripts/emcal_barrel_pions_analysis.cxx
+ 49
− 21
@@ -41,7 +41,12 @@ 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 the barrel region(-1 < eta < 1)
// Taken from : Initial considerations for EMCal of the EIC detector by A. Bazilevsky
// 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)},
@@ -56,7 +61,7 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
ROOT::RDataFrame d0("events", input_fname);
// Sampling Fraction
double samp_frac = 1.0;
double samp_frac = 0.0136;
// Thrown Energy [GeV]
auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
@@ -116,21 +121,22 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
return input[2].pdgID;
};
// Returns the pdgID of the particle daughters(hopefully?)
auto getdau = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
std::vector<int> result;
//for (auto part : input[2].daughters_begin)
// result.push_back(part)
return input[2].daughters_begin;
};
// 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"})
;
// Set sampling Fraction, ideally this will be taken from a json file
auto svec = d1.Mean<double>("fsam");
for (auto ele : svec){cout << ele << endl;}
auto d2 = d1.Define("dE", eResol, {"Esim","Ethr"})
.Define("dE_rel", eResol_rel, {"Esim","Ethr"})
.Define("dau", getdau, {"mcparticles"})
;
// Define Histograms
@@ -139,6 +145,18 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
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 hpid2 = d1.Histo1D({"hpid2", "PID; PID; Count", 100, -220, 220}, "dau");
// 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();
@@ -193,18 +211,28 @@ 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();
TCanvas *c6 = new TCanvas("c6", "c6", 700, 500);
c5->SetLogy(1);
hpid2->GetYaxis()->SetTitleOffset(1.4);
hpid2->SetLineWidth(2);
hpid2->SetLineColor(kBlue);
hpid2->DrawClone();
c6->SaveAs("results/emcal_barrel_pions_pid2.png");
c6->SaveAs("results/emcal_barrel_pions_pid2.pdf");
//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]);
@@ -214,6 +242,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");
@@ -222,8 +252,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");