Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • datasets
  • master
2 results

Target

Select target project
  • EIC/datasets
  • Polakovic/datasets
  • jihee.kim/datasets
  • cpeng/datasets
4 results
Select Git revision
  • datasets
  • master
2 results
Show changes
Commits on Source (18)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 19107 additions and 3164 deletions
data/eic_upsilon.hepmc 0 → 100644
View file @ 86da4437
This diff is collapsed.
data/emcal_electron_0GeVto30GeV_1000kEvt.hepmc 0 → 100644
View file @ 86da4437
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
data/emcal_electron_0GeVto30GeV_100kEvt.hepmc 0 → 100644
View file @ 86da4437
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
data/emcal_electrons.hepmc
View file @ 86da4437
This diff is collapsed.
data/emcal_electrons_upto30GeV_10kevents.hepmc 0 → 100644
View file @ 86da4437
HepMC::Version 3.02.02
HepMC::Asciiv3-START_EVENT_LISTING
data/emcal_pi0_0GeVto30GeV_1000kEvt.hepmc 0 → 100644
View file @ 86da4437
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
data/emcal_pi0_0GeVto30GeV_100kEvt.hepmc 0 → 100644
View file @ 86da4437
File suppressed by a .gitattributes entry or the file's encoding is unsupported.
data/zdc_photons_evt500.hepmc
View file @ 86da4437
Source diff could not be displayed: it is too large. Options to address this: view the blob.
emcal_electrons.cxx
View file @ 86da4437
///////////////////////////
//////////////////////////////////////////////////////////////
// Crystal EMCAL detector
// Electron dataset
// Single Electron dataset
// J.KIM 07/20/2020
///////////////////////////
//
// 10/4/2020
// Changed to have isotropic distribution in momentum sphere
//////////////////////////////////////////////////////////////
#include "HepMC3/GenEvent.h"
#include "HepMC3/ReaderAscii.h"
#include "HepMC3/WriterAscii.h"
......@@ -11,59 +14,77 @@
#include <iostream>
#include<random>
#include<cmath>
#include <math.h>
#include <TMath.h>
using namespace HepMC3;
void emcal_electrons(){
WriterAscii hepmc_output("./data/emcal_electrons.hepmc");
int events_parsed = 0;
GenEvent evt(Units::GEV, Units::MM);
void emcal_electrons(int n_events = 1e6, const char* out_fname = "./data/emcal_electron_0GeVto30GeV_1000kEvt.hepmc")
{
WriterAscii hepmc_output(out_fname);
int events_parsed = 0;
GenEvent evt(Units::GEV, Units::MM);
// Random number pulled from the env variable
std::mt19937 gen( std::stoi(std::getenv("SEED")) );
std::uniform_real_distribution<double> uniform_theta(135.0*TMath::DegToRad(),178.0*TMath::DegToRad()); // 135-178[degree]
std::uniform_real_distribution<double> uniform_phi(0.0,2*TMath::Pi()); // 360[degree]
// Random number generator
TRandom *r1 = new TRandom();
for (events_parsed = 0; events_parsed < 100; events_parsed++) {
// FourVector(px,py,pz,e,pdgid,status)
// type 4 is beam
// pdgid 11 - electron
// pdgid 2212 - proton
GenParticlePtr p1 =
// 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 * (120.0 / 180.0));
double cos_theta_max = std::cos(M_PI);
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);
GenParticlePtr p2 = std::make_shared<GenParticle>(
FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
// Define momentum
// Momentum starting at 0 GeV/c to 30 GeV/c
Double_t p = r1->Uniform(0.0, 30.0);
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);
// Define variables - energy, theta, phi, momentum vectors
double p = 1.0 + events_parsed*0.4; // temp. energy range 1 GeV to 40 GeV
double theta = uniform_theta(gen);
double phi = uniform_phi(gen);
double px = p*sin(theta)*cos(phi);
double py = p*sin(theta)*sin(phi);
double pz = p*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
GenParticlePtr p3 =
std::make_shared<GenParticle>(FourVector(px,py,pz,sqrt((px*px)+(py*py)+(pz*pz)+(0.000511*0.000511))), 11, 1);
// 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);
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);
v1->add_particle_out(p3);
evt.add_vertex(v1);
if (events_parsed == 0) {
std::cout << "First event: " << std::endl;
Print::listing(evt);
}
if (events_parsed == 0) {
std::cout << "First event: " << std::endl;
Print::listing(evt);
}
hepmc_output.write_event(evt);
if (events_parsed % 10 == 0) {
std::cout << "Event: " << events_parsed << std::endl;
}
evt.clear();
}
hepmc_output.close();
std::cout << "Events parsed and written: " << events_parsed << std::endl;
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_electrons_reader.cxx
View file @ 86da4437
......@@ -10,17 +10,35 @@
#include "TH1F.h"
#include <iostream>
#include "TStyle.h"
using namespace HepMC3;
void emcal_electrons_reader(){
ReaderAscii hepmc_input("./data/emcal_electrons.hepmc");
int events_parsed = 0;
GenEvent evt(Units::GEV, Units::MM);
// 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.17);
ReaderAscii hepmc_input("./data/emcal_electron_0GeVto30GeV_1000kEvt.hepmc");
int events_parsed = 0;
GenEvent evt(Units::GEV, Units::MM);
// Histograms
TH1F* h_electrons_energy = new TH1F("electron energy","; E [GeV]",100,0,40);
TH1F* h_electrons_eta = new TH1F("electron #eta","; #eta",20,-10,0);
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
......@@ -31,9 +49,14 @@ void emcal_electrons_reader(){
for(const auto& v : evt.vertices() ) {
for(const auto& p : v->particles_out() ) {
if(p->pid() == 11) {
h_electrons_energy->Fill(p->momentum().e()); // Energy component of momentum
h_electrons_eta->Fill(p->momentum().eta()); // Pseudorapidity
}
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();
......@@ -41,14 +64,63 @@ void emcal_electrons_reader(){
}
std::cout << "Events parsed and written: " << events_parsed << std::endl;
TCanvas* c = new TCanvas();
h_electrons_energy->Draw();
c->SaveAs("results/emcal_electrons_energy_reader.png");
c->SaveAs("results/emcal_electrons_energy_reader.pdf");
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/emcal_electrons_energy_0GeVto30GeV.png");
c->SaveAs("results/emcal_electrons_energy_0GeVto30GeV.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/emcal_electrons_eta_0GeVto30GeV.png");
c1->SaveAs("results/emcal_electrons_eta_0GeVto30GeV.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/emcal_electrons_theta_0GeVto30GeV.png");
c2->SaveAs("results/emcal_electrons_theta_0GeVto30GeV.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()+300.0);
h_electrons_phi->SetLineColor(kBlue);
h_electrons_phi->DrawClone();
c3->SaveAs("results/emcal_electrons_phi_0GeVto30GeV.png");
c3->SaveAs("results/emcal_electrons_phi_0GeVto30GeV.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/emcal_electrons_pzpt_0GeVto30GeV.png");
c4->SaveAs("results/emcal_electrons_pzpt_0GeVto30GeV.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/emcal_electrons_pxpy_0GeVto30GeV.png");
c5->SaveAs("results/emcal_electrons_pxpy_0GeVto30GeV.pdf");
TCanvas* c1 = new TCanvas();
h_electrons_eta->Draw();
c1->SaveAs("results/emcal_electrons_eta_reader.png");
c1->SaveAs("results/emcal_electrons_eta_reader.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/emcal_electrons_p_0GeVto30GeV.png");
c6->SaveAs("results/emcal_electrons_p_0GeVto30GeV.pdf");
}
emcal_pi0.cxx 0 → 100644
View file @ 86da4437
//////////////////////////////////////////////////////////////
// Crystal EMCAL detector
// Single Pion dataset
// J.KIM 10/18/2020
//////////////////////////////////////////////////////////////
#include "HepMC3/GenEvent.h"
#include "HepMC3/ReaderAscii.h"
#include "HepMC3/WriterAscii.h"
#include "HepMC3/Print.h"
#include <iostream>
#include<random>
#include<cmath>
#include <math.h>
#include <TMath.h>
using namespace HepMC3;
void emcal_pi0(int n_events = 1e6, const char* out_fname = "./data/emcal_pi0_0GeVto30GeV_1000kEvt.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 * (120.0 / 180.0));
double cos_theta_max = std::cos(M_PI);
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
// Momentum starting at 0 GeV/c to 30 GeV/c
Double_t p = r1->Uniform(0.0, 30.0);
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 111 - pi0 135 MeV/c^2
GenParticlePtr p3 = std::make_shared<GenParticle>(
FourVector(
px, py, pz,
sqrt(p*p + (0.134976*0.134976))),
111, 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_pi0_reader.cxx 0 → 100644
View file @ 86da4437
//////////////////////////
// Crystal EMCAL detector
// Singe Pion dataset
// J.KIM 10/18/2020
//////////////////////////
#include "HepMC3/GenEvent.h"
#include "HepMC3/ReaderAscii.h"
#include "HepMC3/WriterAscii.h"
#include "HepMC3/Print.h"
#include "TH1F.h"
#include <iostream>
#include "TStyle.h"
using namespace HepMC3;
void emcal_pi0_reader(){
// 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.17);
ReaderAscii hepmc_input("./data/emcal_pi0_0GeVto30GeV_1000kEvt.hepmc");
int events_parsed = 0;
GenEvent evt(Units::GEV, Units::MM);
// Histograms
TH1F* h_pi0_energy = new TH1F("h_pi0_energy", "pi0 energy;E [GeV];Events", 100,-0.5,30.5);
TH1F* h_pi0_eta = new TH1F("h_pi0_eta", "pi0 #eta;#eta;Events", 100,-10.0,10.0);
TH1F* h_pi0_theta = new TH1F("h_pi0_theta", "pi0 #theta;#theta [degree];Events", 100,-0.5,180.5);
TH1F* h_pi0_phi = new TH1F("h_pi0_phi", "pi0 #phi;#phi [degree];Events", 100,-180.5,180.5);
TH2F* h_pi0_pzpt = new TH2F("h_pi0_pzpt", "pi0 pt vs pz;pt [GeV];pz [GeV]", 100,-0.5,30.5,100,-30.5,30.5);
TH2F* h_pi0_pxpy = new TH2F("h_pi0_pxpy", "pi0 px vs py;px [GeV];py [GeV]", 100,-30.5,30.5,100,-30.5,30.5);
TH3F* h_pi0_p = new TH3F("h_pi0_p", "pi0 p;px [GeV];py [GeV];pz [GeV]", 100,-30.5,30.5,100,-30.5,30.5,100,-30.5,30.5);
TH1F* h_pi0_mom = new TH1F("h_pi0_mom", "pi0 momentum;p [GeV];Events", 100,-0.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() == 111) {
h_pi0_energy->Fill(p->momentum().e());
h_pi0_eta->Fill(p->momentum().eta());
h_pi0_theta->Fill(p->momentum().theta()*TMath::RadToDeg());
h_pi0_phi->Fill(p->momentum().phi()*TMath::RadToDeg());
h_pi0_pzpt->Fill(TMath::Sqrt(p->momentum().px()*p->momentum().px() + p->momentum().py()*p->momentum().py()),p->momentum().pz());
h_pi0_pxpy->Fill(p->momentum().px(),p->momentum().py());
h_pi0_p->Fill(p->momentum().px(),p->momentum().py(),p->momentum().pz());
h_pi0_mom->Fill(TMath::Sqrt(p->momentum().px()*p->momentum().px() +
p->momentum().py()*p->momentum().py() +
p->momentum().pz()*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_pi0_energy->GetYaxis()->SetTitleOffset(1.8);
h_pi0_energy->SetLineWidth(2);
h_pi0_energy->SetLineColor(kBlue);
h_pi0_energy->DrawClone();
c->SaveAs("results/emcal_pi0_energy_0GeVto30GeV.png");
c->SaveAs("results/emcal_pi0_energy_0GeVto30GeV.pdf");
TCanvas* c1 = new TCanvas("c1","c1", 500, 500);
h_pi0_eta->GetYaxis()->SetTitleOffset(1.9);
h_pi0_eta->SetLineWidth(2);
h_pi0_eta->SetLineColor(kBlue);
h_pi0_eta->DrawClone();
c1->SaveAs("results/emcal_pi0_eta_0GeVto30GeV.png");
c1->SaveAs("results/emcal_pi0_eta_0GeVto30GeV.pdf");
TCanvas* c2 = new TCanvas("c2","c2", 500, 500);
h_pi0_theta->GetYaxis()->SetTitleOffset(1.8);
h_pi0_theta->SetLineWidth(2);
h_pi0_theta->SetLineColor(kBlue);
h_pi0_theta->DrawClone();
c2->SaveAs("results/emcal_pi0_theta_0GeVto30GeV.png");
c2->SaveAs("results/emcal_pi0_theta_0GeVto30GeV.pdf");
TCanvas* c3 = new TCanvas("c3","c3", 500, 500);
h_pi0_phi->GetYaxis()->SetTitleOffset(1.8);
h_pi0_phi->SetLineWidth(2);
h_pi0_phi->GetYaxis()->SetRangeUser(0.0,h_pi0_phi->GetMaximum()+300.0);
h_pi0_phi->SetLineColor(kBlue);
h_pi0_phi->DrawClone();
c3->SaveAs("results/emcal_pi0_phi_0GeVto30GeV.png");
c3->SaveAs("results/emcal_pi0_phi_0GeVto30GeV.pdf");
TCanvas* c4 = new TCanvas("c4","c4", 500, 500);
h_pi0_pzpt->GetYaxis()->SetTitleOffset(1.4);
h_pi0_pzpt->SetLineWidth(2);
h_pi0_pzpt->SetLineColor(kBlue);
h_pi0_pzpt->DrawClone("COLZ");
c4->SaveAs("results/emcal_pi0_pzpt_0GeVto30GeV.png");
c4->SaveAs("results/emcal_pi0_pzpt_0GeVto30GeV.pdf");
TCanvas* c5 = new TCanvas("c5","c5", 500, 500);
h_pi0_pxpy->GetYaxis()->SetTitleOffset(1.4);
h_pi0_pxpy->SetLineWidth(2);
h_pi0_pxpy->SetLineColor(kBlue);
h_pi0_pxpy->DrawClone("COLZ");
c5->SaveAs("results/emcal_pi0_pxpy_0GeVto30GeV.png");
c5->SaveAs("results/emcal_pi0_pxpy_0GeVto30GeV.pdf");
TCanvas* c6 = new TCanvas("c6","c6", 500, 500);
h_pi0_p->GetYaxis()->SetTitleOffset(1.8);
h_pi0_p->GetXaxis()->SetTitleOffset(1.6);
h_pi0_p->GetZaxis()->SetTitleOffset(1.6);
h_pi0_p->SetLineWidth(2);
h_pi0_p->SetLineColor(kBlue);
h_pi0_p->DrawClone();
c6->SaveAs("results/emcal_pi0_p_0GeVto30GeV.png");
c6->SaveAs("results/emcal_pi0_p_0GeVto30GeV.pdf");
TCanvas* c7 = new TCanvas("c7","c7", 500, 500);
h_pi0_mom->GetYaxis()->SetTitleOffset(1.8);
h_pi0_mom->SetLineWidth(2);
h_pi0_mom->SetLineColor(kBlue);
h_pi0_mom->DrawClone();
c7->SaveAs("results/emcal_pi0_mom_0GeVto30GeV.png");
c7->SaveAs("results/emcal_pi0_mom_0GeVto30GeV.pdf");
}
gen_particles.py 0 → 100644
View file @ 86da4437
import os
from pyHepMC3 import HepMC3 as hm
import numpy as np
import argparse
PARTICLES = {
"pion0": (111, 0.1349766), # pi0
"pion+": (211, 0.13957018), # pi+
"pion-": (-211, 0.13957018), # pi-
"kaon0": (311, 0.497648), # K0
"kaon+": (321, 0.493677), # K+
"kaon-": (-321, 0.493677), # K-
"proton": (2212, 0.938272), # proton
"neutron": (2112, 0.939565), # neutron
"electron": (11, 0.51099895e-3), # electron
"positron": (-11, 0.51099895e-3),# positron
"photon": (22, 0), # photon
}
# p in GeV, angle in degree, vertex in mm
def gen_event(prange=(8, 100), arange=(0, 20), phrange=(0., 360.), parts=[(p, m) for p, m in PARTICLES.values()]):
evt = hm.GenEvent(momentum_unit=hm.Units.MomentumUnit.GEV, length_unit=hm.Units.LengthUnit.MM)
pid, mass = parts[np.random.randint(len(parts))]
# final state
state = 1
# momentum, angles, energy
p = np.random.uniform(*prange)
theta = np.random.uniform(*arange)*np.pi/180.
phi = np.random.uniform(*phrange)*np.pi/180.
e0 = np.sqrt(p*p + mass*mass)
px = np.cos(phi)*np.sin(theta)
py = np.sin(phi)*np.sin(theta)
pz = np.cos(theta)
# beam
pbeam = hm.GenParticle(hm.FourVector(0, 0, 0, 0.938272), 2212, 4)
ebeam = hm.GenParticle(hm.FourVector(0, 0, e0, np.sqrt(e0*e0 + 0.511e-3*0.511e-3)), 11, 4)
hout = hm.GenParticle(hm.FourVector(px*p, py*p, pz*p, e0), pid, state)
# evt.add_particle(part)
vert = hm.GenVertex()
vert.add_particle_in(ebeam)
vert.add_particle_in(pbeam)
vert.add_particle_out(hout)
evt.add_vertex(vert)
return evt
if __name__ == "__main__":
parser = argparse.ArgumentParser('RICH dataset generator')
parser.add_argument('output', help='path to the output file')
parser.add_argument('-n', type=int, default=1000, dest='nev', help='number of events to generate')
parser.add_argument('--parray', type=str, default="", dest='parray', help='an array of momentum in GeV')
parser.add_argument('--pmin', type=float, default=8.0, dest='pmin', help='minimum momentum in GeV')
parser.add_argument('--pmax', type=float, default=100.0, dest='pmax', help='maximum momentum in GeV')
parser.add_argument('--angmin', type=float, default=0.0, dest='angmin', help='minimum angle in degree')
parser.add_argument('--angmax', type=float, default=20.0, dest='angmax', help='maximum angle in degree')
parser.add_argument('--phmin', type=float, default=0.0, dest='phmin', help='minimum angle in degree')
parser.add_argument('--phmax', type=float, default=360.0, dest='phmax', help='maximum angle in degree')
parser.add_argument('--particles', type=str, default='11', dest='particles', help='particles pdg code')
parser.add_argument('-s', type=int, default=-1, dest='seed', help='seed for random generator')
args = parser.parse_args()
# random seed
if args.seed < 0:
args.seed = os.environ.get('SEED', int.from_bytes(os.urandom(4), byteorder='big', signed=False))
np.random.seed(args.seed)
output = hm.WriterAscii(args.output);
if output.failed():
print("Cannot open file \"{}\"".format(args.output))
sys.exit(2)
parts = []
for pid in args.particles.split(','):
pid = pid.strip()
if pid not in PARTICLES.keys():
print('pid {:d} not found in dictionary, ignored.'.format(pid))
continue
parts.append(PARTICLES[pid])
if not args.parray:
count = 0
while count < args.nev:
if (count % 1000 == 0):
print("Generated {} events".format(count), end='\r')
evt = gen_event((args.pmin, args.pmax), (args.angmin, args.angmax), (args.phmin, args.phmax), parts=parts)
output.write_event(evt)
evt.clear()
count += 1
else:
for mo in args.parray.split(','):
p = float(mo.strip())
count = 0
while count < args.nev:
if (count % 1000 == 0):
print("Generated {} events".format(count), end='\r')
evt = gen_event((p, p), (args.angmin, args.angmax), parts=parts)
output.write_event(evt)
evt.clear()
count += 1
print("Generated {} events".format(args.nev))
output.close()
results/emcal_electrons_energy_reader.pdf deleted 100644 → 0
View file @ 05c0966e
File deleted
results/emcal_electrons_energy_reader.png deleted 100644 → 0
View file @ 05c0966e
results/emcal_electrons_energy_reader.png

7.31 KiB

results/emcal_electrons_eta_reader.pdf deleted 100644 → 0
View file @ 05c0966e
File deleted
results/emcal_electrons_eta_reader.png deleted 100644 → 0
View file @ 05c0966e
results/emcal_electrons_eta_reader.png

7.06 KiB

results/zdc_photons_energy_reader.pdf deleted 100644 → 0
View file @ 05c0966e
File deleted
results/zdc_photons_energy_reader.png deleted 100644 → 0
View file @ 05c0966e
results/zdc_photons_energy_reader.png

7.78 KiB

results/zdc_photons_eta_reader.pdf deleted 100644 → 0
View file @ 05c0966e
File deleted