Belle II Software  release-05-02-19
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) || (status == 51) || (status == 52)
90  || (status == 71) ; //hepmc internal status flags. 4 are beam particles, rest could be refined if needed
91  const int pdg_code = (*read_particle)->pdg_id() ;
92  const double mass = (*read_particle)->generated_mass() * mom_conv;
93  auto const mom_tmp = (*read_particle)->momentum();
94 
95  const HepMC::FourVector momentum(
96  mom_tmp.x()*mom_conv * Unit::GeV,
97  mom_tmp.y()*mom_conv * Unit::GeV,
98  mom_tmp.z()*mom_conv * Unit::GeV,
99  mom_tmp.t()*mom_conv * Unit::GeV
100  );
101 
102  B2DEBUG(20, "Read particle: status " << status << " isFinal " << isFinalstate << " isVirtual " << isVirtual << " pdg " << pdg_code
103  << " mass " << mass << " px " << momentum.x() << " py " << momentum.y() << " px " << momentum.z() << " E " << momentum.t());
104  p.addStatus(MCParticle::c_PrimaryParticle); // all particles part of the hepmc file should be set as primary
105  p.setPDG(pdg_code);
106  p.setMomentum(TVector3(momentum.x(), momentum.y(), momentum.z()));
107  p.setEnergy(momentum.t());
108  p.setMass(mass);
109  if (production_vertex) {
110  const auto pos = production_vertex->position();
111  p.setProductionVertex(TVector3(pos.x(), pos.y(), pos.z()) * len_conv * Unit::cm);
112  p.setProductionTime(pos.t() * len_conv * Unit::cm / Const::speedOfLight);
113  p.setValidVertex(true);
114  }
115 
116  if (status == 21) { //removes massless beam particles carried over from MadGraph
117  p.setIgnore(true); //they are just used for internal bookkeeping and serve no other purpose
118  }
119 
120  //assign the parent particle:
121  //for two incoming particles only one of them is assigned as parent
122  if (production_vertex) {
123  if (production_vertex->particles_in_const_begin() != production_vertex->particles_in_const_end()) {
124  auto parent = production_vertex->particles_begin(HepMC::parents);
125  const int parent_index_in_graph = hash_index_map[(*parent)->barcode()];
126  p.comesFrom(graph[event_offset + parent_index_in_graph]);
127  }
128  }
129 
130  //boost particles to lab frame:
131  TLorentzVector p4 = p.get4Vector();
132  if (m_wrongSignPz) { // this means we have to mirror Pz
133  p4.SetPz(-1.0 * p4.Pz());
134  }
135  p4 = m_labboost * p4;
136  p.set4Vector(p4);
137 
139  if (i < m_nVirtual || isVirtual) {
140  p.setVirtual();
141  }
142  ++read_particle;
143  }
144  eventID += 1;
145  B2DEBUG(20, "Returning event id " << eventID);
146  return eventID;
147 }

◆ 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