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 12a78da363d28a2725f022685a8306aa8532ec2f..0000000000000000000000000000000000000000 --- a/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons.cxx +++ /dev/null @@ -1,80 +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 <TMath.h> -#include <cmath> -#include <iostream> -#include <math.h> -#include <random> - -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_analysis.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_analysis.cxx deleted file mode 100644 index 22cb118ce62847bd6ebef8e313537edbd04e3f80..0000000000000000000000000000000000000000 --- a/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_analysis.cxx +++ /dev/null @@ -1,153 +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 "emcal_barrel_common_functions.h" -#include <cctype> - -using ROOT::RDataFrame; -using namespace ROOT::VecOps; - -void emcal_barrel_electrons_analysis(const char* input_fname = "sim_output/sim_emcal_barrel_uniform_electrons.root") -{ - input_fname = "temp_pions_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_fname); - - // 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"); - } - - // Thrown Energy [GeV] - auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) { - auto p = input[2]; - auto energy = TMath::Sqrt(p.psx * p.psx + p.psy * p.psy + p.psz * p.psz + p.mass * p.mass); - return energy; - }; - - // 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; - }; - - // Sampling fraction = Esampling / Ethrown - auto fsam = [](const double sampled, const double thrown) { - return sampled / thrown; - }; - - // Define variables - auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"}) - .Define("nhits", nhits, {"EcalBarrelHits"}) - .Define("Esim", Esim, {"EcalBarrelHits"}) - .Define("fsam", fsam, {"Esim", "Ethr"}); - - // 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", 100, 0.0, 1.0}, - "Esim"); - auto hfsam = d1.Histo1D( - {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1}, - "fsam"); - - addDetectorName(detector_name, hEthr.GetPtr()); - addDetectorName(detector_name, hEsim.GetPtr()); - addDetectorName(detector_name, hfsam.GetPtr()); - - // 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); - auto h = hEthr->DrawCopy(); - //h->GetYaxis()->SetTitleOffset(1.4); - h->SetLineWidth(2); - h->SetLineColor(kBlue); - c1->SaveAs("results/emcal_barrel_electrons_Ethr.png"); - c1->SaveAs("results/emcal_barrel_electrons_Ethr.pdf"); - } - std::cout << "derp1\n"; - - { - TCanvas* c2 = new TCanvas("c2", "c2", 700, 500); - c2->SetLogy(1); - auto h = hNhits->DrawCopy(); - //h->GetYaxis()->SetTitleOffset(1.4); - h->SetLineWidth(2); - h->SetLineColor(kBlue); - c2->SaveAs("results/emcal_barrel_electrons_nhits.png"); - c2->SaveAs("results/emcal_barrel_electrons_nhits.pdf"); - } - - { - TCanvas* c3 = new TCanvas("c3", "c3", 700, 500); - c3->SetLogy(1); - auto h = hEsim->DrawCopy(); - //h->GetYaxis()->SetTitleOffset(1.4); - h->SetLineWidth(2); - h->SetLineColor(kBlue); - c3->SaveAs("results/emcal_barrel_electrons_Esim.png"); - c3->SaveAs("results/emcal_barrel_electrons_Esim.pdf"); - } - - std::cout << "derp4\n"; - { - TCanvas* c4 = new TCanvas("c4", "c4", 700, 500); - c4->SetLogy(1); - auto h = hfsam->DrawCopy(); - //h->GetYaxis()->SetTitleOffset(1.4); - h->SetLineWidth(2); - h->SetLineColor(kBlue); - //h->Fit("gaus", "", "", 0.01, 0.1); - //h->GetFunction("gaus")->SetLineWidth(2); - //h->GetFunction("gaus")->SetLineColor(kRed); - c4->SaveAs("results/emcal_barrel_electrons_fsam.png"); - c4->SaveAs("results/emcal_barrel_electrons_fsam.pdf"); - } -} 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 75c2dae9c9218f036bbca7dd65f65e85bd67eca2..0000000000000000000000000000000000000000 --- a/benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_reader.cxx +++ /dev/null @@ -1,125 +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 "TH1F.h" -#include "TStyle.h" -#include <iostream> - -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"); -} diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_particles_analysis.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_particles_analysis.cxx index 4d52a6e5e8ddd47e0c24747591d529bc50070b3c..8405ff15c41a73838497136968c413427b563beb 100644 --- a/benchmarks/barrel_ecal/scripts/emcal_barrel_particles_analysis.cxx +++ b/benchmarks/barrel_ecal/scripts/emcal_barrel_particles_analysis.cxx @@ -95,7 +95,8 @@ void emcal_barrel_particles_analysis(std::string particle_name = "electron", boo auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"}) .Define("nhits", nhits, {"EcalBarrelHits"}) .Define("Esim", Esim, {"EcalBarrelHits"}) - .Define("fsam", fsam, {"Esim", "Ethr"}); + .Define("EsimScFi", Esim, {"EcalBarrelScFiHits"}) + .Define("fsam", fsam, {"Esim+EsimScFi", "Ethr"}); // Define Histograms auto hEthr = d1.Histo1D(