Skip to content
Snippets Groups Projects
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;
}