estarlight/src/hepmc3writer.cpp

0001 #include "HepMC3/GenVertex.h"
0002 #include "HepMC3/GenVertex_fwd.h"
0003 #include "HepMC3/FourVector.h"
0004 #include "HepMC3/GenEvent.h"
0005 #include "HepMC3/GenParticle.h"
0006 #include "HepMC3/WriterAscii.h"
0007 #include "HepMC3/Print.h"
0008
0009 #include "inputParameters.h"
0010 #include "eXevent.h"
0011 #include "hepmc3writer.h"
0012
0013 using HepMC3::FourVector;
0014 using HepMC3::Print;
0015
0016 hepMC3Writer::hepMC3Writer()
0017 {
0018 }
0019
0020 int hepMC3Writer::initWriter(const inputParameters &param)
0021 {
0022
0023   std::string hepmc3_filename = param.baseFileName() + ".hepmc";
0024
0025   /** need to pass the parameters to HepMC3 **/
0026   initBeamHepMC3( param ); // Inits the beams
0027   _hepmc3_output = new HepMC3::WriterAscii(hepmc3_filename);
0028   return 0;
0029 }
0030
0031 int hepMC3Writer::initBeamHepMC3(const inputParameters &param)
0032 {
0033
0034   //Proton and Electron Mass
0035   double p_mass = 0.93827;
0036   double e_mass = 0.000511;
0037
0038   //Calculate Four Vector
0039   double electronBeam_E  = e_mass*param.electronBeamLorentzGamma();
0040   double targetBeam_E  = p_mass*param.targetBeamLorentzGamma();
0041   double electronBeam_pz = -1.0*sqrt(electronBeam_E*electronBeam_E - e_mass*e_mass);
0042   double targetBeam_pz = sqrt(targetBeam_E*targetBeam_E - p_mass*p_mass);
0043
0044   //Define Four Vector to fill HepMC3
0045   hepmc3_electronBeam_four_vector_ = FourVector(0,0,electronBeam_pz,electronBeam_E);
0046   hepmc3_targetBeam_four_vector_ = FourVector(0,0,targetBeam_pz,targetBeam_E);
0047   electronBeam_pdg_id_ = 11;
0048   int nuclear_pid = param.targetBeamA()*10 + param.targetBeamZ()*10000 + 1000000000;// Form 100ZZZAAAl; l=0
0049   targetBeam_pdg_id_ = ( param.targetBeamA() > 1 ) ? nuclear_pid : 2212; //Everything is a proton right now
0050
0051   return 0;
0052 }
0053
0054 int hepMC3Writer::writeEvent(const eXEvent &event, int eventnumber)
0055 {
0056
0057   /** Make HepMC3 Event and Vertex ... Currently only two Vertices per Event, the beam and gamma  **/
0058   HepMC3::GenEvent hepmc3_evt(HepMC3::Units::GEV , HepMC3::Units::MM);
0059   hepmc3_evt.set_event_number( eventnumber );
0060   //HepMC3::GenVertexPtr hepmc3_beam_vertex_to_write  = std::make_shared<HepMC3::GenVertex>(FourVector(0,0,0,0));
0061   HepMC3::GenVertexPtr hepmc3_gamma_vertex_to_write = std::make_shared<HepMC3::GenVertex>(FourVector(0,0,0,0));
0062   HepMC3::GenVertexPtr hepmc3_ion_vertex_to_write   = std::make_shared<HepMC3::GenVertex>(FourVector(0,0,0,0));
0063
0064   HepMC3::GenParticlePtr hepmc3_electronBeam_particle = std::make_shared<HepMC3::GenParticle>( hepmc3_electronBeam_four_vector_, electronBeam_pdg_id_, 4 );
0065   HepMC3::GenParticlePtr hepmc3_targetBeam_particle = std::make_shared<HepMC3::GenParticle>( hepmc3_targetBeam_four_vector_, targetBeam_pdg_id_, 4 );
0066
0067   /** Add beam particles to beam vertex **/
0068   //hepmc3_beam_vertex_to_write->add_particle_out( hepmc3_electronBeam_particle );
0069   //hepmc3_beam_vertex_to_write->add_particle_out( hepmc3_targetBeam_particle );
0070
0071   hepmc3_gamma_vertex_to_write->add_particle_in( hepmc3_electronBeam_particle );
0072   hepmc3_ion_vertex_to_write->add_particle_in( hepmc3_targetBeam_particle );
0073   // hepmc3_gamma_vertex_to_write->add_particle_in( hepmc3_electronBeam_particle );
0074   // hepmc3_gamma_vertex_to_write->add_particle_in( hepmc3_targetBeam_particle );
0075   //hepmc3_vertex_to_write->add_particle_in( hepmc3_electronBeam_particle );
0076   //hepmc3_vertex_to_write->add_particle_in( hepmc3_targetBeam_particle );
0077
0078
0079   /** Get starlight particle vector **/
0080   const std::vector<starlightParticle> * particle_vector = event.getParticles();
0081 /*
0082   lorentzVector beamElectron_lorentzVec = (*event.getSources())[0];
0083   HepMC3::GenParticlePtr hepmc3_beamelectron = std::make_shared<HepMC3::GenParticle>( hepmc3_electronBeam_four_vector_, 11, 4 );
0084   hepmc3_vertex_to_write->add_particle_in( hepmc3_beamelectron );
0085 */
0086
0087   lorentzVector gamma_lorentzVec = (*event.getGamma())[0];
0088   FourVector hepmc3_gamma_four_vector = FourVector(gamma_lorentzVec.GetPx(),
0089                            gamma_lorentzVec.GetPy(),
0090                            gamma_lorentzVec.GetPz(),
0091                            gamma_lorentzVec.GetE());
0092   HepMC3::GenParticlePtr hepmc3_gamma = std::make_shared<HepMC3::GenParticle>( hepmc3_gamma_four_vector, 22, 13 );
0093   hepmc3_gamma_vertex_to_write->add_particle_out( hepmc3_gamma );
0094   hepmc3_ion_vertex_to_write->add_particle_in( hepmc3_gamma );
0095
0096   /** Takes e_starlight events and converts to hepmc3 format **/
0097   for ( std::vector<starlightParticle>::const_iterator particle_iter = (*particle_vector).begin();
0098     particle_iter != (*particle_vector).end();
0099     ++particle_iter)
0100     {
0101       int hepmc3_pid = (*particle_iter).getPdgCode();
0102       /** pass to HepMC3 FourVector **/
0103       FourVector hepmc3_four_vector = FourVector( (*particle_iter).GetPx(),
0104                           (*particle_iter).GetPy(),
0105                           (*particle_iter).GetPz(),
0106                           (*particle_iter).GetE());
0107
0108       HepMC3::GenParticlePtr hepmc3_particle = std::make_shared<HepMC3::GenParticle>( hepmc3_four_vector, hepmc3_pid, 1 );
0109       hepmc3_ion_vertex_to_write->add_particle_out( hepmc3_particle );
0110     }
0111
0112   lorentzVector electron_lorentzVec = (*event.getSources())[0];
0113   FourVector hepmc3_electron_four_vector = FourVector(electron_lorentzVec.GetPx(),
0114                            electron_lorentzVec.GetPy(),
0115                            electron_lorentzVec.GetPz(),
0116                            electron_lorentzVec.GetE());
0117
0118   HepMC3::GenParticlePtr hepmc3_electron = std::make_shared<HepMC3::GenParticle>( hepmc3_electron_four_vector, 11, 1 );
0119
0120   hepmc3_gamma_vertex_to_write->add_particle_out( hepmc3_electron );
0121
0122   // add code to output ion  SRK August 8 2021
0123   lorentzVector ion_lorentzVec = (*event.getTarget())[0];
0124   FourVector hepmc3_ion_four_vector = FourVector(ion_lorentzVec.GetPx(),
0125                ion_lorentzVec.GetPy(),
0126                ion_lorentzVec.GetPz(),
0127                ion_lorentzVec.GetE());
0128
0129   //  PDG code for protons is 2212; not clear what to do for ions, but use proton ID for now  SRK August 8, 2021
0130   HepMC3::GenParticlePtr hepmc3_ion = std::make_shared<HepMC3::GenParticle>( hepmc3_ion_four_vector, 2212, 1 );
0131
0132   hepmc3_ion_vertex_to_write->add_particle_out( hepmc3_ion );
0133
0134
0135
0136   /** add vertex to HepMC3 Event**/
0137   //hepmc3_evt.add_vertex( hepmc3_beam_vertex_to_write );
0138   hepmc3_evt.add_vertex( hepmc3_gamma_vertex_to_write );
0139   hepmc3_evt.add_vertex( hepmc3_ion_vertex_to_write );
0140   _hepmc3_output->write_event(hepmc3_evt);
0141   //  Print::listing(hepmc3_evt);
0142   //  Print::content(hepmc3_evt);
0143   hepmc3_evt.clear();
0144
0145   return eventnumber;
0146 }
0147