From 3b13cda38a461e226bfc647bf2eb7afc9b2d113b Mon Sep 17 00:00:00 2001
From: Marshall Scott <mbscott@anl.gov>
Date: Wed, 21 Jul 2021 14:50:59 -0500
Subject: [PATCH] removed .cxx analysis.cxx reader.cxx .sh and cleaned up the
config file
---
benchmarks/barrel_ecal/config.yml | 19 --
.../run_emcal_barrel_pions_electrons.sh | 64 ----
.../scripts/emcal_barrel_pions_electrons.cxx | 100 ------
.../emcal_barrel_pions_electrons_analysis.cxx | 308 ------------------
.../emcal_barrel_pions_electrons_reader.cxx | 125 -------
5 files changed, 616 deletions(-)
delete mode 100755 benchmarks/barrel_ecal/run_emcal_barrel_pions_electrons.sh
delete mode 100644 benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons.cxx
delete mode 100644 benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_analysis.cxx
delete mode 100644 benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_reader.cxx
diff --git a/benchmarks/barrel_ecal/config.yml b/benchmarks/barrel_ecal/config.yml
index 074345f9..ff84a74d 100644
--- a/benchmarks/barrel_ecal/config.yml
+++ b/benchmarks/barrel_ecal/config.yml
@@ -25,13 +25,6 @@ sim:emcal_barrel_photons:
- if [[ "$RUN_EXTENDED_BENCHMARK" == "true" ]] ; then bash benchmarks/barrel_ecal/run_emcal_barrel_energy_scan.sh photon ; fi
- 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_electrons.sh
- - bash benchmarks/barrel_ecal/run_emcal_barrel_particles.sh piminus
-
sim:emcal_barrel_pion_rejection:
extends: .det_benchmark
stage: simulate
@@ -88,17 +81,6 @@ bench:emcal_barrel_photons:
- mv sim_output/sim_emcal_barrel_photon.root results/.
- if [[ "$RUN_EXTENDED_BENCHMARK" == "true" ]] ; then root -b -q 'benchmarks/barrel_ecal/scripts/emcal_barrel_energy_scan_analysis.cxx+("photon")' ; fi
-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+
-
bench:emcal_barrel_pion_rejection:
extends: .det_benchmark
stage: benchmarks
@@ -119,7 +101,6 @@ collect_results:barrel_ecal:
- "bench:emcal_barrel_photons"
- "bench:emcal_barrel_pions"
- "bench:emcal_barrel_pi0"
- - "bench:emcal_barrel_pions_electrons"
- "bench:emcal_barrel_pion_rejection"
script:
- ls -lrht
diff --git a/benchmarks/barrel_ecal/run_emcal_barrel_pions_electrons.sh b/benchmarks/barrel_ecal/run_emcal_barrel_pions_electrons.sh
deleted file mode 100755
index aa546843..00000000
--- a/benchmarks/barrel_ecal/run_emcal_barrel_pions_electrons.sh
+++ /dev/null
@@ -1,64 +0,0 @@
-#!/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/
-
diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons.cxx
deleted file mode 100644
index e7c343c1..00000000
--- a/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons.cxx
+++ /dev/null
@@ -1,100 +0,0 @@
-
-//////////////////////////////////////////////////////////////
-// 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;
-}
diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_analysis.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_analysis.cxx
deleted file mode 100644
index e38a5c20..00000000
--- a/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_analysis.cxx
+++ /dev/null
@@ -1,308 +0,0 @@
-////////////////////////////////////////
-// 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"
-#include "emcal_barrel_common_functions.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});
-
- // Environment Variables
- 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");
- }
-
- // Using the detector layers
- 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("EsimImg", Esim, {"EcalBarrelHits"})
- .Define("EsimScFi", Esim, {"EcalBarrelScFiHits"})
- .Define("Esim", "EsimImg+EsimScFi")
- .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");
- addDetectorName(detector_name, hEthr.GetPtr());
- addDetectorName(detector_name, hNhits.GetPtr());
- addDetectorName(detector_name, hfsam.GetPtr());
- addDetectorName(detector_name, hEsim_front.GetPtr());
-
- 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");
- addDetectorName(detector_name, hEsim_ele.GetPtr());
- addDetectorName(detector_name, hEsim_pim_front.GetPtr());
- addDetectorName(detector_name, hEsim_pim.GetPtr());
- addDetectorName(detector_name, hEsim_pim_front.GetPtr());
-
- 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");
- addDetectorName(detector_name, hEsim_ele_front_rej.GetPtr());
- addDetectorName(detector_name, hEsim_pim_front_rej.GetPtr());
-
- 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, 20, 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, 20, 1.0, 7.0}, "EOverP", "Pthr");
- addDetectorName(detector_name, hEpvp_ele.GetPtr());
- addDetectorName(detector_name, hEpvp_pim.GetPtr());
-
-
- TH1D* hElePurity_initial = (TH1D *)hEsim_ele -> Clone();
- hElePurity_initial -> Divide(hEsim.GetPtr());
- hElePurity_initial -> SetTitle("Electron/Pion Rejection");
- addDetectorName(detector_name, hElePurity_initial);
-
- TH1D* hElePurity_ele = (TH1D *)hEsim_ele_front -> Clone();
- hElePurity_ele -> Divide(hEsim_front.GetPtr());
- hElePurity_ele -> SetTitle("Electron/Pion Rejection : Electron");
- addDetectorName(detector_name, hElePurity_ele);
-
- TH1D* hElePurity_pim = (TH1D *)hEsim_pim_front -> Clone();
- hElePurity_pim -> Divide(hEsim_front.GetPtr());
- hElePurity_pim -> SetTitle("Electron/Pion Rejection : Pion Minus");
- addDetectorName(detector_name, hElePurity_pim);
-
- // Rejection
- TH1D* hElePurity_rej = (TH1D *)hEsim_ele_front_rej -> Clone();
- hElePurity_rej -> Add(hEsim_pim_front_rej.GetPtr());
- hElePurity_rej -> SetTitle("Electron/Pion Rejection");
- addDetectorName(detector_name, hElePurity_rej);
-
- TH1D* hElePurity_ele_rej = (TH1D *)hEsim_ele_front_rej -> Clone();
- hElePurity_ele_rej -> Divide(hElePurity_rej);
- hElePurity_ele_rej -> SetTitle("Electron/Pion Rejection : Electron");
- addDetectorName(detector_name, hElePurity_ele_rej);
-
- 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");
- addDetectorName(detector_name, hElePurity_pim_rej);
-
- // 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");
-
-}
diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_reader.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_reader.cxx
deleted file mode 100644
index a5536ad3..00000000
--- a/benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_reader.cxx
+++ /dev/null
@@ -1,125 +0,0 @@
-//////////////////////////////////////////////////////////////
-// 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");
-}
--
GitLab