PYTHIA8/pythia8170/include/BeamParticle.h

   1 // BeamParticle.h is a part of the PYTHIA event generator.
   2 // Copyright (C) 2012 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 information on incoming beams.
   7 // ResolvedParton: an initiator or remnant in beam.
   8 // BeamParticle: contains partons, parton densities, etc.
   9
  10 #ifndef Pythia8_BeamParticle_H
  11 #define Pythia8_BeamParticle_H
  12
  13 #include "Basics.h"
  14 #include "Event.h"
  15 #include "FragmentationFlavZpT.h"
  16 #include "Info.h"
  17 #include "ParticleData.h"
  18 #include "PartonDistributions.h"
  19 #include "PythiaStdlib.h"
  20 #include "Settings.h"
  21
  22 namespace Pythia8 {
  23
  24 //==========================================================================
  25
  26 // This class holds info on a parton resolved inside the incoming beam,
  27 // i.e. either an initiator (part of a hard or a multiparton interaction)
  28 // or a remnant (part of the beam remnant treatment).
  29
  30 // The companion code is -1 from onset and for g, is -2 for an unmatched
  31 // sea quark, is >= 0 for a matched sea quark, with the number giving the
  32 // companion position, and is -3 for a valence quark.
  33
  34 // Rescattering partons properly do not belong here, but bookkeeping is
  35 // simpler with them, so they are stored with companion code -10.
  36
  37 class ResolvedParton {
  38
  39 public:
  40
  41   // Constructor.
  42   ResolvedParton( int iPosIn = 0, int idIn = 0, double xIn = 0.,
  43     int companionIn = -1) : iPosRes(iPosIn), idRes(idIn), xRes(xIn),
  44     companionRes(companionIn), xqCompRes(0.), mRes(0.), factorRes(1.),
  45     colRes(0), acolRes(0) { }
  46
  47   // Set info on initiator or remnant parton.
  48   void iPos( int iPosIn) {iPosRes = iPosIn;}
  49   void id( int idIn) {idRes = idIn;}
  50   void x( double xIn) {xRes = xIn;}
  51   void update( int iPosIn, int idIn, double xIn) {iPosRes = iPosIn;
  52     idRes = idIn; xRes = xIn;}
  53   void companion( int companionIn) {companionRes = companionIn;}
  54   void xqCompanion( double xqCompIn) {xqCompRes = xqCompIn;}
  55   void p(Vec4 pIn) {pRes = pIn;}
  56   void px(double pxIn) {pRes.px(pxIn);}
  57   void py(double pyIn) {pRes.py(pyIn);}
  58   void pz(double pzIn) {pRes.pz(pzIn);}
  59   void e(double eIn) {pRes.e(eIn);}
  60   void m(double mIn) {mRes = mIn;}
  61   void col(int colIn) {colRes = colIn;}
  62   void acol(int acolIn) {acolRes = acolIn;}
  63   void cols(int colIn = 0,int acolIn = 0)
  64     {colRes = colIn; acolRes = acolIn;}
  65   void scalePT( double factorIn) {pRes.px(factorIn * pRes.px());
  66     pRes.py(factorIn * pRes.py()); factorRes *= factorIn;}
  67   void scaleX( double factorIn) {xRes *= factorIn;}
  68
  69   // Get info on initiator or remnant parton.
  70   int    iPos()        const {return iPosRes;}
  71   int    id()          const {return idRes;}
  72   double x()           const {return xRes;}
  73   int    companion()   const {return companionRes;}
  74   bool   isValence()   const {return (companionRes == -3);}
  75   bool   isUnmatched() const {return (companionRes == -2);}
  76   bool   isCompanion() const {return (companionRes >= 0);}
  77   bool   isFromBeam()  const {return (companionRes > -10);}
  78   double xqCompanion() const {return xqCompRes;}
  79   Vec4   p()           const {return pRes;}
  80   double px()          const {return pRes.px();}
  81   double py()          const {return pRes.py();}
  82   double pz()          const {return pRes.pz();}
  83   double e()           const {return pRes.e();}
  84   double m()           const {return mRes;}
  85   double pT()          const {return pRes.pT();}
  86   double mT2()         const {return (mRes >= 0.)
  87     ? mRes*mRes + pRes.pT2() : - mRes*mRes + pRes.pT2();}
  88   double pPos()        const {return pRes.e() +  pRes.pz();}
  89   double pNeg()        const {return pRes.e() -  pRes.pz();}
  90   int    col()         const {return colRes;}
  91   int    acol()        const {return acolRes;}
  92   double pTfactor()    const {return factorRes;}
  93
  94 private:
  95
  96   // Properties of a resolved parton.
  97   int    iPosRes, idRes;
  98   double xRes;
  99   // Companion code and distribution value, if any.
 100   int    companionRes;
 101   double xqCompRes;
 102   // Four-momentum and mass; for remnant kinematics construction.
 103   Vec4   pRes;
 104   double mRes, factorRes;
 105   // Colour codes.
 106   int   colRes, acolRes;
 107
 108 };
 109
 110 //==========================================================================
 111
 112 // This class holds info on a beam particle in the evolution of
 113 // initial-state radiation and multiparton interactions.
 114
 115 class BeamParticle {
 116
 117 public:
 118
 119   // Constructor.
 120   BeamParticle() : nInit(0) {resolved.resize(0); Q2ValFracSav = -1.;}
 121
 122   // Initialize data on a beam particle and save pointers.
 123   void init( int idIn, double pzIn, double eIn, double mIn,
 124     Info* infoPtrIn, Settings& settings, ParticleData* particleDataPtrIn,
 125     Rndm* rndmPtrIn, PDF* pdfInPtr, PDF* pdfHardInPtr, bool isUnresolvedIn,
 126     StringFlav* flavSelPtrIn);
 127
 128   // For mesons like pi0 valence content varies from event to event.
 129   void newValenceContent();
 130
 131   // Set new pZ and E, but keep the rest the same.
 132   void newPzE( double pzIn, double eIn) {pBeam = Vec4( 0., 0., pzIn, eIn);}
 133
 134   // Member functions for output.
 135   int id() const {return idBeam;}
 136   Vec4 p() const {return pBeam;}
 137   double px() const {return pBeam.px();}
 138   double py() const {return pBeam.py();}
 139   double pz() const {return pBeam.pz();}
 140   double e() const {return pBeam.e();}
 141   double m() const {return mBeam;}
 142   bool isLepton() const {return isLeptonBeam;}
 143   bool isUnresolved() const {return isUnresolvedBeam;}
 144   // As hadrons here we only count those we know how to handle remnants for.
 145   bool isHadron() const {return isHadronBeam;}
 146   bool isMeson() const {return isMesonBeam;}
 147   bool isBaryon() const {return isBaryonBeam;}
 148
 149   // Maximum x remaining after previous MPI and ISR, plus safety margin.
 150   double xMax(int iSkip = -1);
 151
 152   // Special hard-process parton distributions (can agree with standard ones).
 153   double xfHard(int idIn, double x, double Q2)
 154     {return pdfHardBeamPtr->xf(idIn, x, Q2);}
 155
 156   // Standard parton distributions.
 157   double xf(int idIn, double x, double Q2)
 158     {return pdfBeamPtr->xf(idIn, x, Q2);}
 159
 160   // Ditto, split into valence and sea parts (where gluon counts as sea).
 161   double xfVal(int idIn, double x, double Q2)
 162     {return pdfBeamPtr->xfVal(idIn, x, Q2);}
 163   double xfSea(int idIn, double x, double Q2)
 164     {return pdfBeamPtr->xfSea(idIn, x, Q2);}
 165
 166   // Rescaled parton distributions, as needed for MPI and ISR.
 167   // For ISR also allow split valence/sea, and only return relevant part.
 168   double xfMPI(int idIn, double x, double Q2)
 169     {return xfModified(-1, idIn, x, Q2);}
 170   double xfISR(int indexMPI, int idIn, double x, double Q2)
 171     {return xfModified( indexMPI, idIn, x, Q2);}
 172
 173   // Decide whether chosen quark is valence, sea or companion.
 174   int pickValSeaComp();
 175
 176   // Initialize kind of incoming beam particle.
 177   void initBeamKind();
 178
 179   // Overload index operator to access a resolved parton from the list.
 180   ResolvedParton& operator[](int i) {return resolved[i];}
 181
 182   // Total number of partons extracted from beam, and initiators only.
 183   int size() const {return resolved.size();}
 184   int sizeInit() const {return nInit;}
 185
 186   // Clear list of resolved partons.
 187   void clear() {resolved.resize(0); nInit = 0;}
 188
 189   // Add a resolved parton to list.
 190   int append( int iPos, int idIn, double x, int companion = -1)
 191     {resolved.push_back( ResolvedParton( iPos, idIn, x, companion) );
 192     return resolved.size() - 1;}
 193
 194   // Print extracted parton list; for debug mainly.
 195   void list(ostream& os = cout) const;
 196
 197   // How many different flavours, and how many quarks of given flavour.
 198   int nValenceKinds() const {return nValKinds;}
 199   int nValence(int idIn) const {for (int i = 0; i < nValKinds; ++i)
 200     if (idIn == idVal[i]) return nVal[i]; return 0;}
 201
 202   // Test whether a lepton is to be considered as unresolved.
 203   bool isUnresolvedLepton();
 204
 205   // Add extra remnant flavours to make valence and sea come out right.
 206   bool remnantFlavours(Event& event);
 207
 208   // Correlate all initiators and remnants to make a colour singlet.
 209   bool remnantColours(Event& event, vector<int>& colFrom,
 210     vector<int>& colTo);
 211
 212   // Pick unrescaled x of remnant parton (valence or sea).
 213   double xRemnant(int i);
 214
 215   // Tell whether a junction has been resolved, and its junction colours.
 216   bool hasJunction() const {return hasJunctionBeam;}
 217   int junctionCol(int i) const {return junCol[i];}
 218   void junctionCol(int i, int col) {junCol[i] = col;}
 219
 220   // For a diffractive system, decide whether to kick out gluon or quark.
 221   bool pickGluon(double mDiff);
 222
 223   // Pick a valence quark at random, and provide the remaining flavour.
 224   int pickValence();
 225   int pickRemnant() const {return idVal2;}
 226
 227   // Share lightcone momentum between two remnants in a diffractive system.
 228   // At the same time generate a relative pT for the two.
 229   double zShare( double mDiff, double m1, double m2);
 230   double pxShare() const {return pxRel;}
 231   double pyShare() const {return pyRel;}
 232
 233 private:
 234
 235   // Constants: could only be changed in the code itself.
 236   static const double XMINUNRESOLVED;
 237
 238   // Pointer to various information on the generation.
 239   Info*         infoPtr;
 240
 241   // Pointer to the particle data table.
 242   ParticleData* particleDataPtr;
 243
 244   // Pointer to the random number generator.
 245   Rndm*         rndmPtr;
 246
 247   // Pointers to PDF sets.
 248   PDF*          pdfBeamPtr;
 249   PDF*          pdfHardBeamPtr;
 250
 251   // Pointer to class for flavour generation.
 252   StringFlav*   flavSelPtr;
 253
 254   // Initialization data, normally only set once.
 255   bool   allowJunction;
 256   int    maxValQuark, companionPower;
 257   double valencePowerMeson, valencePowerUinP, valencePowerDinP,
 258          valenceDiqEnhance, pickQuarkNorm, pickQuarkPower,
 259          diffPrimKTwidth, diffLargeMassSuppress;
 260
 261   // Basic properties of a beam particle.
 262   int    idBeam, idBeamAbs;
 263   Vec4   pBeam;
 264   double mBeam;
 265   // Beam kind. Valence flavour content for hadrons.
 266   bool   isUnresolvedBeam, isLeptonBeam, isHadronBeam, isMesonBeam,
 267          isBaryonBeam;
 268   int    nValKinds, idVal[3], nVal[3];
 269
 270   // Current parton density, by valence, sea and companion.
 271   int    idSave, iSkipSave, nValLeft[3];
 272   double xqgTot, xqVal, xqgSea, xqCompSum;
 273
 274   // The list of resolved partons.
 275   vector<ResolvedParton> resolved;
 276
 277   // Status after all initiators have been accounted for. Junction content.
 278   int    nInit;
 279   bool   hasJunctionBeam;
 280   int    junCol[3];
 281
 282   // Routine to calculate pdf's given previous interactions.
 283   double xfModified( int iSkip, int idIn, double x, double Q2);
 284
 285   // Fraction of hadron momentum sitting in a valence quark distribution.
 286   double xValFrac(int j, double Q2);
 287   double Q2ValFracSav, uValInt, dValInt;
 288
 289   // Fraction of hadron momentum sitting in a companion quark distribution.
 290   double xCompFrac(double xs);
 291
 292   // Value of companion quark PDF, also given the sea quark x.
 293   double xCompDist(double xc, double xs);
 294
 295   // Valence quark subdivision for diffractive systems.
 296   int    idVal1, idVal2, idVal3;
 297   double zRel, pxRel, pyRel;
 298
 299 };
 300
 301 //==========================================================================
 302
 303 } // end namespace Pythia8
 304
 305 #endif // Pythia8_BeamParticle_H