From 579d824c1c11bda7e9cf757972f416b2073578a0 Mon Sep 17 00:00:00 2001
From: Jihee Kim <jihee.kim@anl.gov>
Date: Mon, 23 Nov 2020 10:02:22 -0600
Subject: [PATCH] Added analysis script. WIP
Updated scripts that generate input particles
---
ecal/ecal_config.yml | 2 +-
ecal/emcal_electrons.sh | 6 +
ecal/scripts/emcal_electrons.cxx | 32 +--
ecal/scripts/emcal_electrons_analysis.cxx | 276 ++++++++++++++++++++++
ecal/scripts/emcal_electrons_reader.cxx | 42 ++--
ecal/scripts/emcal_pi0.cxx | 27 ++-
ecal/scripts/emcal_pi0_reader.cxx | 48 ++--
7 files changed, 367 insertions(+), 66 deletions(-)
create mode 100644 ecal/scripts/emcal_electrons_analysis.cxx
diff --git a/ecal/ecal_config.yml b/ecal/ecal_config.yml
index a949bafa..0409b78c 100644
--- a/ecal/ecal_config.yml
+++ b/ecal/ecal_config.yml
@@ -3,7 +3,7 @@ ecal_1_emcal_electrons:
tags:
- silicon
needs: ["configure"]
- timeout: 24 hours
+ timeout: 48 hours
artifacts:
expire_in: 20 weeks
paths:
diff --git a/ecal/emcal_electrons.sh b/ecal/emcal_electrons.sh
index eac56792..8ebe5391 100644
--- a/ecal/emcal_electrons.sh
+++ b/ecal/emcal_electrons.sh
@@ -98,6 +98,12 @@ if [[ "$?" -ne "0" ]] ; then
exit 1
fi
+root -b -q "ecal/scripts/emcal_electrons_analysis.cxx(\"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running analysis script"
+ exit 1
+fi
+
#paste results/sim_${JUGGLER_FILE_NAME_TAG}.txt results/rec_${JUGGLER_FILE_NAME_TAG}.txt > results/eng_${JUGGLER_FILE_NAME_TAG}.txt
#root -b -q "ecal/scripts/read_eng.C(\"results/eng_${JUGGLER_FILE_NAME_TAG}.root\", \"results/eng_${JUGGLER_FILE_NAME_TAG}.txt\")"
#root -b -q "ecal/scripts/cal_eng_res.C(\"results/eng_${JUGGLER_FILE_NAME_TAG}.root\")"
diff --git a/ecal/scripts/emcal_electrons.cxx b/ecal/scripts/emcal_electrons.cxx
index 7707f24a..675556c2 100644
--- a/ecal/scripts/emcal_electrons.cxx
+++ b/ecal/scripts/emcal_electrons.cxx
@@ -19,12 +19,8 @@
using namespace HepMC3;
-void emcal_electrons(int n_events = 1e2, double e_start = 1.0, double e_end = 1.0,
- const char* out_fname = "./data/emcal_electron_0GeVto30GeV_100kEvt.hepmc")
+void emcal_electrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_electron_0GeVto30GeV_1000kEvt.hepmc")
{
- double cos_theta_min = std::cos(M_PI * (155.0 / 180.0));
- double cos_theta_max = std::cos(M_PI);
-
WriterAscii hepmc_output(out_fname);
int events_parsed = 0;
GenEvent evt(Units::GEV, Units::MM);
@@ -32,6 +28,12 @@ void emcal_electrons(int n_events = 1e2, double e_start = 1.0, double e_end = 1.
// Random number generator
TRandom *r1 = new TRandom();
+ // 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);
+
for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
// FourVector(px,py,pz,e,pdgid,status)
// type 4 is beam
@@ -44,18 +46,18 @@ void emcal_electrons(int n_events = 1e2, double e_start = 1.0, double e_end = 1.
FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
// Define momentum
- Double_t p = r1->Uniform(0.0, 30.0);
- Double_t phi = r1->Uniform(0.0, 2.0 * M_PI);
- Double_t costh = r1->Uniform(cos_theta_min, cos_theta_max);
- Double_t th = std::acos(costh);
- Double_t px = p * std::cos(phi) * std::sin(th);
- Double_t py = p * std::sin(phi) * std::sin(th);
- Double_t pz = p * std::cos(th);
+ // Momentum starting at 0 GeV/c to 30 GeV/c
+ Double_t p = r1->Uniform(e_start, e_end);
+ Double_t phi = r1->Uniform(0.0, 2.0 * M_PI);
+ Double_t costheta = r1->Uniform(cos_theta_min, cos_theta_max);
+ Double_t theta = std::acos(costheta);
+ Double_t px = p * std::cos(phi) * std::sin(theta);
+ Double_t py = p * std::sin(phi) * std::sin(theta);
+ Double_t pz = p * std::cos(theta);
+
// Generates random vectors, uniformly distributed over the surface of a
// sphere of given radius, in this case momentum.
- // r1->Sphere(px, py, pz, p);
-
- //std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n";
+ //r1->Sphere(px, py, pz, p);
// type 1 is final state
// pdgid 11 - electron 0.510 MeV/c^2
diff --git a/ecal/scripts/emcal_electrons_analysis.cxx b/ecal/scripts/emcal_electrons_analysis.cxx
new file mode 100644
index 00000000..ce021827
--- /dev/null
+++ b/ecal/scripts/emcal_electrons_analysis.cxx
@@ -0,0 +1,276 @@
+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterData.h"
+
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TMath.h"
+#include "TH1.h"
+#include "TH1D.h"
+
+// If root cannot find a dd4pod header make sure to add the include dir to ROOT_INCLUDE_PATH:
+// export ROOT_INCLUDE_PATH=$HOME/include:$ROOT_INCLUDE_PATH
+// export ROOT_INCLUDE_PATH=/home/jihee/stow/development/include:$ROOT_INCLUDE_PATH
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+void emcal_electrons_analysis(const char* input_fname = "rec_electron_10kEvt.root")
+{
+ // 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);
+
+ ROOT::EnableImplicitMT();
+ ROOT::RDataFrame d0("events", input_fname);
+
+ // Number of Clusters
+ auto ncluster = [] (const std::vector<eic::ClusterData>& evt) {return (int) evt.size(); };
+
+ // Cluster Energy [GeV]
+ auto clusterE = [] (const std::vector<eic::ClusterData>& evt) {
+ std::vector<float> result;
+ for (const auto& i: evt)
+ result.push_back(i.energy);
+ return result;
+ };
+
+ // Scattering Angle [degree]
+ auto theta = [] (const std::vector<eic::ClusterData>& evt) {
+ std::vector<float> result;
+ 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);
+ result.push_back(TMath::ACos(i.position.z/r)*TMath::RadToDeg());
+ }
+ return result;
+ };
+
+ // Pseudo-rapidity
+ auto eta = [] (const std::vector<eic::ClusterData>& evt) {
+ std::vector<float> result;
+ 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();
+ result.push_back(-1.0 * TMath::Log(TMath::Tan((th*TMath::DegToRad())/2.0)));
+ }
+ return result;
+ };
+
+ // Mass [GeV]
+ auto mass2 = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+ std::vector<float> result;
+ result.push_back(input[2].mass*input[2].mass);
+ return result;
+ };
+
+ // Momentum [GeV]
+ auto p_rec = [](const std::vector<float>& energy_term, const std::vector<float>& mass_term) {
+ std::vector<float> result;
+ for (const auto& e : energy_term) {
+ for (const auto& m2 : mass_term) {
+ result.push_back(TMath::Sqrt((e*e) - m2));
+ }
+ }
+ return result;
+ };
+
+ // Thrown Momentum [GeV]
+ auto p_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+ std::vector<float> result;
+ result.push_back(TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz));
+ return result;
+ };
+
+ // Thrown Energy [GeV]
+ auto E_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+ std::vector<float> result;
+ result.push_back(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));
+ return result;
+ };
+
+ // Energy Resolution
+ auto E_res = [] (const std::vector<float>& Erec, const std::vector<float>& Ethr) {
+ std::vector<float> result;
+ for (const auto& E1 : Ethr) {
+ for (const auto& E2 : Erec) {
+ result.push_back((E2-E1)/E1);
+ }
+ }
+ return result;
+ };
+
+ // Momentum Ratio
+ auto p_ratio = [] (const std::vector<float>& Prec, const std::vector<float>& Pthr) {
+ std::vector<float> result;
+ for (const auto& P1 : Pthr) {
+ for (const auto& P2 : Prec) {
+ result.push_back(P2/P1);
+ }
+ }
+ return result;
+ };
+
+ // Define variables
+ auto d1 = d0.Define("ncluster", ncluster, {"EcalClusters"})
+ .Define("clusterE", clusterE, {"EcalClusters"})
+ .Define("theta", theta, {"EcalClusters"})
+ .Define("eta", eta, {"EcalClusters"})
+ .Define("mass2", mass2, {"mcparticles2"})
+ .Define("p_rec", p_rec, {"clusterE","mass2"})
+ .Define("p_thr", p_thr, {"mcparticles2"})
+ .Define("E_thr", E_thr, {"mcparticles2"})
+ .Define("E_res", E_res, {"clusterE","E_thr"})
+ .Define("p_ratio", p_ratio, {"p_rec","p_thr"})
+ ;
+
+ // Define Histograms
+ auto hNCluster = d1.Histo1D({"hNCluster", "Number of Clusters; # of Clusters; Events", 5, -0.5, 4.5}, "ncluster");
+ auto hClusterE = d1.Histo1D({"hClusterE", "Cluster Energy; Cluster Energy [GeV]; Events", 100, -0.5, 30.5}, "clusterE");
+ auto hTheta = d1.Histo1D({"hTheta", "Scattering Angle; #theta [degree]; Events", 100, 130.0, 180.0}, "theta");
+ auto hEta = d1.Histo1D({"hEta", "Pseudo-rapidity; #eta; Events", 100, -5.0, 0.0}, "eta");
+ auto hPrec = d1.Histo1D({"hPrec", "Reconstructed Momentum; p_{rec}[GeV]; Events", 100, -0.5, 30.5}, "p_rec");
+ auto hPthr = d1.Histo1D({"hPthr", "Thrown Momentum; p_{thr}[GeV]; Events", 100, -0.5, 30.5}, "p_thr");
+ auto hEthr = d1.Histo1D({"hEthr", "Thrown Energy; E_{thr}[GeV]; Events", 100, -0.5, 30.5}, "E_thr");
+
+ // Select Events with One Cluster
+ 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 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) {
+ for (const auto& P2 : Prec) {
+ auto Pratio = P2/P1;
+ if (Pratio > 0.70) {
+ return true;
+ }
+ }
+ }
+ 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");
+
+ // Event Counts
+ auto nevents_thrown = d1.Count();
+ auto nevents_cluster1 = d2.Count();
+
+ std::cout << "Number of Events: " << (*nevents_cluster1) << " / " << (*nevents_thrown) << " = single cluster events/all \n";
+
+ // Draw Histograms
+ TCanvas *c1 = new TCanvas("c1", "c1", 500, 500);
+ c1->SetLogy(1);
+ hNCluster->GetYaxis()->SetTitleOffset(1.6);
+ hNCluster->SetLineWidth(2);
+ hNCluster->SetLineColor(kBlue);
+ hNCluster->DrawClone();
+ c1->SaveAs("results/emcal_electrons_ncluster.png");
+ c1->SaveAs("results/emcal_electrons_ncluster.pdf");
+
+ TCanvas *c2 = new TCanvas("c2", "c2", 500, 500);
+ hClusterE->GetYaxis()->SetTitleOffset(1.4);
+ hClusterE->SetLineWidth(2);
+ hClusterE->SetLineColor(kBlue);
+ hClusterE->DrawClone();
+ c2->SaveAs("results/emcal_electrons_clusterE.png");
+ c2->SaveAs("results/emcal_electrons_clusterE.pdf");
+
+ TCanvas *c3 = new TCanvas("c3", "c3", 500, 500);
+ hTheta->GetYaxis()->SetTitleOffset(1.4);
+ hTheta->SetLineWidth(2);
+ hTheta->SetLineColor(kBlue);
+ hTheta->DrawClone();
+ c3->SaveAs("results/emal_electrons_theta.png");
+ c3->SaveAs("results/emal_electrons_theta.pdf");
+
+ TCanvas *c4 = new TCanvas("c4", "c4", 500, 500);
+ hEta->GetYaxis()->SetTitleOffset(1.4);
+ hEta->SetLineWidth(2);
+ hEta->SetLineColor(kBlue);
+ hEta->DrawClone();
+ c4->SaveAs("results/emcal_electrons_eta.png");
+ c4->SaveAs("results/emcal_electrons_eta.pdf");
+
+ TCanvas *c5 = new TCanvas("c5", "c5", 500, 500);
+ hPrec->GetYaxis()->SetTitleOffset(1.4);
+ hPrec->SetLineWidth(2);
+ hPrec->SetLineColor(kBlue);
+ hPrec->DrawClone();
+ c5->SaveAs("results/emcal_electrons_Prec.png");
+ c5->SaveAs("results/emcal_electrons_Prec.pdf");
+
+ TCanvas *c6 = new TCanvas("c6", "c6", 500, 500);
+ hPthr->GetYaxis()->SetTitleOffset(1.4);
+ hPthr->SetLineWidth(2);
+ hPthr->SetLineColor(kBlue);
+ hPthr->DrawClone();
+ c6->SaveAs("results/emcal_electrons_Pthr.png");
+ c6->SaveAs("results/emcal_electrons_Pthr.pdf");
+
+ TCanvas *c7 = new TCanvas("c7", "c7", 500, 500);
+ hEthr->GetYaxis()->SetTitleOffset(1.4);
+ hEthr->SetLineWidth(2);
+ hEthr->SetLineColor(kBlue);
+ hEthr->DrawClone();
+ c7->SaveAs("results/emcal_electrons_Ethr.png");
+ c7->SaveAs("results/emcal_electrons_Ethr.pdf");
+
+ TCanvas *c8 = new TCanvas("c8", "c8", 500, 500);
+ hClusterE->GetYaxis()->SetTitleOffset(1.4);
+ hClusterE->SetLineWidth(2);
+ hClusterE->SetLineColor(kBlue);
+ hClusterE->DrawClone();
+ c8->SaveAs("results/emcal_electrons_clusterE1.png");
+ c8->SaveAs("results/emcal_electrons_clusterE1.pdf");
+
+ TCanvas *c9 = new TCanvas("c9", "c9", 500, 500);
+ hEres->GetYaxis()->SetTitleOffset(1.4);
+ hEres->SetLineWidth(2);
+ hEres->SetLineColor(kBlue);
+ hEres->Fit("gaus");
+ hEres->DrawClone();
+ c9->SaveAs("results/emcal_electrons_Eres.png");
+ c9->SaveAs("results/emcal_electrons_Eres.pdf");
+
+ TCanvas *c10 = new TCanvas("c10", "c10", 500, 500);
+ hPthr_accepted->GetYaxis()->SetTitleOffset(1.4);
+ hPthr_accepted->SetLineWidth(2);
+ hPthr_accepted->SetLineColor(kBlue);
+ hPthr_accepted->DrawClone();
+ c10->SaveAs("results/emcal_electrons_Pthr_accepted.png");
+ c10->SaveAs("results/emcal_electrons_Pthr_accepted.pdf");
+
+ TCanvas *c11 = new TCanvas("c11", "c11", 500, 500);
+ hPratio->GetYaxis()->SetTitleOffset(1.4);
+ hPratio->SetLineWidth(2);
+ hPratio->SetLineColor(kBlue);
+ hPratio->DrawClone();
+ c11->SaveAs("results/emcal_electrons_Pratio.png");
+ c11->SaveAs("results/emcal_electrons_Pratio.pdf");
+
+ TCanvas *c12 = new TCanvas("c12", "c12", 500, 500);
+ TH1D *hPacceptance = new TH1D("hPacceptance","Ratio p_{rec}/p_{thr}", 100, -0.5, 30.5);
+ hPacceptance = (TH1D*) hPthr_accepted->Clone();
+ hPacceptance->Divide(hPthr.GetPtr());
+ hPacceptance->SetTitle("Ratio p_{rec}/p_{thr}");
+ hPacceptance->SetLineColor(kBlue);
+ hPacceptance->SetLineWidth(2);
+ hPacceptance->GetXaxis()->SetTitle("p [GeV]");
+ hPacceptance->GetYaxis()->SetTitle("p_{rec}/p_{thr}");
+ hPacceptance->GetYaxis()->SetTitleOffset(1.4);
+ hPacceptance->DrawClone();
+ c12->SaveAs("results/emcal_electrons_Pacceptance.png");
+ c12->SaveAs("results/emcal_electrons_Pacceptance.pdf");
+}
diff --git a/ecal/scripts/emcal_electrons_reader.cxx b/ecal/scripts/emcal_electrons_reader.cxx
index f916d23b..5b47d016 100644
--- a/ecal/scripts/emcal_electrons_reader.cxx
+++ b/ecal/scripts/emcal_electrons_reader.cxx
@@ -14,7 +14,7 @@
using namespace HepMC3;
-void emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const char* in_fname = "./data/emcal_electron_0GeVto30GeV_100kEvt.hepmc"){
+void emcal_electrons_reader(double e_start = 0.0, double e_end = 30.0, const char* in_fname = "./data/emcal_electron_1000kEvt.hepmc"){
// Setting for graphs
gROOT->SetStyle("Plain");
@@ -25,20 +25,20 @@ void emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const char
gStyle->SetPadGridX(1);
gStyle->SetPadGridY(1);
gStyle->SetPadLeftMargin(0.14);
- gStyle->SetPadRightMargin(0.14);
+ gStyle->SetPadRightMargin(0.17);
ReaderAscii hepmc_input(in_fname);
int events_parsed = 0;
GenEvent evt(Units::GEV, Units::MM);
// Histograms
- TH1F* h_electrons_energy = new TH1F("h_electron_energy", "electron energy;E [GeV];Events", 100,e_start-0.5,e_end+0.5);
+ TH1F* h_electrons_energy = new TH1F("h_electron_energy", "electron energy;E [GeV];Events", 100,-0.5,30.5);
TH1F* h_electrons_eta = new TH1F("h_electron_eta", "electron #eta;#eta;Events", 100,-10.0,10.0);
TH1F* h_electrons_theta = new TH1F("h_electron_theta", "electron #theta;#theta [degree];Events", 100,-0.5,180.5);
TH1F* h_electrons_phi = new TH1F("h_electron_phi", "electron #phi;#phi [degree];Events", 100,-180.5,180.5);
- TH2F* h_electrons_pzpt = new TH2F("h_electrons_pzpt", "electron pt vs pz;pt [GeV];pz [GeV]", 100,e_start-0.5,e_end+0.5,100,-e_end-0.5,e_end+0.5);
- TH2F* h_electrons_pxpy = new TH2F("h_electrons_pxpy", "electron px vs py;px [GeV];py [GeV]", 100,-e_end-0.5,e_end+0.5,100,-e_end-0.5,e_end+0.5);
- TH3F* h_electrons_p = new TH3F("h_electron_p", "electron p;px [GeV];py [GeV];pz [GeV]", 100,-e_end-0.5,e_end+0.5,100,-e_end-0.5,e_end+0.5,100,-e_end-0.5,e_end+0.5);
+ TH2F* h_electrons_pzpt = new TH2F("h_electrons_pzpt", "electron pt vs pz;pt [GeV];pz [GeV]", 100,-0.5,30.5,100,-30.5,30.5);
+ TH2F* h_electrons_pxpy = new TH2F("h_electrons_pxpy", "electron px vs py;px [GeV];py [GeV]", 100,-30.5,30.5,100,-30.5,30.5);
+ TH3F* h_electrons_p = new TH3F("h_electron_p", "electron p;px [GeV];py [GeV];pz [GeV]", 100,-30.5,30.5,100,-30.5,30.5,100,-30.5,30.5);
while(!hepmc_input.failed()) {
// Read event from input file
@@ -69,49 +69,49 @@ void emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const char
h_electrons_energy->SetLineWidth(2);
h_electrons_energy->SetLineColor(kBlue);
h_electrons_energy->DrawClone();
- c->SaveAs("results/input_electrons_energy_0GeVto30GeV.png");
- c->SaveAs("results/input_electrons_energy_0GeVto30GeV.pdf");
+ c->SaveAs("results/input_emcal_electrons_energy.png");
+ c->SaveAs("results/input_emcal_electrons_energy.pdf");
TCanvas* c1 = new TCanvas("c1","c1", 500, 500);
h_electrons_eta->GetYaxis()->SetTitleOffset(1.9);
h_electrons_eta->SetLineWidth(2);
h_electrons_eta->SetLineColor(kBlue);
h_electrons_eta->DrawClone();
- c1->SaveAs("results/input_electrons_eta_0GeVto30GeV.png");
- c1->SaveAs("results/input_electrons_eta_0GeVto30GeV.pdf");
+ c1->SaveAs("results/input_emcal_electrons_eta.png");
+ c1->SaveAs("results/input_emcal_electrons_eta.pdf");
TCanvas* c2 = new TCanvas("c2","c2", 500, 500);
h_electrons_theta->GetYaxis()->SetTitleOffset(1.8);
h_electrons_theta->SetLineWidth(2);
h_electrons_theta->SetLineColor(kBlue);
h_electrons_theta->DrawClone();
- c2->SaveAs("results/input_electrons_theta_0GeVto30GeV.png");
- c2->SaveAs("results/input_electrons_theta_0GeVto30GeV.pdf");
+ c2->SaveAs("results/input_emcal_electrons_theta.png");
+ c2->SaveAs("results/input_emcal_electrons_theta.pdf");
TCanvas* c3 = new TCanvas("c3","c3", 500, 500);
h_electrons_phi->GetYaxis()->SetTitleOffset(1.8);
h_electrons_phi->SetLineWidth(2);
- h_electrons_phi->GetYaxis()->SetRangeUser(0.0,h_electrons_phi->GetMaximum()+100.0);
+ h_electrons_phi->GetYaxis()->SetRangeUser(0.0,h_electrons_phi->GetMaximum()+300.0);
h_electrons_phi->SetLineColor(kBlue);
h_electrons_phi->DrawClone();
- c3->SaveAs("results/input_electrons_phi_0GeVto30GeV.png");
- c3->SaveAs("results/input_electrons_phi_0GeVto30GeV.pdf");
+ c3->SaveAs("results/input_emcal_electrons_phi.png");
+ c3->SaveAs("results/input_emcal_electrons_phi.pdf");
TCanvas* c4 = new TCanvas("c4","c4", 500, 500);
h_electrons_pzpt->GetYaxis()->SetTitleOffset(1.4);
h_electrons_pzpt->SetLineWidth(2);
h_electrons_pzpt->SetLineColor(kBlue);
h_electrons_pzpt->DrawClone("COLZ");
- c4->SaveAs("results/input_electrons_pzpt_0GeVto30GeV.png");
- c4->SaveAs("results/input_electrons_pzpt_0GeVto30GeV.pdf");
+ c4->SaveAs("results/input_emcal_electrons_pzpt.png");
+ c4->SaveAs("results/input_emcal_electrons_pzpt.pdf");
TCanvas* c5 = new TCanvas("c5","c5", 500, 500);
h_electrons_pxpy->GetYaxis()->SetTitleOffset(1.4);
h_electrons_pxpy->SetLineWidth(2);
h_electrons_pxpy->SetLineColor(kBlue);
h_electrons_pxpy->DrawClone("COLZ");
- c5->SaveAs("results/input_electrons_pxpy_0GeVto30GeV.png");
- c5->SaveAs("results/input_electrons_pxpy_0GeVto30GeV.pdf");
+ c5->SaveAs("results/input_emcal_electrons_pxpy.png");
+ c5->SaveAs("results/input_emcal_electrons_pxpy.pdf");
TCanvas* c6 = new TCanvas("c6","c6", 500, 500);
h_electrons_p->GetYaxis()->SetTitleOffset(1.8);
@@ -120,7 +120,7 @@ void emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const char
h_electrons_p->SetLineWidth(2);
h_electrons_p->SetLineColor(kBlue);
h_electrons_p->DrawClone();
- c6->SaveAs("results/input_electrons_p_0GeVto30GeV.png");
- c6->SaveAs("results/input_electrons_p_0GeVto30GeV.pdf");
+ c6->SaveAs("results/input_emcal_electrons_p.png");
+ c6->SaveAs("results/input_emcal_electrons_p.pdf");
}
diff --git a/ecal/scripts/emcal_pi0.cxx b/ecal/scripts/emcal_pi0.cxx
index b9a70d1e..69018e78 100644
--- a/ecal/scripts/emcal_pi0.cxx
+++ b/ecal/scripts/emcal_pi0.cxx
@@ -16,11 +16,8 @@
using namespace HepMC3;
-void emcal_pi0(int n_events = 1e6, const char* out_fname = "./data/emcal_pi0_0GeVto30GeV_100kEvt.hepmc")
+void emcal_pi0(int n_events = 1e6, const char* out_fname = "./data/emcal_pi0_0GeVto30GeV_1000kEvt.hepmc")
{
- double cos_theta_min = std::cos(M_PI*(120.0/180.0));
- double cos_theta_max = std::cos(M_PI);
-
WriterAscii hepmc_output(out_fname);
int events_parsed = 0;
GenEvent evt(Units::GEV, Units::MM);
@@ -28,6 +25,12 @@ void emcal_pi0(int n_events = 1e6, const char* out_fname = "./data/emcal_pi0_0Ge
// Random number generator
TRandom *r1 = new TRandom();
+ // 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 * (120.0 / 180.0));
+ double cos_theta_max = std::cos(M_PI);
+
for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
// FourVector(px,py,pz,e,pdgid,status)
// type 4 is beam
@@ -40,13 +43,15 @@ void emcal_pi0(int n_events = 1e6, const char* out_fname = "./data/emcal_pi0_0Ge
FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
// Define momentum
- Double_t p = r1->Uniform(0.0, 30.0);
- Double_t phi = r1->Uniform(0.0,2.0*M_PI);
- Double_t costh = r1->Uniform(cos_theta_min,cos_theta_max);
- Double_t th = std::acos(costh);
- Double_t px = p * std::cos(phi)*std::sin(th) ;
- Double_t py = p * std::sin(phi)*std::sin(th) ;
- Double_t pz = p * std::cos(th) ;
+ // Momentum starting at 0 GeV/c to 30 GeV/c
+ Double_t p = r1->Uniform(0.0, 30.0);
+ Double_t phi = r1->Uniform(0.0, 2.0 * M_PI);
+ Double_t costheta = r1->Uniform(cos_theta_min, cos_theta_max);
+ Double_t theta = std::acos(costheta);
+ Double_t px = p * std::cos(phi) * std::sin(theta);
+ Double_t py = p * std::sin(phi) * std::sin(theta);
+ Double_t pz = p * std::cos(theta);
+
// Generates random vectors, uniformly distributed over the surface of a
// sphere of given radius, in this case momentum.
//r1->Sphere(px, py, pz, p);
diff --git a/ecal/scripts/emcal_pi0_reader.cxx b/ecal/scripts/emcal_pi0_reader.cxx
index c43958ac..95f631df 100644
--- a/ecal/scripts/emcal_pi0_reader.cxx
+++ b/ecal/scripts/emcal_pi0_reader.cxx
@@ -14,7 +14,7 @@
using namespace HepMC3;
-void emcal_pi0_reader(const char* in_fname = "./data/emcal_electron_0GeVto30GeV_100kEvt.hepmc"){
+void emcal_pi0_reader(const char* in_fname = "./data/emcal_pi0_0GeVto30GeV_1000kEvt.hepmc"){
// Setting for graphs
gROOT->SetStyle("Plain");
@@ -25,7 +25,7 @@ void emcal_pi0_reader(const char* in_fname = "./data/emcal_electron_0GeVto30GeV_
gStyle->SetPadGridX(1);
gStyle->SetPadGridY(1);
gStyle->SetPadLeftMargin(0.14);
- gStyle->SetPadRightMargin(0.14);
+ gStyle->SetPadRightMargin(0.17);
ReaderAscii hepmc_input(in_fname);
int events_parsed = 0;
@@ -39,7 +39,8 @@ void emcal_pi0_reader(const char* in_fname = "./data/emcal_electron_0GeVto30GeV_
TH2F* h_pi0_pzpt = new TH2F("h_pi0_pzpt", "pi0 pt vs pz;pt [GeV];pz [GeV]", 100,-0.5,30.5,100,-30.5,30.5);
TH2F* h_pi0_pxpy = new TH2F("h_pi0_pxpy", "pi0 px vs py;px [GeV];py [GeV]", 100,-30.5,30.5,100,-30.5,30.5);
TH3F* h_pi0_p = new TH3F("h_pi0_p", "pi0 p;px [GeV];py [GeV];pz [GeV]", 100,-30.5,30.5,100,-30.5,30.5,100,-30.5,30.5);
-
+ TH1F* h_pi0_mom = new TH1F("h_pi0_mom", "pi0 momentum;p [GeV];Events", 100,-0.5,30.5);
+
while(!hepmc_input.failed()) {
// Read event from input file
hepmc_input.read_event(evt);
@@ -56,6 +57,9 @@ void emcal_pi0_reader(const char* in_fname = "./data/emcal_electron_0GeVto30GeV_
h_pi0_pzpt->Fill(TMath::Sqrt(p->momentum().px()*p->momentum().px() + p->momentum().py()*p->momentum().py()),p->momentum().pz());
h_pi0_pxpy->Fill(p->momentum().px(),p->momentum().py());
h_pi0_p->Fill(p->momentum().px(),p->momentum().py(),p->momentum().pz());
+ h_pi0_mom->Fill(TMath::Sqrt(p->momentum().px()*p->momentum().px() +
+ p->momentum().py()*p->momentum().py() +
+ p->momentum().pz()*p->momentum().pz()));
}
}
}
@@ -69,49 +73,49 @@ void emcal_pi0_reader(const char* in_fname = "./data/emcal_electron_0GeVto30GeV_
h_pi0_energy->SetLineWidth(2);
h_pi0_energy->SetLineColor(kBlue);
h_pi0_energy->DrawClone();
- c->SaveAs("results/input_pi0_energy_0GeVto30GeV.png");
- c->SaveAs("results/input_pi0_energy_0GeVto30GeV.pdf");
+ c->SaveAs("results/input_emcal_pi0_energy_0GeVto30GeV.png");
+ c->SaveAs("results/input_emcal_pi0_energy_0GeVto30GeV.pdf");
TCanvas* c1 = new TCanvas("c1","c1", 500, 500);
h_pi0_eta->GetYaxis()->SetTitleOffset(1.9);
h_pi0_eta->SetLineWidth(2);
h_pi0_eta->SetLineColor(kBlue);
h_pi0_eta->DrawClone();
- c1->SaveAs("results/input_pi0_eta_0GeVto30GeV.png");
- c1->SaveAs("results/input_pi0_eta_0GeVto30GeV.pdf");
+ c1->SaveAs("results/input_emcal_pi0_eta_0GeVto30GeV.png");
+ c1->SaveAs("results/input_emcal_pi0_eta_0GeVto30GeV.pdf");
TCanvas* c2 = new TCanvas("c2","c2", 500, 500);
h_pi0_theta->GetYaxis()->SetTitleOffset(1.8);
h_pi0_theta->SetLineWidth(2);
h_pi0_theta->SetLineColor(kBlue);
h_pi0_theta->DrawClone();
- c2->SaveAs("results/input_pi0_theta_0GeVto30GeV.png");
- c2->SaveAs("results/input_pi0_theta_0GeVto30GeV.pdf");
+ c2->SaveAs("results/input_emcal_pi0_theta_0GeVto30GeV.png");
+ c2->SaveAs("results/input_emcal_pi0_theta_0GeVto30GeV.pdf");
TCanvas* c3 = new TCanvas("c3","c3", 500, 500);
h_pi0_phi->GetYaxis()->SetTitleOffset(1.8);
h_pi0_phi->SetLineWidth(2);
- h_pi0_phi->GetYaxis()->SetRangeUser(0.0,h_pi0_phi->GetMaximum()+100.0);
+ h_pi0_phi->GetYaxis()->SetRangeUser(0.0,h_pi0_phi->GetMaximum()+300.0);
h_pi0_phi->SetLineColor(kBlue);
h_pi0_phi->DrawClone();
- c3->SaveAs("results/input_pi0_phi_0GeVto30GeV.png");
- c3->SaveAs("results/input_pi0_phi_0GeVto30GeV.pdf");
+ c3->SaveAs("results/input_emcal_pi0_phi_0GeVto30GeV.png");
+ c3->SaveAs("results/input_emcal_pi0_phi_0GeVto30GeV.pdf");
TCanvas* c4 = new TCanvas("c4","c4", 500, 500);
h_pi0_pzpt->GetYaxis()->SetTitleOffset(1.4);
h_pi0_pzpt->SetLineWidth(2);
h_pi0_pzpt->SetLineColor(kBlue);
h_pi0_pzpt->DrawClone("COLZ");
- c4->SaveAs("results/input_pi0_pzpt_0GeVto30GeV.png");
- c4->SaveAs("results/input_pi0_pzpt_0GeVto30GeV.pdf");
+ c4->SaveAs("results/input_emcal_pi0_pzpt_0GeVto30GeV.png");
+ c4->SaveAs("results/input_emcal_pi0_pzpt_0GeVto30GeV.pdf");
TCanvas* c5 = new TCanvas("c5","c5", 500, 500);
h_pi0_pxpy->GetYaxis()->SetTitleOffset(1.4);
h_pi0_pxpy->SetLineWidth(2);
h_pi0_pxpy->SetLineColor(kBlue);
h_pi0_pxpy->DrawClone("COLZ");
- c5->SaveAs("results/input_pi0_pxpy_0GeVto30GeV.png");
- c5->SaveAs("results/input_pi0_pxpy_0GeVto30GeV.pdf");
+ c5->SaveAs("results/input_emcal_pi0_pxpy_0GeVto30GeV.png");
+ c5->SaveAs("results/input_emcal_pi0_pxpy_0GeVto30GeV.pdf");
TCanvas* c6 = new TCanvas("c6","c6", 500, 500);
h_pi0_p->GetYaxis()->SetTitleOffset(1.8);
@@ -120,7 +124,15 @@ void emcal_pi0_reader(const char* in_fname = "./data/emcal_electron_0GeVto30GeV_
h_pi0_p->SetLineWidth(2);
h_pi0_p->SetLineColor(kBlue);
h_pi0_p->DrawClone();
- c6->SaveAs("results/input_pi0_p_0GeVto30GeV.png");
- c6->SaveAs("results/input_pi0_p_0GeVto30GeV.pdf");
+ c6->SaveAs("results/input_emcal_pi0_p_0GeVto30GeV.png");
+ c6->SaveAs("results/input_emcal_pi0_p_0GeVto30GeV.pdf");
+
+ TCanvas* c7 = new TCanvas("c7","c7", 500, 500);
+ h_pi0_mom->GetYaxis()->SetTitleOffset(1.8);
+ h_pi0_mom->SetLineWidth(2);
+ h_pi0_mom->SetLineColor(kBlue);
+ h_pi0_mom->DrawClone();
+ c7->SaveAs("results/input_emcal_pi0_mom_0GeVto30GeV.png");
+ c7->SaveAs("results/input_emcal_pi0_mom_0GeVto30GeV.pdf");
}
--
GitLab