Belle II Software development
HepMCReader Class Reference

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

#include <HepMCReader.h>

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.
 
 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.
 
void closeCurrentInputFile ()
 Closes the current input file to allow opening the next one.
 
int nextValidEvent (HepMC::GenEvent &evt)
 Got to next event in bounds.
 
int getEvent (MCParticleGraph &graph, double &weight)
 Reads the next event and stores the result in the given MCParticle graph.
 
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
 

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 32 of file HepMCReader.h.

Constructor & Destructor Documentation

◆ HepMCReader()

HepMCReader ( const int minEvent = 0,
const int maxEvent = INT_MAX )
inline

Constructor.

Definition at line 56 of file HepMCReader.h.

56 : m_nVirtual(0), m_wrongSignPz(false), m_minEvent(minEvent),
57 m_maxEvent(maxEvent) {}

◆ ~HepMCReader()

~HepMCReader ( )
inline

Destructor.

Definition at line 62 of file HepMCReader.h.

62{ if (m_input) {m_input.close();} }

Member Function Documentation

◆ closeCurrentInputFile()

void closeCurrentInputFile ( )
inline

Closes the current input file to allow opening the next one.

Definition at line 73 of file HepMCReader.h.

73{m_input.close();}

◆ countEvents()

int countEvents ( const std::string & filename)

Count events in file.

Definition at line 171 of file HepMCReader.cc.

172{
173 //different way to read file for consistency check
174 HepMC::IO_GenEvent ascii_in(filename.c_str(), std::ios::in);
175 int count = 0;
176 HepMC::GenEvent* evt = ascii_in.read_next_event();
177 while (evt) {
178 evt = ascii_in.read_next_event();
179 count++;
180 }
181 B2INFO("Counted " << count << " events in " << filename << ".");
182 return count;
183}

◆ 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 number if the event could be read and the number was provided in the file.

Definition at line 44 of file HepMCReader.cc.

45{
46 int eventID;
47
48 HepMC::GenEvent evt;
49 // read event number once
50 eventID = nextValidEvent(evt); // check if this event is in the bounds if not go to next one in bounds
51 if (eventID == -999) { throw HepMCInvalidEventError(); }
52
53 eventWeight = 1; // why is evt.weights() a std::vector?
54 const int nparticles = evt.particles_size();
55
56 B2DEBUG(20, "Found eventID " << eventID << " with " << nparticles << " particles.");
57
58 if (nparticles <= 0) {
59 throw (HepMCInvalidEventError());
60 }
61 const int event_offset = graph.size(); //offset
62 //Make list of particles; Prepare graph
63 for (int i = 0; i < nparticles; i++) {
64 graph.addParticle();
65 }
66
67 std::unordered_map<int, int> hash_index_map;
68 HepMC::GenEvent::particle_iterator tmp_particle = evt.particles_begin();
69 for (int i = 0; i < nparticles; ++i) {
70 //the barcode seems to just do 0++ etc. but im not sure if there are exceptions
71 const int hash = (*tmp_particle)->barcode();
72 hash_index_map[hash] = i;
73 ++tmp_particle;
74 }
75 const double len_conv = HepMC::Units::conversion_factor(evt.length_unit(), HepMC::Units::CM); // from, to
76 const double mom_conv = HepMC::Units::conversion_factor(evt.momentum_unit(), HepMC::Units::GEV); // from, to
77 auto read_particle = evt.particles_begin();
78 //Read particles from file
79 for (int i = 0; i < nparticles; ++i) {
80 auto* decay_vertex = (*read_particle)->end_vertex();
81 auto* production_vertex = (*read_particle)->production_vertex();
82
83 MCParticleGraph::GraphParticle& p = graph[event_offset + i];
84
85 const int status = (*read_particle)->status();
86 const bool isFinalstate = !decay_vertex && status == 1;
87 //the status numbers can be read in more detail here: https://pythia.org/latest-manual/ParticleProperties.html
88 const bool isVirtual = (status == 4) || (status == 21) || (status == 22) || (status == 23)
89 //41 through 44 appear for beam particles when there is an ISR photon involved, which would have been originally status numbers 21/22/23 without the ISR photon
90 || (status == 41) || (status == 42) || (status == 43) || (status == 44)
91 //51 through 73 are numbers which come up for outgoing mediators and particles which are virtual
92 || (status == 51) || (status == 52) || (status == 71) || (status == 73)
93 //83 through 91 are numbers which appear for intermediate outgoing particles and their decay particles, but are not the final state
94 || (status == 83) || (status == 84)
95 || (status == 91); //hepmc internal status flags. 4 are beam particles, rest could be refined if needed
96 const int pdg_code = (*read_particle)->pdg_id() ;
97 const double mass = (*read_particle)->generated_mass() * mom_conv;
98 auto const mom_tmp = (*read_particle)->momentum();
99
100 const HepMC::FourVector momentum(
101 mom_tmp.x()*mom_conv * Unit::GeV,
102 mom_tmp.y()*mom_conv * Unit::GeV,
103 mom_tmp.z()*mom_conv * Unit::GeV,
104 mom_tmp.t()*mom_conv * Unit::GeV
105 );
106
107 B2DEBUG(20, "Read particle: status " << status << " isFinal " << isFinalstate << " isVirtual " << isVirtual << " pdg " << pdg_code
108 << " mass " << mass << " px " << momentum.x() << " py " << momentum.y() << " px " << momentum.z() << " E " << momentum.t());
109 p.addStatus(MCParticle::c_PrimaryParticle); // all particles part of the hepmc file should be set as primary
110 p.setPDG(pdg_code);
111 p.setMomentum(ROOT::Math::XYZVector(momentum.x(), momentum.y(), momentum.z()));
112 p.setEnergy(momentum.t());
113 p.setMass(mass);
114 if (production_vertex) {
115 const auto pos = production_vertex->position();
116 p.setProductionVertex(ROOT::Math::XYZVector(pos.x(), pos.y(), pos.z()) * len_conv * Unit::cm);
117 p.setProductionTime(pos.t() * len_conv * Unit::cm / Const::speedOfLight);
118 p.setValidVertex(true);
119 }
120
121 if ((status == 21) || (status == 41)
122 || (status ==
123 42)) { //removes massless beam particles carried over from MadGraph (41 and 42 are equivalent to 21 but for when there's an ISR photon)
124 p.setIgnore(true); //they are just used for internal bookkeeping and serve no other purpose
125 }
126
127 //assign the parent particle:
128 //for two incoming particles only one of them is assigned as parent
129 if (production_vertex) {
130 if (production_vertex->particles_in_const_begin() != production_vertex->particles_in_const_end()) {
131 auto parent = production_vertex->particles_begin(HepMC::parents);
132 const int parent_index_in_graph = hash_index_map[(*parent)->barcode()];
133 p.comesFrom(graph[event_offset + parent_index_in_graph]);
134 }
135 }
136
137 if (m_wrongSignPz) { // this means we have to mirror Pz
138 ROOT::Math::PxPyPzEVector p4 = p.get4Vector();
139 p4.SetPz(-1.0 * p4.Pz());
140 p.set4Vector(p4);
141 }
142
143
145 if (i < m_nVirtual || isVirtual) {
146 p.setVirtual();
147 }
148 ++read_particle;
149 }
150 eventID += 1;
151 B2DEBUG(20, "Returning event id " << eventID);
152 return eventID;
153}
static const double speedOfLight
[cm/ns]
Definition Const.h:695
bool m_wrongSignPz
Bool to indicate that HER and LER were swapped.
Definition HepMCReader.h:89
int nextValidEvent(HepMC::GenEvent &evt)
Got to next event in bounds.
int m_nVirtual
The number of particles in each event with a set Virtual flag.
Definition HepMCReader.h:88
size_t size() const
Return the number of particles in the graph.
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition MCParticle.h:47
static const double cm
Standard units with the value = 1.
Definition Unit.h:47
static const double GeV
Standard of [energy, momentum, mass].
Definition Unit.h:51
GraphParticle & addParticle()
Add new particle to the graph.

◆ nextValidEvent()

int nextValidEvent ( HepMC::GenEvent & evt)

Got to next event in bounds.

Throw InvalidEventException if we went over the bounds

Definition at line 29 of file HepMCReader.cc.

30{
31 readNextEvent(evt);
32 int eventID = evt.event_number();
33 // if eventID is not in the valid bounds get next event
34 if (eventID < m_minEvent) {
35 while (!(eventID >= m_minEvent)) {
36 readNextEvent(evt);
37 eventID = evt.event_number();
38 }
39 }
40 if (eventID >= m_maxEvent) { return -999; } // went too far; there is no int nan
41 return eventID;
42}
void readNextEvent(HepMC::GenEvent &evt)
read the next event from the IO stream and write into evt
const int m_maxEvent
max events to process
Definition HepMCReader.h:91
const int m_minEvent
min event nr to process
Definition HepMCReader.h:90

◆ 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 22 of file HepMCReader.cc.

23{
24 B2INFO("Reading HEPMC inputfile at " << filename);
25 m_input.open(filename.c_str());
26 if (!m_input) { throw (HepMCCouldNotOpenFileError() << filename); }
27}
std::ifstream m_input
The input stream of the ascii file.
Definition HepMCReader.h:97

◆ readNextEvent()

void readNextEvent ( HepMC::GenEvent & evt)
protected

read the next event from the IO stream and write into evt

Definition at line 155 of file HepMCReader.cc.

156{
157 // walk through the IOStream
158 if (m_input) {
159 evt.read(m_input);
160 if (evt.is_valid()) {
161 const int nparticles = evt.particles_size();
162 B2DEBUG(20, "Found valid event.i N particles " << nparticles);
163 return; //
164 } else {
165 B2DEBUG(20, "The next event was invalid. Will stop reading now.");
166 }
167 }
168 return;
169}

Member Data Documentation

◆ m_input

std::ifstream m_input
protected

The input stream of the ascii file.

Definition at line 97 of file HepMCReader.h.

◆ m_maxEvent

const int m_maxEvent

max events to process

Definition at line 91 of file HepMCReader.h.

◆ m_minEvent

const int m_minEvent

min event nr to process

Definition at line 90 of file HepMCReader.h.

◆ m_nVirtual

int m_nVirtual

The number of particles in each event with a set Virtual flag.

Definition at line 88 of file HepMCReader.h.

◆ m_wrongSignPz

bool m_wrongSignPz

Bool to indicate that HER and LER were swapped.

Definition at line 89 of file HepMCReader.h.


The documentation for this class was generated from the following files: