Skip to content
Snippets Groups Projects
Commit 08045def authored by Dmitry Kalinkin's avatar Dmitry Kalinkin Committed by Wouter Deconinck
Browse files

benchmarks/barrel_ecal: remove unused emcal_barrel_electrons (covered by emcal_barrel_particles)

parent 293dd0af
No related tags found
No related merge requests found
This commit is part of merge request !127. Comments created here will be created in the context of that merge request.
#!/bin/bash
if [[ ! -n "${DETECTOR}" ]] ; then
export DETECTOR="topside"
fi
if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
export JUGGLER_N_EVENTS=100
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_electrons"
export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
echo "DETECTOR = ${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\")"
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\")"
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running script: plotting input events"
exit 1
fi
# Run geant4 simulations
ddsim --runType batch \
-v WARNING \
--part.minimalKineticEnergy 0.5*GeV \
--filter.tracker edep0 \
--numberOfEvents ${JUGGLER_N_EVENTS} \
--compactFile ${DETECTOR_PATH}/${DETECTOR_CONFIG}.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/
//////////////////////////////////////////////////////////////
// EMCAL Barrel detector
// Single Electron dataset
// J.KIM 04/02/2021
//////////////////////////////////////////////////////////////
#include "HepMC3/GenEvent.h"
#include "HepMC3/Print.h"
#include "HepMC3/ReaderAscii.h"
#include "HepMC3/WriterAscii.h"
#include <cmath>
#include <iostream>
#include <math.h>
#include <random>
#include "TMath.h"
#include "TRandom.h"
using namespace HepMC3;
void emcal_barrel_electrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_barrel_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 11 - electron 0.510 MeV/c^2
GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.000511 * 0.000511))), 11, 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;
}
//////////////////////////
// EMCAL Barrel detector
// Electron dataset
// J.KIM 04/02/2021
//////////////////////////
#include "HepMC3/GenEvent.h"
#include "HepMC3/Print.h"
#include "HepMC3/ReaderAscii.h"
#include "HepMC3/WriterAscii.h"
#include <iostream>
#include "TROOT.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TMath.h"
#include "TStyle.h"
#include "TCanvas.h"
using namespace HepMC3;
void emcal_barrel_electrons_reader(double e_start = 0.0, double e_end = 30.0, const char* in_fname = "./data/emcal_barrel_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_electrons_energy = new TH1F("h_electron_energy", "electron energy;E [GeV];Events", 100, -0.5, 30.5);
TH1F* h_electrons_eta = new TH1F("h_electron_eta", "electron #eta;#eta;Events", 100, -10.0, 10.0);
TH1F* h_electrons_theta = new TH1F("h_electron_theta", "electron #theta;#theta [degree];Events", 100, -0.5, 180.5);
TH1F* h_electrons_phi = new TH1F("h_electron_phi", "electron #phi;#phi [degree];Events", 100, -180.5, 180.5);
TH2F* h_electrons_pzpt = new TH2F("h_electrons_pzpt", "electron pt vs pz;pt [GeV];pz [GeV]", 100, -0.5, 30.5, 100, -30.5, 30.5);
TH2F* h_electrons_pxpy = new TH2F("h_electrons_pxpy", "electron px vs py;px [GeV];py [GeV]", 100, -30.5, 30.5, 100, -30.5, 30.5);
TH3F* h_electrons_p = new TH3F("h_electron_p", "electron 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_electrons_energy->Fill(p->momentum().e());
h_electrons_eta->Fill(p->momentum().eta());
h_electrons_theta->Fill(p->momentum().theta() * TMath::RadToDeg());
h_electrons_phi->Fill(p->momentum().phi() * TMath::RadToDeg());
h_electrons_pzpt->Fill(TMath::Sqrt(p->momentum().px() * p->momentum().px() + p->momentum().py() * p->momentum().py()), p->momentum().pz());
h_electrons_pxpy->Fill(p->momentum().px(), p->momentum().py());
h_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_electrons_energy->GetYaxis()->SetTitleOffset(1.8);
h_electrons_energy->SetLineWidth(2);
h_electrons_energy->SetLineColor(kBlue);
h_electrons_energy->DrawClone();
c->SaveAs("results/input_emcal_barrel_electrons_energy.png");
c->SaveAs("results/input_emcal_barrel_electrons_energy.pdf");
TCanvas* c1 = new TCanvas("c1", "c1", 500, 500);
h_electrons_eta->GetYaxis()->SetTitleOffset(1.9);
h_electrons_eta->SetLineWidth(2);
h_electrons_eta->SetLineColor(kBlue);
h_electrons_eta->DrawClone();
c1->SaveAs("results/input_emcal_barrel_electrons_eta.png");
c1->SaveAs("results/input_emcal_barrel_electrons_eta.pdf");
TCanvas* c2 = new TCanvas("c2", "c2", 500, 500);
h_electrons_theta->GetYaxis()->SetTitleOffset(1.8);
h_electrons_theta->SetLineWidth(2);
h_electrons_theta->SetLineColor(kBlue);
h_electrons_theta->DrawClone();
c2->SaveAs("results/input_emcal_barrel_electrons_theta.png");
c2->SaveAs("results/input_emcal_barrel_electrons_theta.pdf");
TCanvas* c3 = new TCanvas("c3", "c3", 500, 500);
h_electrons_phi->GetYaxis()->SetTitleOffset(1.8);
h_electrons_phi->SetLineWidth(2);
h_electrons_phi->GetYaxis()->SetRangeUser(0.0, h_electrons_phi->GetMaximum() + 100.0);
h_electrons_phi->SetLineColor(kBlue);
h_electrons_phi->DrawClone();
c3->SaveAs("results/input_emcal_barrel_electrons_phi.png");
c3->SaveAs("results/input_emcal_barrel_electrons_phi.pdf");
TCanvas* c4 = new TCanvas("c4", "c4", 500, 500);
h_electrons_pzpt->GetYaxis()->SetTitleOffset(1.4);
h_electrons_pzpt->SetLineWidth(2);
h_electrons_pzpt->SetLineColor(kBlue);
h_electrons_pzpt->DrawClone("COLZ");
c4->SaveAs("results/input_emcal_barrel_electrons_pzpt.png");
c4->SaveAs("results/input_emcal_barrel_electrons_pzpt.pdf");
TCanvas* c5 = new TCanvas("c5", "c5", 500, 500);
h_electrons_pxpy->GetYaxis()->SetTitleOffset(1.4);
h_electrons_pxpy->SetLineWidth(2);
h_electrons_pxpy->SetLineColor(kBlue);
h_electrons_pxpy->DrawClone("COLZ");
c5->SaveAs("results/input_emcal_barrel_electrons_pxpy.png");
c5->SaveAs("results/input_emcal_barrel_electrons_pxpy.pdf");
TCanvas* c6 = new TCanvas("c6", "c6", 500, 500);
h_electrons_p->GetYaxis()->SetTitleOffset(1.8);
h_electrons_p->GetXaxis()->SetTitleOffset(1.6);
h_electrons_p->GetZaxis()->SetTitleOffset(1.6);
h_electrons_p->SetLineWidth(2);
h_electrons_p->SetLineColor(kBlue);
h_electrons_p->DrawClone();
c6->SaveAs("results/input_emcal_barrel_electrons_p.png");
c6->SaveAs("results/input_emcal_barrel_electrons_p.pdf");
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment