diff --git a/ecal/options/crystal_calorimeter_reco.py b/ecal/options/crystal_calorimeter_reco.py
index fa85b808b7f964c69c4b4d0122cf10fdc648224d..f42d78f91eea8d4fae998a060d5bc8b9d3f3bf2e 100644
--- a/ecal/options/crystal_calorimeter_reco.py
+++ b/ecal/options/crystal_calorimeter_reco.py
@@ -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",
diff --git a/ecal/scripts/emcal_electrons.cxx b/ecal/scripts/emcal_electrons.cxx
index 675556c28a6c0bbebc903a1aef87f2a1b4c0cefb..8229652199e7ae6f404aae2d7ba988533da8aad4 100644
--- a/ecal/scripts/emcal_electrons.cxx
+++ b/ecal/scripts/emcal_electrons.cxx
@@ -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)
diff --git a/ecal/scripts/emcal_electrons_analysis.cxx b/ecal/scripts/emcal_electrons_analysis.cxx
index 9698b3286b9d64cd2c1c89ef101f6cfe1dd90ed8..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);
@@ -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");
+
}