Skip to content
Snippets Groups Projects
Commit 3b13cda3 authored by Marshall Scott's avatar Marshall Scott
Browse files

removed .cxx analysis.cxx reader.cxx .sh and cleaned up the config file

parent 7d82cd39
Branches
No related tags found
1 merge request!58Resolve "remove pions electrons benchmark"
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 616 deletions
benchmarks/barrel_ecal/config.yml
+ 0
− 19
View file @ 3b13cda3
......@@ -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
......
benchmarks/barrel_ecal/run_emcal_barrel_pions_electrons.sh deleted 100755 → 0
+ 0
− 64
View file @ 7d82cd39
#!/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 deleted 100644 → 0
+ 0
− 100
View file @ 7d82cd39
//////////////////////////////////////////////////////////////
// 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 deleted 100644 → 0
+ 0
− 308
View file @ 7d82cd39
////////////////////////////////////////
// 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");
}
benchmarks/barrel_ecal/scripts/emcal_barrel_pions_electrons_reader.cxx deleted 100644 → 0
+ 0
− 125
View file @ 7d82cd39
//////////////////////////////////////////////////////////////
// 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");
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment