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

Resolve "Broken pi0 ecal barrel analysis script"

parent 011487d8
No related branches found
No related tags found
1 merge request!73Resolve "Broken pi0 ecal barrel analysis script"
Hide whitespace changes
Inline Side-by-side
benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_analysis.cxx
+ 100
− 119
View file @ ed85f282
...@@ -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);
...@@ -97,33 +99,28 @@ void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_b ...@@ -97,33 +99,28 @@ 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"); 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 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, 2.0}, "Esim"); auto hEsim = d1.Histo1D({"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 100, 0.0, 2.0}, "Esim");
auto hfsam = d1.Histo1D({"hfsam", "Sampling Fraction; Sampling Fraction; Events", 150, 0.0, 0.15}, "fsam"); auto hfsam = d1.Histo1D({"hfsam", "Sampling Fraction; Sampling Fraction; Events", 150, 0.0, 0.15}, "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");
auto hdau = d1.Histo1D({"hdau", "Number of Daughters; Number of Daughters; Count", 10, 0, 10}, "dau"); auto hdau = d1.Histo1D({"hdau", "Number of Daughters; Number of Daughters; Count", 10, 0, 10}, "dau");
auto hdE_rel = d1.Histo1D({"hdE_rel", "dE Relative; dE Relative; Events", 100, -3.0, 3.0}, "dE_rel");
// put this line above GetMean below so it can be lazily evaluated with the histgorams. const double meanE = hEthr->GetMean();
auto nevents_thrown = d1.Count();
const double meanE = hEthr->GetMean();
// Event Counts
std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
// Draw Histograms // Draw Histograms
TCanvas *c1 = new TCanvas("c1", "c1", 700, 500); TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
...@@ -138,105 +135,89 @@ void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_b ...@@ -138,105 +135,89 @@ void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_b
fmt::print("1\n"); fmt::print("1\n");
//TCanvas *c2 = new TCanvas("c2", "c2", 700, 500); TCanvas *c2 = new TCanvas("c2", "c2", 700, 500);
//c2->SetLogy(1); c2->SetLogy(1);
//h1 = hNhits->DrawCopy(); h1 = hNhits->DrawCopy();
//h1->GetYaxis()->SetTitleOffset(1.4); h1->GetYaxis()->SetTitleOffset(1.4);
//h1->SetLineWidth(2); h1->SetLineWidth(2);
//h1->SetLineColor(kBlue); h1->SetLineColor(kBlue);
//c2->SaveAs("results/emcal_barrel_pi0_nhits.png"); c2->SaveAs("results/emcal_barrel_pi0_nhits.png");
//c2->SaveAs("results/emcal_barrel_pi0_nhits.pdf"); c2->SaveAs("results/emcal_barrel_pi0_nhits.pdf");
//TCanvas *c3 = new TCanvas("c3", "c3", 700, 500); TCanvas *c3 = new TCanvas("c3", "c3", 700, 500);
//c3->SetLogy(1); c3->SetLogy(1);
//h1 = hEsim->DrawCopy(); h1 = hEsim->DrawCopy();
//h1->GetYaxis()->SetTitleOffset(1.4); h1->GetYaxis()->SetTitleOffset(1.4);
//h1->SetLineWidth(2); h1->SetLineWidth(2);
//h1->SetLineColor(kBlue); h1->SetLineColor(kBlue);
//c3->SaveAs("results/emcal_barrel_pi0_Esim.png"); c3->SaveAs("results/emcal_barrel_pi0_Esim.png");
//c3->SaveAs("results/emcal_barrel_pi0_Esim.pdf"); c3->SaveAs("results/emcal_barrel_pi0_Esim.pdf");
//TCanvas *c4 = new TCanvas("c4", "c4", 700, 500);
//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{ TCanvas *c4 = new TCanvas("c4", "c4", 700, 500);
// {{"name", fmt::format("{}_energy_resolution", test_tag)}, c4->SetLogy(1);
// {"title", "Pi0 Energy resolution"}, h1 = hfsam->DrawCopy();
// {"description", h1->GetYaxis()->SetTitleOffset(1.4);
// fmt::format("Pi0 energy resolution for {}, estimated using a Gaussian fit.", detEle)}, h1->SetLineWidth(2);
// {"quantity", "resolution (in %)"}, h1->SetLineColor(kBlue);
// {"target", std::to_string(resolutionTarget)}} 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 *cdE_rel = new TCanvas("cdE_rel", "cdE_rel", 700, 500);
h1 = hdE_rel->DrawCopy();
h1->GetYaxis()->SetTitleOffset(1.4);
h1->SetLineWidth(2);
h1->SetLineColor(kBlue);
cdE_rel->SaveAs("results/emcal_barrel_pi0_dE_rel.png");
cdE_rel->SaveAs("results/emcal_barrel_pi0_dE_rel.pdf");
h1->Delete();
// 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 //// Histograms and Fitting
//auto hdE = d1.Histo1D({"hdE", "dE; dE[GeV]; Events", 100, -3.0, 3.0}, "dE"); 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 h2 = hdE->DrawCopy();
//auto hdEcopy = hdE->DrawCopy(); h2->Fit("gaus", "", "", -3.0, 3.0);
//hdEcopy->Fit("gaus", "", "", -3.0, 3.0); double* res = h2->GetFunction("gaus")->GetParameters();
//double* res = hdEcopy->GetFunction("gaus")->GetParameters(); double sigmaOverE = res[2] / meanE;
//double sigmaOverE = res[2] / meanE; auto tf = h2->GetFunction("gaus");
//// 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));
std::printf("Energy Resolution is %f\n", res[2]);
//// Energy Resolution Histogram Plotting //// Energy Resolution Histogram Plotting
//auto *cdE = new TCanvas("cdE", "cdE", 700, 500); auto *cdE = new TCanvas("cdE", "cdE", 700, 500);
//cdE->SetLogy(1); cdE->SetLogy(1);
//h1 = hdEcopy->DrawCopy(); //auto h3 = hdE->DrawCopy();
//hdEcopy->GetYaxis()->SetTitleOffset(1.4); h2->GetYaxis()->SetTitleOffset(1.4);
//hdEcopy->SetLineWidth(2); h2->SetLineWidth(2);
//hdEcopy->SetLineColor(kBlue); h2->SetLineColor(kBlue);
//hdEcopy->GetFunction("gaus")->SetLineWidth(2); tf->SetLineWidth(2);
//hdEcopy->GetFunction("gaus")->SetLineColor(kRed); tf->SetLineColor(kRed);
//cdE->SaveAs("results/emcal_barrel_pi0_dE.png"); tf->DrawClone("same");
//cdE->SaveAs("results/emcal_barrel_pi0_dE.pdf"); 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