include/Pythia8/Event.h

0001 // Event.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2024 Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005
0006 // Header file for the Particle and Event classes.
0007 // Particle: information on an instance of a particle.
0008 // Junction: information on a junction between three colours.
0009 // Event: list of particles in the current event.
0010
0011 #ifndef Pythia8_Event_H
0012 #define Pythia8_Event_H
0013
0014 #include "Pythia8/Basics.h"
0015 #include "Pythia8/ParticleData.h"
0016 #include "Pythia8/PythiaStdlib.h"
0017
0018 namespace Pythia8 {
0019
0020 //==========================================================================
0021
0022 // Forward references to ParticleDataEntry and ResonanceWidths classes.
0023 class ParticleDataEntry;
0024 class ResonanceWidths;
0025 class Event;
0026
0027 //==========================================================================
0028
0029 // Particle class.
0030 // This class holds info on a particle in general.
0031
0032 class Particle {
0033
0034 public:
0035
0036   // Constructors.
0037   Particle() : idSave(0), statusSave(0), mother1Save(0), mother2Save(0),
0038     daughter1Save(0), daughter2Save(0), colSave(0), acolSave(0),
0039     pSave(Vec4(0.,0.,0.,0.)), mSave(0.), scaleSave(0.), polSave(9.),
0040     hasVertexSave(false), vProdSave(Vec4(0.,0.,0.,0.)), tauSave(0.),
0041     pdePtr(0), evtPtr(0) { }
0042   Particle(int idIn, int statusIn = 0, int mother1In = 0,
0043     int mother2In = 0, int daughter1In = 0, int daughter2In = 0,
0044     int colIn = 0, int acolIn = 0, double pxIn = 0., double pyIn = 0.,
0045     double pzIn = 0., double eIn = 0., double mIn = 0.,
0046     double scaleIn = 0., double polIn = 9.)
0047     : idSave(idIn), statusSave(statusIn), mother1Save(mother1In),
0048     mother2Save(mother2In), daughter1Save(daughter1In),
0049     daughter2Save(daughter2In), colSave(colIn), acolSave(acolIn),
0050     pSave(Vec4(pxIn, pyIn, pzIn, eIn)), mSave(mIn), scaleSave(scaleIn),
0051     polSave(polIn), hasVertexSave(false), vProdSave(Vec4(0.,0.,0.,0.)),
0052     tauSave(0.), pdePtr(0), evtPtr(0) { }
0053   Particle(int idIn, int statusIn, int mother1In, int mother2In,
0054     int daughter1In, int daughter2In, int colIn, int acolIn,
0055     Vec4 pIn, double mIn = 0., double scaleIn = 0., double polIn = 9.)
0056     : idSave(idIn), statusSave(statusIn), mother1Save(mother1In),
0057     mother2Save(mother2In), daughter1Save(daughter1In),
0058     daughter2Save(daughter2In), colSave(colIn), acolSave(acolIn),
0059     pSave(pIn), mSave(mIn), scaleSave(scaleIn), polSave(polIn),
0060     hasVertexSave(false), vProdSave(Vec4(0.,0.,0.,0.)), tauSave(0.),
0061     pdePtr(0), evtPtr(0) { }
0062   Particle(const Particle& pt) : idSave(pt.idSave),
0063     statusSave(pt.statusSave), mother1Save(pt.mother1Save),
0064     mother2Save(pt.mother2Save), daughter1Save(pt.daughter1Save),
0065     daughter2Save(pt.daughter2Save), colSave(pt.colSave),
0066     acolSave(pt.acolSave), pSave(pt.pSave), mSave(pt.mSave),
0067     scaleSave(pt.scaleSave), polSave(pt.polSave),
0068     hasVertexSave(pt.hasVertexSave), vProdSave(pt.vProdSave),
0069     tauSave(pt.tauSave), pdePtr(pt.pdePtr), evtPtr(pt.evtPtr) { }
0070   Particle& operator=(const Particle& pt) {if (this != &pt) {
0071     idSave = pt.idSave; statusSave = pt.statusSave;
0072     mother1Save = pt.mother1Save; mother2Save = pt.mother2Save;
0073     daughter1Save = pt.daughter1Save; daughter2Save = pt.daughter2Save;
0074     colSave = pt.colSave; acolSave = pt.acolSave; pSave = pt.pSave;
0075     mSave = pt.mSave; scaleSave = pt.scaleSave; polSave = pt.polSave;
0076     hasVertexSave = pt.hasVertexSave; vProdSave = pt.vProdSave;
0077     tauSave = pt.tauSave; pdePtr = pt.pdePtr; evtPtr = pt.evtPtr; }
0078     return *this; }
0079
0080   // Destructor.
0081   virtual ~Particle() {}
0082
0083   // Member functions to set the Event and ParticleDataEntry pointers.
0084   void setEvtPtr(Event* evtPtrIn) { evtPtr = evtPtrIn; setPDEPtr();}
0085   void setPDEPtr(ParticleDataEntryPtr pdePtrIn = nullptr);
0086
0087   // Member functions for input.
0088   void id(int idIn) {idSave = idIn; setPDEPtr();}
0089   void status(int statusIn) {statusSave = statusIn;}
0090   void statusPos() {statusSave = abs(statusSave);}
0091   void statusNeg() {statusSave = -abs(statusSave);}
0092   void statusCode(int statusIn) {statusSave =
0093     (statusSave > 0) ? abs(statusIn) : -abs(statusIn);}
0094   void mother1(int mother1In) {mother1Save = mother1In;}
0095   void mother2(int mother2In) {mother2Save = mother2In;}
0096   void mothers(int mother1In = 0, int mother2In = 0)
0097     {mother1Save = mother1In; mother2Save = mother2In;}
0098   void daughter1(int daughter1In) {daughter1Save = daughter1In;}
0099   void daughter2(int daughter2In) {daughter2Save = daughter2In;}
0100   void daughters(int daughter1In = 0, int daughter2In = 0)
0101     {daughter1Save = daughter1In; daughter2Save = daughter2In;}
0102   void col(int colIn) {colSave = colIn;}
0103   void acol(int acolIn) {acolSave = acolIn;}
0104   void cols(int colIn = 0,int acolIn = 0) {colSave = colIn;
0105     acolSave = acolIn;}
0106   void p(Vec4 pIn) {pSave = pIn;}
0107   void p(double pxIn, double pyIn, double pzIn, double eIn)
0108     {pSave.p(pxIn, pyIn, pzIn, eIn);}
0109   void px(double pxIn) {pSave.px(pxIn);}
0110   void py(double pyIn) {pSave.py(pyIn);}
0111   void pz(double pzIn) {pSave.pz(pzIn);}
0112   void e(double eIn) {pSave.e(eIn);}
0113   void m(double mIn) {mSave = mIn;}
0114   void scale(double scaleIn) {scaleSave = scaleIn;}
0115   void pol(double polIn) {polSave = polIn;}
0116   void vProd(Vec4 vProdIn) {vProdSave = vProdIn; hasVertexSave = true;}
0117   void vProd(double xProdIn, double yProdIn, double zProdIn, double tProdIn)
0118     {vProdSave.p(xProdIn, yProdIn, zProdIn, tProdIn); hasVertexSave = true;}
0119   void xProd(double xProdIn) {vProdSave.px(xProdIn); hasVertexSave = true;}
0120   void yProd(double yProdIn) {vProdSave.py(yProdIn); hasVertexSave = true;}
0121   void zProd(double zProdIn) {vProdSave.pz(zProdIn); hasVertexSave = true;}
0122   void tProd(double tProdIn) {vProdSave.e(tProdIn); hasVertexSave = true;}
0123   void vProdAdd(Vec4 vProdIn) {vProdSave += vProdIn; hasVertexSave = true;}
0124   void tau(double tauIn) {tauSave = tauIn;}
0125
0126   // Member functions for output.
0127   int    id()        const {return idSave;}
0128   int    status()    const {return statusSave;}
0129   int    mother1()   const {return mother1Save;}
0130   int    mother2()   const {return mother2Save;}
0131   int    daughter1() const {return daughter1Save;}
0132   int    daughter2() const {return daughter2Save;}
0133   int    col()       const {return colSave;}
0134   int    acol()      const {return acolSave;}
0135   Vec4   p()         const {return pSave;}
0136   double px()        const {return pSave.px();}
0137   double py()        const {return pSave.py();}
0138   double pz()        const {return pSave.pz();}
0139   double e()         const {return pSave.e();}
0140   double m()         const {return mSave;}
0141   double scale()     const {return scaleSave;}
0142   double pol()       const {return polSave;}
0143   bool   hasVertex() const {return hasVertexSave;}
0144   Vec4   vProd()     const {return vProdSave;}
0145   double xProd()     const {return vProdSave.px();}
0146   double yProd()     const {return vProdSave.py();}
0147   double zProd()     const {return vProdSave.pz();}
0148   double tProd()     const {return vProdSave.e();}
0149   double tau()       const {return tauSave;}
0150
0151   // Member functions for output; derived int and bool quantities.
0152   int    idAbs()     const {return abs(idSave);}
0153   int    statusAbs() const {return abs(statusSave);}
0154   bool   isFinal()   const {return (statusSave > 0);}
0155   int    intPol()    const;
0156   bool   isRescatteredIncoming() const {return
0157     (statusSave == -34 || statusSave == -45 ||
0158      statusSave == -46 || statusSave == -54);}
0159
0160   // Member functions for output; derived double quantities.
0161   double m2()        const {return (mSave >= 0.) ?  mSave*mSave
0162                                                  : -mSave*mSave;}
0163   double mCalc()     const {return pSave.mCalc();}
0164   double m2Calc()    const {return pSave.m2Calc();}
0165   double eCalc()     const {return sqrt(abs(m2() + pSave.pAbs2()));}
0166   double pT()        const {return pSave.pT();}
0167   double pT2()       const {return pSave.pT2();}
0168   double mT()        const {double temp = m2() + pSave.pT2();
0169     return (temp >= 0.) ? sqrt(temp) : -sqrt(-temp);}
0170   double mT2()       const {return m2() + pSave.pT2();}
0171   double pAbs()      const {return pSave.pAbs();}
0172   double pAbs2()     const {return pSave.pAbs2();}
0173   double eT()        const {return pSave.eT();}
0174   double eT2()       const {return pSave.eT2();}
0175   double theta()     const {return pSave.theta();}
0176   double phi()       const {return pSave.phi();}
0177   double thetaXZ()   const {return pSave.thetaXZ();}
0178   double pPos()      const {return pSave.pPos();}
0179   double pNeg()      const {return pSave.pNeg();}
0180   double y()         const;
0181   double eta()       const;
0182   double y(double mCut) const;
0183   double y(double mCut, RotBstMatrix& M) const;
0184   Vec4   vDec()      const {return (tauSave > 0. && mSave > 0.)
0185     ? vProdSave + tauSave * pSave / mSave : vProdSave;}
0186   double xDec()      const {return (tauSave > 0. && mSave > 0.)
0187     ? vProdSave.px() + tauSave * pSave.px() / mSave : vProdSave.px();}
0188   double yDec()      const {return (tauSave > 0. && mSave > 0.)
0189     ? vProdSave.py() + tauSave * pSave.py() / mSave : vProdSave.py();}
0190   double zDec()      const {return (tauSave > 0. && mSave > 0.)
0191     ? vProdSave.pz() + tauSave * pSave.pz() / mSave : vProdSave.pz();}
0192   double tDec()      const {return (tauSave > 0. && mSave > 0.)
0193     ? vProdSave.e()  + tauSave * pSave.e()  / mSave : vProdSave.e();}
0194
0195   // Methods that can refer back to the event the particle belongs to.
0196   virtual int index() const;
0197   int iTopCopy()      const;
0198   int iBotCopy()      const;
0199   int iTopCopyId(bool simplify = false) const;
0200   int iBotCopyId(bool simplify = false) const;
0201   vector<int> motherList()   const;
0202   vector<int> daughterList() const;
0203   vector<int> daughterListRecursive() const;
0204   vector<int> sisterList(bool traceTopBot = false) const;
0205   bool isAncestor(int iAncestor) const;
0206   int statusHepMC()  const;
0207   bool isFinalPartonLevel() const;
0208   bool undoDecay();
0209
0210   // Get and set Hidden Valley colours, referring back to the event.
0211   int  colHV() const;
0212   int  acolHV() const;
0213   void colHV(int colHVin);
0214   void acolHV(int acolHVin);
0215   void colsHV(int colHVin, int acolHVin);
0216
0217   // Further output, based on a pointer to a ParticleDataEntry object.
0218   string name()      const {
0219     return (pdePtr != 0) ? pdePtr->name(idSave) : " ";}
0220   string nameWithStatus(int maxLen = 20) const;
0221   int    spinType()  const {
0222     return (pdePtr != 0) ? pdePtr->spinType() : 0;}
0223   int    chargeType() const {
0224     return (pdePtr != 0) ? pdePtr->chargeType(idSave) : 0;}
0225   double charge()    const {
0226     return (pdePtr != 0) ?  pdePtr->charge(idSave) : 0;}
0227   bool   isCharged() const {
0228     return (pdePtr != 0) ? (pdePtr->chargeType(idSave) != 0) : false;}
0229   bool   isNeutral() const {
0230     return (pdePtr != 0) ? (pdePtr->chargeType(idSave) == 0) : false;}
0231   int    colType()   const {
0232     return (pdePtr != 0) ? pdePtr->colType(idSave) : 0;}
0233   double m0()        const {
0234     return (pdePtr != 0) ? pdePtr->m0() : 0.;}
0235   double mWidth()    const {
0236     return (pdePtr != 0) ? pdePtr->mWidth() : 0.;}
0237   double mMin()      const {
0238     return (pdePtr != 0) ? pdePtr->mMin() : 0.;}
0239   double mMax()      const {
0240     return (pdePtr != 0) ? pdePtr->mMax() : 0.;}
0241   double mSel()      const {
0242     return (pdePtr != 0) ? pdePtr->mSel() : 0.;}
0243   double constituentMass() const {
0244     return (pdePtr != 0) ? pdePtr->constituentMass() : 0.;}
0245   double tau0()      const {
0246     return (pdePtr != 0) ? pdePtr->tau0() : 0.;}
0247   bool   mayDecay()  const {
0248     return (pdePtr != 0) ? pdePtr->mayDecay() : false;}
0249   bool   canDecay()  const {
0250     return (pdePtr != 0) ? pdePtr->canDecay() : false;}
0251   bool   doExternalDecay() const {
0252     return (pdePtr != 0) ? pdePtr->doExternalDecay() : false;}
0253   bool   isResonance() const {
0254     return (pdePtr != 0) ? pdePtr->isResonance() : false;}
0255   bool   isVisible() const {
0256     return (pdePtr != 0) ? pdePtr->isVisible() : false;}
0257   bool   isLepton()  const {
0258     return (pdePtr != 0) ? pdePtr->isLepton() : false;}
0259   bool   isQuark()   const {
0260     return (pdePtr != 0) ? pdePtr->isQuark() : false;}
0261   bool   isGluon()   const {
0262     return (pdePtr != 0) ? pdePtr->isGluon() : false;}
0263   bool   isDiquark()   const {
0264     return (pdePtr != 0) ? pdePtr->isDiquark() : false;}
0265   bool   isParton()   const {
0266     return (pdePtr != 0) ? pdePtr->isParton() : false;}
0267   bool   isHadron()  const {
0268     return (pdePtr != 0) ? pdePtr->isHadron() : false;}
0269   bool   isExotic()  const {
0270     return (pdePtr != 0) ? pdePtr->isExotic() : false;}
0271   ParticleDataEntry& particleDataEntry() const {return *pdePtr;}
0272
0273   // Member functions that perform operations.
0274   void rescale3(double fac) {pSave.rescale3(fac);}
0275   void rescale4(double fac) {pSave.rescale4(fac);}
0276   void rescale5(double fac) {pSave.rescale4(fac); mSave *= fac;}
0277   void rot(double thetaIn, double phiIn) {pSave.rot(thetaIn, phiIn);
0278     if (hasVertexSave) vProdSave.rot(thetaIn, phiIn);}
0279   void bst(double betaX, double betaY, double betaZ) {
0280     pSave.bst(betaX, betaY, betaZ);
0281     if (hasVertexSave) vProdSave.bst(betaX, betaY, betaZ);}
0282   void bst(double betaX, double betaY, double betaZ, double gamma) {
0283     pSave.bst(betaX, betaY, betaZ, gamma);
0284     if (hasVertexSave) vProdSave.bst(betaX, betaY, betaZ, gamma);}
0285   void bst(const Vec4& pBst) {pSave.bst(pBst);
0286     if (hasVertexSave) vProdSave.bst(pBst);}
0287   void bst(const Vec4& pBst, double mBst) {pSave.bst(pBst, mBst);
0288     if (hasVertexSave) vProdSave.bst(pBst, mBst);}
0289   void bstback(const Vec4& pBst) {pSave.bstback(pBst);
0290     if (hasVertexSave) vProdSave.bstback(pBst);}
0291   void bstback(const Vec4& pBst, double mBst) {pSave.bstback(pBst, mBst);
0292     if (hasVertexSave) vProdSave.bstback(pBst, mBst);}
0293   void rotbst(const RotBstMatrix& M, bool boostVertex = true) {pSave.rotbst(M);
0294     if (hasVertexSave && boostVertex) vProdSave.rotbst(M);}
0295   void offsetHistory( int minMother, int addMother, int minDaughter,
0296     int addDaughter);
0297   void offsetCol( int addCol);
0298
0299 protected:
0300
0301   // Constants: could only be changed in the code itself.
0302   static const double TINY;
0303
0304   // Properties of the current particle.
0305   int    idSave, statusSave, mother1Save, mother2Save, daughter1Save,
0306          daughter2Save, colSave, acolSave;
0307   Vec4   pSave;
0308   double mSave, scaleSave, polSave;
0309   bool   hasVertexSave;
0310   Vec4   vProdSave;
0311   double tauSave;
0312
0313   // Pointer to properties of the particle species.
0314   // Should no be saved in a persistent copy of the event record.
0315   // The //! below is ROOT notation that this member should not be saved.
0316   // Event::restorePtrs() can be called to restore the missing information.
0317   ParticleDataEntryPtr pdePtr;  //!
0318
0319   // Pointer to the whole event record to which the particle belongs (if any).
0320   // As above it should not be saved.
0321   Event*             evtPtr;  //!
0322
0323 };
0324
0325 // Particles invariant mass, mass squared, and momentum dot product.
0326 // (Not part of class proper, but tightly linked.)
0327 double m(const Particle& pp1, const Particle& pp2);
0328 double m2(const Particle& pp1, const Particle& pp2);
0329 double m2(const Particle& pp1, const Particle& pp2, const Particle& pp3);
0330 double dot4(const Particle& pp1, const Particle& pp2);
0331
0332 //==========================================================================
0333
0334 // The junction class stores what kind of junction it is, the colour indices
0335 // of the legs at the junction and as far out as legs have been traced,
0336 // and the status codes assigned for fragmentation of each leg.
0337
0338 class Junction {
0339
0340 public:
0341
0342   // Constructors.
0343   Junction() : remainsSave(true), kindSave(0), colSave(), endColSave(),
0344     statusSave() { }
0345
0346   Junction( int kindIn, int col0In, int col1In, int col2In)
0347     : remainsSave(true), kindSave(kindIn), colSave(), endColSave(),
0348     statusSave() {
0349       colSave[0] = col0In; colSave[1] = col1In; colSave[2] = col2In;
0350       for (int j = 0; j < 3; ++j) {
0351       endColSave[j] = colSave[j];  } }
0352   Junction(const Junction& ju) : remainsSave(ju.remainsSave),
0353     kindSave(ju.kindSave), colSave(), endColSave(), statusSave() {
0354     for (int j = 0; j < 3; ++j) {
0355     colSave[j] = ju.colSave[j]; endColSave[j] = ju.endColSave[j];
0356     statusSave[j] = ju.statusSave[j]; } }
0357   Junction& operator=(const Junction& ju) {if (this != &ju) {
0358     remainsSave = ju.remainsSave; kindSave =  ju.kindSave;
0359     for (int j = 0; j < 3; ++j) { colSave[j] = ju.colSave[j];
0360     endColSave[j] = ju.endColSave[j]; statusSave[j] = ju.statusSave[j]; } }
0361     return *this; }
0362
0363   // Set values.
0364   void remains(bool remainsIn) {remainsSave = remainsIn;}
0365   void col(int j, int colIn) {colSave[j] = colIn; endColSave[j] = colIn;}
0366   void cols(int j, int colIn, int endColIn) {colSave[j] = colIn;
0367     endColSave[j] = endColIn;}
0368   void endCol(int j, int endColIn) {endColSave[j] = endColIn;}
0369   void status(int j, int statusIn) {statusSave[j] = statusIn;}
0370
0371   // Read out value.
0372   bool   remains()     const {return remainsSave;}
0373   int    kind()        const {return kindSave;}
0374   int    col(int j)    const {return colSave[j];}
0375   int    endCol(int j) const {return endColSave[j];}
0376   int    status(int j) const {return statusSave[j];}
0377
0378 private:
0379
0380   // Kind, positions of the three ends and their status codes.
0381   bool remainsSave;
0382   int kindSave, colSave[3], endColSave[3], statusSave[3];
0383
0384 };
0385
0386 //==========================================================================
0387
0388 // The HVcols class stores Hidden Valley colours for HV-coloured particles.
0389
0390 class HVcols {
0391
0392 public:
0393
0394   // Default constructor and constructor with standard input.
0395   HVcols() : iHV(0), colHV(0), acolHV(0) { }
0396   HVcols( int iHVin, int colHVin, int acolHVin) : iHV(iHVin),
0397     colHV(colHVin), acolHV(acolHVin) {}
0398
0399   // Index of HV-particle and its HV-colour and HV-anticolour.
0400   int iHV, colHV, acolHV;
0401
0402 };
0403
0404 //==========================================================================
0405
0406 // The StringBreaks class stores information about how a single string has
0407 // been broken up into quark/antiquark and diquark/anti-diquark pairs.
0408
0409 class StringBreaks {
0410
0411 public:
0412
0413   // Set and get string end indices.
0414   void setEnds(int iPosIn, int iNegIn) {iPos = iPosIn; iNeg = iNegIn;}
0415   int posEnd() {return iPos;}
0416   int negEnd() {return iNeg;}
0417
0418   // Clear the breaks.
0419   void clear() {breaks.clear();}
0420
0421   // Add a string break by particle ID.
0422   void add(int id) {breaks[abs(id)] += 1;}
0423
0424   // Returns the number of breaks for a given particle ID, or the the sum of
0425   // all, light quark, or diquark breaks.
0426   const map<int, int>& nId() {return breaks;}
0427   int nId(int id) const {
0428     auto itr = breaks.find(abs(id));
0429     return itr == breaks.end() ? 0 : itr->second;}
0430   int nAll() const {int nSum = 0; for (auto &id : breaks) nSum += id.second;
0431     return nSum;}
0432   int nQuark() const {return nId(1) + nId(2) + nId(3) + nId(4);}
0433   int nDiquark() const {int nSum = 0;
0434     for (const auto &id : breaks)
0435       if (id.first > 1000 && id.first < 10000 && (id.first/10)%10 == 0)
0436         nSum += id.second;
0437     return nSum;}
0438
0439 private:
0440
0441   // Index of string end partons and sums of all, quark, and diquark breaks.
0442   int iPos{0}, iNeg{0};
0443
0444   // Number of breaks mapped by particle ID.
0445   map<int, int> breaks{};
0446
0447 };
0448
0449 //==========================================================================
0450
0451 // The Event class holds all info on the generated event.
0452
0453 class Event {
0454
0455 public:
0456
0457   // Constructors.
0458   Event(int capacity = 100) : startColTag(100), maxColTag(100),
0459     savedSize(0), savedJunctionSize(0), savedHVcolsSize(0),
0460     savedPartonLevelSize(0), scaleSave(0.), scaleSecondSave(0.),
0461     headerList("----------------------------------------"),
0462     particleDataPtr(0) { entry.reserve(capacity); }
0463   Event& operator=(const Event& oldEvent);
0464   Event(const Event& oldEvent) {*this = oldEvent;}
0465
0466   // Initialize header for event listing, particle data table, and colour.
0467   void init( string headerIn = "", ParticleData* particleDataPtrIn = 0,
0468     int startColTagIn = 100) {
0469     headerList.replace(0, headerIn.length() + 2, headerIn + "  ");
0470      particleDataPtr = particleDataPtrIn; startColTag = startColTagIn;}
0471
0472   // Clear event record.
0473   void clear() {entry.resize(0); maxColTag = startColTag;
0474     savedPartonLevelSize = 0; scaleSave = 0.; scaleSecondSave = 0.;
0475     clearJunctions(); clearHV(); clearStringBreaks();}
0476   void free() {vector<Particle>().swap(entry); maxColTag = startColTag;
0477     savedPartonLevelSize = 0; scaleSave = 0.; scaleSecondSave = 0.;
0478     clearJunctions(); clearHV(); clearStringBreaks();}
0479
0480   // Clear event record, and set first particle empty.
0481   void reset() {clear(); append(90, -11, 0, 0, 0., 0., 0., 0., 0.);}
0482
0483   // Overload index operator to access element of event record.
0484   Particle& operator[](int i) {return entry.at(i);}
0485   const Particle& operator[](int i) const {return entry.at(i);}
0486
0487   // Implement standard references to elements in the particle array.
0488   Particle& front()   {return entry.front();}
0489   Particle& at(int i) {return entry.at(i);}
0490   Particle& back()    {return entry.back();}
0491   const Particle& front()   const {return entry.front();}
0492   const Particle& at(int i) const {return entry.at(i);}
0493   const Particle& back()    const {return entry.back();}
0494
0495   // Implement iterators for the particle array.
0496   vector<Pythia8::Particle>::iterator begin() { return entry.begin(); }
0497   vector<Pythia8::Particle>::iterator end() { return entry.end(); }
0498   vector<Pythia8::Particle>::const_iterator begin() const
0499     { return entry.begin(); }
0500   vector<Pythia8::Particle>::const_iterator end() const
0501     { return entry.end(); }
0502
0503   // Event record size.
0504   int size() const {return entry.size();}
0505
0506   // Put a new particle at the end of the event record; return index.
0507   int append(Particle entryIn) {
0508     entry.push_back(entryIn); setEvtPtr();
0509     if (entryIn.col() > maxColTag) maxColTag = entryIn.col();
0510     if (entryIn.acol() > maxColTag) maxColTag = entryIn.acol();
0511     return entry.size() - 1;
0512   }
0513   int append(int id, int status, int mother1, int mother2, int daughter1,
0514     int daughter2, int col, int acol, double px, double py, double pz,
0515     double e, double m = 0., double scaleIn = 0., double polIn = 9.) {
0516     entry.push_back( Particle(id, status, mother1, mother2, daughter1,
0517     daughter2, col, acol, px, py, pz, e, m, scaleIn, polIn) ); setEvtPtr();
0518     if (col > maxColTag) maxColTag = col;
0519     if (acol > maxColTag) maxColTag = acol;
0520     return entry.size() - 1;
0521   }
0522   int append(int id, int status, int mother1, int mother2, int daughter1,
0523     int daughter2, int col, int acol, Vec4 p, double m = 0.,
0524     double scaleIn = 0., double polIn = 9.) {
0525     entry.push_back( Particle(id, status, mother1, mother2, daughter1,
0526     daughter2, col, acol, p, m, scaleIn, polIn) ); setEvtPtr();
0527     if (col > maxColTag) maxColTag = col;
0528     if (acol > maxColTag) maxColTag = acol;
0529     return entry.size() - 1;
0530   }
0531
0532   // Brief versions of append: no mothers and no daughters.
0533   int append(int id, int status, int col, int acol, double px, double py,
0534     double pz, double e, double m = 0., double scaleIn = 0.,
0535     double polIn = 9.) { entry.push_back( Particle(id, status, 0, 0, 0, 0,
0536     col, acol, px, py, pz, e, m, scaleIn, polIn) ); setEvtPtr();
0537     if (col > maxColTag) maxColTag = col;
0538     if (acol > maxColTag) maxColTag = acol;
0539     return entry.size() - 1;
0540   }
0541   int append(int id, int status, int col, int acol, Vec4 p, double m = 0.,
0542     double scaleIn = 0., double polIn = 9.) {entry.push_back( Particle(id,
0543     status, 0, 0, 0, 0, col, acol, p, m, scaleIn, polIn) ); setEvtPtr();
0544     if (col > maxColTag) maxColTag = col;
0545     if (acol > maxColTag) maxColTag = acol;
0546     return entry.size() - 1;
0547   }
0548
0549   // Set pointer to the event for a particle, by default latest one.
0550   void setEvtPtr(int iSet = -1) {if (iSet < 0) iSet = entry.size() - 1;
0551     entry[iSet].setEvtPtr( this);}
0552
0553   // Add a copy of an existing particle at the end of the event record.
0554   int copy(int iCopy, int newStatus = 0);
0555
0556   // List the particles in an event.
0557   void list(bool showScaleAndVertex = false,
0558     bool showMothersAndDaughters = false, int precision = 3) const;
0559
0560   // Remove last n entries.
0561   void popBack(int nRemove = 1) { if (nRemove ==1) entry.pop_back();
0562     else {int newSize = max( 0, size() - nRemove);
0563     entry.resize(newSize);} }
0564
0565   // Remove entries from iFirst to iLast, including endpoints, and fix history.
0566   // (To the extent possible; history pointers in removed range are zeroed.)
0567   void remove(int iFirst, int iLast, bool shiftHistory = true);
0568
0569   // Restore all ParticleDataEntryPtr pointers in the Particle vector.
0570   // Useful when a persistent copy of the event record is read back in.
0571   void restorePtrs() { for (int i = 0; i < size(); ++i) setEvtPtr(i); }
0572
0573   // Save or restore the size of the event record (throwing at the end).
0574   void saveSize() {savedSize = entry.size();}
0575   void restoreSize() {entry.resize(savedSize);}
0576   int  savedSizeValue() {return savedSize;}
0577
0578   // Initialize and access colour tag information.
0579   void initColTag(int colTag = 0) {maxColTag = max( colTag,startColTag);}
0580   int lastColTag() const {return maxColTag;}
0581   int nextColTag() {return ++maxColTag;}
0582
0583   // Access scale for which event as a whole is defined.
0584   void scale( double scaleIn) {scaleSave = scaleIn;}
0585   double scale() const {return scaleSave;}
0586
0587   // Need a second scale if two hard interactions in event.
0588   void scaleSecond( double scaleSecondIn) {scaleSecondSave = scaleSecondIn;}
0589   double scaleSecond() const {return scaleSecondSave;}
0590
0591   // Find complete list of daughters.
0592   // Note: temporarily retained for CMS compatibility. Do not use!
0593   vector<int> daughterList(int i) const {return entry[i].daughterList();}
0594
0595   // Return number of final-state particles, optionally charged only.
0596   int nFinal(bool chargedOnly = false) const {
0597     int nFin = 0;
0598     for (int i = 0; i < size(); ++i)
0599       if (entry[i].isFinal() && (!chargedOnly || entry[i].isCharged()))
0600         ++nFin;
0601     return nFin; }
0602
0603   // Find separation in y, eta, phi or R between two particles.
0604   double dyAbs(int i1, int i2) const {
0605     return abs( entry[i1].y() - entry[i2].y() ); }
0606   double detaAbs(int i1, int i2) const {
0607     return abs( entry[i1].eta() - entry[i2].eta() ); }
0608   double dphiAbs(int i1, int i2) const {
0609     double dPhiTmp = abs( entry[i1].phi() - entry[i2].phi() );
0610     if (dPhiTmp > M_PI) dPhiTmp = 2. * M_PI - dPhiTmp;
0611     return dPhiTmp; }
0612   double RRapPhi(int i1, int i2) const {
0613     return sqrt( pow2(dyAbs(i1, i2)) + pow2(dphiAbs(i1, i2)) ); }
0614   double REtaPhi(int i1, int i2) const {
0615     return sqrt( pow2(detaAbs(i1, i2)) + pow2(dphiAbs(i1, i2)) ); }
0616
0617   // Member functions for rotations and boosts of an event.
0618   void rot(double theta, double phi)
0619     {for (int i = 0; i < size(); ++i) entry[i].rot(theta, phi);}
0620   void bst(double betaX, double betaY, double betaZ)
0621     {for (int i = 0; i < size(); ++i) entry[i].bst(betaX, betaY, betaZ);}
0622   void bst(double betaX, double betaY, double betaZ, double gamma)
0623     {for (int i = 0; i < size(); ++i) entry[i].bst(betaX, betaY, betaZ,
0624     gamma);}
0625   void bst(const Vec4& vec)
0626     {for (int i = 0; i < size(); ++i) entry[i].bst(vec);}
0627   void rotbst(const RotBstMatrix& M, bool boostVertices = true)
0628     {for (int i = 0; i < size(); ++i) entry[i].rotbst(M, boostVertices);}
0629
0630   // Clear the list of junctions.
0631   void clearJunctions() {junction.resize(0);}
0632
0633   // Add a junction to the list, study it or extra input.
0634   int appendJunction( int kind, int col0, int col1, int col2)
0635     { junction.push_back( Junction( kind, col0, col1, col2) );
0636     return junction.size() - 1;}
0637   int appendJunction(Junction junctionIn) {junction.push_back(junctionIn);
0638     return junction.size() - 1;}
0639   int sizeJunction() const {return junction.size();}
0640   bool remainsJunction(int i) const {return junction[i].remains();}
0641   void remainsJunction(int i, bool remainsIn) {junction[i].remains(remainsIn);}
0642   int kindJunction(int i) const {return junction[i].kind();}
0643   int colJunction( int i, int j) const {return junction[i].col(j);}
0644   void colJunction( int i, int j, int colIn) {junction[i].col(j, colIn);}
0645   int endColJunction( int i, int j) const {return junction[i].endCol(j);}
0646   void endColJunction( int i, int j, int endColIn)
0647     {junction[i].endCol(j, endColIn);}
0648   int statusJunction( int i, int j) const {return junction[i].status(j);}
0649   void statusJunction( int i, int j, int statusIn)
0650     {junction[i].status(j, statusIn);}
0651   Junction& getJunction(int i) {return junction[i];}
0652   const Junction& getJunction(int i) const {return junction[i];}
0653   void eraseJunction(int i);
0654
0655   // Save or restore the size of the junction list (throwing at the end).
0656   void saveJunctionSize() {savedJunctionSize = junction.size();}
0657   void restoreJunctionSize() {junction.resize(savedJunctionSize);}
0658
0659   // List any junctions in the event; for debug mainly.
0660   void listJunctions() const;
0661
0662   // Tell whether event has Hidden Valley colours stored.
0663   bool hasHVcols() const {return (hvCols.size() > 0);}
0664
0665   // List any Hidden Valley colours. Find largest HV colour.
0666   void listHVcols() const;
0667   int  maxHVcols()  const;
0668
0669   // Save or restore the size of the HV-coloured list (throwing at the end).
0670   void saveHVcolsSize() {savedHVcolsSize = hvCols.size();}
0671   void restoreHVcolsSize() {hvCols.resize(savedHVcolsSize);}
0672
0673   // Save event record size at Parton Level, i.e. before hadronization.
0674   void savePartonLevelSize() {savedPartonLevelSize = entry.size();}
0675
0676   // Add information about string breaks for a single string.
0677   void addStringBreaks(StringBreaks& sbIn) {stringBreaks.push_back(sbIn);}
0678
0679   // Clear string breaks information.
0680   void clearStringBreaks() {stringBreaks.clear();}
0681
0682   // Get information about all string breaks in the event.
0683   const vector<StringBreaks> &getStringBreaks() const {return stringBreaks;}
0684
0685   // Operator overloading allows to append one event to an existing one.
0686   // Warning: particles should be OK, but some other information unreliable.
0687   Event& operator+=(const Event& addEvent);
0688
0689   // Direct access to the particles via constant pointer.
0690   const vector<Particle>* particles() const {return &entry;}
0691
0692 private:
0693
0694   // The Particle class needs to access particle data.
0695   friend class Particle;
0696
0697   // Constants: could only be changed in the code itself.
0698   static const int IPERLINE;
0699
0700   // Initialization data, normally only set once.
0701   int startColTag;
0702
0703   // The event: a vector containing all particles (entries).
0704   // The explicit use of Pythia8:: qualifier patches a limitation in ROOT.
0705   vector<Pythia8::Particle> entry;
0706
0707   // The list of junctions.
0708   // The explicit use of Pythia8:: qualifier patches a limitation in ROOT.
0709   vector<Pythia8::Junction> junction;
0710
0711   // The list of Hidden-Valley-coloured partons.
0712   // The explicit use of Pythia8:: qualifier patches a limitation in ROOT.
0713   vector<Pythia8::HVcols> hvCols;
0714
0715   // The list of string break information for each string in the event.
0716   vector<Pythia8::StringBreaks> stringBreaks;
0717
0718   // Find index in Hidden Valley list for a particle in the event record.
0719   bool findIndexHV(int iIn) { if (iIn > 0 && iIn == iEventHV) return true;
0720     for (int i = 0; i < int(hvCols.size()); ++i) if (hvCols[i].iHV == iIn)
0721       {iEventHV = iIn; iIndexHV = i; return true; }
0722     return false; }
0723   int iEventHV, iIndexHV;
0724
0725   // Clear the list of Hidden Valley colours. Reset indices.
0726   void clearHV() {hvCols.resize(0); iEventHV = -1; iIndexHV = -1;}
0727
0728   // The maximum colour tag of the event so far.
0729   int maxColTag;
0730
0731   // Saved entry and junction list sizes, for simple restoration.
0732   int savedSize, savedJunctionSize, savedHVcolsSize, savedPartonLevelSize;
0733
0734   // The scale of the event; linear quantity in GeV.
0735   double scaleSave, scaleSecondSave;
0736
0737   // Header specification in event listing (at most 40 characters wide).
0738   string headerList;
0739
0740   // Pointer to the particle data table.
0741   // The //! below is ROOT notation that this member should not be saved.
0742   ParticleData* particleDataPtr;  //!
0743
0744 };
0745
0746 //==========================================================================
0747
0748 } // end namespace Pythia8
0749
0750 #endif // Pythia8_Event_H