Newer
Older
///////////////////////////
// Crystal EMCAL detector
// Electron dataset
// J.KIM 07/20/2020
///////////////////////////
#include "HepMC3/GenEvent.h"
#include "HepMC3/ReaderAscii.h"
#include "HepMC3/WriterAscii.h"
#include "HepMC3/Print.h"
#include <iostream>
#include<random>
#include<cmath>
using namespace HepMC3;
void emcal_electrons(){
WriterAscii hepmc_output("./data/emcal_electrons.hepmc");
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")) );
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
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]
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 =
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 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);
// 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);
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 % 10 == 0) {
std::cout << "Event: " << events_parsed << std::endl;
}
evt.clear();
}
hepmc_output.close();
std::cout << "Events parsed and written: " << events_parsed << std::endl;
}