Skip to content
Snippets Groups Projects

pi0 data set

Files

+ 76
0
 
//////////////////////////////////////////////////////////////
 
// 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_100kEvt.hepmc")
 
{
 
WriterAscii hepmc_output(out_fname);
 
int events_parsed = 0;
 
GenEvent evt(Units::GEV, Units::MM);
 
 
// Random number generator
 
TRandom *r1 = new TRandom();
 
 
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(0.0,30.0);
 
Double_t px;
 
Double_t py;
 
Double_t pz;
 
// 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;
 
}
Loading