11 #ifndef Pythia8_Event_H
12 #define Pythia8_Event_H
15 #include "ParticleData.h"
16 #include "PythiaStdlib.h"
23 class ParticleDataEntry;
24 class ResonanceWidths;
36 Particle() : idSave(0), statusSave(0), mother1Save(0), mother2Save(0),
37 daughter1Save(0), daughter2Save(0), colSave(0), acolSave(0),
38 pSave(
Vec4(0.,0.,0.,0.)), mSave(0.), scaleSave(0.), polSave(9.),
39 hasVertexSave(
false), vProdSave(
Vec4(0.,0.,0.,0.)), tauSave(0.),
40 pdePtr(0), pdtPtr(0) { }
41 Particle(
int idIn,
int statusIn = 0,
int mother1In = 0,
42 int mother2In = 0,
int daughter1In = 0,
int daughter2In = 0,
43 int colIn = 0,
int acolIn = 0,
double pxIn = 0.,
double pyIn = 0.,
44 double pzIn = 0.,
double eIn = 0.,
double mIn = 0.,
45 double scaleIn = 0.,
double polIn = 9.)
46 : idSave(idIn), statusSave(statusIn), mother1Save(mother1In),
47 mother2Save(mother2In), daughter1Save(daughter1In),
48 daughter2Save(daughter2In), colSave(colIn), acolSave(acolIn),
49 pSave(
Vec4(pxIn, pyIn, pzIn, eIn)), mSave(mIn), scaleSave(scaleIn),
50 polSave(polIn), hasVertexSave(
false), vProdSave(
Vec4(0.,0.,0.,0.)),
51 tauSave(0.), pdePtr(0), pdtPtr(0) { }
52 Particle(
int idIn,
int statusIn,
int mother1In,
int mother2In,
53 int daughter1In,
int daughter2In,
int colIn,
int acolIn,
54 Vec4 pIn,
double mIn = 0.,
double scaleIn = 0.,
double polIn = 9.)
55 : idSave(idIn), statusSave(statusIn), mother1Save(mother1In),
56 mother2Save(mother2In), daughter1Save(daughter1In),
57 daughter2Save(daughter2In), colSave(colIn), acolSave(acolIn),
58 pSave(pIn), mSave(mIn), scaleSave(scaleIn), polSave(polIn),
59 hasVertexSave(
false), vProdSave(
Vec4(0.,0.,0.,0.)), tauSave(0.),
60 pdePtr(0), pdtPtr(0) { }
62 statusSave(pt.statusSave), mother1Save(pt.mother1Save),
63 mother2Save(pt.mother2Save), daughter1Save(pt.daughter1Save),
64 daughter2Save(pt.daughter2Save), colSave(pt.colSave),
65 acolSave(pt.acolSave), pSave(pt.pSave), mSave(pt.mSave),
66 scaleSave(pt.scaleSave), polSave(pt.polSave),
67 hasVertexSave(pt.hasVertexSave), vProdSave(pt.vProdSave),
68 tauSave(pt.tauSave), pdePtr(pt.pdePtr), pdtPtr(pt.pdtPtr) { }
70 idSave = pt.idSave; statusSave = pt.statusSave;
71 mother1Save = pt.mother1Save; mother2Save = pt.mother2Save;
72 daughter1Save = pt.daughter1Save; daughter2Save = pt.daughter2Save;
73 colSave = pt.colSave; acolSave = pt.acolSave; pSave = pt.pSave;
74 mSave = pt.mSave; scaleSave = pt.scaleSave; polSave = pt.polSave;
75 hasVertexSave = pt.hasVertexSave; vProdSave = pt.vProdSave;
76 tauSave = pt.tauSave; pdePtr = pt.pdePtr; pdtPtr = pt.pdtPtr; }
80 void setPDTPtr(
ParticleData* pdtPtrIn) { pdtPtr = pdtPtrIn; setPDEPtr();}
81 void setPDEPtr() {pdePtr = (pdtPtr != nullptr )
82 ? pdtPtr->particleDataEntryPtr( idSave) : 0;}
85 void id(
int idIn) {idSave = idIn; setPDEPtr();}
86 void status(
int statusIn) {statusSave = statusIn;}
87 void statusPos() {statusSave = abs(statusSave);}
88 void statusNeg() {statusSave = -abs(statusSave);}
89 void statusCode(
int statusIn) {statusSave =
90 (statusSave > 0) ? abs(statusIn) : -abs(statusIn);}
91 void mother1(
int mother1In) {mother1Save = mother1In;}
92 void mother2(
int mother2In) {mother2Save = mother2In;}
93 void mothers(
int mother1In = 0,
int mother2In = 0)
94 {mother1Save = mother1In; mother2Save = mother2In;}
95 void daughter1(
int daughter1In) {daughter1Save = daughter1In;}
96 void daughter2(
int daughter2In) {daughter2Save = daughter2In;}
97 void daughters(
int daughter1In = 0,
int daughter2In = 0)
98 {daughter1Save = daughter1In; daughter2Save = daughter2In;}
99 void col(
int colIn) {colSave = colIn;}
100 void acol(
int acolIn) {acolSave = acolIn;}
101 void cols(
int colIn = 0,
int acolIn = 0) {colSave = colIn;
103 void p(
Vec4 pIn) {pSave = pIn;}
104 void p(
double pxIn,
double pyIn,
double pzIn,
double eIn)
105 {pSave.p(pxIn, pyIn, pzIn, eIn);}
106 void px(
double pxIn) {pSave.px(pxIn);}
107 void py(
double pyIn) {pSave.py(pyIn);}
108 void pz(
double pzIn) {pSave.pz(pzIn);}
109 void e(
double eIn) {pSave.e(eIn);}
110 void m(
double mIn) {mSave = mIn;}
111 void scale(
double scaleIn) {scaleSave = scaleIn;}
112 void pol(
double polIn) {polSave = polIn;}
113 void vProd(
Vec4 vProdIn) {vProdSave = vProdIn; hasVertexSave =
true;}
114 void vProd(
double xProdIn,
double yProdIn,
double zProdIn,
double tProdIn)
115 {vProdSave.p(xProdIn, yProdIn, zProdIn, tProdIn); hasVertexSave =
true;}
116 void xProd(
double xProdIn) {vProdSave.px(xProdIn); hasVertexSave =
true;}
117 void yProd(
double yProdIn) {vProdSave.py(yProdIn); hasVertexSave =
true;}
118 void zProd(
double zProdIn) {vProdSave.pz(zProdIn); hasVertexSave =
true;}
119 void tProd(
double tProdIn) {vProdSave.e(tProdIn); hasVertexSave =
true;}
120 void tau(
double tauIn) {tauSave = tauIn;}
123 int id()
const {
return idSave;}
124 int status()
const {
return statusSave;}
125 int mother1()
const {
return mother1Save;}
126 int mother2()
const {
return mother2Save;}
127 int daughter1()
const {
return daughter1Save;}
128 int daughter2()
const {
return daughter2Save;}
129 int col()
const {
return colSave;}
130 int acol()
const {
return acolSave;}
131 Vec4 p()
const {
return pSave;}
132 double px()
const {
return pSave.px();}
133 double py()
const {
return pSave.py();}
134 double pz()
const {
return pSave.pz();}
135 double e()
const {
return pSave.e();}
136 double m()
const {
return mSave;}
137 double scale()
const {
return scaleSave;}
138 double pol()
const {
return polSave;}
139 bool hasVertex()
const {
return hasVertexSave;}
140 Vec4 vProd()
const {
return vProdSave;}
141 double xProd()
const {
return vProdSave.px();}
142 double yProd()
const {
return vProdSave.py();}
143 double zProd()
const {
return vProdSave.pz();}
144 double tProd()
const {
return vProdSave.e();}
145 double tau()
const {
return tauSave;}
148 int idAbs()
const {
return abs(idSave);}
149 int statusAbs()
const {
return abs(statusSave);}
150 bool isFinal()
const {
return (statusSave > 0);}
151 bool isRescatteredIncoming()
const {
return
152 (statusSave == -34 || statusSave == -45 ||
153 statusSave == -46 || statusSave == -54);}
156 double m2()
const {
return (mSave >= 0.) ? mSave*mSave
158 double mCalc()
const {
return pSave.mCalc();}
159 double m2Calc()
const {
return pSave.m2Calc();}
160 double eCalc()
const {
return sqrt(abs(m2() + pSave.pAbs2()));}
161 double pT()
const {
return pSave.pT();}
162 double pT2()
const {
return pSave.pT2();}
163 double mT()
const {
double temp = m2() + pSave.pT2();
164 return (temp >= 0.) ? sqrt(temp) : -sqrt(-temp);}
165 double mT2()
const {
return m2() + pSave.pT2();}
166 double pAbs()
const {
return pSave.pAbs();}
167 double pAbs2()
const {
return pSave.pAbs2();}
168 double eT()
const {
return pSave.eT();}
169 double eT2()
const {
return pSave.eT2();}
170 double theta()
const {
return pSave.theta();}
171 double phi()
const {
return pSave.phi();}
172 double thetaXZ()
const {
return pSave.thetaXZ();}
173 double pPos()
const {
return pSave.pPos();}
174 double pNeg()
const {
return pSave.pNeg();}
177 Vec4 vDec()
const {
return (tauSave > 0. && mSave > 0.)
178 ? vProdSave + tauSave * pSave / mSave : vProdSave;}
179 double xDec()
const {
return (tauSave > 0. && mSave > 0.)
180 ? vProdSave.px() + tauSave * pSave.px() / mSave : vProdSave.px();}
181 double yDec()
const {
return (tauSave > 0. && mSave > 0.)
182 ? vProdSave.py() + tauSave * pSave.py() / mSave : vProdSave.py();}
183 double zDec()
const {
return (tauSave > 0. && mSave > 0.)
184 ? vProdSave.pz() + tauSave * pSave.pz() / mSave : vProdSave.pz();}
185 double tDec()
const {
return (tauSave > 0. && mSave > 0.)
186 ? vProdSave.e() + tauSave * pSave.e() / mSave : vProdSave.e();}
189 string name()
const {
190 return (pdePtr !=
nullptr) ? pdePtr->name(idSave) :
" ";}
191 string nameWithStatus(
int maxLen = 20)
const;
192 int spinType()
const {
193 return (pdePtr !=
nullptr) ? pdePtr->spinType() : 0;}
194 int chargeType()
const {
195 return (pdePtr !=
nullptr) ? pdePtr->chargeType(idSave) : 0;}
196 double charge()
const {
197 return (pdePtr !=
nullptr) ? pdePtr->charge(idSave) : 0;}
198 bool isCharged()
const {
199 return (pdePtr !=
nullptr) ? (pdePtr->chargeType(idSave) != 0) :
false;}
200 bool isNeutral()
const {
201 return (pdePtr !=
nullptr) ? (pdePtr->chargeType(idSave) == 0) :
false;}
202 int colType()
const {
203 return (pdePtr !=
nullptr) ? pdePtr->colType(idSave) : 0;}
205 return (pdePtr !=
nullptr) ? pdePtr->m0() : 0.;}
206 double mWidth()
const {
207 return (pdePtr !=
nullptr) ? pdePtr->mWidth() : 0.;}
208 double mMin()
const {
209 return (pdePtr !=
nullptr) ? pdePtr->mMin() : 0.;}
210 double mMax()
const {
211 return (pdePtr !=
nullptr) ? pdePtr->mMax() : 0.;}
212 double mass()
const {
213 return (pdePtr !=
nullptr) ? pdePtr->mass() : 0.;}
214 double constituentMass()
const {
215 return (pdePtr !=
nullptr) ? pdePtr->constituentMass() : 0.;}
216 double tau0()
const {
217 return (pdePtr !=
nullptr) ? pdePtr->tau0() : 0.;}
218 bool mayDecay()
const {
219 return (pdePtr !=
nullptr) ? pdePtr->mayDecay() :
false;}
220 bool canDecay()
const {
221 return (pdePtr !=
nullptr) ? pdePtr->canDecay() :
false;}
222 bool doExternalDecay()
const {
223 return (pdePtr !=
nullptr) ? pdePtr->doExternalDecay() :
false;}
224 bool isResonance()
const {
225 return (pdePtr !=
nullptr) ? pdePtr->isResonance() :
false;}
226 bool isVisible()
const {
227 return (pdePtr !=
nullptr) ? pdePtr->isVisible() :
false;}
228 bool isLepton()
const {
229 return (pdePtr !=
nullptr) ? pdePtr->isLepton() :
false;}
230 bool isQuark()
const {
231 return (pdePtr !=
nullptr) ? pdePtr->isQuark() :
false;}
232 bool isGluon()
const {
233 return (pdePtr !=
nullptr) ? pdePtr->isGluon() :
false;}
234 bool isDiquark()
const {
235 return (pdePtr !=
nullptr) ? pdePtr->isDiquark() :
false;}
236 bool isParton()
const {
237 return (pdePtr !=
nullptr) ? pdePtr->isParton() :
false;}
238 bool isHadron()
const {
239 return (pdePtr !=
nullptr) ? pdePtr->isHadron() :
false;}
243 void rescale3(
double fac) {pSave.rescale3(fac);}
244 void rescale4(
double fac) {pSave.rescale4(fac);}
245 void rescale5(
double fac) {pSave.rescale4(fac); mSave *= fac;}
246 void rot(
double thetaIn,
double phiIn) {pSave.rot(thetaIn, phiIn);
247 if (hasVertexSave) vProdSave.rot(thetaIn, phiIn);}
248 void bst(
double betaX,
double betaY,
double betaZ) {
249 pSave.bst(betaX, betaY, betaZ);
250 if (hasVertexSave) vProdSave.bst(betaX, betaY, betaZ);}
251 void bst(
double betaX,
double betaY,
double betaZ,
double gamma) {
252 pSave.bst(betaX, betaY, betaZ, gamma);
253 if (hasVertexSave) vProdSave.bst(betaX, betaY, betaZ, gamma);}
254 void bst(
const Vec4& pBst) {pSave.bst(pBst);
255 if (hasVertexSave) vProdSave.bst(pBst);}
256 void bst(
const Vec4& pBst,
double mBst) {pSave.bst(pBst, mBst);
257 if (hasVertexSave) vProdSave.bst(pBst, mBst);}
258 void bstback(
const Vec4& pBst) {pSave.bstback(pBst);
259 if (hasVertexSave) vProdSave.bstback(pBst);}
260 void bstback(
const Vec4& pBst,
double mBst) {pSave.bstback(pBst, mBst);
261 if (hasVertexSave) vProdSave.bstback(pBst, mBst);}
263 if (hasVertexSave) vProdSave.rotbst(M);}
264 void offsetHistory(
int minMother,
int addMother,
int minDaughter,
266 void offsetCol(
int addCol);
271 static const double TINY;
274 int idSave, statusSave, mother1Save, mother2Save, daughter1Save,
275 daughter2Save, colSave, acolSave;
277 double mSave, scaleSave, polSave;
310 Junction() : remainsSave(
true), kindSave(0) {
311 for (
int j = 0; j < 3; ++j) {
312 colSave[j] = 0; endColSave[j] = 0; statusSave[j] = 0; } }
313 Junction(
int kindIn,
int col0In,
int col1In,
int col2In)
314 : remainsSave(
true), kindSave(kindIn) {colSave[0] = col0In;
315 colSave[1] = col1In; colSave[2] = col2In;
316 for (
int j = 0; j < 3; ++j) {
317 endColSave[j] = colSave[j]; statusSave[j] = 0; } }
319 kindSave(ju.kindSave) {
for (
int j = 0; j < 3; ++j) {
320 colSave[j] = ju.colSave[j]; endColSave[j] = ju.endColSave[j];
321 statusSave[j] = ju.statusSave[j]; } }
323 remainsSave = ju.remainsSave; kindSave = ju.kindSave;
324 for (
int j = 0; j < 3; ++j) { colSave[j] = ju.colSave[j];
325 endColSave[j] = ju.endColSave[j]; statusSave[j] = ju.statusSave[j]; } }
329 void remains(
bool remainsIn) {remainsSave = remainsIn;}
330 void col(
int j,
int colIn) {colSave[j] = colIn; endColSave[j] = colIn;}
331 void cols(
int j,
int colIn,
int endColIn) {colSave[j] = colIn;
332 endColSave[j] = endColIn;}
333 void endCol(
int j,
int endColIn) {endColSave[j] = endColIn;}
334 void status(
int j,
int statusIn) {statusSave[j] = statusIn;}
337 bool remains()
const {
return remainsSave;}
338 int kind()
const {
return kindSave;}
339 int col(
int j)
const {
return colSave[j];}
340 int endCol(
int j)
const {
return endColSave[j];}
341 int status(
int j)
const {
return statusSave[j];}
347 int kindSave, colSave[3], endColSave[3], statusSave[3];
360 Event(
int capacity = 100) {entry.reserve(capacity); startColTag = 100;
361 headerList =
"----------------------------------------";
362 particleDataPtr = 0;}
366 void init(
string headerIn =
"",
ParticleData* particleDataPtrIn = 0,
367 int startColTagIn = 100) {
368 headerList.replace(0, headerIn.length() + 2, headerIn +
" ");
369 particleDataPtr = particleDataPtrIn; startColTag = startColTagIn;}
372 void clear() {entry.resize(0); maxColTag = startColTag; scaleSave = 0.;
373 scaleSecondSave = 0.; clearJunctions();}
376 void reset() {clear(); append(90, -11, 0, 0, 0., 0., 0., 0., 0.);}
379 Particle& operator[](
int i) {
return entry[i];}
380 const Particle& operator[](
int i)
const {
return entry[i];}
383 Particle& at(
int i) {
return entry.at(i);}
386 int size()
const {
return entry.size();}
390 entry.push_back(entryIn); setPDTPtr();
391 if (entryIn.col() > maxColTag) maxColTag = entryIn.col();
392 if (entryIn.acol() > maxColTag) maxColTag = entryIn.acol();
393 return entry.size() - 1;
395 int append(
int id,
int status,
int mother1,
int mother2,
int daughter1,
396 int daughter2,
int col,
int acol,
double px,
double py,
double pz,
397 double e,
double m = 0.,
double scaleIn = 0.,
double polIn = 9.) {
398 entry.push_back(
Particle(
id, status, mother1, mother2, daughter1,
399 daughter2, col, acol, px, py, pz, e, m, scaleIn, polIn) ); setPDTPtr();
400 if (col > maxColTag) maxColTag = col;
401 if (acol > maxColTag) maxColTag = acol;
402 return entry.size() - 1;
404 int append(
int id,
int status,
int mother1,
int mother2,
int daughter1,
405 int daughter2,
int col,
int acol,
Vec4 p,
double m = 0.,
406 double scaleIn = 0.,
double polIn = 9.) {
407 entry.push_back(
Particle(
id, status, mother1, mother2, daughter1,
408 daughter2, col, acol, p, m, scaleIn, polIn) ); setPDTPtr();
409 if (col > maxColTag) maxColTag = col;
410 if (acol > maxColTag) maxColTag = acol;
411 return entry.size() - 1;
415 int append(
int id,
int status,
int col,
int acol,
double px,
double py,
416 double pz,
double e,
double m = 0.,
double scaleIn = 0.,
417 double polIn = 9.) { entry.push_back(
Particle(
id, status, 0, 0, 0, 0,
418 col, acol, px, py, pz, e, m, scaleIn, polIn) ); setPDTPtr();
419 if (col > maxColTag) maxColTag = col;
420 if (acol > maxColTag) maxColTag = acol;
421 return entry.size() - 1;
423 int append(
int id,
int status,
int col,
int acol,
Vec4 p,
double m = 0.,
424 double scaleIn = 0.,
double polIn = 9.) {entry.push_back(
Particle(
id,
425 status, 0, 0, 0, 0, col, acol, p, m, scaleIn, polIn) ); setPDTPtr();
426 if (col > maxColTag) maxColTag = col;
427 if (acol > maxColTag) maxColTag = acol;
428 return entry.size() - 1;
432 void setPDTPtr(
int iSet = -1) {
if (iSet < 0) iSet = entry.size() - 1;
433 entry[iSet].setPDTPtr( particleDataPtr);}
436 int copy(
int iCopy,
int newStatus = 0);
439 Particle& back() {
return entry.back();}
443 void list(ostream& os)
const;
444 void list(
bool showScaleAndVertex,
bool showMothersAndDaughters =
false)
446 void list(
bool showScaleAndVertex,
bool showMothersAndDaughters,
450 void popBack(
int nRemove = 1) {
if (nRemove ==1) entry.pop_back();
451 else {
int newSize = max( 0, size() - nRemove);
452 entry.resize(newSize);} }
456 void restorePtrs() {
for (
int i = 0; i < size(); ++i) setPDTPtr(i); }
459 void saveSize() {savedSize = entry.size();}
460 void restoreSize() {entry.resize(savedSize);}
463 void initColTag(
int colTag = 0) {maxColTag = max( colTag,startColTag);}
464 int lastColTag()
const {
return maxColTag;}
465 int nextColTag() {
return ++maxColTag;}
468 void scale(
double scaleIn) {scaleSave = scaleIn;}
469 double scale()
const {
return scaleSave;}
472 void scaleSecond(
double scaleSecondIn) {scaleSecondSave = scaleSecondIn;}
473 double scaleSecond()
const {
return scaleSecondSave;}
476 vector<int> motherList(
int i)
const;
477 vector<int> daughterList(
int i)
const;
480 int statusHepMC(
int i)
const;
483 int iTopCopy(
int i)
const;
484 int iBotCopy(
int i)
const;
487 int iTopCopyId(
int i)
const;
488 int iBotCopyId(
int i)
const;
491 vector<int> sisterList(
int i)
const;
492 vector<int> sisterListTopBot(
int i,
bool widenSearch =
true)
const;
495 bool isAncestor(
int i,
int iAncestor)
const;
498 void rot(
double theta,
double phi)
499 {
for (
int i = 0; i < size(); ++i) entry[i].rot(theta, phi);}
500 void bst(
double betaX,
double betaY,
double betaZ)
501 {
for (
int i = 0; i < size(); ++i) entry[i].bst(betaX, betaY, betaZ);}
502 void bst(
double betaX,
double betaY,
double betaZ,
double gamma)
503 {
for (
int i = 0; i < size(); ++i) entry[i].bst(betaX, betaY, betaZ,
505 void bst(
const Vec4& vec)
506 {
for (
int i = 0; i < size(); ++i) entry[i].bst(vec);}
508 {
for (
int i = 0; i < size(); ++i) entry[i].rotbst(M);}
511 void clearJunctions() {junction.resize(0);}
514 void appendJunction(
int kind,
int col0,
int col1,
int col2)
515 { junction.push_back(
Junction( kind, col0, col1, col2) );}
516 void appendJunction(
Junction junctionIn) {junction.push_back(junctionIn);}
517 int sizeJunction()
const {
return junction.size();}
518 bool remainsJunction(
int i)
const {
return junction[i].remains();}
519 void remainsJunction(
int i,
bool remainsIn) {junction[i].remains(remainsIn);}
520 int kindJunction(
int i)
const {
return junction[i].kind();}
521 int colJunction(
int i,
int j)
const {
return junction[i].col(j);}
522 void colJunction(
int i,
int j,
int colIn) {junction[i].col(j, colIn);}
523 int endColJunction(
int i,
int j)
const {
return junction[i].endCol(j);}
524 void endColJunction(
int i,
int j,
int endColIn)
525 {junction[i].endCol(j, endColIn);}
526 int statusJunction(
int i,
int j)
const {
return junction[i].status(j);}
527 void statusJunction(
int i,
int j,
int statusIn)
528 {junction[i].status(j, statusIn);}
529 Junction& getJunction(
int i) {
return junction[i];}
530 const Junction& getJunction(
int i)
const {
return junction[i];}
531 void eraseJunction(
int i);
534 void saveJunctionSize() {savedJunctionSize = junction.size();}
535 void restoreJunctionSize() {junction.resize(savedJunctionSize);}
538 void listJunctions(ostream& os = cout)
const;
547 static const int IPERLINE;
554 vector<Pythia8::Particle> entry;
558 vector<Pythia8::Junction> junction;
564 int savedSize, savedJunctionSize;
567 double scaleSave, scaleSecondSave;
582 #endif // end Pythia8_Event_H