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

Include removed electron scripts

parent cd0f2487
No related branches found
No related tags found
1 merge request!54Resolve "Add information from SciFi layers to benchmarks"
//////////////////////////////////////////////////////////////
// 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;
}
//////////////////////////
// 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");
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment