diff --git a/benchmarks/barrel_ecal/run_emcal_barrel_electrons.sh b/benchmarks/barrel_ecal/run_emcal_barrel_electrons.sh deleted file mode 100755 index 3b73b2cddbc1ee70df1917a20dde9c3b245c3158..0000000000000000000000000000000000000000 --- a/benchmarks/barrel_ecal/run_emcal_barrel_electrons.sh +++ /dev/null @@ -1,60 +0,0 @@ -#!/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/ - diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons.cxx deleted file mode 100644 index edec1d849d04062b9a4f9cd7ba1447f6d57ff804..0000000000000000000000000000000000000000 --- a/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons.cxx +++ /dev/null @@ -1,83 +0,0 @@ -////////////////////////////////////////////////////////////// -// 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; -} - diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_reader.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_reader.cxx deleted file mode 100644 index a1330970200dc162ca5f1a09cffc9161c084bd89..0000000000000000000000000000000000000000 --- a/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_reader.cxx +++ /dev/null @@ -1,132 +0,0 @@ -////////////////////////// -// 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"); -} -