emcal_electrons.cxx

//////////////////////////////////////////////////////////////
// Crystal EMCAL detector 
// 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"
#include "HepMC3/Print.h"

#include <iostream>
#include<random>
#include<cmath>
#include <math.h>
#include <TMath.h>

using namespace HepMC3;

void emcal_electrons(){
        WriterAscii hepmc_output("./data/emcal_electron_0GeVto30GeV_100kEvt.hepmc");
        int        events_parsed = 0;
        GenEvent   evt(Units::GEV, Units::MM);

	// Random number generator
	TRandom *r1 = new TRandom();

        for (events_parsed = 0; events_parsed < 100000; 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 momentum
	Double_t p = 0.0 + events_parsed * 3.e-4;
	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 11 - electron 0.510 MeV/c^2
        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;
}