Skip to content
Snippets Groups Projects
Commit ea61eb14 authored by Dmitry Kalinkin's avatar Dmitry Kalinkin Committed by Dmitry Kalinkin
Browse files

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

parent 4bdc3d8e
No related branches found
No related tags found
1 merge request!127benchmarks/barrel_ecal: remove unused emcal_barrel_electrons (covered by emcal_barrel_particles)
Hide whitespace changes
Inline Side-by-side
benchmarks/barrel_ecal/run_emcal_barrel_electrons.sh deleted 100755 → 0
+ 0
− 60
View file @ 4bdc3d8e
#!/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/
benchmarks/barrel_ecal/scripts/emcal_barrel_electrons.cxx deleted 100644 → 0
+ 0
− 83
View file @ 4bdc3d8e
//////////////////////////////////////////////////////////////
// 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;
}
benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_reader.cxx deleted 100644 → 0
+ 0
− 132
View file @ 4bdc3d8e
//////////////////////////
// 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% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment