PYTHIA  8.312
BeamParticle.h
1 // BeamParticle.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, 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 "Pythia8/Basics.h"
14 #include "Pythia8/Event.h"
15 #include "Pythia8/FragmentationFlavZpT.h"
16 #include "Pythia8/Info.h"
17 #include "Pythia8/ParticleData.h"
18 #include "Pythia8/PartonDistributions.h"
19 #include "Pythia8/PhysicsBase.h"
20 #include "Pythia8/PythiaStdlib.h"
21 #include "Pythia8/Settings.h"
22 
23 namespace Pythia8 {
24 
25 //==========================================================================
26 
27 // This class holds info on a parton resolved inside the incoming beam,
28 // i.e. either an initiator (part of a hard or a multiparton interaction)
29 // or a remnant (part of the beam remnant treatment).
30 
31 // The companion code is -1 from onset and for g, is -2 for an unmatched
32 // sea quark, is >= 0 for a matched sea quark, with the number giving the
33 // companion position, and is -3 for a valence quark.
34 
35 // Rescattering partons properly do not belong here, but bookkeeping is
36 // simpler with them, so they are stored with companion code -10.
37 
39 
40 public:
41 
42  // Constructor.
43  ResolvedParton( int iPosIn = 0, int idIn = 0, double xIn = 0.,
44  int companionIn = -1) : iPosRes(iPosIn), idRes(idIn), xRes(xIn),
45  companionRes(companionIn), xqCompRes(0.), mRes(0.), factorRes(1.),
46  colRes(0), acolRes(0) { }
47 
48  // Set info on initiator or remnant parton.
49  void iPos( int iPosIn) {iPosRes = iPosIn;}
50  void id( int idIn) {idRes = idIn;}
51  void x( double xIn) {xRes = xIn;}
52  void update( int iPosIn, int idIn, double xIn) {iPosRes = iPosIn;
53  idRes = idIn; xRes = xIn;}
54  void companion( int companionIn) {companionRes = companionIn;}
55  void xqCompanion( double xqCompIn) {xqCompRes = xqCompIn;}
56  void p(Vec4 pIn) {pRes = pIn;}
57  void px(double pxIn) {pRes.px(pxIn);}
58  void py(double pyIn) {pRes.py(pyIn);}
59  void pz(double pzIn) {pRes.pz(pzIn);}
60  void e(double eIn) {pRes.e(eIn);}
61  void m(double mIn) {mRes = mIn;}
62  void col(int colIn) {colRes = colIn;}
63  void acol(int acolIn) {acolRes = acolIn;}
64  void cols(int colIn = 0,int acolIn = 0)
65  {colRes = colIn; acolRes = acolIn;}
66  void scalePT( double factorIn) {pRes.px(factorIn * pRes.px());
67  pRes.py(factorIn * pRes.py()); factorRes *= factorIn;}
68  void scaleX( double factorIn) {xRes *= factorIn;}
69 
70  // Get info on initiator or remnant parton.
71  int iPos() const {return iPosRes;}
72  int id() const {return idRes;}
73  double x() const {return xRes;}
74  int companion() const {return companionRes;}
75  bool isValence() const {return (companionRes == -3);}
76  bool isUnmatched() const {return (companionRes == -2);}
77  bool isCompanion() const {return (companionRes >= 0);}
78  bool isFromBeam() const {return (companionRes > -10);}
79  double xqCompanion() const {return xqCompRes;}
80  Vec4 p() const {return pRes;}
81  double px() const {return pRes.px();}
82  double py() const {return pRes.py();}
83  double pz() const {return pRes.pz();}
84  double e() const {return pRes.e();}
85  double m() const {return mRes;}
86  double pT() const {return pRes.pT();}
87  double mT2() const {return (mRes >= 0.)
88  ? mRes*mRes + pRes.pT2() : - mRes*mRes + pRes.pT2();}
89  double pPos() const {return pRes.e() + pRes.pz();}
90  double pNeg() const {return pRes.e() - pRes.pz();}
91  int col() const {return colRes;}
92  int acol() const {return acolRes;}
93  double pTfactor() const {return factorRes;}
94  bool hasCol() const {return (idRes == 21 || (idRes > 0 && idRes < 9)
95  || (-idRes > 1000 && -idRes < 10000 && (-idRes/10)%10 == 0));}
96  bool hasAcol() const {return (idRes == 21 || (-idRes > 0 && -idRes < 9)
97  || (idRes > 1000 && idRes < 10000 && (idRes/10)%10 == 0));}
98 
99 private:
100 
101  // Properties of a resolved parton.
102  int iPosRes, idRes;
103  double xRes;
104  // Companion code and distribution value, if any.
105  int companionRes;
106  double xqCompRes;
107  // Four-momentum and mass; for remnant kinematics construction.
108  Vec4 pRes;
109  double mRes, factorRes;
110  // Colour codes.
111  int colRes, acolRes;
112 
113 };
114 
115 //==========================================================================
116 
117 // The xfModPrepData struct saves information on the amount of x already used,
118 // to speed up PDF evaluation.
119 
120 typedef struct {
121  double xValTot;
122  double xValLeft;
123  double xLeft;
124  double xCompAdded;
125  double rescaleGS;
126 } xfModPrepData;
127 
128 //==========================================================================
129 
130 // This class holds info on a beam particle in the evolution of
131 // initial-state radiation and multiparton interactions.
132 
133 class BeamParticle : public PhysicsBase {
134 
135 public:
136 
137  // Constructor.
138  BeamParticle() : pdfBeamPtr(), pdfHardBeamPtr(), pdfUnresBeamPtr(),
139  pdfBeamPtrSave(), pdfHardBeamPtrSave(), pdfSavePtrs(),
140  pdfSaveIdx(-1), flavSelPtr(), allowJunction(), beamJunction(),
141  maxValQuark(), companionPower(), valencePowerMeson(), valencePowerUinP(),
142  valencePowerDinP(), valenceDiqEnhance(), pickQuarkNorm(), pickQuarkPower(),
143  diffPrimKTwidth(), diffLargeMassSuppress(), beamSat(), gluonPower(),
144  xGluonCutoff(), heavyQuarkEnhance(), idBeam(), idBeamAbs(), idVMDBeam(),
145  mBeam(), mVMDBeam(), scaleVMDBeam(), isUnresolvedBeam(), isLeptonBeam(),
146  isHadronBeam(), isMesonBeam(), isBaryonBeam(), isGammaBeam(), nValKinds(),
147  idVal(), nVal(), idSave(), iSkipSave(), nValLeft(), xqgTot(), xqVal(),
148  xqgSea(), xqCompSum(), doISR(), doMPI(), doND(), isResolvedGamma(),
149  hasResGammaInBeam(), isResUnres(), hasVMDstateInBeam(), initGammaBeam(),
150  pTminISR(), pTminMPI(), pT2gm2qqbar(), iGamVal(), iPosVal(), gammaMode(),
151  xGm(), Q2gm(), kTgamma(), phiGamma(), cPowerCache(-100), xsCache(-1),
152  resCache(), resolved(), nInit(0), hasJunctionBeam(), junCol(), nJuncs(),
153  nAjuncs(), nDiffJuncs(), allowBeamJunctions(), Q2ValFracSav(-1.),
154  uValInt(), dValInt(), idVal1(), idVal2(), idVal3(), zRel(), pxRel(),
155  pyRel() { }
156 
157  // Initialize data on a beam particle and save pointers.
158  void init( int idIn, double pzIn, double eIn, double mIn,
159  PDFPtr pdfInPtr, PDFPtr pdfHardInPtr, bool isUnresolvedIn,
160  StringFlav* flavSelPtrIn);
161 
162  // Initialize only the id.
163  void initID( int idIn) { idBeam = idIn; initBeamKind();}
164 
165  // Initialize only the two pdf pointers.
166  void initPDFPtr(PDFPtr pdfInPtr, PDFPtr pdfHardInPtr) {
167  pdfBeamPtr = pdfInPtr; pdfHardBeamPtr = pdfHardInPtr; }
168 
169  // Initialize additional PDF pointer for unresolved beam.
170  void initUnres(PDFPtr pdfUnresInPtr);
171 
172  // Initialize array of PDFs for switching between them.
173  void initSwitchID( const vector<PDFPtr>& pdfSavePtrsIn) {
174  pdfSavePtrs = pdfSavePtrsIn; }
175 
176  // For mesons like pi0 valence content varies from event to event.
177  void newValenceContent();
178  void setValenceContent(int idq1, int idq2 = 0, int idq3 = 0);
179 
180  // Switch to new beam particle identity; for similar hadrons only.
181  void setBeamID( int idIn, int iPDFin = -1) {idBeam = idIn;
182  if ( iPDFin >= 0 && iPDFin < int(pdfSavePtrs.size())
183  && iPDFin != pdfSaveIdx ) {
184  pdfBeamPtr = pdfSavePtrs[iPDFin]; pdfHardBeamPtr = pdfBeamPtr;
185  pdfSaveIdx = iPDFin; }
186  mBeam = particleDataPtr->m0(idIn); pdfBeamPtr->setBeamID(idIn); }
187 
188  // Set new pZ and E, but keep the rest the same.
189  void newPzE( double pzIn, double eIn) {pBeam = Vec4( 0., 0., pzIn, eIn);}
190 
191  // Set new mass. Used with photons when virtuality is sampled.
192  void newM( double mIn) { mBeam = mIn; }
193 
194  // Member functions for output.
195  int id() const {return idBeam;}
196  int idVMD() const {return idVMDBeam;}
197  Vec4 p() const {return pBeam;}
198  double px() const {return pBeam.px();}
199  double py() const {return pBeam.py();}
200  double pz() const {return pBeam.pz();}
201  double e() const {return pBeam.e();}
202  double m() const {return mBeam;}
203  double mVMD() const {return mVMDBeam;}
204  double scaleVMD() const {return scaleVMDBeam;}
205  bool isLepton() const {return isLeptonBeam;}
206  bool isUnresolved() const {return isUnresolvedBeam;}
207  // As hadrons here we only count those we know how to handle remnants for.
208  bool isHadron() const {return isHadronBeam;}
209  bool isMeson() const {return isMesonBeam;}
210  bool isBaryon() const {return isBaryonBeam;}
211  bool isGamma() const {return isGammaBeam;}
212  bool hasResGamma() const {return hasResGammaInBeam;}
213  bool hasVMDstate() const {return hasVMDstateInBeam;}
214 
215  // Maximum x remaining after previous MPI and ISR, plus safety margin.
216  double xMax(int iSkip = -1);
217 
218  // Special hard-process parton distributions (can agree with standard ones).
219  double xfHard(int idIn, double x, double Q2)
220  {return pdfHardBeamPtr->xf(idIn, x, Q2);}
221 
222  // Overestimate for PDFs. Same as normal except photons inside leptons.
223  double xfMax(int idIn, double x, double Q2)
224  {return pdfHardBeamPtr->xfMax(idIn, x, Q2);}
225 
226  // Accurate and approximated photon flux and PDFs.
227  double xfFlux(int idIn, double x, double Q2)
228  {return pdfHardBeamPtr->xfFlux(idIn, x, Q2);}
229  double xfApprox(int idIn, double x, double Q2)
230  {return pdfHardBeamPtr->xfApprox(idIn, x, Q2);}
231  double xfGamma(int idIn, double x, double Q2)
232  {return pdfHardBeamPtr->xfGamma(idIn, x, Q2);}
233 
234  // Do not sample the x_gamma value to get correct cross section with
235  // possible second call.
236  double xfSame(int idIn, double x, double Q2)
237  {return pdfHardBeamPtr->xfSame(idIn, x, Q2);}
238 
239  // Standard parton distributions.
240  double xf(int idIn, double x, double Q2)
241  {return pdfBeamPtr->xf(idIn, x, Q2);}
242 
243  // Ditto, split into valence and sea parts (where gluon counts as sea).
244  double xfVal(int idIn, double x, double Q2)
245  {return pdfBeamPtr->xfVal(idIn, x, Q2);}
246  double xfSea(int idIn, double x, double Q2)
247  {return pdfBeamPtr->xfSea(idIn, x, Q2);}
248 
249  // Rescaled parton distributions, as needed for MPI and ISR.
250  // For ISR also allow split valence/sea, and only return relevant part.
251  double xfMPI(int idIn, double x, double Q2, xfModPrepData& xfData)
252  {return xfISR(-1, idIn, x, Q2, xfData);}
253  double xfMPI(int idIn, double x, double Q2)
254  {xfModPrepData xfData = xfModPrep(-1, Q2);
255  return xfISR(-1, idIn, x, Q2, xfData);}
256  double xfISR(int indexMPI, int idIn, double x, double Q2,
257  xfModPrepData& xfData) {return xfModified( indexMPI, idIn, x, Q2, xfData);}
258  double xfISR(int indexMPI, int idIn, double x, double Q2)
259  {xfModPrepData xfData= xfModPrep(indexMPI, Q2);
260  return xfModified( indexMPI, idIn, x, Q2, xfData);}
261 
262  // Check whether x and Q2 values fall inside the fit bounds (LHAPDF6 only).
263  bool insideBounds(double x, double Q2)
264  {return pdfBeamPtr->insideBounds(x,Q2);}
265 
266  // Access the running alpha_s of a PDF set (LHAPDF6 only).
267  double alphaS(double Q2) {return pdfBeamPtr->alphaS(Q2);}
268 
269  // Return quark masses used in the PDF fit (LHAPDF6 only).
270  double mQuarkPDF(int idIn) {return pdfBeamPtr->mQuarkPDF(idIn);}
271 
272  // Return number of members in PDF family (LHAPDF6 only).
273  int nMembers() {return pdfBeamPtr->nMembers();}
274 
275  // Calculate envelope of PDF predictions
276  void calcPDFEnvelope(int idNow, double xNow, double Q2Now, int valSea) {
277  pdfBeamPtr->calcPDFEnvelope(idNow,xNow,Q2Now,valSea);}
278  void calcPDFEnvelope(pair<int,int> idNows, pair<double,double> xNows,
279  double Q2Now, int valSea) {
280  pdfBeamPtr->calcPDFEnvelope(idNows,xNows,Q2Now,valSea);}
281  PDF::PDFEnvelope getPDFEnvelope() { return pdfBeamPtr->getPDFEnvelope(); }
282 
283  // Decide whether chosen quark is valence, sea or companion.
284  int pickValSeaComp();
285 
286  // Initialize kind of incoming beam particle.
287  void initBeamKind();
288 
289  // Overload index operator to access a resolved parton from the list.
290  ResolvedParton& operator[](int i) {return resolved[i];}
291  const ResolvedParton& operator[](int i) const {return resolved[i];}
292 
293  // Total number of partons extracted from beam, and initiators only.
294  int size() const {return resolved.size();}
295  int sizeInit() const {return nInit;}
296 
297  // Clear list of resolved partons.
298  void clear() {resolved.resize(0); nInit = 0;}
299 
300  // Reset variables related to photon beam.
301  void resetGamma() {iGamVal = -1; iPosVal = -1; pT2gm2qqbar = 0.;
302  isResolvedGamma = (gammaMode == 1) ? true : false;}
303 
304  // Reset variables related to photon beam inside a lepton.
305  void resetGammaInLepton() {xGm = 1.; kTgamma = 0.; phiGamma = 0.;}
306 
307  // Add a resolved parton to list.
308  int append( int iPos, int idIn, double x, int companion = -1)
309  {resolved.push_back( ResolvedParton( iPos, idIn, x, companion) );
310  return resolved.size() - 1;}
311 
312  // Remove the last particle from the beam. Reset companion code if needed.
313  void popBack() { int iComp = resolved.back().companion();
314  resolved.pop_back(); if ( iComp >= 0 ) { iSkipSave = iComp;
315  idSave = resolved[iComp].id(); pickValSeaComp(); } }
316 
317  // Print extracted parton list; for debug mainly.
318  void list() const;
319 
320  // How many different flavours, and how many quarks of given flavour.
321  int nValenceKinds() const {return nValKinds;}
322  int nValence(int idIn) const {for (int i = 0; i < nValKinds; ++i)
323  if (idIn == idVal[i]) return nVal[i];
324  return 0;}
325 
326  // Test whether a lepton is to be considered as unresolved.
327  bool isUnresolvedLepton();
328 
329  // Add extra remnant flavours to make valence and sea come out right.
330  bool remnantFlavours(Event& event, bool isDIS = false);
331 
332  // Correlate all initiators and remnants to make a colour singlet.
333  bool remnantColours(Event& event, vector<int>& colFrom,
334  vector<int>& colTo);
335 
336  // Pick unrescaled x of remnant parton (valence or sea).
337  double xRemnant(int i);
338 
339  // Tell whether a junction has been resolved, and its junction colours.
340  bool hasJunction() const {return hasJunctionBeam;}
341  int junctionCol(int i) const {return junCol[i];}
342  void junctionCol(int i, int col) {junCol[i] = col;}
343 
344  // For a diffractive system, decide whether to kick out gluon or quark.
345  bool pickGluon(double mDiff);
346 
347  // Pick a valence quark at random, and provide the remaining flavour.
348  int pickValence();
349  int pickRemnant() const {return idVal2;}
350 
351  // Share lightcone momentum between two remnants in a diffractive system.
352  // At the same time generate a relative pT for the two.
353  double zShare( double mDiff, double m1, double m2);
354  double pxShare() const {return pxRel;}
355  double pyShare() const {return pyRel;}
356 
357  // Add extra remnant flavours to make valence and sea come out right.
358  bool remnantFlavoursNew(Event& event);
359 
360  // Find the colour setup of the removed partons from the scatterings.
361  void findColSetup(Event& event);
362 
363  // Set initial colours.
364  void setInitialCol(Event & event);
365 
366  // Update colours.
367  void updateCol(vector<pair<int,int> > colourChanges);
368 
369  vector<pair <int,int> > getColUpdates() {return colUpdates;}
370 
371  // Set valence content for photon beams and position of first valence quark.
372  bool gammaInitiatorIsVal(int iResolved, int id, double x, double Q2);
373  bool gammaInitiatorIsVal(int iResolved, double Q2);
374  int getGammaValFlavour() { return abs(idVal[0]); }
375  int gammaValSeaComp(int iResolved);
376  void posVal(int iPosValIn) { iPosVal = iPosValIn; }
377  void gamVal(int iGamValIn) { iGamVal = iGamValIn; }
378  int gamVal() { return iGamVal; }
379 
380  // Set and get the state (resolved and/or unresolved) of photon beam.
381  void resolvedGamma(bool isResolved) { isResolvedGamma = isResolved; }
382  bool resolvedGamma() const { return isResolvedGamma; }
383  void setGammaMode(int gammaModeIn);
384  int getGammaMode() const { return gammaMode; }
385  bool isResolvedUnresolved() const { return isResUnres; }
386  void initGammaInBeam() { initGammaBeam = true; }
387  bool gammaInBeam() const { return initGammaBeam; }
388 
389  // Set state of VMD inside gamma.
390  void setVMDstate(bool isVMDIn, int idIn, double mIn, double scaleIn,
391  bool reassignState = false) {
392  hasVMDstateInBeam = isVMDIn;
393  idVMDBeam = idIn;
394  mVMDBeam = mIn;
395  scaleVMDBeam = scaleIn;
396  if (reassignState) {
397  idBeam = idVMDBeam;
398  mBeam = mVMDBeam;
399  pdfBeamPtr->setVMDscale(scaleVMDBeam);
400  }
401  }
402 
403  // Store the pT2 value of gamma->qqbar splitting.
404  void pT2gamma2qqbar(double pT2in) { pT2gm2qqbar = pT2in; }
405  double pT2gamma2qqbar() { return pT2gm2qqbar; }
406 
407  // Store the pT value for the latest MPI.
408  void pTMPI(double pTminMPIin) { pTminMPI = pTminMPIin; }
409 
410  // Check whether room for beam remnants.
411  bool roomFor1Remnant(double eCM);
412  bool roomFor1Remnant(int id1, double x1, double eCM);
413  bool roomFor2Remnants(int id1, double x1, double eCM);
414  bool roomForRemnants(BeamParticle beamOther);
415 
416  // Evaluate the remnant mass with initiator idIn.
417  double remnantMass(int idIn);
418 
419  // Functions to approximate pdfs for ISR.
420  double gammaPDFxDependence(int flavour, double x)
421  { return pdfBeamPtr->gammaPDFxDependence(flavour, x); }
422  double gammaPDFRefScale(int flavour)
423  { return pdfBeamPtr->gammaPDFRefScale(flavour); }
424  double xIntegratedPDFs(double Q2)
425  { return pdfBeamPtr->xfIntegratedTotal(Q2); }
426 
427  // Save the x_gamma value after latest PDF call or set it later if ND.
428  void xGammaPDF() { xGm = pdfHardBeamPtr->xGamma(); }
429  void xGamma(double xGmIn) { xGm = xGmIn; }
430  void Q2Gamma(double Q2GmIn) { Q2gm = Q2GmIn; }
431  void newGammaKTPhi(double kTIn, double phiIn)
432  { kTgamma = kTIn; phiGamma = phiIn; }
433 
434  // Get the kinematic limits for photons emitted by the beam.
435  double xGammaMin() { return pdfHardBeamPtr->getXmin(); }
436  double xGammaHadr() { return pdfHardBeamPtr->getXhadr(); }
437  double gammaFluxIntApprox() { return pdfHardBeamPtr->intFluxApprox(); }
438 
439  // Do photon flux use an approximation for sampling.
440  bool hasApproxGammaFlux() { return pdfHardBeamPtr->hasApproxGammaFlux(); }
441 
442  // Get the kinematics related photons form lepton beams.
443  double xGamma() const { return xGm; }
444  double Q2Gamma() const { return Q2gm; }
445  double gammaKTx() const { return kTgamma*cos(phiGamma); }
446  double gammaKTy() const { return kTgamma*sin(phiGamma); }
447  double gammaKT() const { return kTgamma; }
448  double gammaPhi() const { return phiGamma; }
449 
450  // Keep track of pomeron momentum fraction.
451  void xPom(double xpom = -1.0)
452  { if ( pdfBeamPtr ) pdfBeamPtr->xPom(xpom); }
453 
454  // Sample x and Q2 for emitted photons according to flux.
455  double sampleXgamma(double xMinIn)
456  { xGm = pdfHardBeamPtr->sampleXgamma(xMinIn); return xGm; }
457  double sampleQ2gamma(double Q2min)
458  { Q2gm = pdfHardBeamPtr->sampleQ2gamma(Q2min); return Q2gm;}
459 
460  // Prepare data on how much x has been used, for speedup of PDF evaluation.
461  xfModPrepData xfModPrep( int iSkip, double Q2);
462 
463 private:
464 
465  // Constants: could only be changed in the code itself.
466  static const double XMINUNRESOLVED, POMERONMASS, XMAXCOMPANION, TINYZREL;
467  static const int NMAX, NRANDOMTRIES;
468 
469  // Pointers to PDF sets.
470  PDFPtr pdfBeamPtr;
471  PDFPtr pdfHardBeamPtr;
472 
473  // Pointer to unresolved PDF and two others to save the resolved ptrs.
474  PDFPtr pdfUnresBeamPtr;
475  PDFPtr pdfBeamPtrSave;
476  PDFPtr pdfHardBeamPtrSave;
477 
478  // Array of PDFs to be used when idA can be changed between events.
479  vector<PDFPtr> pdfSavePtrs;
480  int pdfSaveIdx;
481 
482  // Pointer to class for flavour generation.
483  StringFlav* flavSelPtr;
484 
485  // Initialization data, normally only set once.
486  bool allowJunction, beamJunction;
487  int maxValQuark, companionPower;
488  double valencePowerMeson, valencePowerUinP, valencePowerDinP,
489  valenceDiqEnhance, pickQuarkNorm, pickQuarkPower,
490  diffPrimKTwidth, diffLargeMassSuppress, beamSat, gluonPower,
491  xGluonCutoff, heavyQuarkEnhance[6];
492 
493  // Basic properties of a beam particle.
494  int idBeam, idBeamAbs, idVMDBeam;
495  Vec4 pBeam;
496  double mBeam, mVMDBeam, scaleVMDBeam;
497 
498  // Beam kind. Valence flavour content for hadrons.
499  bool isUnresolvedBeam, isLeptonBeam, isHadronBeam, isMesonBeam,
500  isBaryonBeam, isGammaBeam;
501  int nValKinds, idVal[3], nVal[3];
502 
503  // Current parton density, by valence, sea and companion.
504  int idSave, iSkipSave, nValLeft[3];
505  double xqgTot, xqVal, xqgSea, xqCompSum;
506 
507  // Variables related to photon beams (also inside lepton).
508  bool doISR, doMPI, doND, isResolvedGamma, hasResGammaInBeam,
509  isResUnres, hasVMDstateInBeam, initGammaBeam;
510  double pTminISR, pTminMPI, pT2gm2qqbar;
511  int iGamVal, iPosVal, gammaMode;
512 
513  // Variables for photon from lepton.
514  double xGm, Q2gm, kTgamma, phiGamma;
515 
516  // Cache values for xCompFrac method.
517  int cPowerCache;
518  double xsCache, resCache;
519 
520  // The list of resolved partons.
521  vector<ResolvedParton> resolved;
522 
523  // Status after all initiators have been accounted for. Junction content.
524  int nInit;
525  bool hasJunctionBeam;
526  int junCol[3];
527 
528  // Variables for new colour reconnection;
529  pair <int,int> colSetup;
530  vector<int> acols, cols;
531  vector<bool> usedCol,usedAcol;
532  vector< pair<int,int> > colUpdates;
533  int nJuncs, nAjuncs, nDiffJuncs;
534  bool allowBeamJunctions;
535 
536  // Routine to calculate PDFs given previous interactions.
537  double xfModified( int iSkip, int idIn, double x, double Q2) {
538  xfModPrepData xfData = xfModPrep(iSkip, Q2);
539  return xfModified(iSkip, idIn, x, Q2, xfData);}
540  double xfModified( int iSkip, int idIn, double x, double Q2,
541  xfModPrepData& data);
542  // In case of resolved.size() === 0.
543  double xfModified0( int iSkip, int idIn, double x, double Q2);
544 
545  // Fraction of hadron momentum sitting in a valence quark distribution.
546  double xValFrac(int j, double Q2);
547  double Q2ValFracSav, uValInt, dValInt;
548 
549  // Fraction of hadron momentum sitting in a companion quark distribution.
550  double xCompFrac(double xs);
551 
552  // Value of companion quark PDF, also given the sea quark x.
553  double xCompDist(double xc, double xs);
554 
555  // Valence quark subdivision for diffractive systems.
556  int idVal1, idVal2, idVal3;
557  double zRel, pxRel, pyRel;
558 
559  // Update a single (anti-) colour of the event.
560  void updateSingleCol(int oldCol, int newCol);
561 
562  // Find a single (anti-) colour in the beam,
563  // if a beam remnant is set the new colour.
564  int findSingleCol(Event& event, bool isAcol, bool useHardScatters);
565 
566 };
567 
568 //==========================================================================
569 
570 } // end namespace Pythia8
571 
572 #endif // Pythia8_BeamParticle_H
ResolvedParton(int iPosIn=0, int idIn=0, double xIn=0., int companionIn=-1)
Constructor.
Definition: BeamParticle.h:43
double xfMax(int idIn, double x, double Q2)
Overestimate for PDFs. Same as normal except photons inside leptons.
Definition: BeamParticle.h:223
void setBeamID(int idIn, int iPDFin=-1)
Switch to new beam particle identity; for similar hadrons only.
Definition: BeamParticle.h:181
Definition: BeamParticle.h:120
void resolvedGamma(bool isResolved)
Set and get the state (resolved and/or unresolved) of photon beam.
Definition: BeamParticle.h:381
double xfFlux(int idIn, double x, double Q2)
Accurate and approximated photon flux and PDFs.
Definition: BeamParticle.h:227
Definition: PhysicsBase.h:27
The Event class holds all info on the generated event.
Definition: Event.h:453
void newM(double mIn)
Set new mass. Used with photons when virtuality is sampled.
Definition: BeamParticle.h:192
bool hasJunction() const
Tell whether a junction has been resolved, and its junction colours.
Definition: BeamParticle.h:340
Definition: BeamParticle.h:133
double xfHard(int idIn, double x, double Q2)
Special hard-process parton distributions (can agree with standard ones).
Definition: BeamParticle.h:219
void clear()
Clear list of resolved partons.
Definition: BeamParticle.h:298
BeamParticle()
Constructor.
Definition: BeamParticle.h:138
void iPos(int iPosIn)
Set info on initiator or remnant parton.
Definition: BeamParticle.h:49
void setVMDstate(bool isVMDIn, int idIn, double mIn, double scaleIn, bool reassignState=false)
Set state of VMD inside gamma.
Definition: BeamParticle.h:390
void newPzE(double pzIn, double eIn)
Set new pZ and E, but keep the rest the same.
Definition: BeamParticle.h:189
void initPDFPtr(PDFPtr pdfInPtr, PDFPtr pdfHardInPtr)
Initialize only the two pdf pointers.
Definition: BeamParticle.h:166
int nValenceKinds() const
How many different flavours, and how many quarks of given flavour.
Definition: BeamParticle.h:321
bool hasApproxGammaFlux()
Do photon flux use an approximation for sampling.
Definition: BeamParticle.h:440
void resetGammaInLepton()
Reset variables related to photon beam inside a lepton.
Definition: BeamParticle.h:305
double xGamma() const
Get the kinematics related photons form lepton beams.
Definition: BeamParticle.h:443
int size() const
Total number of partons extracted from beam, and initiators only.
Definition: BeamParticle.h:294
void xGammaPDF()
Save the x_gamma value after latest PDF call or set it later if ND.
Definition: BeamParticle.h:428
void pTMPI(double pTminMPIin)
Store the pT value for the latest MPI.
Definition: BeamParticle.h:408
Error envelope from PDF uncertainty.
Definition: PartonDistributions.h:102
double xfMPI(int idIn, double x, double Q2, xfModPrepData &xfData)
Definition: BeamParticle.h:251
bool isHadron() const
As hadrons here we only count those we know how to handle remnants for.
Definition: BeamParticle.h:208
void pT2gamma2qqbar(double pT2in)
Store the pT2 value of gamma->qqbar splitting.
Definition: BeamParticle.h:404
double xf(int idIn, double x, double Q2)
Standard parton distributions.
Definition: BeamParticle.h:240
void popBack()
Remove the last particle from the beam. Reset companion code if needed.
Definition: BeamParticle.h:313
void calcPDFEnvelope(int idNow, double xNow, double Q2Now, int valSea)
Calculate envelope of PDF predictions.
Definition: BeamParticle.h:276
ResolvedParton & operator[](int i)
Overload index operator to access a resolved parton from the list.
Definition: BeamParticle.h:290
double mQuarkPDF(int idIn)
Return quark masses used in the PDF fit (LHAPDF6 only).
Definition: BeamParticle.h:270
bool insideBounds(double x, double Q2)
Check whether x and Q2 values fall inside the fit bounds (LHAPDF6 only).
Definition: BeamParticle.h:263
double gammaPDFxDependence(int flavour, double x)
Functions to approximate pdfs for ISR.
Definition: BeamParticle.h:420
void initSwitchID(const vector< PDFPtr > &pdfSavePtrsIn)
Initialize array of PDFs for switching between them.
Definition: BeamParticle.h:173
double alphaS(double Q2)
Access the running alpha_s of a PDF set (LHAPDF6 only).
Definition: BeamParticle.h:267
int iPos() const
Get info on initiator or remnant parton.
Definition: BeamParticle.h:71
int id() const
Member functions for output.
Definition: BeamParticle.h:195
Definition: BeamParticle.h:38
double m2(const Vec4 &v1)
The squared invariant mass of one or more four-vectors.
Definition: Basics.cc:605
int nMembers()
Return number of members in PDF family (LHAPDF6 only).
Definition: BeamParticle.h:273
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
double m(const Vec4 &v1)
Invariant mass and its square.
Definition: Basics.cc:595
The StringFlav class is used to select quark and hadron flavours.
Definition: FragmentationFlavZpT.h:84
void initID(int idIn)
Initialize only the id.
Definition: BeamParticle.h:163
void resetGamma()
Reset variables related to photon beam.
Definition: BeamParticle.h:301
double xfSame(int idIn, double x, double Q2)
Definition: BeamParticle.h:236
double sampleXgamma(double xMinIn)
Sample x and Q2 for emitted photons according to flux.
Definition: BeamParticle.h:455
double xfVal(int idIn, double x, double Q2)
Ditto, split into valence and sea parts (where gluon counts as sea).
Definition: BeamParticle.h:244
double xGammaMin()
Get the kinematic limits for photons emitted by the beam.
Definition: BeamParticle.h:435
void xPom(double xpom=-1.0)
Keep track of pomeron momentum fraction.
Definition: BeamParticle.h:451
Definition: Basics.h:32
int append(int iPos, int idIn, double x, int companion=-1)
Add a resolved parton to list.
Definition: BeamParticle.h:308