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.

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

~HepMCReader()

~HepMCReader ( )

inline

Destructor.

Definition at line 62 of file HepMCReader.h.

{ 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.

{m_input.close();}

countEvents()

int countEvents

(

const std::string &

filename

)

Count events in file.

Definition at line 171 of file HepMCReader.cc.

{
  //different way to read file for consistency check
  HepMC::IO_GenEvent ascii_in(filename.c_str(), std::ios::in);
  int count = 0;
  HepMC::GenEvent* evt = ascii_in.read_next_event();
  while (evt) {
    evt = ascii_in.read_next_event();
    count++;
  }
  B2INFO("Counted " << count  << " events in  " << filename << ".");
  return count;
}

getEvent()

int getEvent	(	MCParticleGraph &	graph,
		double &	weight )

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

Parameters

graph	Reference to the graph which should be filled with the information from the Hepevt file.
weight	Reference 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.

{
  int eventID;
 
  HepMC::GenEvent evt;
  // read event number once
  eventID = nextValidEvent(evt); // check if this event is in the bounds if not go to next one in bounds
  if (eventID == -999) { throw HepMCInvalidEventError(); }
 
  eventWeight = 1; // why is evt.weights() a std::vector?
  const int nparticles = evt.particles_size();
 
  B2DEBUG(20, "Found eventID " << eventID << " with " << nparticles << " particles.");
 
  if (nparticles <= 0) {
    throw (HepMCInvalidEventError());
  }
  const int event_offset = graph.size(); //offset
  //Make list of particles; Prepare graph
  for (int i = 0; i < nparticles; i++) {
    graph.addParticle();
  }
 
  std::unordered_map<int, int> hash_index_map;
  HepMC::GenEvent::particle_iterator tmp_particle = evt.particles_begin();
  for (int i = 0; i < nparticles; ++i) {
    //the barcode seems to just do 0++ etc. but im not sure if there are exceptions
    const int hash = (*tmp_particle)->barcode();
    hash_index_map[hash] = i;
    ++tmp_particle;
  }
  const double len_conv = HepMC::Units::conversion_factor(evt.length_unit(), HepMC::Units::CM); // from, to
  const double mom_conv = HepMC::Units::conversion_factor(evt.momentum_unit(), HepMC::Units::GEV); // from, to
  auto read_particle = evt.particles_begin();
  //Read particles from file
  for (int i = 0; i < nparticles; ++i) {
    auto* decay_vertex = (*read_particle)->end_vertex();
    auto* production_vertex = (*read_particle)->production_vertex();
 
    MCParticleGraph::GraphParticle& p = graph[event_offset + i];
 
    const int status = (*read_particle)->status();
    const bool isFinalstate =  !decay_vertex && status == 1;
    //the status numbers can be read in more detail here: https://pythia.org/latest-manual/ParticleProperties.html
    const bool isVirtual  = (status == 4)  || (status == 21) || (status == 22) || (status == 23)
                            //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
                            || (status == 41) || (status == 42) || (status == 43) || (status == 44)
                            //51 through 73 are numbers which come up for outgoing mediators and particles which are virtual
                            || (status == 51) || (status == 52) || (status == 71) || (status == 73)
                            //83 through 91 are numbers which appear for intermediate outgoing particles and their decay particles, but are not the final state
                            || (status == 83) || (status == 84)
                            || (status == 91); //hepmc internal status flags. 4 are beam particles, rest could be refined if needed
    const int pdg_code = (*read_particle)->pdg_id() ;
    const double mass = (*read_particle)->generated_mass() * mom_conv;
    auto const mom_tmp = (*read_particle)->momentum();
 
    const HepMC::FourVector momentum(
      mom_tmp.x()*mom_conv * Unit::GeV,
      mom_tmp.y()*mom_conv * Unit::GeV,
      mom_tmp.z()*mom_conv * Unit::GeV,
      mom_tmp.t()*mom_conv * Unit::GeV
    );
 
    B2DEBUG(20, "Read particle: status " << status << " isFinal " << isFinalstate << " isVirtual " << isVirtual << " pdg " << pdg_code
            << " mass " << mass << " px " << momentum.x() << " py " << momentum.y() << " px " << momentum.z() << " E " << momentum.t());
    p.addStatus(MCParticle::c_PrimaryParticle);  // all particles part of the hepmc file should be set as primary
    p.setPDG(pdg_code);
    p.setMomentum(ROOT::Math::XYZVector(momentum.x(), momentum.y(), momentum.z()));
    p.setEnergy(momentum.t());
    p.setMass(mass);
    if (production_vertex) {
      const auto pos = production_vertex->position();
      p.setProductionVertex(ROOT::Math::XYZVector(pos.x(), pos.y(), pos.z()) * len_conv * Unit::cm);
      p.setProductionTime(pos.t() * len_conv * Unit::cm / Const::speedOfLight);
      p.setValidVertex(true);
    }
 
    if ((status == 21) || (status == 41)
        || (status ==
            42))  {      //removes massless beam particles carried over from MadGraph (41 and 42 are equivalent to 21 but for when there's an ISR photon)
      p.setIgnore(true);    //they are just used for internal bookkeeping and serve no other purpose
    }
 
    //assign the parent particle:
    //for two incoming particles only one of them is assigned as parent
    if (production_vertex) {
      if (production_vertex->particles_in_const_begin() != production_vertex->particles_in_const_end()) {
        auto parent = production_vertex->particles_begin(HepMC::parents);
        const int parent_index_in_graph = hash_index_map[(*parent)->barcode()];
        p.comesFrom(graph[event_offset + parent_index_in_graph]);
      }
    }
 
    if (m_wrongSignPz) { // this means we have to mirror Pz
      ROOT::Math::PxPyPzEVector p4 = p.get4Vector();
      p4.SetPz(-1.0 * p4.Pz());
      p.set4Vector(p4);
    }
 

    if (i < m_nVirtual || isVirtual) {
      p.setVirtual();
    }
    ++read_particle;
  }
  eventID += 1;
  B2DEBUG(20, "Returning event id " << eventID);
  return eventID;
}

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.

{
  readNextEvent(evt);
  int eventID = evt.event_number();
  // if eventID is not in the valid bounds get next event
  if (eventID < m_minEvent) {
    while (!(eventID >= m_minEvent)) {
      readNextEvent(evt);
      eventID = evt.event_number();
    }
  }
  if (eventID >= m_maxEvent) { return -999; } // went too far; there is no int nan
  return eventID;
}

open()

void open

(

const std::string &

filename

)

Opens an ascii file and prepares it for reading.

Parameters

filename

The filename of the HepMC2 ascii file which should be read.

Definition at line 22 of file HepMCReader.cc.

{
  B2INFO("Reading HEPMC inputfile at " << filename);
  m_input.open(filename.c_str());
  if (!m_input) { throw (HepMCCouldNotOpenFileError() << filename); }
}

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.

{
  // walk through the IOStream
  if (m_input) {
    evt.read(m_input);
    if (evt.is_valid()) {
      const int nparticles = evt.particles_size();
      B2DEBUG(20, "Found valid event.i N particles " << nparticles);
      return; //
    } else {
      B2DEBUG(20, "The next event was invalid. Will stop reading now.");
    }
  }
  return;
}

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:

• generators/hepmc/include/HepMCReader.h
• generators/hepmc/src/HepMCReader.cc