88918fae5cef2d249dbe7dd6e99f66db9453c475
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / include / Event.h
1 // Event.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2008 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5
6 // Header file for the Particle and Event classes.
7 // Particle: information on an instance of a particle.
8 // Junction: information on a junction between three colours.
9 // Event: list of particles in the current event.
10
11 #ifndef Pythia8_Event_H
12 #define Pythia8_Event_H
13
14 #include "Basics.h"
15 #include "ParticleData.h"
16 #include "PythiaStdlib.h"
17 #include "Settings.h"
18
19 namespace Pythia8 {
20
21 //**************************************************************************
22
23 // Forward references to ParticleData and ResonanceWidths classes.
24 class ParticleDataEntry;
25 class ResonanceWidths;
26
27 //**************************************************************************
28
29 // Particle class.
30 // This class holds info on a particle in general.
31
32 class Particle {
33
34 public:
35
36   // Constructors.
37   Particle() : idSave(0), statusSave(0), mother1Save(0), mother2Save(0), 
38     daughter1Save(0), daughter2Save(0), colSave(0), acolSave(0), 
39     pSave(Vec4(0.,0.,0.,0.)), mSave(0.), scaleSave(0.), 
40     hasVertexSave(false), vProdSave(Vec4(0.,0.,0.,0.)), tauSave(0.), 
41     particlePtr(0) { }  
42   Particle(int idIn, int statusIn = 0, int mother1In = 0, 
43     int mother2In = 0, int daughter1In = 0, int daughter2In = 0,
44     int colIn = 0, int acolIn = 0, double pxIn = 0., double pyIn = 0., 
45     double pzIn = 0., double eIn = 0., double mIn = 0., double scaleIn = 0.) 
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     hasVertexSave(false), vProdSave(Vec4(0.,0.,0.,0.)), tauSave(0.) 
51     {setParticlePtr();}  
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.) 
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), hasVertexSave(false), 
59     vProdSave(Vec4(0.,0.,0.,0.)), tauSave(0.) {setParticlePtr();}  
60   Particle(const Particle& pt) : idSave(pt.idSave), 
61     statusSave(pt.statusSave), mother1Save(pt.mother1Save), 
62     mother2Save(pt.mother2Save), daughter1Save(pt.daughter1Save), 
63     daughter2Save(pt.daughter2Save), colSave(pt.colSave), 
64     acolSave(pt.acolSave), pSave(pt.pSave), mSave(pt.mSave), 
65     scaleSave(pt.scaleSave), hasVertexSave(pt.hasVertexSave), 
66     vProdSave(pt.vProdSave), tauSave(pt.tauSave), 
67     particlePtr(pt.particlePtr) { } 
68   Particle& operator=(const Particle& pt) {if (this != &pt) {
69     idSave = pt.idSave; statusSave = pt.statusSave; 
70     mother1Save = pt.mother1Save; mother2Save = pt.mother2Save; 
71     daughter1Save = pt.daughter1Save; daughter2Save = pt.daughter2Save; 
72     colSave = pt.colSave; acolSave = pt.acolSave; pSave = pt.pSave; 
73     mSave = pt.mSave; scaleSave = pt.scaleSave; 
74     hasVertexSave = pt.hasVertexSave; vProdSave = pt.vProdSave; 
75     tauSave = pt.tauSave; particlePtr = pt.particlePtr; } return *this; } 
76       
77   // Member functions for input.
78   void id(int idIn) {idSave = idIn; setParticlePtr();}
79   void status(int statusIn) {statusSave = statusIn;}
80   void statusPos() {statusSave = abs(statusSave);}
81   void statusNeg() {statusSave = -abs(statusSave);}
82   void statusCode(int statusIn) {statusSave = 
83     (statusSave > 0) ? abs(statusIn) : -abs(statusIn);}
84   void mother1(int mother1In) {mother1Save = mother1In;}
85   void mother2(int mother2In) {mother2Save = mother2In;}
86   void mothers(int mother1In = 0, int mother2In = 0) 
87     {mother1Save = mother1In; mother2Save = mother2In;}
88   void daughter1(int daughter1In) {daughter1Save = daughter1In;}
89   void daughter2(int daughter2In) {daughter2Save = daughter2In;}
90   void daughters(int daughter1In = 0, int daughter2In = 0) 
91     {daughter1Save = daughter1In; daughter2Save = daughter2In;}  
92   void col(int colIn) {colSave = colIn;}
93   void acol(int acolIn) {acolSave = acolIn;}
94   void cols(int colIn = 0,int acolIn = 0) {colSave = colIn; 
95     acolSave = acolIn;}  
96   void p(Vec4 pIn) {pSave = pIn;}
97   void p(double pxIn, double pyIn, double pzIn, double eIn) 
98     {pSave.p(pxIn, pyIn, pzIn, eIn);}
99   void px(double pxIn) {pSave.px(pxIn);}
100   void py(double pyIn) {pSave.py(pyIn);}
101   void pz(double pzIn) {pSave.pz(pzIn);}
102   void e(double eIn) {pSave.e(eIn);}
103   void m(double mIn) {mSave = mIn;}
104   void scale(double scaleIn) {scaleSave = scaleIn;}
105   void vProd(Vec4 vProdIn) {vProdSave = vProdIn; hasVertexSave = true;}
106   void vProd(double xProdIn, double yProdIn, double zProdIn, double tProdIn)
107     {vProdSave.p(xProdIn, yProdIn, zProdIn, tProdIn); hasVertexSave = true;}
108   void xProd(double xProdIn) {vProdSave.px(xProdIn); hasVertexSave = true;} 
109   void yProd(double yProdIn) {vProdSave.py(yProdIn); hasVertexSave = true;} 
110   void zProd(double zProdIn) {vProdSave.pz(zProdIn); hasVertexSave = true;} 
111   void tProd(double tProdIn) {vProdSave.e(tProdIn); hasVertexSave = true;} 
112   void tau(double tauIn) {tauSave = tauIn;} 
113
114   // Member functions for output.
115   int    id()        const {return idSave;}
116   int    status()    const {return statusSave;}
117   int    mother1()   const {return mother1Save;}
118   int    mother2()   const {return mother2Save;}
119   int    daughter1() const {return daughter1Save;}
120   int    daughter2() const {return daughter2Save;}
121   int    col()       const {return colSave;}
122   int    acol()      const {return acolSave;}
123   Vec4   p()         const {return pSave;}
124   double px()        const {return pSave.px();}
125   double py()        const {return pSave.py();}
126   double pz()        const {return pSave.pz();}
127   double e()         const {return pSave.e();}
128   double m()         const {return mSave;}
129   double scale()     const {return scaleSave;}
130   bool   hasVertex() const {return hasVertexSave;}
131   Vec4   vProd()     const {return vProdSave;}
132   double xProd()     const {return vProdSave.px();}
133   double yProd()     const {return vProdSave.py();}
134   double zProd()     const {return vProdSave.pz();}
135   double tProd()     const {return vProdSave.e();}
136   double tau()       const {return tauSave;}
137
138   // Member functions for output; derived int and bool quantities.
139   int    idAbs()     const {return abs(idSave);}
140   int    statusAbs() const {return abs(statusSave);}
141   bool   isFinal()   const {return (statusSave > 0);}
142
143   // Member functions for output; derived double quantities.
144   double m2()        const {return mSave*mSave;}
145   double mCalc()     const {return pSave.mCalc();}
146   double m2Calc()    const {return pSave.m2Calc();}
147   double eCalc()     const {return sqrt(mSave*mSave + pSave.pAbs2());}
148   double pT()        const {return pSave.pT();}
149   double pT2()       const {return pSave.pT2();}
150   double mT()        const {return sqrt(mSave*mSave + pSave.pT2());}
151   double mT2()       const {return mSave*mSave + pSave.pT2();}
152   double pAbs()      const {return pSave.pAbs();}
153   double pAbs2()     const {return pSave.pAbs2();}
154   double theta()     const {return pSave.theta();}
155   double phi()       const {return pSave.phi();}
156   double thetaXZ()   const {return pSave.thetaXZ();}
157   double pPlus()     const {return pSave.pPlus();}
158   double pMinus()    const {return pSave.pMinus();}
159   double y()         const;
160   double eta()       const; 
161   Vec4   vDec()      const {return (tauSave > 0. && mSave > 0.) 
162     ? vProdSave + tauSave * pSave / mSave : vProdSave;}
163   double xDec()      const {return (tauSave > 0. && mSave > 0.) 
164     ? vProdSave.px() + tauSave * pSave.px() / mSave : vProdSave.px();}
165   double yDec()      const {return (tauSave > 0. && mSave > 0.)  
166     ? vProdSave.py() + tauSave * pSave.py() / mSave : vProdSave.py();}
167   double zDec()      const {return (tauSave > 0. && mSave > 0.)  
168     ? vProdSave.pz() + tauSave * pSave.pz() / mSave : vProdSave.pz();}
169   double tDec()      const {return (tauSave > 0. && mSave > 0.)  
170     ? vProdSave.e()  + tauSave * pSave.e()  / mSave : vProdSave.e();}
171
172   // Further output, based on a pointer to a ParticleDataEntry object.
173   string name()      const {return particlePtr->name(idSave);}
174   string nameWithStatus(int maxLen = 20) const;
175   int    spinType()  const {return particlePtr->spinType();}
176   int    chargeType() const {return particlePtr->chargeType(idSave);}
177   double charge()    const {return  particlePtr->charge(idSave);}
178   bool   isCharged() const {return (particlePtr->chargeType(idSave) != 0);}
179   bool   isNeutral() const {return (particlePtr->chargeType(idSave) == 0);}
180   int    colType()   const {return particlePtr->colType(idSave);}
181   double m0()        const {return particlePtr->m0();}
182   double mWidth()    const {return particlePtr->mWidth();}
183   double mMin()      const {return particlePtr->mMin();}
184   double mMax()      const {return particlePtr->mMax();}
185   double mass()      const {return particlePtr->mass();}
186   double constituentMass() const {return particlePtr->constituentMass();}
187   double tau0()      const {return particlePtr->tau0();}
188   bool   mayDecay()  const {return particlePtr->mayDecay();}
189   bool   canDecay()  const {return particlePtr->canDecay();}
190   bool   doExternalDecay() const {return particlePtr->doExternalDecay();}
191   bool   isResonance() const {return particlePtr->isResonance();}
192   bool   isVisible() const {return particlePtr->isVisible();}
193   bool   isLepton()  const {return particlePtr->isLepton();}
194   bool   isQuark()   const {return particlePtr->isQuark();}
195   bool   isGluon()   const {return particlePtr->isGluon();}
196   bool   isHadron()  const {return particlePtr->isHadron();}
197   ParticleDataEntry& particleData() const {return *particlePtr;}
198
199   // Member functions that perform operations.
200   void rescale3(double fac) {pSave.rescale3(fac);}
201   void rescale4(double fac) {pSave.rescale4(fac);}
202   void rescale5(double fac) {pSave.rescale4(fac); mSave *= fac;}
203   void rot(double thetaIn, double phiIn) {pSave.rot(thetaIn, phiIn);
204     if (hasVertexSave) vProdSave.rot(thetaIn, phiIn);} 
205   void bst(double betaX, double betaY, double betaZ) {
206     pSave.bst(betaX, betaY, betaZ);
207     if (hasVertexSave) vProdSave.bst(betaX, betaY, betaZ);}
208   void bst(double betaX, double betaY, double betaZ, double gamma) {
209     pSave.bst(betaX, betaY, betaZ, gamma);
210     if (hasVertexSave) vProdSave.bst(betaX, betaY, betaZ, gamma);}
211   void bst(const Vec4& pBst) {pSave.bst(pBst);
212     if (hasVertexSave) vProdSave.bst(pBst);}
213   void bst(const Vec4& pBst, double mBst) {pSave.bst(pBst, mBst);
214     if (hasVertexSave) vProdSave.bst(pBst, mBst);}
215   void bstback(const Vec4& pBst) {pSave.bstback(pBst);
216     if (hasVertexSave) vProdSave.bstback(pBst);}
217   void bstback(const Vec4& pBst, double mBst) {pSave.bstback(pBst, mBst);
218     if (hasVertexSave) vProdSave.bstback(pBst, mBst);}
219   void rotbst(const RotBstMatrix& M) {pSave.rotbst(M);
220     if (hasVertexSave) vProdSave.rotbst(M);} 
221   void offsetHistory( int minMother, int addMother, int minDaughter, 
222     int addDaughter);
223   void offsetCol( int addCol); 
224
225   // Member function to set the ParticleDataEntry pointer, using idSave.
226   void setParticlePtr();
227
228 private:
229
230   // Constants: could only be changed in the code itself.
231   static const double TINY;
232
233   // Properties of the current particle.
234   int idSave, statusSave, mother1Save, mother2Save, daughter1Save, 
235     daughter2Save, colSave, acolSave;
236   Vec4 pSave;
237   double mSave, scaleSave;
238   bool hasVertexSave;
239   Vec4 vProdSave;
240   double tauSave;
241
242   // Pointer to properties of the particle species.
243   // Should no be saved in a persistent copy of the event record.
244   // The //! below is ROOT notation that this member should not be saved.
245   // Event::restorePtrs() can be called to restore the missing information. 
246   ParticleDataEntry* particlePtr;  //!
247
248 };
249
250 // Invariant mass of a pair and its square.
251 // (Not part of class proper, but tightly linked.)
252 double m(const Particle&, const Particle&); 
253 double m2(const Particle&, const Particle&); 
254
255 //**************************************************************************
256
257 // The juction class stores what kind of junction it is, the colour indices 
258 // of the legs at the junction and as far out as legs have been traced,
259 // and the status codes assigned for fragmentation of each leg.
260
261 class Junction {
262
263 public:
264
265   // Constructors.
266   Junction() : remainsSave(true), kindSave(0) { 
267     for (int j = 0; j < 3; ++j) {
268     colSave[j] = 0; endColSave[j] = 0; statusSave[j] = 0; } }
269   Junction( int kindIn, int col0In, int col1In, int col2In) 
270     : remainsSave(true), kindSave(kindIn) {colSave[0] = col0In; 
271     colSave[1] = col1In; colSave[2] = col2In; 
272     for (int j = 0; j < 3; ++j) {
273     endColSave[j] = colSave[j]; statusSave[j] = 0; } }
274   Junction(const Junction& ju) : remainsSave(ju.remainsSave), 
275     kindSave(ju.kindSave) { for (int j = 0; j < 3; ++j) {
276     colSave[j] = ju.colSave[j]; endColSave[j] = ju.endColSave[j]; 
277     statusSave[j] = ju.statusSave[j]; } }
278   Junction& operator=(const Junction& ju) {if (this != &ju) { 
279     remainsSave = ju.remainsSave; kindSave =  ju.kindSave; 
280     for (int j = 0; j < 3; ++j) { colSave[j] = ju.colSave[j]; 
281     endColSave[j] = ju.endColSave[j]; statusSave[j] = ju.statusSave[j]; } } 
282     return *this; }  
283
284   // Set values.
285   void remains(bool remainsIn) {remainsSave = remainsIn;}
286   void col(int j, int colIn) {colSave[j] = colIn; endColSave[j] = colIn;}
287   void cols(int j, int colIn, int endColIn) {colSave[j] = colIn; 
288     endColSave[j] = endColIn;}
289   void endCol(int j, int endColIn) {endColSave[j] = endColIn;}
290   void status(int j, int statusIn) {statusSave[j] = statusIn;}
291
292   // Read out value.
293   bool   remains()     const {return remainsSave;}
294   int    kind()        const {return kindSave;}
295   int    col(int j)    const {return colSave[j];}
296   int    endCol(int j) const {return endColSave[j];}
297   int    status(int j) const {return statusSave[j];}
298  
299 private:
300
301   // Kind, positions of the three ends and their status codes.
302   bool remainsSave;
303   int kindSave, colSave[3], endColSave[3], statusSave[3];
304
305 };
306
307 //**************************************************************************
308
309 // The Event class holds all info on the generated event.
310
311 class Event {
312     
313 public:
314
315   // Constructors.
316   Event(int capacity = 100) {entry.reserve(capacity); startColTag = 100; 
317     headerList = "----------------------------------------";}
318   Event& operator=(const Event& oldEvent);
319
320   // Initialize colour and header specification for event listing.
321   void init( string headerIn = "");
322
323   // Clear event record.
324   void clear() {entry.resize(0); maxColTag = startColTag; 
325     clearJunctions(); clearSystems();}
326
327   // Clear event record, and set first particle empty.
328   void reset() {clear(); append(90, -11, 0, 0, 0., 0., 0., 0., 0.);}
329
330   // Overload index operator to access element of event record.
331   Particle& operator[](int i) {return entry[i];}
332   const Particle& operator[](int i) const {return entry[i];}
333
334   // Event record size.
335   int size() const {return entry.size();}
336
337   // Put a new particle at the end of the event record; return index.
338   int append(Particle entryIn) {    
339     entry.push_back(entryIn); 
340     if (entryIn.col() > maxColTag) maxColTag = entryIn.col();   
341     if (entryIn.acol() > maxColTag) maxColTag = entryIn.acol();
342     return entry.size() - 1;
343   }
344   int append(int id, int status, int mother1, int mother2, int daughter1, 
345     int daughter2, int col, int acol, double px, double py, double pz, 
346     double e, double m = 0., double scaleIn = 0.) {entry.push_back( 
347     Particle(id, status, mother1, mother2, daughter1, daughter2, col, acol, 
348     px, py, pz, e, m, scaleIn) ); 
349     if (col > maxColTag) maxColTag = col;   
350     if (acol > maxColTag) maxColTag = acol;
351     return entry.size() - 1;
352   }
353   int append(int id, int status, int mother1, int mother2, int daughter1, 
354     int daughter2, int col, int acol, Vec4 p, double m = 0., 
355     double scaleIn = 0.) {entry.push_back( Particle(id, status, mother1, 
356     mother2, daughter1, daughter2, col, acol, p, m, scaleIn) ); 
357     if (col > maxColTag) maxColTag = col;   
358     if (acol > maxColTag) maxColTag = acol;
359     return entry.size() - 1;
360   }
361
362   // Brief versions of append: no mothers and no daughters.
363   int append(int id, int status, int col, int acol, double px, double py, 
364     double pz, double e, double m = 0.) {entry.push_back( Particle(id, 
365     status, 0, 0, 0, 0, col, acol, px, py, pz, e, m, 0.) ); 
366     if (col > maxColTag) maxColTag = col;   
367     if (acol > maxColTag) maxColTag = acol;
368     return entry.size() - 1;
369   }
370   int append(int id, int status, int col, int acol, Vec4 p, double m = 0.) 
371     {entry.push_back( Particle(id, status, 0, 0, 0, 0, col, acol, p, m, 0.) ); 
372     if (col > maxColTag) maxColTag = col;   
373     if (acol > maxColTag) maxColTag = acol;
374     return entry.size() - 1;
375   }
376
377   // Add a copy of an existing particle at the end of the event record.
378   int copy(int iCopy, int newStatus = 0);
379
380   // Implement reference "back" to access last element.
381   Particle& back() {return entry.back();}
382
383   // List the particles in an event.
384   void list(ostream& os = cout) {list(false, false, os);}  
385   void list(bool showScaleAndVertex, bool showMothersAndDaughters = false, 
386     ostream& os = cout);  
387
388   // Remove last n entries.
389   void popBack(int nRemove = 1) { if (nRemove ==1) entry.pop_back();
390     else {int newSize = max( 0, size() - nRemove); 
391     entry.resize(newSize);} } 
392
393   // Restore all ParticleDataEntry* pointers in the Particle vector.
394   // Useful when a persistent copy of the event record is read back in.
395   void restorePtrs() { 
396     for (int i = 0; i < size(); ++i) entry[i].setParticlePtr(); } 
397
398   // Save or restore the size of the event record (throwing at the end).
399   void saveSize() {savedSize = entry.size();}
400   void restoreSize() {entry.resize(savedSize);}   
401
402   // Initialize and access colour tag information.
403   void initColTag(int colTag = 0) {maxColTag = max( colTag,startColTag);}
404   int lastColTag() const {return maxColTag;}
405   int nextColTag() {return ++maxColTag;}
406
407   // Access scale for which event as a whole is defined.
408   void scale( double scaleIn) {scaleSave = scaleIn;}
409   double scale() const {return scaleSave;}
410
411   // Need a second scale if two hard interactions in event.
412   void scaleSecond( double scaleSecondIn) {scaleSecondSave = scaleSecondIn;}
413   double scaleSecond() const {return scaleSecondSave;}
414
415   // Find complete list of daughters and mothers.
416   vector<int> motherList(int i) const;
417   vector<int> daughterList(int i) const;
418  
419   // Trace the first and last copy of one and the same particle.
420   int iTopCopy(int i) const;
421   int iBotCopy(int i) const;
422
423   // Trace the first and last copy of a particle, using flavour match.
424   int iTopCopyId(int i) const;
425   int iBotCopyId(int i) const;
426
427   // Find list of sisters, also tracking up and down identical copies.
428   vector<int> sisterList(int i) const;
429   vector<int> sisterListTopBot(int i, bool widenSearch = true) const;
430
431   // Check whether two particles have a direct mother-daughter relation.
432   bool isAncestor(int i, int iAncestor) const;
433
434   // Member functions for rotations and boosts of an event.
435   void rot(double theta, double phi) 
436     {for (int i = 0; i < size(); ++i) entry[i].rot(theta, phi);} 
437   void bst(double betaX, double betaY, double betaZ) 
438     {for (int i = 0; i < size(); ++i) entry[i].bst(betaX, betaY, betaZ);}
439   void bst(double betaX, double betaY, double betaZ, double gamma) 
440     {for (int i = 0; i < size(); ++i) entry[i].bst(betaX, betaY, betaZ, 
441     gamma);}
442   void bst(const Vec4& vec) 
443     {for (int i = 0; i < size(); ++i) entry[i].bst(vec);}
444   void rotbst(const RotBstMatrix& M) 
445     {for (int i = 0; i < size(); ++i) entry[i].rotbst(M);}
446
447   // Clear the list of junctions.
448   void clearJunctions() {junction.resize(0);}
449  
450   // Add a junction to the list, study it or extra input.
451   void appendJunction( int kind, int col0, int col1, int col2)  
452     { junction.push_back( Junction( kind, col0, col1, col2) );} 
453   void appendJunction(Junction junctionIn) {junction.push_back(junctionIn);} 
454   int sizeJunction() const {return junction.size();}
455   bool remainsJunction(int i) const {return junction[i].remains();}
456   void remainsJunction(int i, bool remainsIn) {junction[i].remains(remainsIn);}
457   int kindJunction(int i) const {return junction[i].kind();}
458   int colJunction( int i, int j) const {return junction[i].col(j);}
459   void colJunction( int i, int j, int colIn) {junction[i].col(j, colIn);}
460   int endColJunction( int i, int j) const {return junction[i].endCol(j);}
461   void endColJunction( int i, int j, int endColIn) 
462     {junction[i].endCol(j, endColIn);}
463   int statusJunction( int i, int j) const {return junction[i].status(j);}
464   void statusJunction( int i, int j, int statusIn) 
465     {junction[i].status(j, statusIn);}
466   Junction& getJunction(int i) {return junction[i];}
467   const Junction& getJunction(int i) const {return junction[i];}
468   void eraseJunction(int i);
469
470   // Save or restore the size of the junction list (throwing at the end).
471   void saveJunctionSize() {savedJunctionSize = junction.size();}
472   void restoreJunctionSize() {junction.resize(savedJunctionSize);}   
473
474   // List any junctions in the event; for debug mainly.
475   void listJunctions(ostream& os = cout) const;
476
477   // Operations with grouped systems of partons for internal use only.
478   // (Used by combined MI, ISR, FSR and BR machinery in PartonLevel.)
479
480   // Reset all systems and system number to empty.
481   void clearSystems() {beginSys.resize(0); sizeSys.resize(0); 
482     memberSys.resize(0);}
483   
484   // Get number of systems or number of members in a system. 
485   int sizeSystems() const {return beginSys.size();}
486   int sizeSystem(int iSys) const {return sizeSys[iSys];}
487
488   // New system or new parton in system.
489   int newSystem() {beginSys.push_back(memberSys.size()); 
490     sizeSys.push_back(0); return (beginSys.size() - 1);}
491   void addToSystem(int iSys, int iPos);
492
493   // Get or set value of given member in given system. Replace value by new.
494   int getInSystem(int iSys, int iMem) const {
495     return memberSys[beginSys[iSys] + iMem];}
496   void setInSystem(int iSys, int iMem, int iPos) {
497     memberSys[beginSys[iSys] + iMem] = iPos;}
498   void replaceInSystem(int iSys, int iPosOld, int iPosNew);
499
500   // List members in systems; for debug mainly.
501   void listSystems(ostream& os = cout) const;
502
503   // Operator overloading allows to append one event to an existing one.
504   // Warning: particles should be OK, but some other information unreliable.
505   Event& operator+=(const Event& addEvent);
506
507 private: 
508
509   // Constants: could only be changed in the code itself.
510   static const int IPERLINE;
511
512   // Initialization data, normally only set once.
513   int startColTag;
514
515   // The event: a vector containing all particles (entries).
516   vector<Particle> entry;
517
518   // The list of junctions.
519   vector<Junction> junction;
520
521   // The maximum colour tag of the event so far.
522   int maxColTag;
523
524   // Saved entry and junction list sizes, for simple restoration.
525   int savedSize, savedJunctionSize;
526
527   // The scale of the event; linear quantity in GeV.
528   double scaleSave, scaleSecondSave;
529
530   // Header specification in event listing (at most 40 characters wide).
531   string headerList;
532
533   // Offsets, sizes and values of systems.
534   vector<int> beginSys, sizeSys, memberSys;
535   
536 };
537
538 //**************************************************************************
539
540 } // end namespace Pythia8
541
542 #endif // end Pythia8_Event_H