dox/html/Pythia8__2__35_2include_2Pythia8_2BeamParticle_8h_source.html

 // BeamParticle.h is a part of the PYTHIA event generator.

 // Copyright (C) 2018 Torbjorn Sjostrand.

 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.

 // Please respect the MCnet Guidelines, see GUIDELINES for details.


 // Header file for information on incoming beams.

 // ResolvedParton: an initiator or remnant in beam.

 // BeamParticle: contains partons, parton densities, etc.


 #ifndef Pythia8_BeamParticle_H

 #define Pythia8_BeamParticle_H


 #include "Pythia8/Basics.h"

 #include "Pythia8/Event.h"

 #include "Pythia8/FragmentationFlavZpT.h"

 #include "Pythia8/Info.h"

 #include "Pythia8/ParticleData.h"

 #include "Pythia8/PartonDistributions.h"

 #include "Pythia8/PythiaStdlib.h"

 #include "Pythia8/Settings.h"


 namespace Pythia8 {


 //==========================================================================


 // This class holds info on a parton resolved inside the incoming beam,

 // i.e. either an initiator (part of a hard or a multiparton interaction)

 // or a remnant (part of the beam remnant treatment).


 // The companion code is -1 from onset and for g, is -2 for an unmatched

 // sea quark, is >= 0 for a matched sea quark, with the number giving the

 // companion position, and is -3 for a valence quark.


 // Rescattering partons properly do not belong here, but bookkeeping is

 // simpler with them, so they are stored with companion code -10.


 class ResolvedParton {


 public:


   // Constructor.

   ResolvedParton( int iPosIn = 0, int idIn = 0, double xIn = 0.,

     int companionIn = -1) : iPosRes(iPosIn), idRes(idIn), xRes(xIn),

     companionRes(companionIn), xqCompRes(0.), mRes(0.), factorRes(1.),

     colRes(0), acolRes(0) { }


   // Set info on initiator or remnant parton.

   void iPos( int iPosIn) {iPosRes = iPosIn;}

   void id( int idIn) {idRes = idIn;}

   void x( double xIn) {xRes = xIn;}

   void update( int iPosIn, int idIn, double xIn) {iPosRes = iPosIn;

     idRes = idIn; xRes = xIn;}

   void companion( int companionIn) {companionRes = companionIn;}

   void xqCompanion( double xqCompIn) {xqCompRes = xqCompIn;}

   void p(Vec4 pIn) {pRes = pIn;}

   void px(double pxIn) {pRes.px(pxIn);}

   void py(double pyIn) {pRes.py(pyIn);}

   void pz(double pzIn) {pRes.pz(pzIn);}

   void e(double eIn) {pRes.e(eIn);}

   void m(double mIn) {mRes = mIn;}

   void col(int colIn) {colRes = colIn;}

   void acol(int acolIn) {acolRes = acolIn;}

   void cols(int colIn = 0,int acolIn = 0)

     {colRes = colIn; acolRes = acolIn;}

   void scalePT( double factorIn) {pRes.px(factorIn * pRes.px());

     pRes.py(factorIn * pRes.py()); factorRes *= factorIn;}

   void scaleX( double factorIn) {xRes *= factorIn;}


   // Get info on initiator or remnant parton.

   int    iPos()        const {return iPosRes;}

   int    id()          const {return idRes;}

   double x()           const {return xRes;}

   int    companion()   const {return companionRes;}

   bool   isValence()   const {return (companionRes == -3);}

   bool   isUnmatched() const {return (companionRes == -2);}

   bool   isCompanion() const {return (companionRes >= 0);}

   bool   isFromBeam()  const {return (companionRes > -10);}

   double xqCompanion() const {return xqCompRes;}

   Vec4   p()           const {return pRes;}

   double px()          const {return pRes.px();}

   double py()          const {return pRes.py();}

   double pz()          const {return pRes.pz();}

   double e()           const {return pRes.e();}

   double m()           const {return mRes;}

   double pT()          const {return pRes.pT();}

   double mT2()         const {return (mRes >= 0.)

     ? mRes*mRes + pRes.pT2() : - mRes*mRes + pRes.pT2();}

   double pPos()        const {return pRes.e() +  pRes.pz();}

   double pNeg()        const {return pRes.e() -  pRes.pz();}

   int    col()         const {return colRes;}

   int    acol()        const {return acolRes;}

   double pTfactor()    const {return factorRes;}

   bool hasCol()        const {return (idRes == 21 || (idRes > 0 && idRes < 9)

     || (-idRes > 1000 && -idRes < 10000 && (-idRes/10)%10 == 0));}

   bool hasAcol()       const {return (idRes == 21 || (-idRes > 0 && -idRes < 9)

     || (idRes > 1000 && idRes < 10000 && (idRes/10)%10 == 0));}


 private:


   // Properties of a resolved parton.

   int    iPosRes, idRes;

   double xRes;

   // Companion code and distribution value, if any.

   int    companionRes;

   double xqCompRes;

   // Four-momentum and mass; for remnant kinematics construction.

   Vec4   pRes;

   double mRes, factorRes;

   // Colour codes.

   int   colRes, acolRes;


 };


 //==========================================================================


 // This class holds info on a beam particle in the evolution of

 // initial-state radiation and multiparton interactions.


 class BeamParticle {


 public:


   // Constructor.

   BeamParticle() : nInit(0) {resolved.resize(0); Q2ValFracSav = -1.;}


   // Initialize data on a beam particle and save pointers.

   void init( int idIn, double pzIn, double eIn, double mIn,

     Info* infoPtrIn, Settings& settings, ParticleData* particleDataPtrIn,

     Rndm* rndmPtrIn, PDF* pdfInPtr, PDF* pdfHardInPtr, bool isUnresolvedIn,

     StringFlav* flavSelPtrIn);


   // Initialize only the two pdf pointers.

   void initPDFPtr(PDF* pdfInPtr, PDF* pdfHardInPtr) {

     pdfBeamPtr = pdfInPtr; pdfHardBeamPtr = pdfHardInPtr; }


   // Initialize additional PDF pointer for unresolved beam.

   void initUnres(PDF* pdfUnresInPtr);


   // For mesons like pi0 valence content varies from event to event.

   void newValenceContent();


   // Set new pZ and E, but keep the rest the same.

   void newPzE( double pzIn, double eIn) {pBeam = Vec4( 0., 0., pzIn, eIn);}


   // Set new mass. Used with photons when virtuality is sampled.

   void newM( double mIn) { mBeam = mIn; }


   // Member functions for output.

   int id()            const {return idBeam;}

   int  idVMD()        const {return idVMDBeam;}

   Vec4 p()            const {return pBeam;}

   double px()         const {return pBeam.px();}

   double py()         const {return pBeam.py();}

   double pz()         const {return pBeam.pz();}

   double e()          const {return pBeam.e();}

   double m()          const {return mBeam;}

   double mVMD()       const {return mVMDBeam;}

   double scaleVMD()   const {return scaleVMDBeam;}

   bool isLepton()     const {return isLeptonBeam;}

   bool isUnresolved() const {return isUnresolvedBeam;}

   // As hadrons here we only count those we know how to handle remnants for.

   bool isHadron()     const {return isHadronBeam;}

   bool isMeson()      const {return isMesonBeam;}

   bool isBaryon()     const {return isBaryonBeam;}

   bool isGamma()      const {return isGammaBeam;}

   bool hasResGamma()  const {return hasResGammaInBeam;}

   bool hasVMDstate()  const {return hasVMDstateInBeam;}


   // Maximum x remaining after previous MPI and ISR, plus safety margin.

   double xMax(int iSkip = -1);


   // Special hard-process parton distributions (can agree with standard ones).

   double xfHard(int idIn, double x, double Q2)

     {return pdfHardBeamPtr->xf(idIn, x, Q2);}


   // Overestimate for PDFs. Same as normal except photons inside leptons.

   double xfMax(int idIn, double x, double Q2)

     {return pdfHardBeamPtr->xfMax(idIn, x, Q2);}


   // Accurate and approximated photon flux and PDFs.

   double xfFlux(int idIn, double x, double Q2)

     {return pdfHardBeamPtr->xfFlux(idIn, x, Q2);}

   double xfApprox(int idIn, double x, double Q2)

     {return pdfHardBeamPtr->xfApprox(idIn, x, Q2);}

   double xfGamma(int idIn, double x, double Q2)

     {return pdfHardBeamPtr->xfGamma(idIn, x, Q2);}


   // Do not sample the x_gamma value to get correct cross section with

   // possible second call.

   double xfSame(int idIn, double x, double Q2)

     {return pdfHardBeamPtr->xfSame(idIn, x, Q2);}


   // Standard parton distributions.

   double xf(int idIn, double x, double Q2)

     {return pdfBeamPtr->xf(idIn, x, Q2);}


   // Ditto, split into valence and sea parts (where gluon counts as sea).

   double xfVal(int idIn, double x, double Q2)

     {return pdfBeamPtr->xfVal(idIn, x, Q2);}

   double xfSea(int idIn, double x, double Q2)

     {return pdfBeamPtr->xfSea(idIn, x, Q2);}


   // Rescaled parton distributions, as needed for MPI and ISR.

   // For ISR also allow split valence/sea, and only return relevant part.

   double xfMPI(int idIn, double x, double Q2)

     {return xfModified(-1, idIn, x, Q2);}

   double xfISR(int indexMPI, int idIn, double x, double Q2)

     {return xfModified( indexMPI, idIn, x, Q2);}


   // Check whether x and Q2 values fall inside the fit bounds (LHAPDF6 only).

   bool insideBounds(double x, double Q2)

     {return pdfBeamPtr->insideBounds(x,Q2);}


   // Access the running alpha_s of a PDF set (LHAPDF6 only).

   double alphaS(double Q2) {return pdfBeamPtr->alphaS(Q2);}


   // Return quark masses used in the PDF fit (LHAPDF6 only).

   double mQuarkPDF(int idIn) {return pdfBeamPtr->mQuarkPDF(idIn);}


   // Return number of members in PDF family (LHAPDF6 only).

   int nMembers() {return pdfBeamPtr->nMembers();}


   // Calculate envelope of PDF predictions

   void calcPDFEnvelope(int idNow, double xNow, double Q2Now, int valSea) {

     pdfBeamPtr->calcPDFEnvelope(idNow,xNow,Q2Now,valSea);}

   void calcPDFEnvelope(pair<int,int> idNows, pair<double,double> xNows,

     double Q2Now, int valSea) {

     pdfBeamPtr->calcPDFEnvelope(idNows,xNows,Q2Now,valSea);}

   PDF::PDFEnvelope getPDFEnvelope() { return pdfBeamPtr->getPDFEnvelope(); }


   // Decide whether chosen quark is valence, sea or companion.

   int pickValSeaComp();


   // Initialize kind of incoming beam particle.

   void initBeamKind();


   // Overload index operator to access a resolved parton from the list.

   ResolvedParton& operator[](int i) {return resolved[i];}

   const ResolvedParton& operator[](int i) const {return resolved[i];}


   // Total number of partons extracted from beam, and initiators only.

   int size() const {return resolved.size();}

   int sizeInit() const {return nInit;}


   // Clear list of resolved partons.

   void clear() {resolved.resize(0); nInit = 0;}


   // Reset variables related to photon beam.

   void resetGamma() {iGamVal = -1; iPosVal = -1; pT2gm2qqbar = 0.;

     isResolvedGamma = (gammaMode == 1) ? true : false;}


   // Reset variables related to photon beam inside a lepton.

   void resetGammaInLepton() {xGm = 1.; kTgamma = 0.; phiGamma = 0.;}


   // Add a resolved parton to list.

   int append( int iPos, int idIn, double x, int companion = -1)

     {resolved.push_back( ResolvedParton( iPos, idIn, x, companion) );

     return resolved.size() - 1;}


   // Remove the last particle from the beam. Reset companion code if needed.

   void popBack() { int iComp = resolved.back().companion();

     resolved.pop_back(); if ( iComp >= 0 ) { iSkipSave = iComp;

       idSave = resolved[iComp].id(); pickValSeaComp(); } }


   // Print extracted parton list; for debug mainly.

   void list() const;


   // How many different flavours, and how many quarks of given flavour.

   int nValenceKinds() const {return nValKinds;}

   int nValence(int idIn) const {for (int i = 0; i < nValKinds; ++i)

       if (idIn == idVal[i]) return nVal[i];

     return 0;}


   // Test whether a lepton is to be considered as unresolved.

   bool isUnresolvedLepton();


   // Add extra remnant flavours to make valence and sea come out right.

   bool remnantFlavours(Event& event, bool isDIS = false);


   // Correlate all initiators and remnants to make a colour singlet.

   bool remnantColours(Event& event, vector<int>& colFrom,

     vector<int>& colTo);


   // Pick unrescaled x of remnant parton (valence or sea).

   double xRemnant(int i);


   // Tell whether a junction has been resolved, and its junction colours.

   bool hasJunction() const {return hasJunctionBeam;}

   int junctionCol(int i) const {return junCol[i];}

   void junctionCol(int i, int col) {junCol[i] = col;}


   // For a diffractive system, decide whether to kick out gluon or quark.

   bool pickGluon(double mDiff);


   // Pick a valence quark at random, and provide the remaining flavour.

   int pickValence();

   int pickRemnant() const {return idVal2;}


   // Share lightcone momentum between two remnants in a diffractive system.

   // At the same time generate a relative pT for the two.

   double zShare( double mDiff, double m1, double m2);

   double pxShare() const {return pxRel;}

   double pyShare() const {return pyRel;}


   // Add extra remnant flavours to make valence and sea come out right.

   bool remnantFlavoursNew(Event& event);


   // Find the colour setup of the removed partons from the scatterings.

   void findColSetup(Event& event);


   // Set initial colours.

   void setInitialCol(Event & event);


   // Update colours.

   void updateCol(vector<pair<int,int> > colourChanges);


   vector<pair <int,int> > getColUpdates() {return colUpdates;}


   // Set valence content for photon beams and position of first valence quark.

   bool gammaInitiatorIsVal(int iResolved, int id, double x, double Q2);

   bool gammaInitiatorIsVal(int iResolved, double Q2);

   int  getGammaValFlavour() { return abs(idVal[0]); }

   int  gammaValSeaComp(int iResolved);

   void posVal(int iPosValIn)          { iPosVal = iPosValIn; }

   void gamVal(int iGamValIn)          { iGamVal = iGamValIn; }

   int  gamVal()                       { return iGamVal; }


   // Set and get the state (resolved and/or unresolved) of photon beam.

   void resolvedGamma(bool isResolved) { isResolvedGamma = isResolved; }

   bool resolvedGamma()                { return isResolvedGamma; }

   void setGammaMode(int gammaModeIn);

   int  getGammaMode()                 { return gammaMode; }

   bool isResolvedUnresolved()         { return isResUnres; }


   // Set state of VMD inside gamma.

   void setVMDstate(bool isVMDIn, int idIn, double mIn, double scaleIn,

     bool reassignState = false) {

     hasVMDstateInBeam = isVMDIn;

     idVMDBeam         = idIn;

     mVMDBeam          = mIn;

     scaleVMDBeam      = scaleIn;

     if (reassignState) {

       idBeam = idVMDBeam;

       mBeam  = mVMDBeam;

       pdfBeamPtr->setVMDscale(scaleVMDBeam);

     }

   }


   // Store the pT2 value of gamma->qqbar splitting.

   void   pT2gamma2qqbar(double pT2in) { pT2gm2qqbar = pT2in; }

   double pT2gamma2qqbar()             { return pT2gm2qqbar; }


   // Store the pT value for the latest MPI.

   void   pTMPI(double pTminMPIin)     { pTminMPI = pTminMPIin; }


   // Check whether room for beam remnants.

   bool roomFor1Remnant(double eCM);

   bool roomFor1Remnant(int id1, double x1, double eCM);

   bool roomFor2Remnants(int id1, double x1, double eCM);

   bool roomForRemnants(BeamParticle beamOther);


   // Evaluate the remnant mass with initiator idIn.

   double remnantMass(int idIn);


   // Functions to approximate pdfs for ISR.

   double gammaPDFxDependence(int flavour, double x)

     { return pdfBeamPtr->gammaPDFxDependence(flavour, x); }

   double gammaPDFRefScale(int flavour)

     { return pdfBeamPtr->gammaPDFRefScale(flavour); }

   double xIntegratedPDFs(double Q2)

     { return pdfBeamPtr->xfIntegratedTotal(Q2); }


   // Save the x_gamma value after latest PDF call or set it later if ND.

   void xGammaPDF()            { xGm = pdfHardBeamPtr->xGamma(); }

   void xGamma(double xGmIn)   { xGm = xGmIn; }

   void Q2Gamma(double Q2GmIn) { Q2gm = Q2GmIn; }

   void newGammaKTPhi(double kTIn, double phiIn)

     { kTgamma = kTIn; phiGamma = phiIn; }


   // Get the kinematic limits for photons emitted by the beam.

   double Q2minPDF()      { return pdfHardBeamPtr->getQ2min(); }

   double xGammaMin()     { return pdfHardBeamPtr->getXmin(); }

   double xGammaHadr()    { return pdfHardBeamPtr->getXhadr(); }

   double gammaFluxNorm() { return pdfHardBeamPtr->getGammaFluxNorm(); }


   // Get the kinematics related photons form lepton beams.

   double xGamma()   const { return xGm; }

   double Q2Gamma()  const { return Q2gm; }

   double gammaKTx() const { return kTgamma*cos(phiGamma); }

   double gammaKTy() const { return kTgamma*sin(phiGamma); }

   double gammaKT()  const { return kTgamma; }

   double gammaPhi() const { return phiGamma; }


   // Keep track of pomeron momentum fraction.

   void xPom(double xpom = -1.0)

     { if ( pdfBeamPtr ) pdfBeamPtr->xPom(xpom); }


   // Sample x and Q2 for emitted photons according to flux.

   double sampleXgamma(double xMinIn)

     { xGm = pdfHardBeamPtr->sampleXgamma(xMinIn); return xGm; }

   double sampleQ2gamma(double Q2min)

     { Q2gm = pdfHardBeamPtr->sampleQ2gamma(Q2min); return Q2gm;}


 private:


   // Constants: could only be changed in the code itself.

   static const double XMINUNRESOLVED, POMERONMASS, XMAXCOMPANION, TINYZREL;

   static const int NMAX, NRANDOMTRIES;


   // Pointer to various information on the generation.

   Info*         infoPtr;


   // Pointer to the particle data table.

   ParticleData* particleDataPtr;


   // Pointer to the random number generator.

   Rndm*         rndmPtr;


   // Pointers to PDF sets.

   PDF*          pdfBeamPtr;

   PDF*          pdfHardBeamPtr;


   // Pointer to unresolved PDF and two others to save the resolved ptrs.

   PDF*          pdfUnresBeamPtr;

   PDF*          pdfBeamPtrSave;

   PDF*          pdfHardBeamPtrSave;


   // Pointer to class for flavour generation.

   StringFlav*   flavSelPtr;


   // Initialization data, normally only set once.

   bool   allowJunction, beamJunction;

   int    maxValQuark, companionPower;

   double valencePowerMeson, valencePowerUinP, valencePowerDinP,

          valenceDiqEnhance, pickQuarkNorm, pickQuarkPower,

          diffPrimKTwidth, diffLargeMassSuppress, beamSat, gluonPower,

          xGluonCutoff;


   // Basic properties of a beam particle.

   int    idBeam, idBeamAbs, idVMDBeam;

   Vec4   pBeam;

   double mBeam, mVMDBeam, scaleVMDBeam;


   // Beam kind. Valence flavour content for hadrons.

   bool   isUnresolvedBeam, isLeptonBeam, isHadronBeam, isMesonBeam,

          isBaryonBeam, isGammaBeam;

   int    nValKinds, idVal[3], nVal[3];


   // Current parton density, by valence, sea and companion.

   int    idSave, iSkipSave, nValLeft[3];

   double xqgTot, xqVal, xqgSea, xqCompSum;


   // Variables related to photon beams (also inside lepton).

   bool   doISR, doMPI, doND, isResolvedGamma, hasResGammaInBeam,

          isResUnres, hasVMDstateInBeam;

   double pTminISR, pTminMPI, pT2gm2qqbar;

   int    iGamVal, iPosVal, gammaMode;


   // Variables for photon from lepton.

   double xGm, Q2gm, kTgamma, phiGamma;


   // The list of resolved partons.

   vector<ResolvedParton> resolved;


   // Status after all initiators have been accounted for. Junction content.

   int    nInit;

   bool   hasJunctionBeam;

   int    junCol[3];


   // Variables for new colour reconnection;

   pair <int,int> colSetup;

   vector<int> acols, cols;

   vector<bool> usedCol,usedAcol;

   vector< pair<int,int> > colUpdates;

   int nJuncs, nAjuncs, nDiffJuncs;

   bool allowBeamJunctions;


   // Routine to calculate pdf's given previous interactions.

   double xfModified( int iSkip, int idIn, double x, double Q2);


   // Fraction of hadron momentum sitting in a valence quark distribution.

   double xValFrac(int j, double Q2);

   double Q2ValFracSav, uValInt, dValInt;


   // Fraction of hadron momentum sitting in a companion quark distribution.

   double xCompFrac(double xs);


   // Value of companion quark PDF, also given the sea quark x.

   double xCompDist(double xc, double xs);


   // Valence quark subdivision for diffractive systems.

   int    idVal1, idVal2, idVal3;

   double zRel, pxRel, pyRel;


   // Update a single (anti-) colour of the event.

   void updateSingleCol(int oldCol, int newCol);


   // Find a single (anti-) colour in the beam,

   // if a beam remnant is set the new colour.

   int findSingleCol(Event& event, bool isAcol, bool useHardScatters);


 };


 //==========================================================================


 } // end namespace Pythia8


 #endif // Pythia8_BeamParticle_H

Event
Definition: AgUStep.h:26