Skip to content
Snippets Groups Projects
Commit b132d1ef authored by Jihee Kim's avatar Jihee Kim
Browse files

Update Event Display Script

parent 14155fff
No related branches found
No related tags found
1 merge request!54Update Event Display Script
Hide whitespace changes
Inline Side-by-side
ecal/options/crystal_calorimeter_reco.py
+ 1
− 1
View file @ b132d1ef
......@@ -59,7 +59,7 @@ emcalreco = CrystalEndcapsReco("ecal_reco",
emcalcluster = IslandCluster("emcal_cluster",
inputHitCollection="RecoEcalHits",
outputClusterCollection="EcalClusters",
minClusterCenterEdep=0.00001*units.MeV,
minClusterCenterEdep=50.0*units.MeV,
groupRange=2.0)
clusterreco = RecoCoG("cluster_reco",
......
ecal/scripts/emcal_electrons.cxx
+ 2
− 2
View file @ b132d1ef
......@@ -31,8 +31,8 @@ void emcal_electrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30
// Constraining the solid angle, but larger than that subtended by the detector
// https://indico.bnl.gov/event/7449/contributions/35966/attachments/27177/41430/EIC-DWG-Calo-03192020.pdf
// See a figure on slide 26
double cos_theta_min = std::cos(M_PI * (155.0 / 180.0));
double cos_theta_max = std::cos(M_PI);
double cos_theta_min = std::cos(M_PI * (160.0 / 180.0));
double cos_theta_max = std::cos(M_PI * (175.0 / 180.0));
for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
// FourVector(px,py,pz,e,pdgid,status)
......
ecal/scripts/emcal_electrons_analysis.cxx
+ 178
− 1
View file @ b132d1ef
......@@ -23,7 +23,7 @@
using ROOT::RDataFrame;
using namespace ROOT::VecOps;
void emcal_electrons_analysis(const char* input_fname = "rec_electron_10kEvt.root")
void emcal_electrons_analysis(const char* input_fname = "rec_emcal_uniform_electrons.root")
{
// Setting for graphs
gROOT->SetStyle("Plain");
......@@ -34,6 +34,7 @@ void emcal_electrons_analysis(const char* input_fname = "rec_electron_10kEvt.roo
gStyle->SetPadGridX(1);
gStyle->SetPadGridY(1);
gStyle->SetPadLeftMargin(0.14);
gStyle->SetPadRightMargin(0.14);
ROOT::EnableImplicitMT();
ROOT::RDataFrame d0("events", input_fname);
......@@ -150,6 +151,38 @@ void emcal_electrons_analysis(const char* input_fname = "rec_electron_10kEvt.roo
auto d2 = d1.Filter("ncluster==1");
auto hClusterE1 = d2.Histo1D({"hClusterE1", "One Cluster Energy; Cluster Energy [GeV]; Events", 100, -0.5, 30.5}, "clusterE");
auto hEres = d2.Histo1D({"hEres", "Energy Resolution; #DeltaE/E; Events", 100, -1.0, 1.0}, "E_res");
auto hEres_lmom = d2.Filter([=](std::vector<dd4pod::Geant4ParticleData> const& input) {
auto Pthr = TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass);
if (Pthr <= 15.0)
return true;
return false;}, {"mcparticles2"})
.Histo1D({"hEres_lmom", "Energy Resolution; #DeltaE/E; Events", 100, -1.0, 1.0}, "E_res");
auto hEres_hmom = d2.Filter([=](std::vector<dd4pod::Geant4ParticleData> const& input) {
auto Pthr = TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass);
if (Pthr > 15.0)
return true;
return false;}, {"mcparticles2"})
.Histo1D({"hEres_hmom", "Energy Resolution; #DeltaE/E; Events", 100, -1.0, 1.0}, "E_res");
auto hEres_leta = d2.Filter([=](const std::vector<eic::ClusterData>& evt) {
for(const auto& i: evt) {
auto r = TMath::Sqrt(i.position.x*i.position.x + i.position.y*i.position.y + i.position.z*i.position.z);
auto th = TMath::ACos(i.position.z/r)*TMath::RadToDeg();
auto eta = -1.0 * TMath::Log(TMath::Tan((th*TMath::DegToRad())/2.0));
if (eta <= -2.25)
return true;
}
return false;}, {"EcalClusters"})
.Histo1D({"hEres_leta", "Energy Resolution; #DeltaE/E; Events", 100, -1.0, 1.0}, "E_res");
auto hEres_heta = d2.Filter([=](const std::vector<eic::ClusterData>& evt) {
for(const auto& i: evt) {
auto r = TMath::Sqrt(i.position.x*i.position.x + i.position.y*i.position.y + i.position.z*i.position.z);
auto th = TMath::ACos(i.position.z/r)*TMath::RadToDeg();
auto eta = -1.0 * TMath::Log(TMath::Tan((th*TMath::DegToRad())/2.0));
if (eta > -2.25)
return true;
}
return false;}, {"EcalClusters"})
.Histo1D({"hEres_heta", "Energy Resolution; #DeltaE/E; Events", 100, -1.0, 1.0}, "E_res");
auto hPratio = d2.Histo1D({"hPratio", "Momentum ratio; p_{rec}/p_{thr}; Events", 100, 0.0, 1.0}, "p_ratio");
auto hPthr_accepted = d2.Filter([=] (const std::vector<float>& Prec, const std::vector<float>& Pthr) {
for (const auto& P1 : Pthr) {
......@@ -162,6 +195,9 @@ void emcal_electrons_analysis(const char* input_fname = "rec_electron_10kEvt.roo
}
return false;}, {"p_rec","p_thr"})
.Histo1D({"hPthr_accepted", "Thrown momentum for reconstructed particle; p_{thr} [GeV]; Events", 100, -0.5, 30.5}, "p_thr");
// 2D histograms
auto hPrec_Pthr = d2.Histo2D({"hPrec_Pthr", "p_{rec} vs p_{thr}; p_{rec}[GeV]; p_{thr[GeV]}", 100, -0.5, 30.5, 100, -0.5, 30.5}, "p_rec", "p_thr");
auto hEres_Ethr = d2.Histo2D({"hEres_Ethr", "Energy Resolution vs Ethr; #DeltaE/E; E_{thr}[GeV]", 100, -1.0, 1.0, 100, -0.5, 30.5}, "E_res", "E_thr");
// Cut on Radial Distance
auto d3 = d2.Filter([=] (const std::vector<eic::ClusterData>& evt) {
......@@ -174,9 +210,45 @@ void emcal_electrons_analysis(const char* input_fname = "rec_electron_10kEvt.roo
}
return false;}, {"EcalClusters"});
auto hEres_cut = d3.Histo1D({"hEres_cut", "Energy Resolution; #DeltaE/E; Events", 100, -1.0, 1.0}, "E_res");
auto hEres_cut_lmom = d3.Filter([=](std::vector<dd4pod::Geant4ParticleData> const& input) {
auto Pthr = TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass);
if (Pthr <= 15.0)
return true;
return false;}, {"mcparticles2"})
.Histo1D({"hEres_cut_lmom", "Energy Resolution; #DeltaE/E; Events", 100, -1.0, 1.0}, "E_res");
auto hEres_cut_hmom = d3.Filter([=](std::vector<dd4pod::Geant4ParticleData> const& input) {
auto Pthr = TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass);
if (Pthr > 15.0)
return true;
return false;}, {"mcparticles2"})
.Histo1D({"hEres_cut_hmom", "Energy Resolution; #DeltaE/E; Events", 100, -1.0, 1.0}, "E_res");
auto hEres_cut_leta = d3.Filter([=](const std::vector<eic::ClusterData>& evt) {
for(const auto& i: evt) {
auto r = TMath::Sqrt(i.position.x*i.position.x + i.position.y*i.position.y + i.position.z*i.position.z);
auto th = TMath::ACos(i.position.z/r)*TMath::RadToDeg();
auto eta = -1.0 * TMath::Log(TMath::Tan((th*TMath::DegToRad())/2.0));
if (eta <= -2.25)
return true;
}
return false;}, {"EcalClusters"})
.Histo1D({"hEres_cut_leta", "Energy Resolution; #DeltaE/E; Events", 100, -1.0, 1.0}, "E_res");
auto hEres_cut_heta = d3.Filter([=](const std::vector<eic::ClusterData>& evt) {
for(const auto& i: evt) {
auto r = TMath::Sqrt(i.position.x*i.position.x + i.position.y*i.position.y + i.position.z*i.position.z);
auto th = TMath::ACos(i.position.z/r)*TMath::RadToDeg();
auto eta = -1.0 * TMath::Log(TMath::Tan((th*TMath::DegToRad())/2.0));
if (eta > -2.25)
return true;
}
return false;}, {"EcalClusters"})
.Histo1D({"hEres_cut_heta", "Energy Resolution; #DeltaE/E; Events", 100, -1.0, 1.0}, "E_res");
auto hTheta_cut = d3.Histo1D({"hTheta_cut", "Scattering Angle; #theta [degree]; Events", 100, 130.0, 180.0}, "theta");
auto hEta_cut = d3.Histo1D({"hEta_cut", "Pseudo-rapidity; #eta; Events", 100, -5.0, 0.0}, "eta");
// 2D histograms
auto hPrec_Pthr_cut = d3.Histo2D({"hPrec_Pthr_cut", "p_{rec} vs p_{thr}; p_{rec}[GeV]; p_{thr[GeV]}", 100, -0.5, 30.5, 100, -0.5, 30.5}, "p_rec", "p_thr");
auto hEres_Ethr_cut = d3.Histo2D({"hEres_Ethr_cut", "Energy Resolution vs Ethr; #DeltaE/E; E_{thr}[GeV]", 100, -1.0, 1.0, 100, -0.5, 30.5}, "E_res", "E_thr");
// Event Counts
auto nevents_thrown = d1.Count();
auto nevents_cluster1 = d2.Count();
......@@ -312,4 +384,109 @@ void emcal_electrons_analysis(const char* input_fname = "rec_electron_10kEvt.roo
hEta_cut->DrawClone();
c15->SaveAs("results/emcal_electrons_eta_cut.png");
c15->SaveAs("results/emcal_electrons_eta_cut.pdf");
TCanvas *c16 = new TCanvas("c16", "c16", 500, 500);
hEres_lmom->GetYaxis()->SetTitleOffset(1.4);
hEres_lmom->SetLineWidth(2);
hEres_lmom->SetLineColor(kBlue);
hEres_lmom->Fit("gaus");
hEres_lmom->DrawClone();
c16->SaveAs("results/emcal_electrons_Eres_lmom.png");
c16->SaveAs("results/emcal_electrons_Eres_lmom.pdf");
TCanvas *c17 = new TCanvas("c17", "c17", 500, 500);
hEres_hmom->GetYaxis()->SetTitleOffset(1.4);
hEres_hmom->SetLineWidth(2);
hEres_hmom->SetLineColor(kBlue);
hEres_hmom->Fit("gaus");
hEres_hmom->DrawClone();
c17->SaveAs("results/emcal_electrons_Eres_hmom.png");
c17->SaveAs("results/emcal_electrons_Eres_hmom.pdf");
TCanvas *c18 = new TCanvas("c18", "c18", 500, 500);
hEres_cut_lmom->GetYaxis()->SetTitleOffset(1.4);
hEres_cut_lmom->SetLineWidth(2);
hEres_cut_lmom->SetLineColor(kBlue);
hEres_cut_lmom->Fit("gaus");
hEres_cut_lmom->DrawClone();
c18->SaveAs("results/emcal_electrons_Eres_cut_lmom.png");
c18->SaveAs("results/emcal_electrons_Eres_cut_lmom.pdf");
TCanvas *c19 = new TCanvas("c19", "c19", 500, 500);
hEres_cut_hmom->GetYaxis()->SetTitleOffset(1.4);
hEres_cut_hmom->SetLineWidth(2);
hEres_cut_hmom->SetLineColor(kBlue);
hEres_cut_hmom->Fit("gaus");
hEres_cut_hmom->DrawClone();
c19->SaveAs("results/emcal_electrons_Eres_cut_hmom.png");
c19->SaveAs("results/emcal_electrons_Eres_cut_hmom.pdf");
TCanvas *c20 = new TCanvas("c20", "c20", 500, 500);
hPrec_Pthr->GetYaxis()->SetTitleOffset(1.4);
hPrec_Pthr->SetLineWidth(2);
hPrec_Pthr->SetLineColor(kBlue);
hPrec_Pthr->DrawClone("colz");
c20->SaveAs("results/emcal_electrons_Prec_Pthr.png");
c20->SaveAs("results/emcal_electrons_Prec_Pthr.pdf");
TCanvas *c21 = new TCanvas("c21", "c21", 500, 500);
hEres_Ethr->GetYaxis()->SetTitleOffset(1.4);
hEres_Ethr->SetLineWidth(2);
hEres_Ethr->SetLineColor(kBlue);
hEres_Ethr->DrawClone("colz");
c21->SaveAs("results/emcal_electrons_Eres_Ethr.png");
c21->SaveAs("results/emcal_electrons_Eres_Ethr.pdf");
TCanvas *c22 = new TCanvas("c22", "c22", 500, 500);
hEres_leta->GetYaxis()->SetTitleOffset(1.4);
hEres_leta->SetLineWidth(2);
hEres_leta->SetLineColor(kBlue);
hEres_leta->Fit("gaus");
hEres_leta->DrawClone();
c22->SaveAs("results/emcal_electrons_Eres_leta.png");
c22->SaveAs("results/emcal_electrons_Eres_leta.pdf");
TCanvas *c23 = new TCanvas("c23", "c23", 500, 500);
hEres_heta->GetYaxis()->SetTitleOffset(1.4);
hEres_heta->SetLineWidth(2);
hEres_heta->SetLineColor(kBlue);
hEres_heta->Fit("gaus");
hEres_heta->DrawClone();
c23->SaveAs("results/emcal_electrons_Eres_heta.png");
c23->SaveAs("results/emcal_electrons_Eres_heta.pdf");
TCanvas *c24 = new TCanvas("c24", "c24", 500, 500);
hEres_cut_leta->GetYaxis()->SetTitleOffset(1.4);
hEres_cut_leta->SetLineWidth(2);
hEres_cut_leta->SetLineColor(kBlue);
hEres_cut_leta->Fit("gaus");
hEres_cut_leta->DrawClone();
c24->SaveAs("results/emcal_electrons_Eres_cut_leta.png");
c24->SaveAs("results/emcal_electrons_Eres_cut_leta.pdf");
TCanvas *c25 = new TCanvas("c25", "c25", 500, 500);
hEres_cut_heta->GetYaxis()->SetTitleOffset(1.4);
hEres_cut_heta->SetLineWidth(2);
hEres_cut_heta->SetLineColor(kBlue);
hEres_cut_heta->Fit("gaus");
hEres_cut_heta->DrawClone();
c25->SaveAs("results/emcal_electrons_Eres_cut_heta.png");
c25->SaveAs("results/emcal_electrons_Eres_cut_heta.pdf");
TCanvas *c26 = new TCanvas("c26", "c26", 500, 500);
hPrec_Pthr_cut->GetYaxis()->SetTitleOffset(1.4);
hPrec_Pthr_cut->SetLineWidth(2);
hPrec_Pthr_cut->SetLineColor(kBlue);
hPrec_Pthr_cut->DrawClone("colz");
c26->SaveAs("results/emcal_electrons_Prec_Pthr_cut.png");
c26->SaveAs("results/emcal_electrons_Prec_Pthr_cut.pdf");
TCanvas *c27 = new TCanvas("c27", "c27", 500, 500);
hEres_Ethr_cut->GetYaxis()->SetTitleOffset(1.4);
hEres_Ethr_cut->SetLineWidth(2);
hEres_Ethr_cut->SetLineColor(kBlue);
hEres_Ethr_cut->DrawClone("colz");
c27->SaveAs("results/emcal_electrons_Eres_Ethr_cut.png");
c27->SaveAs("results/emcal_electrons_Eres_Ethr_cut.pdf");
}
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