Resolve "Pion rejection plot"

benchmarks/barrel_ecal/config.yml

@@ -25,6 +25,14 @@ sim:emcal_barrel_photons:
     - bash benchmarks/barrel_ecal/run_emcal_barrel_energy_scan.sh photon
     - bash benchmarks/barrel_ecal/run_emcal_barrel_particles.sh photon
 
+sim:emcal_barrel_pions_electrons:
+  extends: .det_benchmark
+  stage: simulate
+  script:
+    #- bash benchmarks/barrel_ecal/run_emcal_barrel_particles.sh electron
+    - bash benchmarks/barrel_ecal/run_emcal_barrel_electrons.sh
+    - bash benchmarks/barrel_ecal/run_emcal_barrel_particles.sh piminus
+
 bench:emcal_barrel_pions:
   extends: .det_benchmark
   stage: benchmarks
@@ -63,11 +71,22 @@ bench:emcal_barrel_photons:
     - root -b -q 'benchmarks/barrel_ecal/scripts/emcal_barrel_particles_analysis.cxx+("photon")'
     - root -b -q 'benchmarks/barrel_ecal/scripts/emcal_barrel_energy_scan_analysis.cxx+("photon")'
 
+bench:emcal_barrel_pions_electrons:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs:
+    - ["sim:emcal_barrel_pions_electrons"]
+  script:
+    - ls -lhtR sim_output/
+    - rootls -t sim_output/sim_emcal_barrel_piminus.root
+    - rootls -t sim_output/sim_emcal_barrel_uniform_electrons.root
+    - root -b -q benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_analysis.cxx+
+
 collect_results:barrel_ecal:
   extends: .det_benchmark
   stage: collect
   needs: 
-    - ["bench:emcal_barrel_electrons", "bench:emcal_barrel_photons", "bench:emcal_barrel_pions", "bench:emcal_barrel_pi0"]
+    - ["bench:emcal_barrel_electrons", "bench:emcal_barrel_photons", "bench:emcal_barrel_pions", "bench:emcal_barrel_pi0", "bench:emcal_barrel_pions_electrons"]
   script:
     - ls -lrht
     - echo " FIX ME" 

benchmarks/barrel_ecal/run_emcal_barrel_electrons.sh

@@ -8,12 +8,12 @@ if [[ ! -n  "${JUGGLER_N_EVENTS}" ]] ; then
   export JUGGLER_N_EVENTS=100
 fi
 
-if [[ ! -n  "${E_start}" ]] ; then
-  export E_start=5.0
+if [[ ! -n  "${E_START}" ]] ; then
+  export E_START=5.0
 fi
 
-if [[ ! -n  "${E_end}" ]] ; then
-  export E_end=5.0
+if [[ ! -n  "${E_END}" ]] ; then
+  export E_END=5.0
 fi
 export JUGGLER_FILE_NAME_TAG="emcal_barrel_uniform_electrons"
 export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
@@ -25,13 +25,13 @@ 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_electrons.cxx(${JUGGLER_N_EVENTS}, ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+root -b -q "benchmarks/barrel_ecal/scripts/emcal_barrel_electrons.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_electrons_reader.cxx(${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+root -b -q "benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_reader.cxx(${E_START}, ${E_END}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
 if [[ "$?" -ne "0" ]] ; then
   echo "ERROR running script: plotting input events"
   exit 1

benchmarks/barrel_ecal/run_emcal_barrel_pions_electrons.sh

+#!/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_pions_electrons"
+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_pions_electrons.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_pions_electrons_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/
+

benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons.cxx

+
+//////////////////////////////////////////////////////////////
+// EMCAL Barrel detector
+// Electron & Pion Minus 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_pions_electrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_barrel_pions_electrons.hepmc") {
+	/*
+	n_events = 1000; 
+	e_start = 5;
+	e_end = 5; 
+	out_fname = "temp_pions_electrons.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 -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);
+
+    p        = r1->Uniform(e_start, e_end);
+    phi      = r1->Uniform(0.0, 2.0 * M_PI);
+    costheta = r1->Uniform(cos_theta_min, cos_theta_max);
+    theta    = std::acos(costheta);
+    px       = p * std::cos(phi) * std::sin(theta);
+    py       = p * std::sin(phi) * std::sin(theta);
+    pz       = p * std::cos(theta);
+    GenParticlePtr p4 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.000511 * 0.000511))), 11, 1);
+
+    GenVertexPtr v1 = std::make_shared<GenVertex>();
+    GenVertexPtr v2 = std::make_shared<GenVertex>();
+    v1->add_particle_in(p1);
+    v1->add_particle_in(p2);
+
+    if (r1 -> Uniform(0,1) <= 0.5) {v1->add_particle_out(p3);}
+    else {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;
+}

benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_analysis.cxx

+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/CalorimeterHitCollection.h"
+
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TMath.h"
+#include "TH1.h"
+#include "TF1.h"
+#include "TH1D.h"
+#include "TFitResult.h"
+#include "TLegend.h"
+#include "TString.h"
+
+R__LOAD_LIBRARY(libfmt.so)
+#include "fmt/core.h"
+#include "DD4hep/Detector.h"
+#include "DDG4/Geant4Data.h"
+#include "DDRec/CellIDPositionConverter.h"
+
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+void emcal_barrel_pions_electrons_analysis(const char* input_fname1 = "sim_output/sim_emcal_barrel_piminus.root", const char* input_fname2 = "sim_output/sim_emcal_barrel_uniform_electrons.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);
+
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame d0("events", {input_fname1, input_fname2});
+
+  //Using the detector layers 
+  std::string detector_path = "";
+  std::string detector_name = "topside";
+  if(std::getenv("DETECTOR_PATH")) {
+    detector_path = std::getenv("DETECTOR_PATH");
+  }
+  if(std::getenv("JUGGLER_DETECTOR")) {
+    detector_name = std::getenv("JUGGLER_DETECTOR");
+  }
+
+  dd4hep::Detector& detector = dd4hep::Detector::getInstance();
+  detector.fromCompact(fmt::format("{}/{}.xml", detector_path,detector_name));
+  //dd4hep::rec::CellIDPositionConverter cellid_converter(detector);
+
+  auto decoder = detector.readout("EcalBarrelHits").idSpec().decoder();
+  fmt::print("{}\n", decoder->fieldDescription());
+  auto layer_index = decoder->index("layer");
+  fmt::print(" layer index is {}.\n", layer_index);
+
+  // Rejection Value [GeV] based upon Energy deposit in the first 4 layers.
+  // Currently constructed from looking at tthe pion and electron Energy deposit plots
+  // should eventually grab a value from a json file
+  // Suggested cut at 0.005
+  double rejectCut = 0.0;
+
+  // Thrown Energy [GeV]
+  auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+    return 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);
+  };
+
+  // Thrown Momentum [GeV]
+  auto Pthr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+    return TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz);
+  };
+
+  // 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) {
+    auto total_edep = 0.0;
+    for (const auto& i: evt)
+      total_edep += i.energyDeposit;
+  return total_edep;
+  };
+
+  // Energy deposititon [GeV] in the first 4 layers
+  auto Esim_front = [=](const std::vector<dd4pod::CalorimeterHitData>& evt) {
+    auto total_edep = 0.0;
+    for (const auto& i: evt) {
+      //fmt::print("cell id {}, layer {}\n",i.cellID, decoder->get(i.cellID, layer_index));
+      if( decoder->get(i.cellID, layer_index) < 4 ){
+        total_edep += i.energyDeposit;
+      }
+    }
+    return total_edep;
+  };
+
+  // Sampling fraction = Esampling / Ethrown
+  auto fsam = [](const double& sampled, const double& thrown) {
+    return sampled / thrown;
+  };
+
+  // E_front / p
+  auto fEp = [](const double& E_front, const double& mom) {
+    return E_front / mom;
+  };
+
+  // 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;
+  };
+
+  // Filter function to get electrons
+  auto is_electron = [](std::vector<dd4pod::Geant4ParticleData> const& input){
+    if (input[2].pdgID == 11){return true;}
+    else {return false;}
+  };
+
+  // Filter function to get just negative pions
+  auto is_piMinus = [](std::vector<dd4pod::Geant4ParticleData> const& input){
+    if (input[2].pdgID == -211){return true;}
+    else {return false;}
+  };
+
+  // Define variables
+  auto d1 = d0.Define("Ethr",       Ethr,       {"mcparticles"})
+              .Define("Pthr",       Pthr,       {"mcparticles"})
+              .Define("nhits",      nhits,      {"EcalBarrelHits"})
+              .Define("Esim",       Esim,       {"EcalBarrelHits"})
+              .Define("fsam",       fsam,       {"Esim","Ethr"})
+              .Define("pid",        getpid,     {"mcparticles"})
+              .Define("dau",        getdau,     {"mcparticles"})
+              .Define("Esim_front", Esim_front, {"EcalBarrelHits"})
+              .Define("EOverP",     fEp,        {"Esim_front", "Pthr"})
+              ;
+
+  // Particle Filters
+  auto d_ele = d1.Filter(is_electron, {"mcparticles"});
+  auto d_pim = d1.Filter(is_piMinus,  {"mcparticles"});
+  
+  // Energy Deposit Filters
+  auto d_ele_rej = d_ele.Filter("Esim_front > " + to_string(rejectCut));
+  auto d_pim_rej = d_pim.Filter("Esim_front > " + to_string(rejectCut));
+
+  // 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",                           10,  0.0,   0.05}, "Esim");
+  auto hfsam       = d1.Histo1D({"hfsam",  "Sampling Fraction; Sampling Fraction; Events",                          100,  0.0,    0.1}, "fsam");
+  auto hEsim_front = d1.Histo1D({"hEsim_front",  "Energy Deposit Front; Energy Deposit [GeV]; Events",               10,  0.0,   0.05}, "Esim_front");
+
+  auto hEsim_ele        = d_ele.Histo1D({"hEsim_ele",        "Energy Deposit Electron; Energy Deposit [GeV]; Events",            10,  0.0,    0.05}, "Esim");
+  auto hEsim_ele_front  = d_ele.Histo1D({"hEsim_ele_front",  "Energy Deposit Front Electron; Energy Deposit [GeV]; Events",      10,  0.0,    0.05}, "Esim_front");
+  auto hEsim_pim        = d_pim.Histo1D({"hEsim_pim",        "Energy Deposit Electron; Energy Deposit [GeV]; Events",            10,  0.0,    0.05}, "Esim");
+  auto hEsim_pim_front  = d_pim.Histo1D({"hEsim_pim_front",  "Energy Deposit Front Pion-; Energy Deposit [GeV]; Events",         10,  0.0,    0.05}, "Esim_front");
+
+  auto hEsim_ele_front_rej  = d_ele_rej.Histo1D({"hEsim_ele_front_rej",  "Energy Deposit Front Electron; Energy Deposit [GeV]; Events",      10,  0.0,    0.05}, "Esim_front");
+  auto hEsim_pim_front_rej  = d_pim_rej.Histo1D({"hEsim_pim_front_rej",  "Energy Deposit Front Pion-; Energy Deposit [GeV]; Events",         10,  0.0,    0.05}, "Esim_front");
+
+  auto hEpvp_ele = d_ele.Histo2D({"hEpvp_ele", "Energy Deposit 1st 4 Layers/P vs P; E/P; P [GeV]", 20, 0.0, 0.01, 10, 1.0, 7.0}, "EOverP", "Pthr");
+  auto hEpvp_pim = d_pim.Histo2D({"hEpvp_pim", "Energy Deposit 1st 4 Layers/P vs P; E/P; P [GeV]", 20, 0.0, 0.01, 10, 1.0, 7.0}, "EOverP", "Pthr");
+
+
+  TH1D* hElePurity_initial = (TH1D *)hEsim_ele -> Clone();
+  hElePurity_initial -> Divide(hEsim.GetPtr());
+  hElePurity_initial -> SetTitle("Electron/Pion Rejection");
+
+  TH1D* hElePurity_ele = (TH1D *)hEsim_ele_front -> Clone();
+  hElePurity_ele -> Divide(hEsim_front.GetPtr());
+  hElePurity_ele -> SetTitle("Electron/Pion Rejection : Electron");
+
+  TH1D* hElePurity_pim = (TH1D *)hEsim_pim_front -> Clone();
+  hElePurity_pim -> Divide(hEsim_front.GetPtr());
+  hElePurity_pim -> SetTitle("Electron/Pion Rejection : Pion Minus");
+
+  // Rejection
+  TH1D* hElePurity_rej = (TH1D *)hEsim_ele_front_rej -> Clone();
+  hElePurity_rej -> Add(hEsim_pim_front_rej.GetPtr());
+  hElePurity_rej -> SetTitle("Electron/Pion Rejection");
+
+  TH1D* hElePurity_ele_rej = (TH1D *)hEsim_ele_front_rej -> Clone();
+  hElePurity_ele_rej -> Divide(hElePurity_rej);
+  hElePurity_ele_rej -> SetTitle("Electron/Pion Rejection : Electron");
+
+  TH1D* hElePurity_pim_rej = (TH1D *)hEsim_pim_front_rej -> Clone();
+  hElePurity_pim_rej -> Divide(hElePurity_rej);
+  hElePurity_pim_rej -> SetTitle("Electron/Pion Rejection : Pi Minus");
+
+  // 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_pions_electrons_Ethr.png");
+  c1->SaveAs("results/emcal_barrel_pions_electrons_Ethr.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_pions_electrons_Esim.png"); 
+  c3->SaveAs("results/emcal_barrel_pions_electrons_Esim.pdf");
+
+  TCanvas *c4 = new TCanvas("c4", "c4", 700, 500);
+  c4->SetLogy(1);
+  hEsim_ele->GetYaxis()->SetTitleOffset(1.4);
+  hEsim_ele->SetLineWidth(2);
+  hEsim_ele->SetLineColor(kBlue);
+  hEsim_ele->DrawClone();
+  c4->SaveAs("results/emcal_barrel_pions_electrons_Esim_ele.png");
+  c4->SaveAs("results/emcal_barrel_pions_electrons_Esim_ele.pdf");
+
+  TCanvas *c5 = new TCanvas("c5", "c5", 700, 500);
+  //c5->SetLogy(1);
+  hElePurity_initial->GetYaxis()->SetTitleOffset(1.4);
+  hElePurity_initial->SetLineWidth(2);
+  hElePurity_initial->SetLineColor(kBlue);
+  hElePurity_initial->DrawClone();
+  c5->SaveAs("results/emcal_barrel_pions_electrons_rejection_initial.png");
+  c5->SaveAs("results/emcal_barrel_pions_electrons_rejection_initial.pdf");
+
+  TCanvas *c6 = new TCanvas("c6", "c6", 700, 500);
+  //c6->SetLogy(1);
+  hElePurity_ele->GetYaxis()->SetTitleOffset(1.4);
+  hElePurity_ele->SetLineWidth(2);
+  hElePurity_ele->SetLineColor(kBlue);
+  hElePurity_ele->DrawClone();
+  hElePurity_ele_rej->SetLineWidth(2);
+  hElePurity_ele_rej->SetLineColor(kRed);
+  hElePurity_ele_rej->SetLineStyle(10);
+  hElePurity_ele_rej->DrawClone("Same");
+  c6->SaveAs("results/emcal_barrel_pions_electrons_rejection_ele.png");
+  c6->SaveAs("results/emcal_barrel_pions_electrons_rejection_ele.pdf");
+
+  auto leg = new TLegend(0.7, 0.8, 0.8, 0.9);
+  leg->AddEntry(hElePurity_initial, "Initial", "l");
+  leg->AddEntry(hElePurity_ele, "Final", "l");
+
+  TCanvas *c7 = new TCanvas("c7", "c7", 700, 500);
+  //c6->SetLogy(1);
+  hElePurity_pim->GetYaxis()->SetTitleOffset(1.4);
+  hElePurity_pim->SetLineWidth(2);
+  hElePurity_pim->SetLineColor(kBlue);
+  hElePurity_pim->DrawClone();
+  hElePurity_pim_rej->SetLineWidth(2);
+  hElePurity_pim_rej->SetLineColor(kRed);
+  hElePurity_pim_rej->SetLineStyle(10);
+  hElePurity_pim_rej->DrawClone("Same");
+  c7->SaveAs("results/emcal_barrel_pions_electrons_rejection_pim.png");
+  c7->SaveAs("results/emcal_barrel_pions_electrons_rejection_pim.pdf");
+
+  TCanvas *c8 = new TCanvas("c8", "c8", 700, 500);
+  //c6->SetLogy(1);
+  hEpvp_ele->GetYaxis()->SetTitleOffset(1.4);
+  hEpvp_ele->DrawClone("COLZ");
+  c8->SaveAs("results/emcal_barrel_pions_electrons_Epvp_ele.png");
+  c8->SaveAs("results/emcal_barrel_pions_electrons_Epvp_ele.pdf");
+
+  TCanvas *c9 = new TCanvas("c9", "c9", 700, 500);
+  //c6->SetLogy(1);
+  hEpvp_pim->GetYaxis()->SetTitleOffset(1.4);
+  hEpvp_pim->DrawClone("COLZ");
+  c9->SaveAs("results/emcal_barrel_pions_electrons_Epvp_pim.png");
+  c9->SaveAs("results/emcal_barrel_pions_electrons_Epvp_pim.pdf");
+
+}

benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_reader.cxx

+//////////////////////////////////////////////////////////////
+// EMCAL Barrel detector
+// Electron & Pion Minus 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_pions_electrons_reader(double e_start = 0.0, double e_end = 30.0, const char* in_fname = "./data/emcal_barrel_uniform_pions_electrons.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_pions_electrons_energy = new TH1F("h_pions_electrons_energy", "pion energy;E [GeV];Events",         100, -0.5, 30.5);
+  TH1F* h_pions_electrons_eta    = new TH1F("h_pions_electrons_eta",    "pion #eta;#eta;Events",              100, -10.0, 10.0);
+  TH1F* h_pions_electrons_theta  = new TH1F("h_pions_electrons_theta",  "pion #theta;#theta [degree];Events", 100, -0.5, 180.5);
+  TH1F* h_pions_electrons_phi    = new TH1F("h_pions_electrons_phi",    "pion #phi;#phi [degree];Events",     100, -180.5, 180.5);
+  TH2F* h_pions_electrons_pzpt   = new TH2F("h_pions_electrons_pzpt",   "pion pt vs pz;pt [GeV];pz [GeV]",    100, -0.5, 30.5, 100, -30.5, 30.5);
+  TH2F* h_pions_electrons_pxpy   = new TH2F("h_pions_electrons_pxpy",   "pion px vs py;px [GeV];py [GeV]",    100, -30.5, 30.5, 100, -30.5, 30.5);
+  TH3F* h_pions_electrons_p      = new TH3F("h_pions_electrons_p",      "pion 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_pions_electrons_energy->Fill(p->momentum().e());
+          h_pions_electrons_eta->Fill(p->momentum().eta());
+          h_pions_electrons_theta->Fill(p->momentum().theta() * TMath::RadToDeg());
+          h_pions_electrons_phi->Fill(p->momentum().phi() * TMath::RadToDeg());
+          h_pions_electrons_pzpt->Fill(TMath::Sqrt(p->momentum().px() * p->momentum().px() + p->momentum().py() * p->momentum().py()), p->momentum().pz());
+          h_pions_electrons_pxpy->Fill(p->momentum().px(), p->momentum().py());
+          h_pions_electrons_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_pions_electrons_energy->GetYaxis()->SetTitleOffset(1.8);
+  h_pions_electrons_energy->SetLineWidth(2);
+  h_pions_electrons_energy->SetLineColor(kBlue);
+  h_pions_electrons_energy->DrawClone();
+  c->SaveAs("results/input_emcal_barrel_pions_electrons_energy.png");
+  c->SaveAs("results/input_emcal_barrel_pions_electrons_energy.pdf");
+
+  TCanvas* c1 = new TCanvas("c1", "c1", 500, 500);
+  h_pions_electrons_eta->GetYaxis()->SetTitleOffset(1.9);
+  h_pions_electrons_eta->SetLineWidth(2);
+  h_pions_electrons_eta->SetLineColor(kBlue);
+  h_pions_electrons_eta->DrawClone();
+  c1->SaveAs("results/input_emcal_barrel_pions_electrons_eta.png");
+  c1->SaveAs("results/input_emcal_barrel_pions_electrons_eta.pdf");
+
+  TCanvas* c2 = new TCanvas("c2", "c2", 500, 500);
+  h_pions_electrons_theta->GetYaxis()->SetTitleOffset(1.8);
+  h_pions_electrons_theta->SetLineWidth(2);
+  h_pions_electrons_theta->SetLineColor(kBlue);
+  h_pions_electrons_theta->DrawClone();
+  c2->SaveAs("results/input_emcal_barrel_pions_electrons_theta.png");
+  c2->SaveAs("results/input_emcal_barrel_pions_electrons_theta.pdf");
+
+  TCanvas* c3 = new TCanvas("c3", "c3", 500, 500);
+  h_pions_electrons_phi->GetYaxis()->SetTitleOffset(1.8);
+  h_pions_electrons_phi->SetLineWidth(2);
+  h_pions_electrons_phi->GetYaxis()->SetRangeUser(0.0, h_pions_electrons_phi->GetMaximum() + 100.0);
+  h_pions_electrons_phi->SetLineColor(kBlue);
+  h_pions_electrons_phi->DrawClone();
+  c3->SaveAs("results/input_emcal_barrel_pions_electrons_phi.png");
+  c3->SaveAs("results/input_emcal_barrel_pions_electrons_phi.pdf");
+
+  TCanvas* c4 = new TCanvas("c4", "c4", 500, 500);
+  h_pions_electrons_pzpt->GetYaxis()->SetTitleOffset(1.4);
+  h_pions_electrons_pzpt->SetLineWidth(2);
+  h_pions_electrons_pzpt->SetLineColor(kBlue);
+  h_pions_electrons_pzpt->DrawClone("COLZ");
+  c4->SaveAs("results/input_emcal_barrel_pions_electrons_pzpt.png");
+  c4->SaveAs("results/input_emcal_barrel_pions_electrons_pzpt.pdf");
+
+  TCanvas* c5 = new TCanvas("c5", "c5", 500, 500);
+  h_pions_electrons_pxpy->GetYaxis()->SetTitleOffset(1.4);
+  h_pions_electrons_pxpy->SetLineWidth(2);
+  h_pions_electrons_pxpy->SetLineColor(kBlue);
+  h_pions_electrons_pxpy->DrawClone("COLZ");
+  c5->SaveAs("results/input_emcal_barrel_pions_electrons_pxpy.png");
+  c5->SaveAs("results/input_emcal_barrel_pions_electrons_pxpy.pdf");
+
+  TCanvas* c6 = new TCanvas("c6", "c6", 500, 500);
+  h_pions_electrons_p->GetYaxis()->SetTitleOffset(1.8);
+  h_pions_electrons_p->GetXaxis()->SetTitleOffset(1.6);
+  h_pions_electrons_p->GetZaxis()->SetTitleOffset(1.6);
+  h_pions_electrons_p->SetLineWidth(2);
+  h_pions_electrons_p->SetLineColor(kBlue);
+  h_pions_electrons_p->DrawClone();
+  c6->SaveAs("results/input_emcal_barrel_pions_electrons_p.png");
+  c6->SaveAs("results/input_emcal_barrel_pions_electrons_p.pdf");
+}