Belle II Software  release-05-01-25
HepMCReader Class Reference

Class to read HepMC files and store the content in a MCParticle graph. More...

#include <HepMCReader.h>

Collaboration diagram for HepMCReader:

Public Member Functions

 BELLE2_DEFINE_EXCEPTION (HepMCCouldNotOpenFileError, "Could not open file %1% !")
 Exception is thrown if the HepMC file could not be opened.
 
 BELLE2_DEFINE_EXCEPTION (HepMCInvalidDaughterIndicesError, "Event %1%: Invalid daughter indices d1=%2%, d2=%3%, N=%4% (0<=d1<=d2<=N required)")
 Exception is thrown if the given indices of the daughters are not valid.
 
 BELLE2_DEFINE_EXCEPTION (HepMCHeaderNotValidError, "Event %1%: Event header not understood: %2%")
 Exception is thrown if the header specifying the event header could not be parsed.
 
 BELLE2_DEFINE_EXCEPTION (HepMCConvertFieldError, "Event %1%: Could not convert field %2%: %3%")
 Exception is thrown if a field in the HepMC file could not be converted to a number.
 
 BELLE2_DEFINE_EXCEPTION (HepMCParticleFormatError, "Event %1%: Particle format not understood, got %2% fields !")
 Exception is thrown if the format of a line of the HepMC file could not be parsed.
 
 BELLE2_DEFINE_EXCEPTION (HepMCInvalidEventError, "Event is invalid.")
 Exception is thrown if the format of a line of the HepMC file could not be parsed.
 
 BELLE2_DEFINE_EXCEPTION (HepMCEmptyEventError, "Event %1%: Number of particles in event is %2% ! (This could mean EOF is reached.) ")
 Exception is thrown if the number of particles for this event is 0 or less. More...
 
 HepMCReader (const int minEvent=0, const int maxEvent=INT_MAX)
 Constructor.
 
 ~HepMCReader ()
 Destructor.
 
void open (const std::string &filename)
 Opens an ascii file and prepares it for reading. More...
 
void closeCurrentInputFile ()
 Closes the current input file to allow opening the next one.
 
int nextValidEvent (HepMC::GenEvent &evt)
 Got to next event in bounds. More...
 
int getEvent (MCParticleGraph &graph, double &weight)
 Reads the next event and stores the result in the given MCParticle graph. More...
 
int countEvents (const std::string &filename)
 Count events in file.
 

Public Attributes

int m_nVirtual
 The number of particles in each event with a set Virtual flag.
 
bool m_wrongSignPz
 Bool to indicate that HER and LER were swapped.
 
const int m_minEvent
 min event nr to process
 
const int m_maxEvent
 max events to process
 
TLorentzRotation m_labboost
 Boost&rotation vector for boost from CM to LAB.
 

Protected Member Functions

void readNextEvent (HepMC::GenEvent &evt)
 read the next event from the IO stream and write into evt
 

Protected Attributes

std::ifstream m_input
 The input stream of the ascii file.
 

Detailed Description

Class to read HepMC files and store the content in a MCParticle graph.

The reader supports retrieving the HepMC2 information from an ascii text file as provided by for example pythia8.

the events are read sequentially with the option of skipping a certain number of events. Random access of events is not possible.

Definition at line 44 of file HepMCReader.h.

Member Function Documentation

◆ BELLE2_DEFINE_EXCEPTION()

BELLE2_DEFINE_EXCEPTION ( HepMCEmptyEventError  ,
"Event %1%: Number of particles in event is %2% ! (This could mean EOF is reached.) "   
)

Exception is thrown if the number of particles for this event is 0 or less.


◆ getEvent()

int getEvent ( MCParticleGraph graph,
double &  weight 
)

Reads the next event and stores the result in the given MCParticle graph.

Parameters
graphReference to the graph which should be filled with the information from the Hepevt file.
weightReference to the event weight which can be filled from the file.
Returns
event numer if the event could be read and the number was provided in the file.

Definition at line 46 of file HepMCReader.cc.

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 }

◆ nextValidEvent()

int nextValidEvent ( HepMC::GenEvent &  evt)

Got to next event in bounds.

Throw InvalidEventException if we went over the bounds

Definition at line 31 of file HepMCReader.cc.

◆ open()

void open ( const std::string &  filename)

Opens an ascii file and prepares it for reading.

Parameters
filenameThe filename of the HepMC2 ascii file which should be read.

Definition at line 24 of file HepMCReader.cc.


The documentation for this class was generated from the following files:
Belle2::Unit::cm
static const double cm
Standard units with the value = 1.
Definition: Unit.h:57
Belle2::MCParticleGraph::size
size_t size() const
Return the number of particles in the graph.
Definition: MCParticleGraph.h:253
Belle2::HepMCReader::m_wrongSignPz
bool m_wrongSignPz
Bool to indicate that HER and LER were swapped.
Definition: HepMCReader.h:101
Belle2::Const::speedOfLight
static const double speedOfLight
[cm/ns]
Definition: Const.h:568
Belle2::HepMCReader::m_nVirtual
int m_nVirtual
The number of particles in each event with a set Virtual flag.
Definition: HepMCReader.h:100
Belle2::HepMCReader::nextValidEvent
int nextValidEvent(HepMC::GenEvent &evt)
Got to next event in bounds.
Definition: HepMCReader.cc:31
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:305
Belle2::HepMCReader::m_labboost
TLorentzRotation m_labboost
Boost&rotation vector for boost from CM to LAB.
Definition: HepMCReader.h:104
Belle2::Unit::GeV
static const double GeV
Standard of [energy, momentum, mass].
Definition: Unit.h:61
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86
Belle2::MCParticle::c_PrimaryParticle
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:58