-
Tom Polakovic authoredTom Polakovic authored
forward_ions.cxx 1.78 KiB
#include "HepMC3/GenEvent.h"
#include "HepMC3/ReaderAscii.h"
#include "HepMC3/WriterAscii.h"
#include "HepMC3/Print.h"
#include <iostream>
#include <random>
using namespace HepMC3;
void forward_ions(){
//-------------------------------------
WriterAscii hepmc_output("data/forward_ions.hepmc");
int events_parsed = 0;
GenEvent evt(Units::GEV, Units::MM);
std::mt19937 gen(1337);
std::uniform_real_distribution<double> phi_distr(0,6.2831853072);
std::uniform_real_distribution<double> theta_distr(0,0.01745329252);
for (events_parsed = 0; events_parsed < 100; events_parsed++) {
// type 4 is beam
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);
double p = 80 + 1.95*events_parsed;
double phi = phi_distr(gen);
double theta = theta_distr(gen);
double px = p * std::sin(theta) * std::cos(phi);
double py = p * std::sin(theta) * std::sin(phi);
double pz = p * std::cos(theta);
auto pvec = FourVector(px, py, pz, std::sqrt(px*px + py*py + pz*pz + 0.938*0.938));
// type 1 is final state
GenParticlePtr p3 =
std::make_shared<GenParticle>(pvec, 2212, 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;
}