9#include <framework/gearbox/Unit.h>
10#include <framework/gearbox/Const.h>
11#include <framework/logging/Logger.h>
12#include <generators/lhe/LHEReader.h>
13#include <mdst/dataobjects/MCParticle.h>
16#include <boost/lexical_cast.hpp>
17#include <boost/algorithm/string.hpp>
18#include <boost/format.hpp>
31 if (!
m_input)
throw(LHECouldNotOpenFileError() << filename);
37 if (nparticles <= 0) {
38 throw (LHEEmptyEventError() <<
m_lineNr << nparticles);
41 int first = graph.
size();
43 for (
int i = 0; i < nparticles; i++) {
47 double r, x = 0, y = 0, z = 0, t = 0;
49 for (
int i = 0; i < nparticles; ++i) {
61 ROOT::Math::PxPyPzEVector p4 = p.get4Vector();
62 p4.SetPz(-1.0 * p4.Pz());
63 p4.SetPx(-1.0 * p4.Px());
70 TF1 fr(
"fr",
"exp(-x/[0])", 0, 1000000);
71 ROOT::Math::PxPyPzEVector p4 = p.get4Vector();
75 x = r * p4.Px() / p4.P();
76 y = r * p4.Py() / p4.P();
77 z = r * p4.Pz() / p4.P();
78 p.setDecayVertex(x, y, z);
81 p.setValidVertex(
true);
86 p.setProductionVertex(x, y, z);
87 p.setProductionTime(t);
88 p.setValidVertex(
true);
113 for (
int i = 0; i < n; i++) {
116 if (nparticles < 0)
return false;
117 for (
int j = 0; j < nparticles; j++)
getLine();
132 size_t commentPos = line.find_first_of(
'#');
133 if (commentPos != string::npos) {
134 line = line.substr(0, commentPos);
138 }
while (line ==
"" && !
m_input.eof());
151 size_t commentPos = line2.find_first_of(
'#');
152 if (commentPos != string::npos) {
153 line2 = line2.substr(0, commentPos);
157 }
while (line2 !=
"<event>" && !
m_input.eof());
163 if (line ==
"" ||
m_input.eof())
return -1;
165 vector<double> fields;
171 for (
const string& tok : tokens) {
174 fields.push_back(boost::lexical_cast<double>(tok));
175 }
catch (boost::bad_lexical_cast& e) {
176 throw (LHEConvertFieldError() <<
m_lineNr << index << tok);
180 if (fields.size() < 3) {
181 nparticles =
static_cast<int>(fields[0]);
184 nparticles =
static_cast<int>(fields[0]);
185 eventWeight =
static_cast<double>(fields[2]);
195 vector<double> fields;
201 for (
const string& tok : tokens) {
204 fields.push_back(boost::lexical_cast<double>(tok));
205 }
catch (boost::bad_lexical_cast& e) {
206 throw (LHEConvertFieldError() <<
m_lineNr << index << tok);
210 switch (fields.size()) {
213 particle.setPDG(
static_cast<int>(fields[0]));
214 mother =
static_cast<int>(fields[2]);
215 particle.setMomentum(ROOT::Math::XYZVector(fields[6], fields[7], fields[8]));
216 particle.setEnergy(fields[9]);
217 particle.setMass(fields[10]);
220 throw (LHEParticleFormatError() <<
m_lineNr << fields.size());
static const double speedOfLight
[cm/ns]
int readParticle(MCParticleGraph::GraphParticle &particle)
Reads the information for a single particle from the LHE file.
bool skipEvents(int n)
Skips a given number of events.
std::string getLine()
Returns the current line from the LHE ascii file.
double m_Rmin
Minimum of vertex distance to IP.
int m_pdgDisplaced
PDG code of the displaced particle being studied.
static const boost::char_separator< char > sep
The characters at which to split, defaults to ",; \t".
int m_lineNr
The current line number within the ascii file.
int m_indexVirtual
Maximum index of particles in each event that must be set as c_IsVirtual (1-based).
double m_Rmax
Maximum of vertex distance to IP.
void open(const std::string &filename)
Opens an ascii file and prepares it for reading.
std::ifstream m_input
The input stream of the ascii file.
boost::tokenizer< boost::char_separator< char > > tokenizer
Just a typedef for simple use of the boost::tokenizer to split the lines.
int m_indexInitial
Maximum index of particles in each event that must be set as c_Initial (1-based).
int getEvent(MCParticleGraph &graph, double &weight)
Reads the next event and stores the result in the given MCParticle graph.
bool m_wrongSignPz
Bool to indicate that HER and LER were swapped.
int readEventHeader(double &eventWeight)
Reads the event header from the hepevt file.
double m_meanDecayLength
Mean lifetime*c of displaced particle.
Class to represent Particle data in graph.
Class to build, validate and sort a particle decay chain.
size_t size() const
Return the number of particles in the graph.
@ c_Initial
bit 5: Particle is initial such as e+ or e- and not going to Geant4
@ c_PrimaryParticle
bit 0: Particle is primary particle.
@ c_IsVirtual
bit 4: Particle is virtual and not going to Geant4.
GraphParticle & addParticle()
Add new particle to the graph.
Abstract base class for different kinds of events.