diff --git a/benchmarks/barrel_ecal/config.yml b/benchmarks/barrel_ecal/config.yml
index 9ab153840d29151bf98e19e48d3b18ea8d45d2bb..209524de86225f9b30996120e9ed5048b2b00ba5 100644
--- a/benchmarks/barrel_ecal/config.yml
+++ b/benchmarks/barrel_ecal/config.yml
@@ -4,6 +4,12 @@ sim:emcal_barrel_pions:
script:
- bash benchmarks/barrel_ecal/run_emcal_barrel_pions.sh
+sim:emcal_barrel_pi0:
+ extends: .det_benchmark
+ stage: simulate
+ script:
+ - bash benchmarks/barrel_ecal/run_emcal_barrel_pi0.sh
+
sim:emcal_barrel_electrons:
extends: .det_benchmark
stage: simulate
@@ -26,6 +32,14 @@ bench:emcal_barrel_pions:
script:
- root -b -q benchmarks/barrel_ecal/scripts/emcal_barrel_pions_analysis.cxx+
+bench:emcal_barrel_pi0:
+ extends: .det_benchmark
+ stage: benchmarks
+ needs:
+ - ["sim:emcal_barrel_pi0"]
+ script:
+ - root -b -q benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_analysis.cxx+
+
bench:emcal_barrel_electrons:
extends: .det_benchmark
stage: benchmarks
@@ -52,7 +66,7 @@ collect_results:barrel_ecal:
extends: .det_benchmark
stage: collect
needs:
- - ["bench:emcal_barrel_electrons", "bench:emcal_barrel_photons"]
+ - ["bench:emcal_barrel_electrons", "bench:emcal_barrel_photons", "bench:emcal_barrel_pions", "bench:emcal_barrel_pi0"]
script:
- ls -lrht
- echo " FIX ME"
diff --git a/benchmarks/barrel_ecal/run_emcal_barrel_pi0.sh b/benchmarks/barrel_ecal/run_emcal_barrel_pi0.sh
new file mode 100755
index 0000000000000000000000000000000000000000..62384d07db67f4da02cd6211697a84bfc2bd94ad
--- /dev/null
+++ b/benchmarks/barrel_ecal/run_emcal_barrel_pi0.sh
@@ -0,0 +1,64 @@
+#!/bin/bash
+
+if [[ ! -n "${JUGGLER_DETECTOR}" ]] ; then
+ export JUGGLER_DETECTOR="topside"
+fi
+
+if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
+ export JUGGLER_N_EVENTS=1000
+fi
+
+if [[ ! -n "${JUGGLER_INSTALL_PREFIX}" ]] ; then
+ export JUGGLER_INSTALL_PREFIX="/usr/local"
+fi
+
+if [[ ! -n "${E_start}" ]] ; then
+ export E_start=5.0
+fi
+
+if [[ ! -n "${E_end}" ]] ; then
+ export E_end=5.0
+fi
+
+export JUGGLER_FILE_NAME_TAG="emcal_barrel_uniform_pi0"
+export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
+
+export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
+export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
+
+echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
+echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
+
+# Generate the input events
+root -b -q "benchmarks/barrel_ecal/scripts/emcal_barrel_pi0.cxx(${JUGGLER_N_EVENTS}, ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running script: generating input events"
+ exit 1
+fi
+# Plot the input events
+root -b -q "benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_reader.cxx(${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running script: plotting input events"
+ exit 1
+fi
+
+# Run geant4 simulations
+npsim --runType batch \
+ -v WARNING \
+ --part.minimalKineticEnergy 0.5*GeV \
+ --numberOfEvents ${JUGGLER_N_EVENTS} \
+ --compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml \
+ --inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \
+ --outputFile sim_output/${JUGGLER_SIM_FILE}
+
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running npdet"
+ exit 1
+fi
+
+# Directory for plots
+mkdir -p results
+
+# Move ROOT output file
+#mv ${JUGGLER_REC_FILE} sim_output/
+
diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_pi0.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_pi0.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..bd5211ab02664cef9a3bca071a9d0fa09c6a14f3
--- /dev/null
+++ b/benchmarks/barrel_ecal/scripts/emcal_barrel_pi0.cxx
@@ -0,0 +1,83 @@
+//////////////////////////////////////////////////////////////
+// EMCAL Barrel detector
+// Single Pi0 dataset
+// M. Scott 05/2021
+//////////////////////////////////////////////////////////////
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/Print.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+
+#include <TMath.h>
+#include <cmath>
+#include <iostream>
+#include <math.h>
+#include <random>
+
+using namespace HepMC3;
+
+void emcal_barrel_pi0(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_barrel_pi0.hepmc") {
+ WriterAscii hepmc_output(out_fname);
+ int events_parsed = 0;
+ GenEvent evt(Units::GEV, Units::MM);
+
+ // 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 * (45.0 / 180.0));
+ double cos_theta_max = std::cos(M_PI * (135.0 / 180.0));
+
+ for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
+ // FourVector(px,py,pz,e,pdgid,status)
+ // type 4 is beam
+ // pdgid 11 - electron
+ // pdgid 111 - pi0
+ // pdgid 2212 - proton
+ 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);
+
+ // Define momentum
+ 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);
+
+ // type 1 is final state
+ // pdgid 211 - pion+ 139.570 MeV/c^2
+ // 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>();
+ v1->add_particle_in(p1);
+ v1->add_particle_in(p2);
+
+ //v1->add_particle_out(p3);
+ v1->add_particle_out(p4);
+ evt.add_vertex(v1);
+
+ if (events_parsed == 0) {
+ std::cout << "First event: " << std::endl;
+ Print::listing(evt);
+ }
+
+ hepmc_output.write_event(evt);
+ if (events_parsed % 10000 == 0) {
+ std::cout << "Event: " << events_parsed << std::endl;
+ }
+ evt.clear();
+ }
+ hepmc_output.close();
+ std::cout << "Events parsed and written: " << events_parsed << std::endl;
+}
diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_analysis.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_analysis.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..e7171383f905a3ef18dc2b55afcac178d0d7e482
--- /dev/null
+++ b/benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_analysis.cxx
@@ -0,0 +1,253 @@
+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/CalorimeterHitCollection.h"
+
+#include "benchmark.h"
+#include "mt.h"
+#include "util.h"
+
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TMath.h"
+#include "TH1.h"
+#include "TF1.h"
+#include "TH1D.h"
+#include "TFitResult.h"
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+void emcal_barrel_pi0_analysis(const char* input_fname = "sim_output/sim_emcal_barrel_uniform_pi0.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);
+ 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::RDataFrame d0("events", input_fname);
+
+ // Sampling Fraction
+ double samp_frac = 0.0136;
+
+ // Thrown Energy [GeV]
+ auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+ std::vector<double> 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;
+ };
+
+ // Number of hits
+ auto nhits = [] (const std::vector<dd4pod::CalorimeterHitData>& evt) {return (int) evt.size(); };
+
+ // Energy deposition [GeV]
+ auto Esim = [](const std::vector<dd4pod::CalorimeterHitData>& evt) {
+ std::vector<double> result;
+ auto total_edep = 0.0;
+ for (const auto& i: evt)
+ total_edep += i.energyDeposit;
+ result.push_back(total_edep);
+ return result;
+ };
+
+ // Sampling fraction = Esampling / Ethrown
+ auto fsam = [](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 / *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;
+ };
+
+ // 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
+ auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
+ .Define("nhits", nhits, {"EcalBarrelHits"})
+ .Define("Esim", Esim, {"EcalBarrelHits"})
+ .Define("fsam", fsam, {"Esim","Ethr"})
+ .Define("pid", getpid, {"mcparticles"})
+ .Define("dau", getdau, {"mcparticles"})
+ ;
+
+ // Define Histograms
+ 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 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 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
+ auto nevents_thrown = d1.Count();
+ std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
+
+ // Draw Histograms
+ TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
+ c1->SetLogy(1);
+ hEthr->GetYaxis()->SetTitleOffset(1.4);
+ hEthr->SetLineWidth(2);
+ hEthr->SetLineColor(kBlue);
+ hEthr->DrawClone();
+ c1->SaveAs("results/emcal_barrel_pi0_Ethr.png");
+ c1->SaveAs("results/emcal_barrel_pi0_Ethr.pdf");
+
+ TCanvas *c2 = new TCanvas("c2", "c2", 700, 500);
+ c2->SetLogy(1);
+ hNhits->GetYaxis()->SetTitleOffset(1.4);
+ hNhits->SetLineWidth(2);
+ hNhits->SetLineColor(kBlue);
+ hNhits->DrawClone();
+ c2->SaveAs("results/emcal_barrel_pi0_nhits.png");
+ c2->SaveAs("results/emcal_barrel_pi0_nhits.pdf");
+
+ TCanvas *c3 = new TCanvas("c3", "c3", 700, 500);
+ c3->SetLogy(1);
+ hEsim->GetYaxis()->SetTitleOffset(1.4);
+ hEsim->SetLineWidth(2);
+ hEsim->SetLineColor(kBlue);
+ hEsim->DrawClone();
+ c3->SaveAs("results/emcal_barrel_pi0_Esim.png");
+ c3->SaveAs("results/emcal_barrel_pi0_Esim.pdf");
+
+ TCanvas *c4 = new TCanvas("c4", "c4", 700, 500);
+ c4->SetLogy(1);
+ hfsam->GetYaxis()->SetTitleOffset(1.4);
+ hfsam->SetLineWidth(2);
+ hfsam->SetLineColor(kBlue);
+ hfsam->Fit("gaus","","",0.005,0.1);
+ hfsam->GetFunction("gaus")->SetLineWidth(2);
+ hfsam->GetFunction("gaus")->SetLineColor(kRed);
+ hfsam->DrawClone();
+ c4->SaveAs("results/emcal_barrel_pi0_fsam.png");
+ c4->SaveAs("results/emcal_barrel_pi0_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_pi0_pid.png");
+ c5->SaveAs("results/emcal_barrel_pi0_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_pi0_dau.png");
+ c6->SaveAs("results/emcal_barrel_pi0_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("{}_pi0.json", detector));
+
+
+}
diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_reader.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_reader.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..597f090fa1927ddd52c0a3b9e509ade997bfc272
--- /dev/null
+++ b/benchmarks/barrel_ecal/scripts/emcal_barrel_pi0_reader.cxx
@@ -0,0 +1,125 @@
+//////////////////////////
+// EMCAL Barrel detector
+// Pi0 dataset
+// M. Scott 05/2021
+//////////////////////////
+#include "HepMC3/GenEvent.h"
+#include "HepMC3/Print.h"
+#include "HepMC3/ReaderAscii.h"
+#include "HepMC3/WriterAscii.h"
+
+#include "TH1F.h"
+#include "TStyle.h"
+#include <iostream>
+
+using namespace HepMC3;
+
+void emcal_barrel_pi0_reader(double e_start = 0.0, double e_end = 30.0, const char* in_fname = "./data/emcal_barrel_pi0.hepmc") {
+ // Setting for graphs
+ gROOT->SetStyle("Plain");
+ gStyle->SetOptFit(1);
+ gStyle->SetLineWidth(1);
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetPadGridX(1);
+ gStyle->SetPadGridY(1);
+ gStyle->SetPadLeftMargin(0.14);
+ gStyle->SetPadRightMargin(0.17);
+
+ ReaderAscii hepmc_input(in_fname);
+ int events_parsed = 0;
+ GenEvent evt(Units::GEV, Units::MM);
+
+ // Histograms
+ TH1F* h_pi0_energy = new TH1F("h_pi0_energy", "pi0 energy;E [GeV];Events", 100, -0.5, 30.5);
+ TH1F* h_pi0_eta = new TH1F("h_pi0_eta", "pi0 #eta;#eta;Events", 100, -10.0, 10.0);
+ TH1F* h_pi0_theta = new TH1F("h_pi0_theta", "pi0 #theta;#theta [degree];Events", 100, -0.5, 180.5);
+ TH1F* h_pi0_phi = new TH1F("h_pi0_phi", "pi0 #phi;#phi [degree];Events", 100, -180.5, 180.5);
+ 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);
+
+ while (!hepmc_input.failed()) {
+ // Read event from input file
+ hepmc_input.read_event(evt);
+ // If reading failed - exit loop
+ if (hepmc_input.failed())
+ break;
+
+ for (const auto& v : evt.vertices()) {
+ for (const auto& p : v->particles_out()) {
+ if (p->pid() == 11) {
+ h_pi0_energy->Fill(p->momentum().e());
+ h_pi0_eta->Fill(p->momentum().eta());
+ h_pi0_theta->Fill(p->momentum().theta() * TMath::RadToDeg());
+ h_pi0_phi->Fill(p->momentum().phi() * TMath::RadToDeg());
+ 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());
+ }
+ }
+ }
+ evt.clear();
+ events_parsed++;
+ }
+ std::cout << "Events parsed and written: " << events_parsed << std::endl;
+
+ TCanvas* c = new TCanvas("c", "c", 500, 500);
+ h_pi0_energy->GetYaxis()->SetTitleOffset(1.8);
+ h_pi0_energy->SetLineWidth(2);
+ h_pi0_energy->SetLineColor(kBlue);
+ h_pi0_energy->DrawClone();
+ c->SaveAs("results/input_emcal_barrel_pi0_energy.png");
+ c->SaveAs("results/input_emcal_barrel_pi0_energy.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_emcal_barrel_pi0_eta.png");
+ c1->SaveAs("results/input_emcal_barrel_pi0_eta.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_emcal_barrel_pi0_theta.png");
+ c2->SaveAs("results/input_emcal_barrel_pi0_theta.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->SetLineColor(kBlue);
+ h_pi0_phi->DrawClone();
+ c3->SaveAs("results/input_emcal_barrel_pi0_phi.png");
+ c3->SaveAs("results/input_emcal_barrel_pi0_phi.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_emcal_barrel_pi0_pzpt.png");
+ c4->SaveAs("results/input_emcal_barrel_pi0_pzpt.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_emcal_barrel_pi0_pxpy.png");
+ c5->SaveAs("results/input_emcal_barrel_pi0_pxpy.pdf");
+
+ TCanvas* c6 = new TCanvas("c6", "c6", 500, 500);
+ h_pi0_p->GetYaxis()->SetTitleOffset(1.8);
+ h_pi0_p->GetXaxis()->SetTitleOffset(1.6);
+ h_pi0_p->GetZaxis()->SetTitleOffset(1.6);
+ h_pi0_p->SetLineWidth(2);
+ h_pi0_p->SetLineColor(kBlue);
+ h_pi0_p->DrawClone();
+ c6->SaveAs("results/input_emcal_barrel_pi0_p.png");
+ c6->SaveAs("results/input_emcal_barrel_pi0_p.pdf");
+}