PYTHIA8/pythia8175/include/Event.h

   1 // Event.h is a part of the PYTHIA event generator.
   2 // Copyright (C) 2013 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
  18 namespace Pythia8 {
  19
  20 //==========================================================================
  21
  22 // Forward references to ParticleDataEntry and ResonanceWidths classes.
  23 class ParticleDataEntry;
  24 class ResonanceWidths;
  25
  26 //==========================================================================
  27
  28 // Particle class.
  29 // This class holds info on a particle in general.
  30
  31 class Particle {
  32
  33 public:
  34
  35   // Constructors.
  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) { }
  61   Particle(const Particle& pt) : idSave(pt.idSave),
  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) { }
  69   Particle& operator=(const Particle& pt) {if (this != &pt) {
  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; }
  77     return *this; }
  78
  79   // Member functions to set the ParticleData and ParticleDataEntry pointers.
  80   void setPDTPtr(ParticleData* pdtPtrIn) { pdtPtr = pdtPtrIn; setPDEPtr();}
  81   void setPDEPtr() {pdePtr = (pdtPtr != 0)
  82     ? pdtPtr->particleDataEntryPtr( idSave) : 0;}
  83
  84   // Member functions for input.
  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;
 102     acolSave = acolIn;}
 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;}
 121
 122   // Member functions for output.
 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;}
 146
 147   // Member functions for output; derived int and bool quantities.
 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);}
 154
 155   // Member functions for output; derived double quantities.
 156   double m2()        const {return (mSave >= 0.) ?  mSave*mSave
 157                                                  : -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();}
 175   double y()         const;
 176   double eta()       const;
 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();}
 187
 188   // Further output, based on a pointer to a ParticleDataEntry object.
 189   string name()      const {
 190     return (pdePtr != 0) ? pdePtr->name(idSave) : " ";}
 191   string nameWithStatus(int maxLen = 20) const;
 192   int    spinType()  const {
 193     return (pdePtr != 0) ? pdePtr->spinType() : 0;}
 194   int    chargeType() const {
 195     return (pdePtr != 0) ? pdePtr->chargeType(idSave) : 0;}
 196   double charge()    const {
 197     return (pdePtr != 0) ?  pdePtr->charge(idSave) : 0;}
 198   bool   isCharged() const {
 199     return (pdePtr != 0) ? (pdePtr->chargeType(idSave) != 0) : false;}
 200   bool   isNeutral() const {
 201     return (pdePtr != 0) ? (pdePtr->chargeType(idSave) == 0) : false;}
 202   int    colType()   const {
 203     return (pdePtr != 0) ? pdePtr->colType(idSave) : 0;}
 204   double m0()        const {
 205     return (pdePtr != 0) ? pdePtr->m0() : 0.;}
 206   double mWidth()    const {
 207     return (pdePtr != 0) ? pdePtr->mWidth() : 0.;}
 208   double mMin()      const {
 209     return (pdePtr != 0) ? pdePtr->mMin() : 0.;}
 210   double mMax()      const {
 211     return (pdePtr != 0) ? pdePtr->mMax() : 0.;}
 212   double mass()      const {
 213     return (pdePtr != 0) ? pdePtr->mass() : 0.;}
 214   double constituentMass() const {
 215     return (pdePtr != 0) ? pdePtr->constituentMass() : 0.;}
 216   double tau0()      const {
 217     return (pdePtr != 0) ? pdePtr->tau0() : 0.;}
 218   bool   mayDecay()  const {
 219     return (pdePtr != 0) ? pdePtr->mayDecay() : false;}
 220   bool   canDecay()  const {
 221     return (pdePtr != 0) ? pdePtr->canDecay() : false;}
 222   bool   doExternalDecay() const {
 223     return (pdePtr != 0) ? pdePtr->doExternalDecay() : false;}
 224   bool   isResonance() const {
 225     return (pdePtr != 0) ? pdePtr->isResonance() : false;}
 226   bool   isVisible() const {
 227     return (pdePtr != 0) ? pdePtr->isVisible() : false;}
 228   bool   isLepton()  const {
 229     return (pdePtr != 0) ? pdePtr->isLepton() : false;}
 230   bool   isQuark()   const {
 231     return (pdePtr != 0) ? pdePtr->isQuark() : false;}
 232   bool   isGluon()   const {
 233     return (pdePtr != 0) ? pdePtr->isGluon() : false;}
 234   bool   isDiquark()   const {
 235     return (pdePtr != 0) ? pdePtr->isDiquark() : false;}
 236   bool   isParton()   const {
 237     return (pdePtr != 0) ? pdePtr->isParton() : false;}
 238   bool   isHadron()  const {
 239     return (pdePtr != 0) ? pdePtr->isHadron() : false;}
 240   ParticleDataEntry& particleDataEntry() const {return *pdePtr;}
 241
 242   // Member functions that perform operations.
 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);}
 262   void rotbst(const RotBstMatrix& M) {pSave.rotbst(M);
 263     if (hasVertexSave) vProdSave.rotbst(M);}
 264   void offsetHistory( int minMother, int addMother, int minDaughter,
 265     int addDaughter);
 266   void offsetCol( int addCol);
 267
 268 private:
 269
 270   // Constants: could only be changed in the code itself.
 271   static const double TINY;
 272
 273   // Properties of the current particle.
 274   int    idSave, statusSave, mother1Save, mother2Save, daughter1Save,
 275          daughter2Save, colSave, acolSave;
 276   Vec4   pSave;
 277   double mSave, scaleSave, polSave;
 278   bool   hasVertexSave;
 279   Vec4   vProdSave;
 280   double tauSave;
 281
 282   // Pointer to properties of the particle species.
 283   // Should no be saved in a persistent copy of the event record.
 284   // The //! below is ROOT notation that this member should not be saved.
 285   // Event::restorePtrs() can be called to restore the missing information.
 286   ParticleDataEntry* pdePtr;  //!
 287
 288   // Pointer to the whole particle data table. Used to update the above
 289   // pointer when id(...) changes identity. As above it should not be saved.
 290   ParticleData*      pdtPtr;  //!
 291
 292 };
 293
 294 // Invariant mass of a pair and its square.
 295 // (Not part of class proper, but tightly linked.)
 296 double m(const Particle&, const Particle&);
 297 double m2(const Particle&, const Particle&);
 298
 299 //==========================================================================
 300
 301 // The junction class stores what kind of junction it is, the colour indices
 302 // of the legs at the junction and as far out as legs have been traced,
 303 // and the status codes assigned for fragmentation of each leg.
 304
 305 class Junction {
 306
 307 public:
 308
 309   // Constructors.
 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; } }
 318   Junction(const Junction& ju) : remainsSave(ju.remainsSave),
 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]; } }
 322   Junction& operator=(const Junction& ju) {if (this != &ju) {
 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]; } }
 326     return *this; }
 327
 328   // Set values.
 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;}
 335
 336   // Read out value.
 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];}
 342
 343 private:
 344
 345   // Kind, positions of the three ends and their status codes.
 346   bool remainsSave;
 347   int kindSave, colSave[3], endColSave[3], statusSave[3];
 348
 349 };
 350
 351 //==========================================================================
 352
 353 // The Event class holds all info on the generated event.
 354
 355 class Event {
 356
 357 public:
 358
 359   // Constructors.
 360   Event(int capacity = 100) : startColTag(100), maxColTag(100),
 361     savedSize(0), savedJunctionSize(0), scaleSave(0.), scaleSecondSave(0.),
 362     headerList("----------------------------------------"),
 363     particleDataPtr(0) { entry.reserve(capacity); }
 364   Event& operator=(const Event& oldEvent);
 365
 366   // Initialize header for event listing, particle data table, and colour.
 367   void init( string headerIn = "", ParticleData* particleDataPtrIn = 0,
 368     int startColTagIn = 100) {
 369     headerList.replace(0, headerIn.length() + 2, headerIn + "  ");
 370      particleDataPtr = particleDataPtrIn; startColTag = startColTagIn;}
 371
 372   // Clear event record.
 373   void clear() {entry.resize(0); maxColTag = startColTag; scaleSave = 0.;
 374     scaleSecondSave = 0.; clearJunctions();}
 375
 376   // Clear event record, and set first particle empty.
 377   void reset() {clear(); append(90, -11, 0, 0, 0., 0., 0., 0., 0.);}
 378
 379   // Overload index operator to access element of event record.
 380   Particle& operator[](int i) {return entry[i];}
 381   const Particle& operator[](int i) const {return entry[i];}
 382
 383   // Other way to access an element of event record.
 384   Particle& at(int i) {return entry.at(i);}
 385
 386   // Event record size.
 387   int size() const {return entry.size();}
 388
 389   // Put a new particle at the end of the event record; return index.
 390   int append(Particle entryIn) {
 391     entry.push_back(entryIn); setPDTPtr();
 392     if (entryIn.col() > maxColTag) maxColTag = entryIn.col();
 393     if (entryIn.acol() > maxColTag) maxColTag = entryIn.acol();
 394     return entry.size() - 1;
 395   }
 396   int append(int id, int status, int mother1, int mother2, int daughter1,
 397     int daughter2, int col, int acol, double px, double py, double pz,
 398     double e, double m = 0., double scaleIn = 0., double polIn = 9.) {
 399     entry.push_back( Particle(id, status, mother1, mother2, daughter1,
 400     daughter2, col, acol, px, py, pz, e, m, scaleIn, polIn) ); setPDTPtr();
 401     if (col > maxColTag) maxColTag = col;
 402     if (acol > maxColTag) maxColTag = acol;
 403     return entry.size() - 1;
 404   }
 405   int append(int id, int status, int mother1, int mother2, int daughter1,
 406     int daughter2, int col, int acol, Vec4 p, double m = 0.,
 407     double scaleIn = 0., double polIn = 9.) {
 408     entry.push_back( Particle(id, status, mother1, mother2, daughter1,
 409     daughter2, col, acol, p, m, scaleIn, polIn) ); setPDTPtr();
 410     if (col > maxColTag) maxColTag = col;
 411     if (acol > maxColTag) maxColTag = acol;
 412     return entry.size() - 1;
 413   }
 414
 415   // Brief versions of append: no mothers and no daughters.
 416   int append(int id, int status, int col, int acol, double px, double py,
 417     double pz, double e, double m = 0., double scaleIn = 0.,
 418     double polIn = 9.) { entry.push_back( Particle(id, status, 0, 0, 0, 0,
 419     col, acol, px, py, pz, e, m, scaleIn, polIn) ); setPDTPtr();
 420     if (col > maxColTag) maxColTag = col;
 421     if (acol > maxColTag) maxColTag = acol;
 422     return entry.size() - 1;
 423   }
 424   int append(int id, int status, int col, int acol, Vec4 p, double m = 0.,
 425     double scaleIn = 0., double polIn = 9.) {entry.push_back( Particle(id,
 426     status, 0, 0, 0, 0, col, acol, p, m, scaleIn, polIn) ); setPDTPtr();
 427     if (col > maxColTag) maxColTag = col;
 428     if (acol > maxColTag) maxColTag = acol;
 429     return entry.size() - 1;
 430   }
 431
 432   // Set pointer to ParticleData for a particle, by default latest one.
 433   void setPDTPtr(int iSet = -1) {if (iSet < 0) iSet = entry.size() - 1;
 434     entry[iSet].setPDTPtr( particleDataPtr);}
 435
 436   // Add a copy of an existing particle at the end of the event record.
 437   int copy(int iCopy, int newStatus = 0);
 438
 439   // Implement reference "back" to access last element.
 440   Particle& back() {return entry.back();}
 441
 442   // List the particles in an event.
 443   void list() const;
 444   void list(ostream& os) const;
 445   void list(bool showScaleAndVertex, bool showMothersAndDaughters = false)
 446     const;
 447   void list(bool showScaleAndVertex, bool showMothersAndDaughters,
 448     ostream& os) const;
 449
 450   // Remove last n entries.
 451   void popBack(int nRemove = 1) { if (nRemove ==1) entry.pop_back();
 452     else {int newSize = max( 0, size() - nRemove);
 453     entry.resize(newSize);} }
 454
 455   // Undo the decay of a single particle (where daughters well-defined).
 456   bool undoDecay(int i);
 457
 458   // Restore all ParticleDataEntry* pointers in the Particle vector.
 459   // Useful when a persistent copy of the event record is read back in.
 460   void restorePtrs() { for (int i = 0; i < size(); ++i) setPDTPtr(i); }
 461
 462   // Save or restore the size of the event record (throwing at the end).
 463   void saveSize() {savedSize = entry.size();}
 464   void restoreSize() {entry.resize(savedSize);}
 465   int  savedSizeValue() {return savedSize;}
 466
 467   // Initialize and access colour tag information.
 468   void initColTag(int colTag = 0) {maxColTag = max( colTag,startColTag);}
 469   int lastColTag() const {return maxColTag;}
 470   int nextColTag() {return ++maxColTag;}
 471
 472   // Access scale for which event as a whole is defined.
 473   void scale( double scaleIn) {scaleSave = scaleIn;}
 474   double scale() const {return scaleSave;}
 475
 476   // Need a second scale if two hard interactions in event.
 477   void scaleSecond( double scaleSecondIn) {scaleSecondSave = scaleSecondIn;}
 478   double scaleSecond() const {return scaleSecondSave;}
 479
 480   // Find complete list of daughters and mothers.
 481   vector<int> motherList(int i) const;
 482   vector<int> daughterList(int i) const;
 483
 484   // Convert to HepMC status code conventions.
 485   int statusHepMC(int i) const;
 486
 487   // Trace the first and last copy of one and the same particle.
 488   int iTopCopy(int i) const;
 489   int iBotCopy(int i) const;
 490
 491   // Trace the first and last copy of a particle, using flavour match.
 492   int iTopCopyId(int i) const;
 493   int iBotCopyId(int i) const;
 494
 495   // Find list of sisters, also tracking up and down identical copies.
 496   vector<int> sisterList(int i) const;
 497   vector<int> sisterListTopBot(int i, bool widenSearch = true) const;
 498
 499   // Check whether two particles have a direct mother-daughter relation.
 500   bool isAncestor(int i, int iAncestor) const;
 501
 502   // Member functions for rotations and boosts of an event.
 503   void rot(double theta, double phi)
 504     {for (int i = 0; i < size(); ++i) entry[i].rot(theta, phi);}
 505   void bst(double betaX, double betaY, double betaZ)
 506     {for (int i = 0; i < size(); ++i) entry[i].bst(betaX, betaY, betaZ);}
 507   void bst(double betaX, double betaY, double betaZ, double gamma)
 508     {for (int i = 0; i < size(); ++i) entry[i].bst(betaX, betaY, betaZ,
 509     gamma);}
 510   void bst(const Vec4& vec)
 511     {for (int i = 0; i < size(); ++i) entry[i].bst(vec);}
 512   void rotbst(const RotBstMatrix& M)
 513     {for (int i = 0; i < size(); ++i) entry[i].rotbst(M);}
 514
 515   // Clear the list of junctions.
 516   void clearJunctions() {junction.resize(0);}
 517
 518   // Add a junction to the list, study it or extra input.
 519   void appendJunction( int kind, int col0, int col1, int col2)
 520     { junction.push_back( Junction( kind, col0, col1, col2) );}
 521   void appendJunction(Junction junctionIn) {junction.push_back(junctionIn);}
 522   int sizeJunction() const {return junction.size();}
 523   bool remainsJunction(int i) const {return junction[i].remains();}
 524   void remainsJunction(int i, bool remainsIn) {junction[i].remains(remainsIn);}
 525   int kindJunction(int i) const {return junction[i].kind();}
 526   int colJunction( int i, int j) const {return junction[i].col(j);}
 527   void colJunction( int i, int j, int colIn) {junction[i].col(j, colIn);}
 528   int endColJunction( int i, int j) const {return junction[i].endCol(j);}
 529   void endColJunction( int i, int j, int endColIn)
 530     {junction[i].endCol(j, endColIn);}
 531   int statusJunction( int i, int j) const {return junction[i].status(j);}
 532   void statusJunction( int i, int j, int statusIn)
 533     {junction[i].status(j, statusIn);}
 534   Junction& getJunction(int i) {return junction[i];}
 535   const Junction& getJunction(int i) const {return junction[i];}
 536   void eraseJunction(int i);
 537
 538   // Save or restore the size of the junction list (throwing at the end).
 539   void saveJunctionSize() {savedJunctionSize = junction.size();}
 540   void restoreJunctionSize() {junction.resize(savedJunctionSize);}
 541
 542   // List any junctions in the event; for debug mainly.
 543   void listJunctions(ostream& os = cout) const;
 544
 545   // Operator overloading allows to append one event to an existing one.
 546   // Warning: particles should be OK, but some other information unreliable.
 547   Event& operator+=(const Event& addEvent);
 548
 549 private:
 550
 551   // Constants: could only be changed in the code itself.
 552   static const int IPERLINE;
 553
 554   // Initialization data, normally only set once.
 555   int startColTag;
 556
 557   // The event: a vector containing all particles (entries).
 558   // The explicit use of Pythia8:: qualifier patches a limitation in ROOT.
 559   vector<Pythia8::Particle> entry;
 560
 561   // The list of junctions.
 562   // The explicit use of Pythia8:: qualifier patches a limitation in ROOT.
 563   vector<Pythia8::Junction> junction;
 564
 565   // The maximum colour tag of the event so far.
 566   int maxColTag;
 567
 568   // Saved entry and junction list sizes, for simple restoration.
 569   int savedSize, savedJunctionSize;
 570
 571   // The scale of the event; linear quantity in GeV.
 572   double scaleSave, scaleSecondSave;
 573
 574   // Header specification in event listing (at most 40 characters wide).
 575   string headerList;
 576
 577   // Pointer to the particle data table.
 578   // The //! below is ROOT notation that this member should not be saved.
 579   ParticleData* particleDataPtr;  //!
 580
 581 };
 582
 583 //==========================================================================
 584
 585 } // end namespace Pythia8
 586
 587 #endif // end Pythia8_Event_H