diff --git a/ecal/scripts/emcal_electrons_analysis.cxx b/ecal/scripts/emcal_electrons_analysis.cxx
index d70a412f9591663c38415a058a2331a046d1ccdd..f680b81ddd487d2b2e040081586f2e79a088795d 100644
--- a/ecal/scripts/emcal_electrons_analysis.cxx
+++ b/ecal/scripts/emcal_electrons_analysis.cxx
@@ -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);
@@ -151,11 +152,37 @@ void emcal_electrons_analysis(const char* input_fname = "rec_electron_10kEvt.roo
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 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) {
@@ -168,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) {
@@ -180,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();
@@ -328,4 +394,99 @@ void emcal_electrons_analysis(const char* input_fname = "rec_electron_10kEvt.roo
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");
+
}