Skip to content
Snippets Groups Projects
forward_ions.cxx 1.81 KiB
Newer Older
  • Learn to ignore specific revisions
  • Tom Polakovic's avatar
    Tom Polakovic committed
    #include "HepMC3/GenEvent.h"
    #include "HepMC3/ReaderAscii.h"
    #include "HepMC3/WriterAscii.h"
    #include "HepMC3/Print.h"
    
    #include <iostream>
    #include <random>
    using namespace HepMC3;
    
    
    Tom Polakovic's avatar
    Tom Polakovic committed
    void forward_ions(){
    
    Tom Polakovic's avatar
    Tom Polakovic committed
    
      //-------------------------------------
      WriterAscii hepmc_output("data/forward_ions.hepmc");
      int        events_parsed = 0;
      GenEvent   evt(Units::GEV, Units::MM);
    
      std::mt19937 gen( std::stoi(std::getenv("SEED")) );
    
    Tom Polakovic's avatar
    Tom Polakovic committed
      std::uniform_real_distribution<double> phi_distr(0,6.2831853072);
    
      std::uniform_real_distribution<double> theta_distr(0,0.00872664626);
    
    Tom Polakovic's avatar
    Tom Polakovic committed
    
    
      for (events_parsed = 0; events_parsed < 300; events_parsed++) {
    
    Tom Polakovic's avatar
    Tom Polakovic committed
        // 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 + 0.65*events_parsed;
    
    Tom Polakovic's avatar
    Tom Polakovic committed
        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;
    }