11 #include <generators/hepmc/HepMCReader.h>
12 #include <framework/gearbox/Unit.h>
13 #include <framework/gearbox/Const.h>
14 #include <framework/logging/Logger.h>
16 #include <HepMC/IO_GenEvent.h>
18 #include <TLorentzVector.h>
24 void HepMCReader::open(
const string& filename)
26 B2INFO(
"Reading HEPMC inputfile at " << filename);
27 m_input.open(filename.c_str());
28 if (!m_input) {
throw (HepMCCouldNotOpenFileError() << filename); }
31 int HepMCReader::nextValidEvent(HepMC::GenEvent& evt)
34 int eventID = evt.event_number();
36 if (eventID < m_minEvent) {
37 while (!(eventID >= m_minEvent)) {
39 eventID = evt.event_number();
42 if (eventID >= m_maxEvent) {
return -999; }
52 eventID = nextValidEvent(evt);
53 if (eventID == -999) {
throw HepMCInvalidEventError(); }
56 const int nparticles = evt.particles_size();
58 B2DEBUG(20,
"Found eventID " << eventID <<
" with " << nparticles <<
" particles.");
60 if (nparticles <= 0) {
61 throw (HepMCInvalidEventError());
63 const int event_offset = graph.
size();
65 for (
int i = 0; i < nparticles; i++) {
69 std::unordered_map<int, int> hash_index_map;
70 HepMC::GenEvent::particle_iterator tmp_particle = evt.particles_begin();
71 for (
int i = 0; i < nparticles; ++i) {
73 const int hash = (*tmp_particle)->barcode();
74 hash_index_map[hash] = i;
77 const double len_conv = HepMC::Units::conversion_factor(evt.length_unit(), HepMC::Units::CM);
78 const double mom_conv = HepMC::Units::conversion_factor(evt.momentum_unit(), HepMC::Units::GEV);
79 auto read_particle = evt.particles_begin();
81 for (
int i = 0; i < nparticles; ++i) {
82 auto* decay_vertex = (*read_particle)->end_vertex();
83 auto* production_vertex = (*read_particle)->production_vertex();
87 const int status = (*read_particle)->status();
88 const bool isFinalstate = !decay_vertex && status == 1;
89 const bool isVirtual = (status == 4) || (status == 21) || (status == 22) || (status == 23);
90 const int pdg_code = (*read_particle)->pdg_id() ;
91 const double mass = (*read_particle)->generated_mass() * mom_conv;
92 auto const mom_tmp = (*read_particle)->momentum();
94 const HepMC::FourVector momentum(
95 mom_tmp.x()*mom_conv * Unit::GeV,
96 mom_tmp.y()*mom_conv * Unit::GeV,
97 mom_tmp.z()*mom_conv * Unit::GeV,
98 mom_tmp.t()*mom_conv * Unit::GeV
101 B2DEBUG(20,
"Read particle: status " << status <<
" isFinal " << isFinalstate <<
" isVirtual " << isVirtual <<
" pdg " << pdg_code
102 <<
" mass " << mass <<
" px " << momentum.x() <<
" py " << momentum.y() <<
" px " << momentum.z() <<
" E " << momentum.t());
103 p.addStatus(MCParticle::c_PrimaryParticle);
105 p.setMomentum(TVector3(momentum.x(), momentum.y(), momentum.z()));
106 p.setEnergy(momentum.t());
108 if (production_vertex) {
109 const auto pos = production_vertex->position();
110 p.setProductionVertex(TVector3(pos.x(), pos.y(), pos.z()) * len_conv * Unit::cm);
111 p.setProductionTime(pos.t() * len_conv * Unit::cm / Const::speedOfLight);
112 p.setValidVertex(
true);
116 int daughter_idx = 0;
117 for (
auto daughter = decay_vertex->particles_begin(HepMC::children);
118 daughter != decay_vertex->particles_end(HepMC::children); ++daughter) {
119 const int daughter_index_in_graph = hash_index_map[(*daughter)->barcode()];
120 p.decaysInto(graph[event_offset + daughter_index_in_graph]);
125 TLorentzVector p4 = p.get4Vector();
127 p4.SetPz(-1.0 * p4.Pz());
129 p4 = m_labboost * p4;
133 if (i < m_nVirtual || isVirtual) {
139 B2DEBUG(20,
"Returning event id " << eventID);
143 void HepMCReader::readNextEvent(HepMC::GenEvent& evt)
148 if (evt.is_valid()) {
149 const int nparticles = evt.particles_size();
150 B2DEBUG(20,
"Found valid event.i N particles " << nparticles);
153 B2DEBUG(20,
"The next event was invalid. Will stop reading now.");
159 int HepMCReader::countEvents(
const std::string& filename)
162 HepMC::IO_GenEvent ascii_in(filename.c_str(), std::ios::in);
164 HepMC::GenEvent* evt = ascii_in.read_next_event();
166 evt = ascii_in.read_next_event();
169 B2INFO(
"Counted " << count <<
" events in " << filename <<
".");