Skip to content
Snippets Groups Projects
Commit 578054fb authored by Marshall Scott's avatar Marshall Scott Committed by Whitney Armstrong
Browse files

Resolve "pi0 resolution for ECal barrel"

parent a116f39e
No related branches found
No related tags found
1 merge request!24Resolve "pi0 resolution for ECal barrel"
{
// top-level include-dir
gROOT->ProcessLine(".include include");
// Setting for Graphs
gROOT->SetStyle("Plain");
gStyle->SetOptFit(1);
gStyle->SetLineWidth(2);
gStyle->SetPadTickX(1);
gStyle->SetPadTickY(1);
gStyle->SetPadGridX(1);
gStyle->SetPadGridY(1);
gStyle->SetPadLeftMargin(0.14);
}
...@@ -34,4 +34,5 @@ collect_results:barrel_ecal: ...@@ -34,4 +34,5 @@ collect_results:barrel_ecal:
- ["bench:emcal_barrel_pions","bench:emcal_barrel_electrons"] - ["bench:emcal_barrel_pions","bench:emcal_barrel_electrons"]
script: script:
- ls -lrht - ls -lrht
- echo " FIX ME"
{ {
/*
// top-level include-dir
gROOT->ProcessLine(".include detector_benchmarks/include");
// Setting for Graphs // Setting for Graphs
gROOT->SetStyle("Plain"); gROOT->SetStyle("Plain");
gStyle->SetOptFit(1); gStyle->SetOptFit(1);
...@@ -8,4 +12,6 @@ ...@@ -8,4 +12,6 @@
gStyle->SetPadGridX(1); gStyle->SetPadGridX(1);
gStyle->SetPadGridY(1); gStyle->SetPadGridY(1);
gStyle->SetPadLeftMargin(0.14); gStyle->SetPadLeftMargin(0.14);
*/
} }
...@@ -35,7 +35,7 @@ void emcal_barrel_pions(int n_events = 1e6, double e_start = 0.0, double e_end = ...@@ -35,7 +35,7 @@ void emcal_barrel_pions(int n_events = 1e6, double e_start = 0.0, double e_end =
// FourVector(px,py,pz,e,pdgid,status) // FourVector(px,py,pz,e,pdgid,status)
// type 4 is beam // type 4 is beam
// pdgid 11 - electron // pdgid 11 - electron
// pdgid 111 - pi0 // pdgid 111 - pi0
// pdgid 2212 - proton // pdgid 2212 - proton
GenParticlePtr p1 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4); GenParticlePtr p1 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
GenParticlePtr p2 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4); GenParticlePtr p2 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
...@@ -55,13 +55,16 @@ void emcal_barrel_pions(int n_events = 1e6, double e_start = 0.0, double e_end = ...@@ -55,13 +55,16 @@ void emcal_barrel_pions(int n_events = 1e6, double e_start = 0.0, double e_end =
// type 1 is final state // type 1 is final state
// pdgid 211 - pion+ 139.570 MeV/c^2 // pdgid 211 - pion+ 139.570 MeV/c^2
GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.139570 * 0.139570))), 211, 1); // pdgid 111 - pion0 134.977 MeV/c^2
//GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.139570 * 0.139570))), 211, 1);
GenParticlePtr p4 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.134977 * 0.134977))), 111, 1);
GenVertexPtr v1 = std::make_shared<GenVertex>(); GenVertexPtr v1 = std::make_shared<GenVertex>();
v1->add_particle_in(p1); v1->add_particle_in(p1);
v1->add_particle_in(p2); v1->add_particle_in(p2);
v1->add_particle_out(p3); //v1->add_particle_out(p3);
v1->add_particle_out(p4);
evt.add_vertex(v1); evt.add_vertex(v1);
if (events_parsed == 0) { if (events_parsed == 0) {
......
...@@ -9,12 +9,17 @@ ...@@ -9,12 +9,17 @@
#include "dd4pod/Geant4ParticleCollection.h" #include "dd4pod/Geant4ParticleCollection.h"
#include "dd4pod/CalorimeterHitCollection.h" #include "dd4pod/CalorimeterHitCollection.h"
#include "benchmark.h"
#include "mt.h"
#include "util.h"
#include "TCanvas.h" #include "TCanvas.h"
#include "TStyle.h" #include "TStyle.h"
#include "TMath.h" #include "TMath.h"
#include "TH1.h" #include "TH1.h"
#include "TF1.h" #include "TF1.h"
#include "TH1D.h" #include "TH1D.h"
#include "TFitResult.h"
using ROOT::RDataFrame; using ROOT::RDataFrame;
using namespace ROOT::VecOps; using namespace ROOT::VecOps;
...@@ -32,9 +37,32 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal ...@@ -32,9 +37,32 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
gStyle->SetPadLeftMargin(0.14); gStyle->SetPadLeftMargin(0.14);
gStyle->SetPadRightMargin(0.14); gStyle->SetPadRightMargin(0.14);
//Tests
std::string test_tag = "Barrel_emcal_pi0";
//TODO: Change test_tag to something else
std:string detector = "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 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)},
{"title", "Pion0 Energy resolution"},
{"description",
fmt::format("Pion0 energy resolution with {}, estimated using a Gaussian fit.", detector)},
{"quantity", "resolution (in %)"},
{"target", std::to_string(resolutionTarget)}}};
ROOT::EnableImplicitMT(); ROOT::EnableImplicitMT();
ROOT::RDataFrame d0("events", input_fname); ROOT::RDataFrame d0("events", input_fname);
// Sampling Fraction
double samp_frac = 0.0136;
// Thrown Energy [GeV] // Thrown Energy [GeV]
auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) { auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
std::vector<double> result; std::vector<double> result;
...@@ -58,25 +86,69 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal ...@@ -58,25 +86,69 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
// Sampling fraction = Esampling / Ethrown // Sampling fraction = Esampling / Ethrown
auto fsam = [](const std::vector<double>& sampled, const std::vector<double>& thrown) { auto fsam = [](const std::vector<double>& sampled, const std::vector<double>& thrown) {
std::vector<double> result; std::vector<double> result;
for (const auto& E1 : thrown) { auto it_sam = sampled.cbegin();
for (const auto& E2 : sampled) auto it_thr = thrown.cbegin();
result.push_back(E2/E1); for (; it_sam != sampled.end() && it_thr != thrown.end(); ++it_sam, ++it_thr) {
result.push_back(*it_sam / *it_thr);
}
return result;
};
// 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 / samp_frac - *it_thr);
} }
return result; return result;
}; };
// 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 / samp_frac - *it_thr) / *it_thr);
}
return result;
};
// Returns the pdgID of the particle
auto getpid = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
return input[2].pdgID;
};
// Returns number of particle daughters
auto getdau = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
return input[2].daughters_begin;
};
// 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("Esim", Esim, {"EcalBarrelHits"}) .Define("Esim", Esim, {"EcalBarrelHits"})
.Define("fsam", fsam, {"Esim","Ethr"}) .Define("fsam", fsam, {"Esim","Ethr"})
; .Define("pid", getpid, {"mcparticles"})
.Define("dau", getdau, {"mcparticles"})
;
// 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, 1.0}, "Esim"); 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 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 hdau = d1.Histo1D({"hdau", "Number of Daughters; Number of Daughters; Count", 10, 0, 10}, "dau");
// Set sampling Fraction, ideally this will be taken from a json file
samp_frac = hfsam -> GetMean();
auto d2 = d1.Define("dE", eResol, {"Esim","Ethr"})
.Define("dE_rel", eResol_rel, {"Esim","Ethr"})
;
// Event Counts // Event Counts
auto nevents_thrown = d1.Count(); auto nevents_thrown = d1.Count();
...@@ -121,4 +193,61 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal ...@@ -121,4 +193,61 @@ void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/sim_emcal
hfsam->DrawClone(); hfsam->DrawClone();
c4->SaveAs("results/emcal_barrel_pions_fsam.png"); c4->SaveAs("results/emcal_barrel_pions_fsam.png");
c4->SaveAs("results/emcal_barrel_pions_fsam.pdf"); c4->SaveAs("results/emcal_barrel_pions_fsam.pdf");
TCanvas *c5 = new TCanvas("c5", "c5", 700, 500);
c5->SetLogy(1);
hpid->GetYaxis()->SetTitleOffset(1.4);
hpid->SetLineWidth(2);
hpid->SetLineColor(kBlue);
hpid->DrawClone();
c5->SaveAs("results/emcal_barrel_pions_pid.png");
c5->SaveAs("results/emcal_barrel_pions_pid.pdf");
TCanvas *c6 = new TCanvas("c6", "c6", 700, 500);
c5->SetLogy(1);
hdau->GetYaxis()->SetTitleOffset(1.4);
hdau->SetLineWidth(2);
hdau->SetLineColor(kBlue);
hdau->DrawClone();
c6->SaveAs("results/emcal_barrel_pions_dau.png");
c6->SaveAs("results/emcal_barrel_pions_dau.pdf");
//Energy Resolution Calculation
auto hdE = d2.Histo1D({"hdE", "dE; dE[GeV]; Events", 100, -3.0, 3.0}, "dE");
auto hdE_rel = d2.Histo1D({"hdE_rel", "dE Relative; dE Relative; Events", 100, -3.0, 3.0}, "dE_rel");
hdE->Fit("gaus", "", "", -3.0, 3.0);
double* res = hdE->GetFunction("gaus")->GetParameters();
double sigmaOverE = res[2] / thrown_energy;
//Pass/Fail
if (sigmaOverE <= resolutionTarget) {
pi0_energy_resolution.pass(sigmaOverE);
} else {
pi0_energy_resolution.fail(sigmaOverE);
}
//std::printf("Energy Resolution is %f\n", res[2]);
//Energy Resolution Histogram Plotting
auto *cdE = new TCanvas("cdE", "cdE", 700, 500);
cdE->SetLogy(1);
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");
auto *cdE_rel = new TCanvas("cdE_rel", "cdE_rel", 700, 500);
hdE_rel->GetYaxis()->SetTitleOffset(1.4);
hdE_rel->SetLineWidth(2);
hdE_rel->SetLineColor(kBlue);
hdE_rel->DrawClone();
cdE_rel->SaveAs("results/emcal_barrel_pi0_dE_rel.png");
cdE_rel->SaveAs("results/emcal_barrel_pi0_dE_rel.pdf");
eic::util::write_test({pi0_energy_resolution}, fmt::format("{}_pions.json", detector));
} }
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment