Resolve "Pion CI Workflow for EM Barrel"

.gitlab-ci.yml

@@ -31,25 +31,6 @@ get_data:
     - mkdir -p results
     - mkdir -p sim_output
 
-#generate_config:
-#  stage: ci_gen
-#  tags:
-#    - silicon
-#  script:
-#    - ./bin/gen_ci_config -p track_test_ -i trackers > results/trackers_config.yml
-#    - ./bin/gen_ci_config -p cal_test_   -i calorimeters > results/calorimeters_config.yml
-#    - ./bin/gen_ci_config -p pid_test_   -i pid > results/pid_config.yml
-
-      #cal_test_3_zdc_neutrons_reader:
-      #stage: benchmarks
-      #tags:
-      #- silicon
-      #needs: 
-      #- ["zdc_simulation"]
-      #script:
-      #- root -b -q calorimeters/zdc_neutrons_reader.cxx
-      #allow_failure: true
-
 roman_pot_simu:
   stage: simulate
   needs: 
@@ -81,83 +62,13 @@ roman_pot_eta:
     - root -b -q trackers/roman_pot_hit_eta.cxx+
   allow_failure: true
 
-    #zdc_simulation:
-    #stage: simulate
-    #needs: 
-    #- ["get_data"]
-    #tags:
-    #- silicon
-    #script:
-    #- cp NPDet/src/detectors/calorimeters/compact/elements.xml ./.
-    #- cp NPDet/src/detectors/calorimeters/compact/materials.xml ./.
-    #- bash calorimeters/run_simulation_zdc.sh
-
-    #zdc_benchmark:
-    #stage: benchmarks
-    #tags:
-    #- silicon
-    #needs: 
-    #- ["zdc_simulation"]
-    #dependencies:
-    #- zdc_simulation
-    #script:
-    #- ls -lrth sim_output
-    #- root -b -q calorimeters/simple_checking.cxx+
-    #allow_failure: true
-
-    #zdc_benchmark_info_histogram:
-    #stage: benchmarks
-    #needs: 
-    #- ["zdc_simulation"]
-    #tags:
-    #- silicon
-    #dependencies:
-    #- zdc_simulation
-    #script:
-    #- cp NPDet/src/detectors/calorimeters/compact/elements.xml calorimeters/
-    #- cp NPDet/src/detectors/calorimeters/compact/materials.xml calorimeters/
-    #- root -b -q calorimeters/simple_info_plot_histograms.cxx+
-    #allow_failure: true
-
-crystal_emcal_simulation:
-  stage: simulate
-  needs: 
-    - ["get_data"]
-  tags:
-    - silicon
-  script:
-    - cp NPDet/src/detectors/calorimeters/compact/elements.xml ./.
-    - cp NPDet/src/detectors/calorimeters/compact/materials.xml ./.
-    - bash calorimeters/run_simulation_crystal.sh
-
-crystal_benchmark:
-  stage: benchmarks
-  tags:
-    - silicon
-  needs:
-    - ["crystal_emcal_simulation"]
-  script:
-    - ls -lrth sim_output
-    - root -b -q calorimeters/simple_checking_crystal.cxx+
-  allow_failure: true
-
-crystal_pion_simulation:
-  stage: simulate
-  needs:
-    - ["get_data"]
-  tags:
-    - silicon
-  script:
-    - cd topside  && ls -l
-    - npsim --runType batch --numberOfEvents 100 --compactFile topside.xml --inputFiles  ../data/emcal_electrons.hepmc  --outputFile  ../sim_output/output_emcal_electrons.root
-
 include:
   - local: 'calorimeters/calorimeters_config.yml'
 
 deploy_results:
   stage: deploy
-    #needs:
-    #- ["zdc_benchmark","zdc_benchmark_info_histogram"]
+  needs:
+    - ["ben_emcal_barrel_electrons"]
   tags:
     - silicon
   script:
@@ -177,4 +88,3 @@ pages:
     paths:
       - public
 
-

calorimeters/calorimeters_config.yml

+#####################
+# Simulations
+# - Generate datasets
+# - Run Geant4
+# - Run Juggler
+#####################
+
+sim_emcal_barrel_pions:
+  image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
+  stage: simulate
+  tags:
+    - silicon
+  artifacts:
+    expire_in: 20 weeks
+    paths:
+      - results/
+      - sim_output/
+  script:
+    - bash calorimeters/run_emcal_barrel_pions.sh
+  allow_failure: true
+
 sim_emcal_barrel_electrons:
   image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
   stage: simulate
   tags:
-    - sodium
+    - silicon
   artifacts:
     expire_in: 20 weeks
     paths:
@@ -12,8 +33,57 @@ sim_emcal_barrel_electrons:
     - bash calorimeters/run_emcal_barrel_electrons.sh
   allow_failure: true
 
+crystal_emcal_simulation:
+  stage: simulate
+  needs: 
+    - ["get_data"]
+  tags:
+    - silicon
+  script:
+    - cp NPDet/src/detectors/calorimeters/compact/elements.xml ./.
+    - cp NPDet/src/detectors/calorimeters/compact/materials.xml ./.
+    - bash calorimeters/run_simulation_crystal.sh
+
+crystal_pion_simulation:
+  stage: simulate
+  needs:
+    - ["get_data"]
+  tags:
+    - silicon
+  script:
+    - cd topside  && ls -l
+    - npsim --runType batch --numberOfEvents 100 --compactFile topside.xml --inputFiles  ../data/emcal_electrons.hepmc  --outputFile  ../sim_output/output_emcal_electrons.root
+
+    #zdc_simulation:
+    #stage: simulate
+    #needs: 
+    #- ["get_data"]
+    #tags:
+    #- silicon
+    #script:
+    #- cp NPDet/src/detectors/calorimeters/compact/elements.xml ./.
+    #- cp NPDet/src/detectors/calorimeters/compact/materials.xml ./.
+    #- bash calorimeters/run_simulation_zdc.sh
+
+###################
+# Benchmarks
+###################
+
+ben_emcal_barrel_pions:
+  stage: benchmarks
+  tags:
+    - silicon
+  needs:
+    - ["sim_emcal_barrel_pions"]
+  script:
+    - ls -lrth sim_output
+    - root -b -q calorimeters/scripts/emcal_barrel_pions_analysis.cxx+
+  allow_failure: true
+
 ben_emcal_barrel_electrons:
   stage: benchmarks
+  tags:
+    - silicon
   needs:
     - ["sim_emcal_barrel_electrons"]
   script:
@@ -21,6 +91,17 @@ ben_emcal_barrel_electrons:
     - root -b -q calorimeters/scripts/emcal_barrel_electrons_analysis.cxx+
   allow_failure: true
 
+crystal_benchmark:
+  stage: benchmarks
+  tags:
+    - silicon
+  needs:
+    - ["crystal_emcal_simulation"]
+  script:
+    - ls -lrth sim_output
+    - root -b -q calorimeters/simple_checking_crystal.cxx+
+  allow_failure: true
+
     #cal_test_3_zdc_neutrons_reader:
     #stage: benchmarks
     #tags:
@@ -32,3 +113,30 @@ ben_emcal_barrel_electrons:
       #    - results/
       #allow_failure: true
 
+    #zdc_benchmark:
+    #stage: benchmarks
+    #tags:
+    #- silicon
+    #needs: 
+    #- ["zdc_simulation"]
+    #dependencies:
+    #- zdc_simulation
+    #script:
+    #- ls -lrth sim_output
+    #- root -b -q calorimeters/simple_checking.cxx+
+    #allow_failure: true
+
+    #zdc_benchmark_info_histogram:
+    #stage: benchmarks
+    #needs: 
+    #- ["zdc_simulation"]
+    #tags:
+    #- silicon
+    #dependencies:
+    #- zdc_simulation
+    #script:
+    #- cp NPDet/src/detectors/calorimeters/compact/elements.xml calorimeters/
+    #- cp NPDet/src/detectors/calorimeters/compact/materials.xml calorimeters/
+    #- root -b -q calorimeters/simple_info_plot_histograms.cxx+
+    #allow_failure: true
+

calorimeters/run_emcal_barrel_pions.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"
+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 "calorimeters/scripts/emcal_barrel_pions.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 "calorimeters/scripts/emcal_barrel_pions_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 topside/${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
+
+# Run Juggler
+xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+        gaudirun.py calorimeters/options/emcal_barrel_reco.py
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running juggler"
+  exit 1
+fi
+
+# Directory for plots
+mkdir -p results
+
+# Move ROOT output file
+mv ${JUGGLER_REC_FILE} sim_output/
+

calorimeters/scripts/emcal_barrel_pions.cxx

+//////////////////////////////////////////////////////////////
+// EMCAL Barrel detector
+// Single Pion dataset
+// J.KIM 04/04/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(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_barrel_pions.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
+    GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.139570 * 0.139570))), 211, 1);
+
+    GenVertexPtr v1 = std::make_shared<GenVertex>();
+    v1->add_particle_in(p1);
+    v1->add_particle_in(p2);
+
+    v1->add_particle_out(p3);
+    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;
+}

calorimeters/scripts/emcal_barrel_pions_analysis.cxx

+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/CalorimeterHitCollection.h"
+#include "dd4pod/TrackerHitCollection.h"
+#include "eicd/RawCalorimeterHitCollection.h"
+#include "eicd/RawCalorimeterHitData.h"
+#include "eicd/CalorimeterHitCollection.h"
+#include "eicd/CalorimeterHitData.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterData.h"
+
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TMath.h"
+#include "TH1.h"
+#include "TF1.h"
+#include "TH1D.h"
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+void emcal_barrel_pions_analysis(const char* input_fname = "sim_output/rec_emcal_barrel_uniform_pions.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_fname);
+
+  // 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;
+  };
+
+  // Reconstructed Energy [GeV] in XY merger
+  auto ErecXY = [] (const std::vector<eic::CalorimeterHitData> & evt) {
+    std::vector<double> result;
+    auto total_eng = 0.0;
+    for (const auto& i: evt)
+      total_eng += i.energy;
+    result.push_back(total_eng / 1.e+3);
+    return result;
+  };
+
+  // Reconstructed Energy [GeV] in Z merger
+  auto ErecZ = [] (const std::vector<eic::CalorimeterHitData> & evt) {
+    std::vector<double> result;
+    auto total_eng = 0.0;
+    for (const auto& i: evt)
+      total_eng += i.energy;
+    result.push_back(total_eng / 1.e+3);
+    return result;
+  };
+
+  // Number of Clusters
+  auto ncluster = [] (const std::vector<eic::ClusterData>& evt) {return (int) evt.size(); };
+
+  // Cluster Energy [GeV]
+  auto Ecluster = [] (const std::vector<eic::ClusterData>& evt) {
+    std::vector<double> result;
+    for (const auto& i: evt)
+      result.push_back(i.energy / 1.e+3);
+    return result;
+  };
+
+  // Sampling fraction = Esampling / Ethrown
+  auto fsam = [](const std::vector<double>& sampled, const std::vector<double>& thrown) {
+    std::vector<double> result;
+    for (const auto& E1 : thrown) {
+      for (const auto& E2 : sampled)
+        result.push_back(E2/E1);
+    }
+    return result;
+  };
+
+  // Define variables
+  auto d1 = d0.Define("Ethr",      Ethr,       {"mcparticles2"})
+	      .Define("ErecXY",    ErecXY,     {"RecoEcalBarrelHitsXY"})
+	      .Define("ErecZ",     ErecZ,      {"RecoEcalBarrelHitsZ"})
+	      .Define("ncluster",  ncluster,   {"EcalBarrelClusters"})
+	      .Define("Ecluster",  Ecluster,   {"EcalBarrelClusters"})
+	      .Define("fsam",      fsam,       {"Ecluster","Ethr"})
+	      ;
+
+  // Define Histograms
+  auto hEthr     = d1.Histo1D({"hEthr",     "Thrown Energy; Thrown Energy [GeV]; Events",                            100, -0.5, 10.5}, "Ethr");
+  auto hErecXY   = d1.Histo1D({"hErecXY",   "Reconstructed Energy in XY merger; Reconstructed Energy [GeV]; Events", 100, -0.5, 10.5}, "ErecXY");
+  auto hErecZ    = d1.Histo1D({"hErecZ",    "Reconstructed Energy in Z merger; Reconstructed Energy [GeV]; Events",  100, -0.5, 10.5}, "ErecZ");
+  auto hNCluster = d1.Histo1D({"hNCluster", "Number of Clusters; # of Clusters; Events",                              20, -0.5, 20.5}, "ncluster");
+  auto hEcluster = d1.Histo1D({"hEcluster", "Cluster Energy; Cluster Energy [GeV]; Events",                          100, -0.5, 10.5}, "Ecluster");
+  auto hfsam     = d1.Histo1D({"hfsam",     "Sampling Fraction; Sampling Fraction; Events",                          100,  0.0,  1.0}, "fsam");
+
+  // 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_pions_Ethr.png");
+  c1->SaveAs("results/emcal_pions_Ethr.pdf");
+
+  TCanvas *c2 = new TCanvas("c2", "c2", 700, 500);
+  c2->SetLogy(1);
+  hErecXY->GetYaxis()->SetTitleOffset(1.4);
+  hErecXY->SetLineWidth(2);
+  hErecXY->SetLineColor(kBlue);
+  hErecXY->DrawClone();
+  c2->SaveAs("results/emcal_pions_ErecXY.png");
+  c2->SaveAs("results/emcal_pions_ErecXY.pdf");
+
+  TCanvas *c3 = new TCanvas("c3", "c3", 700, 500);
+  c3->SetLogy(1);
+  hErecZ->GetYaxis()->SetTitleOffset(1.4);
+  hErecZ->SetLineWidth(2);
+  hErecZ->SetLineColor(kBlue);
+  hErecZ->DrawClone();
+  c3->SaveAs("results/emal_pions_ErecZ.png"); 
+  c3->SaveAs("results/emal_pions_ErecZ.pdf");
+
+  TCanvas *c4 = new TCanvas("c4", "c4", 700, 500);
+  c4->SetLogy(1); 
+  hNCluster->GetYaxis()->SetTitleOffset(1.6);
+  hNCluster->SetLineWidth(2);
+  hNCluster->SetLineColor(kBlue);
+  hNCluster->DrawClone();
+  c4->SaveAs("results/emcal_pions_ncluster.png");
+  c4->SaveAs("results/emcal_pions_ncluster.pdf");
+
+  TCanvas *c5 = new TCanvas("c5", "c5", 700, 500); 
+  c5->SetLogy(1);
+  hEcluster->GetYaxis()->SetTitleOffset(1.4);
+  hEcluster->SetLineWidth(2);
+  hEcluster->SetLineColor(kBlue);
+  hEcluster->DrawClone();
+  c5->SaveAs("results/emcal_pions_Ecluster.png");
+  c5->SaveAs("results/emcal_pions_Ecluster.pdf");
+
+  TCanvas *c6 = new TCanvas("c6", "c6", 700, 500);
+  c6->SetLogy(1);
+  hfsam->GetYaxis()->SetTitleOffset(1.4);
+  hfsam->SetLineWidth(2);
+  hfsam->SetLineColor(kBlue);
+  hfsam->Fit("gaus","","",0.1,1.0);
+  hfsam->GetFunction("gaus")->SetLineWidth(2);
+  hfsam->GetFunction("gaus")->SetLineColor(kRed);
+  hfsam->DrawClone();
+  c6->SaveAs("results/emcal_pions_fsam.png");
+  c6->SaveAs("results/emcal_pions_fsam.pdf");
+}

calorimeters/scripts/emcal_barrel_pions_reader.cxx

+//////////////////////////
+// EMCAL Barrel detector
+// Pion dataset
+// J.KIM 04/04/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_reader(double e_start = 0.0, double e_end = 30.0, const char* in_fname = "./data/emcal_barrel_pions.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_energy = new TH1F("h_pions_energy", "pion energy;E [GeV];Events",         100, -0.5, 30.5);
+  TH1F* h_pions_eta    = new TH1F("h_pions_eta",    "pion #eta;#eta;Events",              100, -10.0, 10.0);
+  TH1F* h_pions_theta  = new TH1F("h_pions_theta",  "pion #theta;#theta [degree];Events", 100, -0.5, 180.5);
+  TH1F* h_pions_phi    = new TH1F("h_pions_phi",    "pion #phi;#phi [degree];Events",     100, -180.5, 180.5);
+  TH2F* h_pions_pzpt   = new TH2F("h_pions_pzpt",   "pion pt vs pz;pt [GeV];pz [GeV]",    100, -0.5, 30.5, 100, -30.5, 30.5);
+  TH2F* h_pions_pxpy   = new TH2F("h_pions_pxpy",   "pion px vs py;px [GeV];py [GeV]",    100, -30.5, 30.5, 100, -30.5, 30.5);
+  TH3F* h_pions_p      = new TH3F("h_pions_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_energy->Fill(p->momentum().e());
+          h_pions_eta->Fill(p->momentum().eta());
+          h_pions_theta->Fill(p->momentum().theta() * TMath::RadToDeg());
+          h_pions_phi->Fill(p->momentum().phi() * TMath::RadToDeg());
+          h_pions_pzpt->Fill(TMath::Sqrt(p->momentum().px() * p->momentum().px() + p->momentum().py() * p->momentum().py()), p->momentum().pz());
+          h_pions_pxpy->Fill(p->momentum().px(), p->momentum().py());
+          h_pions_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_energy->GetYaxis()->SetTitleOffset(1.8);
+  h_pions_energy->SetLineWidth(2);
+  h_pions_energy->SetLineColor(kBlue);
+  h_pions_energy->DrawClone();
+  c->SaveAs("results/input_emcal_barrel_pions_energy.png");
+  c->SaveAs("results/input_emcal_barrel_pions_energy.pdf");
+
+  TCanvas* c1 = new TCanvas("c1", "c1", 500, 500);
+  h_pions_eta->GetYaxis()->SetTitleOffset(1.9);
+  h_pions_eta->SetLineWidth(2);
+  h_pions_eta->SetLineColor(kBlue);
+  h_pions_eta->DrawClone();
+  c1->SaveAs("results/input_emcal_barrel_pions_eta.png");
+  c1->SaveAs("results/input_emcal_barrel_pions_eta.pdf");
+
+  TCanvas* c2 = new TCanvas("c2", "c2", 500, 500);
+  h_pions_theta->GetYaxis()->SetTitleOffset(1.8);
+  h_pions_theta->SetLineWidth(2);
+  h_pions_theta->SetLineColor(kBlue);
+  h_pions_theta->DrawClone();
+  c2->SaveAs("results/input_emcal_barrel_pions_theta.png");
+  c2->SaveAs("results/input_emcal_barrel_pions_theta.pdf");
+
+  TCanvas* c3 = new TCanvas("c3", "c3", 500, 500);
+  h_pions_phi->GetYaxis()->SetTitleOffset(1.8);
+  h_pions_phi->SetLineWidth(2);
+  h_pions_phi->GetYaxis()->SetRangeUser(0.0, h_pions_phi->GetMaximum() + 100.0);
+  h_pions_phi->SetLineColor(kBlue);
+  h_pions_phi->DrawClone();
+  c3->SaveAs("results/input_emcal_barrel_pions_phi.png");
+  c3->SaveAs("results/input_emcal_barrel_pions_phi.pdf");
+
+  TCanvas* c4 = new TCanvas("c4", "c4", 500, 500);
+  h_pions_pzpt->GetYaxis()->SetTitleOffset(1.4);
+  h_pions_pzpt->SetLineWidth(2);
+  h_pions_pzpt->SetLineColor(kBlue);
+  h_pions_pzpt->DrawClone("COLZ");
+  c4->SaveAs("results/input_emcal_barrel_pions_pzpt.png");
+  c4->SaveAs("results/input_emcal_barrel_pions_pzpt.pdf");
+
+  TCanvas* c5 = new TCanvas("c5", "c5", 500, 500);
+  h_pions_pxpy->GetYaxis()->SetTitleOffset(1.4);
+  h_pions_pxpy->SetLineWidth(2);
+  h_pions_pxpy->SetLineColor(kBlue);
+  h_pions_pxpy->DrawClone("COLZ");
+  c5->SaveAs("results/input_emcal_barrel_pions_pxpy.png");
+  c5->SaveAs("results/input_emcal_barrel_pions_pxpy.pdf");
+
+  TCanvas* c6 = new TCanvas("c6", "c6", 500, 500);
+  h_pions_p->GetYaxis()->SetTitleOffset(1.8);
+  h_pions_p->GetXaxis()->SetTitleOffset(1.6);
+  h_pions_p->GetZaxis()->SetTitleOffset(1.6);
+  h_pions_p->SetLineWidth(2);
+  h_pions_p->SetLineColor(kBlue);
+  h_pions_p->DrawClone();
+  c6->SaveAs("results/input_emcal_barrel_pions_p.png");
+  c6->SaveAs("results/input_emcal_barrel_pions_p.pdf");
+}