Skip to content
Snippets Groups Projects
Commit 27542e8c authored by Marshall Scott's avatar Marshall Scott Committed by Maria Zurek
Browse files

Resolve "Broken pi0 ecal barrel analysis script"

parent 3a27c116
No related branches found
No related tags found
1 merge request!78Resolve "Broken pi0 ecal barrel analysis script"
Hide whitespace changes
Inline Side-by-side
benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_analysis.cxx
+ 103
− 137
View file @ 27542e8c
...@@ -29,9 +29,11 @@ using ROOT::RDataFrame; ...@@ -29,9 +29,11 @@ using ROOT::RDataFrame;
using namespace ROOT::VecOps; using namespace ROOT::VecOps;
using json = nlohmann::json; using json = nlohmann::json;
void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_barrel_pi0.root") void emcal_barrel_pi0_analysis(
const char* input_fname = "sim_output/sim_emcal_barrel_pi0.root"
//const char* input_fname = "../sim_output/sim_emcal_barrel_uniform_pi0.root"
)
{ {
//input_fname = "../sim_output/sim_emcal_barrel_pi0.root";
// Setting for graphs // Setting for graphs
gROOT->SetStyle("Plain"); gROOT->SetStyle("Plain");
gStyle->SetOptFit(1); gStyle->SetOptFit(1);
...@@ -54,7 +56,6 @@ void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_b ...@@ -54,7 +56,6 @@ void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_b
double samp_frac = j["electron"]["sampling_fraction"]; double samp_frac = j["electron"]["sampling_fraction"];
// Thrown Energy [GeV] // Thrown Energy [GeV]
//double meanE = 5; // Calculated later
auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) { auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
return TMath::Sqrt(input[2].ps.x*input[2].ps.x + input[2].ps.y*input[2].ps.y + input[2].ps.z*input[2].ps.z + input[2].mass*input[2].mass); return TMath::Sqrt(input[2].ps.x*input[2].ps.x + input[2].ps.y*input[2].ps.y + input[2].ps.z*input[2].ps.z + input[2].mass*input[2].mass);
}; };
...@@ -97,146 +98,111 @@ void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_b ...@@ -97,146 +98,111 @@ void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_b
}; };
// Define variables // Define variables
auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"}) auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
.Define("nhits", nhits, {"EcalBarrelHits"}) .Define("nhits", nhits, {"EcalBarrelHits"})
.Define("EsimImg", Esim, {"EcalBarrelHits"}) .Define("EsimImg", Esim, {"EcalBarrelHits"})
.Define("EsimScFi", Esim, {"EcalBarrelScFiHits"}) .Define("EsimScFi", Esim, {"EcalBarrelScFiHits"})
.Define("Esim", "EsimImg+EsimScFi") .Define("Esim", "EsimImg+EsimScFi")
.Define("fsam", fsam, {"Esim","Ethr"}) .Define("fsam", fsam, {"Esim","Ethr"})
.Define("pid", getpid, {"mcparticles"}) .Define("pid", getpid, {"mcparticles"})
.Define("dau", getdau, {"mcparticles"}) .Define("dau", getdau, {"mcparticles"})
.Define("dE", eResol, {"Esim","Ethr"}) .Define("dE", eResol, {"Esim","Ethr"})
.Define("dE_rel", eResol_rel, {"Esim","Ethr"}) .Define("dE_rel", eResol_rel, {"Esim","Ethr"})
; ;
// Define Histograms // Define Histograms
auto hEthr = d1.Histo1D({"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 7.5}, "Ethr"); std::vector <std::string> titleStr = {
auto hNhits = d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events", 100, 0.0, 2000.0}, "nhits"); "Thrown Energy; Thrown Energy [GeV]; Events",
auto hEsim = d1.Histo1D({"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 100, 0.0, 2.0}, "Esim"); "Number of hits per events; Number of hits; Events",
auto hfsam = d1.Histo1D({"hfsam", "Sampling Fraction; Sampling Fraction; Events", 150, 0.0, 0.15}, "fsam"); "Energy Deposit; Energy Deposit [GeV]; Events",
auto hpid = d1.Histo1D({"hpid", "PID; PID; Count", 100, -220, 220}, "pid"); "dE Relative; dE Relative; Events"
auto hdau = d1.Histo1D({"hdau", "Number of Daughters; Number of Daughters; Count", 10, 0, 10}, "dau"); };
// put this line above GetMean below so it can be lazily evaluated with the histgorams.
auto nevents_thrown = d1.Count();
const double meanE = hEthr->GetMean();
// Event Counts std::vector<std::vector<double>> range = {{0, 7.5}, {0, 2000}, {0, 2}, {-3, 3}};
std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n"; std::vector<std::string> col = {"Ethr", "nhits", "Esim", "dE_rel"};
// Draw Histograms double meanE = 5;
TCanvas *c1 = new TCanvas("c1", "c1", 700, 500); int nCol = range.size();
for (int i = 0; i < nCol; i++){
int binNum = 100;
auto h = d1.Histo1D({"hist", titleStr[i].c_str(), binNum, range[i][0], range[i][1]}, col[i].c_str());
if (col[i] == "Ethr"){
meanE = h->GetMean();
}
auto *c = new TCanvas("c", "c", 700, 500);
c->SetLogy(1);
auto h1 = h->DrawCopy();
h1->GetYaxis()->SetTitleOffset(1.4);
h1->SetLineWidth(2);
h1->SetLineColor(kBlue);
c->SaveAs((fmt::format("results/emcal_barrel_pi0_{}.png", col[i])).c_str());
c->SaveAs((fmt::format("results/emcal_barrel_pi0_{}.pdf", col[i])).c_str());
//std::printf("Generic %d\n", i);
}
// Resolution Plots
titleStr = {
"Sampling Fraction; Sampling Fraction; Events",
"dE; dE[GeV]; Events"
};
range = {{0,0.15}, {-3, 3}};
col = {"fsam", "dE"};
nCol = range.size();
std::printf("Here %d\n", 10);
std::vector<std::vector<double>> fitRange = {{0.005, 0.1}, {-3, 3}};
double sigmaOverE = 0;
auto hr = d1.Histo1D({"histr", titleStr[0].c_str(), 150, range[0][0], range[0][1]}, col[0].c_str());
auto *c = new TCanvas("c", "c", 700, 500);
c->SetLogy(1);
auto h2 = hr->DrawCopy();
h2->GetYaxis()->SetTitleOffset(1.4);
h2->SetLineWidth(2);
h2->SetLineColor(kBlue);
h2->Fit("gaus","","", fitRange[0][0], fitRange[0][1]);
h2->GetFunction("gaus")->SetLineWidth(2);
h2->GetFunction("gaus")->SetLineColor(kRed);
c->SaveAs((fmt::format("results/emcal_barrel_pi0_{}.png", col[0])).c_str());
c->SaveAs((fmt::format("results/emcal_barrel_pi0_{}.pdf", col[0])).c_str());
std::printf("Resolution %d\n", 0);
auto hs = d1.Histo1D({"hists", titleStr[1].c_str(), 100, range[1][0], range[1][1]}, col[1].c_str());
auto *c1 = new TCanvas("c1", "c1", 700, 500);
c1->SetLogy(1); c1->SetLogy(1);
fmt::print("0\n"); auto h3 = hs->DrawCopy();
auto h1 = hEthr->DrawCopy(); h3->GetYaxis()->SetTitleOffset(1.4);
//h1->GetYaxis()->SetTitleOffset(1.4); h3->SetLineWidth(2);
h1->SetLineWidth(2); h3->SetLineColor(kBlue);
h1->SetLineColor(kBlue); auto fit = h3->Fit("gaus","","", fitRange[1][0], fitRange[1][1]);
c1->SaveAs("results/emcal_barrel_pi0_Ethr.png"); double* res = h3->GetFunction("gaus")->GetParameters();
c1->SaveAs("results/emcal_barrel_pi0_Ethr.pdf"); sigmaOverE = res[2] / meanE;
fmt::print("1\n"); c1->SaveAs((fmt::format("results/emcal_barrel_pi0_{}.png", col[1])).c_str());
c1->SaveAs((fmt::format("results/emcal_barrel_pi0_{}.pdf", col[1])).c_str());
//TCanvas *c2 = new TCanvas("c2", "c2", 700, 500); std::printf("Resolution %d\n", 1);
//c2->SetLogy(1);
//h1 = hNhits->DrawCopy(); // Energy Resolution Calculation
//h1->GetYaxis()->SetTitleOffset(1.4); std::string test_tag = "Barrel_emcal_pi0";// TODO: Change test_tag to something else
//h1->SetLineWidth(2); std:string detEle = "Barrel_emcal";
//h1->SetLineColor(kBlue);
//c2->SaveAs("results/emcal_barrel_pi0_nhits.png"); // Energy resolution in the barrel region (-1 < eta < 1)
//c2->SaveAs("results/emcal_barrel_pi0_nhits.pdf"); // Taken from : Initial considerations for EMCal of the EIC detector by A. Bazilevsky
// sigma_E / E = 12% / E^0.5 convoluted with 2%
//TCanvas *c3 = new TCanvas("c3", "c3", 700, 500); // sigma_E / E = [ (0.12/E^0.5)^2 + 0.02^2]^0.5, with E in [GeV]
//c3->SetLogy(1); double resolutionTarget = TMath::Sqrt(0.12 * 0.12 / meanE + 0.02 * 0.02);
//h1 = hEsim->DrawCopy();
//h1->GetYaxis()->SetTitleOffset(1.4); common_bench::Test pi0_energy_resolution{
//h1->SetLineWidth(2); {{"name", fmt::format("{}_energy_resolution", test_tag)},
//h1->SetLineColor(kBlue); {"title", "Pi0 Energy resolution"},
//c3->SaveAs("results/emcal_barrel_pi0_Esim.png"); {"description",
//c3->SaveAs("results/emcal_barrel_pi0_Esim.pdf"); fmt::format("Pi0 energy resolution for {}, estimated using a Gaussian fit.", detEle)},
{"quantity", "resolution (in %)"},
//TCanvas *c4 = new TCanvas("c4", "c4", 700, 500); {"target", std::to_string(resolutionTarget)}}
//c4->SetLogy(1); };
//h1 = hfsam->DrawCopy();
//h1->GetYaxis()->SetTitleOffset(1.4);
//h1->SetLineWidth(2);
//h1->SetLineColor(kBlue);
//h1->Fit("gaus","","",0.005,0.1);
//h1->GetFunction("gaus")->SetLineWidth(2);
//h1->GetFunction("gaus")->SetLineColor(kRed);
//c4->SaveAs("results/emcal_barrel_pi0_fsam.png");
//c4->SaveAs("results/emcal_barrel_pi0_fsam.pdf");
//TCanvas *c5 = new TCanvas("c5", "c5", 700, 500);
//c5->SetLogy(1);
//h1 = hpid->DrawCopy();
//h1->GetYaxis()->SetTitleOffset(1.4);
//h1->SetLineWidth(2);
//h1->SetLineColor(kBlue);
//c5->SaveAs("results/emcal_barrel_pi0_pid.png");
//c5->SaveAs("results/emcal_barrel_pi0_pid.pdf");
//TCanvas *c6 = new TCanvas("c6", "c6", 700, 500);
//c5->SetLogy(1);
//h1 = hdau->DrawCopy();
//h1->GetYaxis()->SetTitleOffset(1.4);
//h1->SetLineWidth(2);
//h1->SetLineColor(kBlue);
//c6->SaveAs("results/emcal_barrel_pi0_dau.png");
//c6->SaveAs("results/emcal_barrel_pi0_dau.pdf");
//// Energy Resolution Calculation
//std::string test_tag = "Barrel_emcal_pi0";// TODO: Change test_tag to something else
//std:string detEle = "Barrel_emcal";
//// 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 resolutionTarget = TMath::Sqrt(0.12 * 0.12 / meanE + 0.02 * 0.02);
//common_bench::Test pi0_energy_resolution{
// {{"name", fmt::format("{}_energy_resolution", test_tag)},
// {"title", "Pi0 Energy resolution"},
// {"description",
// fmt::format("Pi0 energy resolution for {}, estimated using a Gaussian fit.", detEle)},
// {"quantity", "resolution (in %)"},
// {"target", std::to_string(resolutionTarget)}}
//};
//// Histograms and Fitting
//auto hdE = d1.Histo1D({"hdE", "dE; dE[GeV]; Events", 100, -3.0, 3.0}, "dE");
//auto hdE_rel = d1.Histo1D({"hdE_rel", "dE Relative; dE Relative; Events", 100, -3.0, 3.0}, "dE_rel");
//auto hdEcopy = hdE->DrawCopy();
//hdEcopy->Fit("gaus", "", "", -3.0, 3.0);
//double* res = hdEcopy->GetFunction("gaus")->GetParameters();
//double sigmaOverE = res[2] / meanE;
//// Pass/Fail //// Pass/Fail
//sigmaOverE <= resolutionTarget ? pi0_energy_resolution.pass(sigmaOverE) : pi0_energy_resolution.fail(sigmaOverE); sigmaOverE <= resolutionTarget ? pi0_energy_resolution.pass(sigmaOverE) : pi0_energy_resolution.fail(sigmaOverE);
//std::printf("Energy Resolution is %f\n", res[2]); common_bench::write_test({pi0_energy_resolution}, fmt::format("results/{}_pi0.json", detEle));
//// Energy Resolution Histogram Plotting
//auto *cdE = new TCanvas("cdE", "cdE", 700, 500);
//cdE->SetLogy(1);
//h1 = hdEcopy->DrawCopy();
//hdEcopy->GetYaxis()->SetTitleOffset(1.4);
//hdEcopy->SetLineWidth(2);
//hdEcopy->SetLineColor(kBlue);
//hdEcopy->GetFunction("gaus")->SetLineWidth(2);
//hdEcopy->GetFunction("gaus")->SetLineColor(kRed);
//cdE->SaveAs("results/emcal_barrel_pi0_dE.png");
//cdE->SaveAs("results/emcal_barrel_pi0_dE.pdf");
//auto *cdE_rel = new TCanvas("cdE_rel", "cdE_rel", 700, 500);
//h1 = hdE_rel->DrawCopy();
//hdE_rel->GetYaxis()->SetTitleOffset(1.4);
//hdE_rel->SetLineWidth(2);
//hdE_rel->SetLineColor(kBlue);
//cdE_rel->SaveAs("results/emcal_barrel_pi0_dE_rel.png");
//cdE_rel->SaveAs("results/emcal_barrel_pi0_dE_rel.pdf");
//common_bench::write_test({pi0_energy_resolution}, fmt::format("results/{}_pi0.json", detEle));
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment