Skip to content
Snippets Groups Projects
Commit 13318567 authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

Revert "Resolve "Broken pi0 ecal barrel analysis script"" 

This reverts commit ed85f282.
parent ed85f282
No related branches found
No related tags found
1 merge request!76Fix broken pipeline
Hide whitespace changes
Inline Side-by-side
benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_analysis.cxx
+ 119
− 100
View file @ 13318567
...@@ -29,11 +29,9 @@ using ROOT::RDataFrame; ...@@ -29,11 +29,9 @@ using ROOT::RDataFrame;
using namespace ROOT::VecOps; using namespace ROOT::VecOps;
using json = nlohmann::json; using json = nlohmann::json;
void emcal_barrel_pi0_analysis( 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_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);
...@@ -99,28 +97,33 @@ void emcal_barrel_pi0_analysis( ...@@ -99,28 +97,33 @@ void emcal_barrel_pi0_analysis(
}; };
// 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");
const double meanE = hEthr->GetMean(); // 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::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);
...@@ -135,89 +138,105 @@ void emcal_barrel_pi0_analysis( ...@@ -135,89 +138,105 @@ void emcal_barrel_pi0_analysis(
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); //TCanvas *c4 = new TCanvas("c4", "c4", 700, 500);
c4->SetLogy(1); //c4->SetLogy(1);
h1 = hfsam->DrawCopy(); //h1 = hfsam->DrawCopy();
h1->GetYaxis()->SetTitleOffset(1.4); //h1->GetYaxis()->SetTitleOffset(1.4);
h1->SetLineWidth(2); //h1->SetLineWidth(2);
h1->SetLineColor(kBlue); //h1->SetLineColor(kBlue);
h1->Fit("gaus","","",0.005,0.1); //h1->Fit("gaus","","",0.005,0.1);
h1->GetFunction("gaus")->SetLineWidth(2); //h1->GetFunction("gaus")->SetLineWidth(2);
h1->GetFunction("gaus")->SetLineColor(kRed); //h1->GetFunction("gaus")->SetLineColor(kRed);
c4->SaveAs("results/emcal_barrel_pi0_fsam.png"); //c4->SaveAs("results/emcal_barrel_pi0_fsam.png");
c4->SaveAs("results/emcal_barrel_pi0_fsam.pdf"); //c4->SaveAs("results/emcal_barrel_pi0_fsam.pdf");
TCanvas *cdE_rel = new TCanvas("cdE_rel", "cdE_rel", 700, 500); //TCanvas *c5 = new TCanvas("c5", "c5", 700, 500);
h1 = hdE_rel->DrawCopy(); //c5->SetLogy(1);
h1->GetYaxis()->SetTitleOffset(1.4); //h1 = hpid->DrawCopy();
h1->SetLineWidth(2); //h1->GetYaxis()->SetTitleOffset(1.4);
h1->SetLineColor(kBlue); //h1->SetLineWidth(2);
cdE_rel->SaveAs("results/emcal_barrel_pi0_dE_rel.png"); //h1->SetLineColor(kBlue);
cdE_rel->SaveAs("results/emcal_barrel_pi0_dE_rel.pdf"); //c5->SaveAs("results/emcal_barrel_pi0_pid.png");
//c5->SaveAs("results/emcal_barrel_pi0_pid.pdf");
h1->Delete();
//TCanvas *c6 = new TCanvas("c6", "c6", 700, 500);
// Energy Resolution Calculation //c5->SetLogy(1);
std::string test_tag = "Barrel_emcal_pi0";// TODO: Change test_tag to something else //h1 = hdau->DrawCopy();
std:string detEle = "Barrel_emcal"; //h1->GetYaxis()->SetTitleOffset(1.4);
//h1->SetLineWidth(2);
// Energy resolution in the barrel region (-1 < eta < 1) //h1->SetLineColor(kBlue);
// Taken from : Initial considerations for EMCal of the EIC detector by A. Bazilevsky //c6->SaveAs("results/emcal_barrel_pi0_dau.png");
// sigma_E / E = 12% / E^0.5 convoluted with 2% //c6->SaveAs("results/emcal_barrel_pi0_dau.pdf");
// sigma_E / E = [ (0.12/E^0.5)^2 + 0.02^2]^0.5, with E in [GeV]
//// 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); //double resolutionTarget = TMath::Sqrt(0.12 * 0.12 / meanE + 0.02 * 0.02);
common_bench::Test pi0_energy_resolution{ //common_bench::Test pi0_energy_resolution{
{{"name", fmt::format("{}_energy_resolution", test_tag)}, // {{"name", fmt::format("{}_energy_resolution", test_tag)},
{"title", "Pi0 Energy resolution"}, // {"title", "Pi0 Energy resolution"},
{"description", // {"description",
fmt::format("Pi0 energy resolution for {}, estimated using a Gaussian fit.", detEle)}, // fmt::format("Pi0 energy resolution for {}, estimated using a Gaussian fit.", detEle)},
{"quantity", "resolution (in %)"}, // {"quantity", "resolution (in %)"},
{"target", std::to_string(resolutionTarget)}} // {"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 h2 = hdE->DrawCopy(); //auto hdE_rel = d1.Histo1D({"hdE_rel", "dE Relative; dE Relative; Events", 100, -3.0, 3.0}, "dE_rel");
h2->Fit("gaus", "", "", -3.0, 3.0); //auto hdEcopy = hdE->DrawCopy();
double* res = h2->GetFunction("gaus")->GetParameters(); //hdEcopy->Fit("gaus", "", "", -3.0, 3.0);
double sigmaOverE = res[2] / meanE; //double* res = hdEcopy->GetFunction("gaus")->GetParameters();
auto tf = h2->GetFunction("gaus"); //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);
common_bench::write_test({pi0_energy_resolution}, fmt::format("results/{}_pi0.json", detEle)); //std::printf("Energy Resolution is %f\n", res[2]);
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);
//auto h3 = hdE->DrawCopy(); //h1 = hdEcopy->DrawCopy();
h2->GetYaxis()->SetTitleOffset(1.4); //hdEcopy->GetYaxis()->SetTitleOffset(1.4);
h2->SetLineWidth(2); //hdEcopy->SetLineWidth(2);
h2->SetLineColor(kBlue); //hdEcopy->SetLineColor(kBlue);
tf->SetLineWidth(2); //hdEcopy->GetFunction("gaus")->SetLineWidth(2);
tf->SetLineColor(kRed); //hdEcopy->GetFunction("gaus")->SetLineColor(kRed);
tf->DrawClone("same"); //cdE->SaveAs("results/emcal_barrel_pi0_dE.png");
cdE->SaveAs("results/emcal_barrel_pi0_dE.png"); //cdE->SaveAs("results/emcal_barrel_pi0_dE.pdf");
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