Belle II Software development
LHEReader.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
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>
14
15#include <string>
16#include <boost/lexical_cast.hpp>
17#include <boost/algorithm/string.hpp>
18#include <boost/format.hpp>
19
20#include <TF1.h>
21
22using namespace std;
23using namespace Belle2;
24
25const boost::char_separator<char> LHEReader::sep(",; \t");
26
27void LHEReader::open(const string& filename)
28{
29 m_lineNr = 0;
30 m_input.open(filename.c_str());
31 if (!m_input) throw(LHECouldNotOpenFileError() << filename);
32}
33
34int LHEReader::getEvent(MCParticleGraph& graph, double& eventWeight)
35{
36 int nparticles = readEventHeader(eventWeight);
37 if (nparticles <= 0) {
38 throw (LHEEmptyEventError() << m_lineNr << nparticles);
39 }
40
41 int first = graph.size();
42 //Make list of particles
43 for (int i = 0; i < nparticles; i++) {
44 graph.addParticle();
45 }
46
47 double r, x = 0, y = 0, z = 0, t = 0;
48 //Read particles from file
49 for (int i = 0; i < nparticles; ++i) {
50 MCParticleGraph::GraphParticle& p = graph[first + i];
51 int mother = readParticle(p);
52
53 // add the mother
54 if (mother > 0) {
55 MCParticleGraph::GraphParticle* q = &graph[mother - 1];
56 p.comesFrom(*q);
57 }
58
59 // if positron goes to positive z-direction, we have to rotate along y-axis by 180deg
60 if (m_wrongSignPz) {
61 ROOT::Math::PxPyPzEVector p4 = p.get4Vector();
62 p4.SetPz(-1.0 * p4.Pz());
63 p4.SetPx(-1.0 * p4.Px());
64 p.set4Vector(p4);
65 }
66
67 //move vertex position of selected particle and its daughters
68 if (m_meanDecayLength > 0) {
69 if (p.getPDG() == m_pdgDisplaced) {
70 TF1 fr("fr", "exp(-x/[0])", 0, 1000000);
71 ROOT::Math::PxPyPzEVector p4 = p.get4Vector();
72 fr.SetRange(m_Rmin, m_Rmax);
73 fr.SetParameter(0, m_meanDecayLength * p4.Gamma());
74 r = fr.GetRandom();
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);
79 t = (r / Const::speedOfLight) * (p4.E() / p4.P());
80 p.setDecayTime(t);
81 p.setValidVertex(true);
82 }
83
84 if (mother > 0) {
85 if (graph[mother - 1].getPDG() == m_pdgDisplaced) {
86 p.setProductionVertex(x, y, z);
87 p.setProductionTime(t);
88 p.setValidVertex(true);
89 }
90 }
91 }
92
93 // initial 2 (e+/e-), virtual 3 (Z/gamma*)
94 // check if particle should be made virtual according to steering options:
95 if (i < m_indexVirtual && i >= m_indexInitial)
96 p.addStatus(MCParticle::c_IsVirtual);
97
98 if (i < m_indexInitial)
99 p.addStatus(MCParticle::c_Initial);
100
101 if (m_indexVirtual < m_indexInitial) B2WARNING("IsVirtual particle requested but is overwritten by Initial");
102
103 }
104// return eventID;
105 return -1;
106}
107
108
110{
111// int eventID;
112 double weight;
113 for (int i = 0; i < n; i++) {
114// int nparticles = readEventHeader(eventID, weight);
115 int nparticles = readEventHeader(weight);
116 if (nparticles < 0) return false;
117 for (int j = 0; j < nparticles; j++) getLine();
118 }
119 return true;
120}
121
122//===================================================================
123// Protected methods
124//===================================================================
125
127{
128 std::string line;
129 do {
130 getline(m_input, line);
131 m_lineNr++;
132 size_t commentPos = line.find_first_of('#');
133 if (commentPos != string::npos) {
134 line = line.substr(0, commentPos);
135 }
136 boost::trim(line);
137
138 } while (line == "" && !m_input.eof());
139
140 return line;
141}
142
143int LHEReader::readEventHeader(double& eventWeight)
144{
145
146 // Search for next <event>
147 std::string line2;
148 do {
149 getline(m_input, line2);
150 m_lineNr++;
151 size_t commentPos = line2.find_first_of('#');
152 if (commentPos != string::npos) {
153 line2 = line2.substr(0, commentPos);
154 }
155 boost::trim(line2);
156
157 } while (line2 != "<event>" && !m_input.eof());
158
159 //Get number of particles from file
160 int nparticles = -1;
161 string line = getLine();
162
163 if (line == "" || m_input.eof()) return -1;
164
165 vector<double> fields;
166 fields.reserve(15);
167
168 tokenizer tokens(line, sep);
169 int index(0);
170
171 for (const string& tok : tokens) {
172 ++index;
173 try {
174 fields.push_back(boost::lexical_cast<double>(tok));
175 } catch (boost::bad_lexical_cast& e) {
176 throw (LHEConvertFieldError() << m_lineNr << index << tok);
177 }
178 }
179
180 if (fields.size() < 3) {
181 nparticles = static_cast<int>(fields[0]);
182 eventWeight = 1.0;
183 } else {
184 nparticles = static_cast<int>(fields[0]);
185 eventWeight = static_cast<double>(fields[2]);
186 }
187 return nparticles;
188}
189
191{
192 int mother = -1;
193
194 string line = getLine();
195 vector<double> fields;
196 fields.reserve(15);
197
198 tokenizer tokens(line, sep);
199 int index(0);
200
201 for (const string& tok : tokens) {
202 ++index;
203 try {
204 fields.push_back(boost::lexical_cast<double>(tok));
205 } catch (boost::bad_lexical_cast& e) {
206 throw (LHEConvertFieldError() << m_lineNr << index << tok);
207 }
208 }
209
210 switch (fields.size()) {
211 case 13:
212 particle.addStatus(MCParticle::c_PrimaryParticle);
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]);
218 break;
219 default:
220 throw (LHEParticleFormatError() << m_lineNr << fields.size());
221 }
222
223 return mother;
224}
static const double speedOfLight
[cm/ns]
Definition Const.h:695
int readParticle(MCParticleGraph::GraphParticle &particle)
Reads the information for a single particle from the LHE file.
Definition LHEReader.cc:190
bool skipEvents(int n)
Skips a given number of events.
Definition LHEReader.cc:109
std::string getLine()
Returns the current line from the LHE ascii file.
Definition LHEReader.cc:126
double m_Rmin
Minimum of vertex distance to IP.
Definition LHEReader.h:104
int m_pdgDisplaced
PDG code of the displaced particle being studied.
Definition LHEReader.h:106
static const boost::char_separator< char > sep
The characters at which to split, defaults to ",; \t".
Definition LHEReader.h:112
int m_lineNr
The current line number within the ascii file.
Definition LHEReader.h:114
int m_indexVirtual
Maximum index of particles in each event that must be set as c_IsVirtual (1-based).
Definition LHEReader.h:142
double m_Rmax
Maximum of vertex distance to IP.
Definition LHEReader.h:105
void open(const std::string &filename)
Opens an ascii file and prepares it for reading.
Definition LHEReader.cc:27
std::ifstream m_input
The input stream of the ascii file.
Definition LHEReader.h:115
boost::tokenizer< boost::char_separator< char > > tokenizer
Just a typedef for simple use of the boost::tokenizer to split the lines.
Definition LHEReader.h:110
int m_indexInitial
Maximum index of particles in each event that must be set as c_Initial (1-based).
Definition LHEReader.h:140
int getEvent(MCParticleGraph &graph, double &weight)
Reads the next event and stores the result in the given MCParticle graph.
Definition LHEReader.cc:34
bool m_wrongSignPz
Bool to indicate that HER and LER were swapped.
Definition LHEReader.h:102
int readEventHeader(double &eventWeight)
Reads the event header from the hepevt file.
Definition LHEReader.cc:143
double m_meanDecayLength
Mean lifetime*c of displaced particle.
Definition LHEReader.h:103
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
Definition MCParticle.h:57
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition MCParticle.h:47
@ c_IsVirtual
bit 4: Particle is virtual and not going to Geant4.
Definition MCParticle.h:55
GraphParticle & addParticle()
Add new particle to the graph.
Abstract base class for different kinds of events.
STL namespace.