Belle II Software  release-05-01-25
HepMCReader.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Jo-Frederik Krohn *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <generators/hepmc/HepMCReader.h>
12 #include <framework/gearbox/Unit.h>
13 #include <framework/gearbox/Const.h>
14 #include <framework/logging/Logger.h>
15 
16 #include <HepMC/IO_GenEvent.h>
17 
18 #include <TLorentzVector.h>
19 
20 using namespace std;
21 using namespace Belle2;
22 
23 
24 void HepMCReader::open(const string& filename)
25 {
26  B2INFO("Reading HEPMC inputfile at " << filename);
27  m_input.open(filename.c_str());
28  if (!m_input) { throw (HepMCCouldNotOpenFileError() << filename); }
29 }
30 
31 int HepMCReader::nextValidEvent(HepMC::GenEvent& evt)
32 {
33  readNextEvent(evt);
34  int eventID = evt.event_number();
35  // if eventID is not in the valid bounds get next event
36  if (eventID < m_minEvent) {
37  while (!(eventID >= m_minEvent)) {
38  readNextEvent(evt);
39  eventID = evt.event_number();
40  }
41  }
42  if (eventID >= m_maxEvent) { return -999; } // went too far; there is no int nan
43  return eventID;
44 }
45 
46 int HepMCReader::getEvent(MCParticleGraph& graph, double& eventWeight)
47 {
48  int eventID;
49 
50  HepMC::GenEvent evt;
51  // read event number once
52  eventID = nextValidEvent(evt); // check if this event is in the bounds if not go to next one in bounds
53  if (eventID == -999) { throw HepMCInvalidEventError(); }
54 
55  eventWeight = 1; // why is evt.weights() a std::vector?
56  const int nparticles = evt.particles_size();
57 
58  B2DEBUG(20, "Found eventID " << eventID << " with " << nparticles << " particles.");
59 
60  if (nparticles <= 0) {
61  throw (HepMCInvalidEventError());
62  }
63  const int event_offset = graph.size(); //offset
64  //Make list of particles; Prepare graph
65  for (int i = 0; i < nparticles; i++) {
66  graph.addParticle();
67  }
68 
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) {
72  //the barcode seems to just do 0++ etc. but im not sure if there are exceptions
73  const int hash = (*tmp_particle)->barcode();
74  hash_index_map[hash] = i;
75  ++tmp_particle;
76  }
77  const double len_conv = HepMC::Units::conversion_factor(evt.length_unit(), HepMC::Units::CM); // from, to
78  const double mom_conv = HepMC::Units::conversion_factor(evt.momentum_unit(), HepMC::Units::GEV); // from, to
79  auto read_particle = evt.particles_begin();
80  //Read particles from file
81  for (int i = 0; i < nparticles; ++i) {
82  auto* decay_vertex = (*read_particle)->end_vertex();
83  auto* production_vertex = (*read_particle)->production_vertex();
84 
85  MCParticleGraph::GraphParticle& p = graph[event_offset + i];
86 
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); //TODO refine
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();
93  //whatever genius wrote this vector class did not implement an operator for multiplication with a scalar or even access via []
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
99  );
100 
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); // bit needs to be set in order for the GeneratedVertexDisplacer to work
104  p.setPDG(pdg_code);
105  p.setMomentum(TVector3(momentum.x(), momentum.y(), momentum.z()));
106  p.setEnergy(momentum.t());
107  p.setMass(mass);
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);
113  }
114 
115  if (decay_vertex) {
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]);
121  daughter_idx++;
122  }
123  }
124  //boost particles to lab frame:
125  TLorentzVector p4 = p.get4Vector();
126  if (m_wrongSignPz) { // this means we have to mirror Pz
127  p4.SetPz(-1.0 * p4.Pz());
128  }
129  p4 = m_labboost * p4;
130  p.set4Vector(p4);
131 
133  if (i < m_nVirtual || isVirtual) {
134  p.setVirtual();
135  }
136  ++read_particle;
137  }
138  eventID += 1;
139  B2DEBUG(20, "Returning event id " << eventID);
140  return eventID;
141 }
142 
143 void HepMCReader::readNextEvent(HepMC::GenEvent& evt)
144 {
145  // walk through the IOStream
146  if (m_input) {
147  evt.read(m_input);
148  if (evt.is_valid()) {
149  const int nparticles = evt.particles_size();
150  B2DEBUG(20, "Found valid event.i N particles " << nparticles);
151  return; //
152  } else {
153  B2DEBUG(20, "The next event was invalid. Will stop reading now.");
154  }
155  }
156  return;
157 }
158 
159 int HepMCReader::countEvents(const std::string& filename)
160 {
161  //different way to read file for consitency check
162  HepMC::IO_GenEvent ascii_in(filename.c_str(), std::ios::in);
163  int count = 0;
164  HepMC::GenEvent* evt = ascii_in.read_next_event();
165  while (evt) {
166  evt = ascii_in.read_next_event();
167  count++;
168  }
169  B2INFO("Counted " << count << " events in " << filename << ".");
170  return count;
171 }
172 
173 
Belle2::MCParticleGraph::size
size_t size() const
Return the number of particles in the graph.
Definition: MCParticleGraph.h:253
Belle2::MCParticleGraph
Class to build, validate and sort a particle decay chain.
Definition: MCParticleGraph.h:48
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:305
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86