Skip to content
Snippets Groups Projects
Commit cd0f2487 authored by Maria Zurek's avatar Maria Zurek
Browse files

Add Energy from SciFi layers

parent dab6bbab
No related branches found
No related tags found
1 merge request!54Resolve "Add information from SciFi layers to benchmarks"
This commit is part of merge request !54. Comments created here will be created in the context of that merge request.
Hide whitespace changes
Inline Side-by-side
benchmarks/barrel_ecal/scripts/emcal_barrel_electrons.cxx deleted 100644 → 0
+ 0
− 80
View file @ dab6bbab
//////////////////////////////////////////////////////////////
// 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;
}
benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_analysis.cxx deleted 100644 → 0
+ 0
− 153
View file @ dab6bbab
////////////////////////////////////////
// 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");
}
}
benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_reader.cxx deleted 100644 → 0
+ 0
− 125
View file @ dab6bbab
//////////////////////////
// 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");
}
benchmarks/barrel_ecal/scripts/emcal_barrel_particles_analysis.cxx
+ 2
− 1
View file @ cd0f2487
...@@ -95,7 +95,8 @@ void emcal_barrel_particles_analysis(std::string particle_name = "electron", boo ...@@ -95,7 +95,8 @@ void emcal_barrel_particles_analysis(std::string particle_name = "electron", boo
auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"}) auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
.Define("nhits", nhits, {"EcalBarrelHits"}) .Define("nhits", nhits, {"EcalBarrelHits"})
.Define("Esim", Esim, {"EcalBarrelHits"}) .Define("Esim", Esim, {"EcalBarrelHits"})
.Define("fsam", fsam, {"Esim", "Ethr"}); .Define("EsimScFi", Esim, {"EcalBarrelScFiHits"})
.Define("fsam", fsam, {"Esim+EsimScFi", "Ethr"});
// Define Histograms // Define Histograms
auto hEthr = d1.Histo1D( auto hEthr = d1.Histo1D(
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment