PYTHIA  8.312
MultipartonInteractions.h
1 // MultipartonInteractions.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 // This file contains the main classes for multiparton interactions physics.
7 // SigmaMultiparton stores allowed processes by in-flavour combination.
8 // MultipartonInteractions: generates multiparton parton-parton interactions.
9 
10 #ifndef Pythia8_MultipartonInteractions_H
11 #define Pythia8_MultipartonInteractions_H
12 
13 #include "Pythia8/Basics.h"
14 #include "Pythia8/BeamParticle.h"
15 #include "Pythia8/Event.h"
16 #include "Pythia8/Info.h"
17 #include "Pythia8/PartonSystems.h"
18 #include "Pythia8/PartonVertex.h"
19 #include "Pythia8/PhysicsBase.h"
20 #include "Pythia8/PythiaStdlib.h"
21 #include "Pythia8/Settings.h"
22 #include "Pythia8/SigmaTotal.h"
23 #include "Pythia8/SigmaProcess.h"
24 #include "Pythia8/StandardModel.h"
25 #include "Pythia8/UserHooks.h"
26 
27 namespace Pythia8 {
28 
29 //==========================================================================
30 
31 // SigmaMultiparton is a helper class to MultipartonInteractions.
32 // It packs pointers to the allowed processes for different
33 // flavour combinations and levels of ambition.
34 
36 
37 public:
38 
39  // Constructor.
40  SigmaMultiparton() : nChan(), needMasses(), useNarrowBW3(), useNarrowBW4(),
41  m3Fix(), m4Fix(), sHatMin(), sigmaT(), sigmaU(), sigmaTval(), sigmaUval(),
42  sigmaTsum(), sigmaUsum(), pickOther(), pickedU(), particleDataPtr(),
43  rndmPtr() {}
44 
45  // Initialize list of processes.
46  bool init(int inState, int processLevel, Info* infoPtr,
47  BeamParticle* beamAPtr, BeamParticle* beamBPtr);
48 
49  // Switch to new beam particle identities.
50  void updateBeamIDs() {
51  for (int i = 0; i < nChan; ++i) sigmaT[i]->updateBeamIDs();
52  for (int i = 0; i < nChan; ++i) sigmaU[i]->updateBeamIDs(); }
53 
54  // Calculate cross section summed over possibilities.
55  double sigma( int id1, int id2, double x1, double x2, double sHat,
56  double tHat, double uHat, double alpS, double alpEM,
57  bool restore = false, bool pickOtherIn = false);
58 
59  // Return whether the other, rare processes were selected.
60  bool pickedOther() {return pickOther;}
61 
62  // Return one subprocess, picked according to relative cross sections.
63  SigmaProcessPtr sigmaSel();
64  bool swapTU() {return pickedU;}
65 
66  // Return code or name of a specified process, for statistics table.
67  int nProc() const {return nChan;}
68  int codeProc(int iProc) const {return sigmaT[iProc]->code();}
69  string nameProc(int iProc) const {return sigmaT[iProc]->name();}
70 
71 private:
72 
73  // Constants: could only be changed in the code itself.
74  static const double MASSMARGIN, OTHERFRAC;
75 
76  // Number of processes. Some use massive matrix elements.
77  int nChan;
78  vector<bool> needMasses, useNarrowBW3, useNarrowBW4;
79  vector<double> m3Fix, m4Fix, sHatMin;
80 
81  // Vector of process list, one for t-channel and one for u-channel.
82  vector<SigmaProcessPtr> sigmaT, sigmaU;
83 
84  // Values of cross sections in process list above.
85  vector<double> sigmaTval, sigmaUval;
86  double sigmaTsum, sigmaUsum;
87  bool pickOther, pickedU;
88 
89  // Pointers to the particle data and random number generator.
90  ParticleData* particleDataPtr;
91  Rndm* rndmPtr;
92 
93 };
94 
95 //==========================================================================
96 
97 // The MultipartonInteractions class contains the main methods for the
98 // generation of multiparton parton-parton interactions in hadronic collisions.
99 
101 
102 public:
103 
104  // Constructor.
105  MultipartonInteractions() : allowRescatter(), allowDoubleRes(), canVetoMPI(),
106  doPartonVertex(), doVarEcm(), setAntiSame(), setAntiSameNow(),
107  pTmaxMatch(), alphaSorder(),
108  alphaEMorder(), alphaSnfmax(), bProfile(), processLevel(), bSelScale(),
109  rescatterMode(), nQuarkIn(), nSample(), enhanceScreening(), pT0paramMode(),
110  alphaSvalue(), Kfactor(), pT0Ref(), ecmRef(), ecmPow(), pTmin(),
111  coreRadius(), coreFraction(), expPow(), ySepResc(), deltaYResc(),
112  sigmaPomP(), mPomP(), pPomP(), mMaxPertDiff(), mMinPertDiff(), a1(),
113  a0now(), a02now(), bstepNow(), a2max(), b2now(), enhanceBmax(),
114  enhanceBnow(), id1Save(), id2Save(), pT2Save(), x1Save(), x2Save(),
115  sHatSave(), tHatSave(), uHatSave(), alpSsave(), alpEMsave(), pT2FacSave(),
116  pT2RenSave(), xPDF1nowSave(), xPDF2nowSave(), scaleLimitPTsave(),
117  dSigmaDtSelSave(), vsc1(), vsc2(), hasPomeronBeams(), hasLowPow(),
118  globalRecoilFSR(), iDiffSys(), nMaxGlobalRecoilFSR(), bSelHard(), eCM(),
119  sCM(), pT0(), pT20(), pT2min(), pTmax(), pT2max(), pT20R(), pT20minR(),
120  pT20maxR(), pT20min0maxR(), pT2maxmin(), sigmaND(), pT4dSigmaMax(),
121  pT4dProbMax(), dSigmaApprox(), sigmaInt(), sudExpPT(), zeroIntCorr(),
122  normOverlap(), nAvg(), kNow(), normPi(), bAvg(), bDiv(), probLowB(),
123  radius2B(), radius2C(), fracA(), fracB(), fracC(), fracAhigh(),
124  fracBhigh(), fracChigh(), fracABChigh(), expRev(), cDiv(), cMax(),
125  enhanceBavg(), bIsSet(false), bSetInFirst(), isAtLowB(), pickOtherSel(),
126  id1(), id2(), i1Sel(), i2Sel(), id1Sel(), id2Sel(), iPDFA(), nPDFA(1),
127  bNow(), enhanceB(), pT2(), pT2shift(), pT2Ren(), pT2Fac(), x1(),
128  x2(), xT(), xT2(), tau(), y(), sHat(), tHat(), uHat(), alpS(), alpEM(),
129  xPDF1now(), xPDF2now(), dSigmaSum(), x1Sel(), x2Sel(), sHatSel(),
130  tHatSel(), uHatSel(), iPDFAsave(), nStep(), iStepFrom(), iStepTo(),
131  eCMsave(), eStepMin(), eStepMax(), eStepSize(), eStepMix(), eStepFrom(),
132  eStepTo(), beamOffset(), mGmGmMin(), mGmGmMax(), hasGamma(),
133  isGammaGamma(), isGammaHadron(), isHadronGamma(), partonVertexPtr(),
134  sigma2Sel(), dSigmaDtSel() {}
135 
136  // Initialize the generation process for given beams.
137  bool init( bool doMPIinit, int iDiffSysIn,
138  BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
139  PartonVertexPtr partonVertexPtrIn, bool hasGammaIn = false);
140 
141  // Special setup to allow switching between beam PDFs for MPI handling.
142  void initSwitchID( const vector<int>& idAListIn) {
143  idAList = idAListIn; nPDFA = idAList.size();
144  mpis = vector<MPIInterpolationInfo>(nPDFA);}
145 
146  // Switch to new beam particle identities, and possibly PDFs.
147  void setBeamID(int iPDFAin) { iPDFA = iPDFAin;
148  sigma2gg.updateBeamIDs(); sigma2qg.updateBeamIDs();
149  sigma2qqbarSame.updateBeamIDs(); sigma2qq.updateBeamIDs();
150  setAntiSameNow = setAntiSame && particleDataPtr->hasAnti(infoPtr->idA())
151  && particleDataPtr->hasAnti(infoPtr->idB());}
152 
153  // Reset impact parameter choice and update the CM energy.
154  void reset();
155 
156  // Select first = hardest pT in nondiffractive process.
157  void pTfirst();
158 
159  // Set up kinematics for first = hardest pT in nondiffractive process.
160  void setupFirstSys( Event& process);
161 
162  // Find whether to limit maximum scale of emissions.
163  // Provide sum pT / 2 as potential limit where relevant.
164  bool limitPTmax( Event& event);
165  double scaleLimitPT() const {return scaleLimitPTsave;}
166 
167  // Prepare system for evolution.
168  void prepare(Event& event, double pTscale = 1000., bool rehashB = false) {
169  if (!bSetInFirst) overlapNext(event, pTscale, rehashB); }
170 
171  // Select next pT in downwards evolution.
172  double pTnext( double pTbegAll, double pTendAll, Event& event);
173 
174  // Set up kinematics of acceptable interaction.
175  bool scatter( Event& event);
176 
177  // Set "empty" values to avoid query of undefined quantities.
178  void setEmpty() {pT2Ren = alpS = alpEM = x1 = x2 = pT2Fac
179  = xPDF1now = xPDF2now = 0.; bIsSet = false;}
180 
181  // Get some information on current interaction.
182  double Q2Ren() const {return pT2Ren;}
183  double alphaSH() const {return alpS;}
184  double alphaEMH() const {return alpEM;}
185  double x1H() const {return x1;}
186  double x2H() const {return x2;}
187  double Q2Fac() const {return pT2Fac;}
188  double pdf1() const {return (id1 == 21 ? 4./9. : 1.) * xPDF1now;}
189  double pdf2() const {return (id2 == 21 ? 4./9. : 1.) * xPDF2now;}
190  double bMPI() const {return (bIsSet) ? bNow : 0.;}
191  double enhanceMPI() const {return (bIsSet) ? enhanceB / zeroIntCorr : 1.;}
192  double enhanceMPIavg() const {return enhanceBavg;}
193 
194  // For x-dependent matter profile, return incoming valence/sea
195  // decision from trial interactions.
196  int getVSC1() const {return vsc1;}
197  int getVSC2() const {return vsc2;}
198 
199  // Set the offset wrt. to normal beam particle positions for hard diffraction
200  // and for photon beam from lepton.
201  int getBeamOffset() const {return beamOffset;}
202  void setBeamOffset(int offsetIn) {beamOffset = offsetIn;}
203 
204  // Update and print statistics on number of processes.
205  // Note: currently only valid for nondiffractive systems, not diffraction??
206  void accumulate() { int iBeg = (infoPtr->isNonDiffractive()) ? 0 : 1;
207  for (int i = iBeg; i < infoPtr->nMPI(); ++i)
208  ++nGen[ infoPtr->codeMPI(i) ];}
209  void statistics(bool resetStat = false);
210  void resetStatistics() { for ( map<int, int>::iterator iter = nGen.begin();
211  iter != nGen.end(); ++iter) iter->second = 0; }
212 
213 private:
214 
215  // Constants: could only be changed in the code itself.
216  static const bool SHIFTFACSCALE, PREPICKRESCATTER;
217  static const double SIGMAFUDGE, RPT20, PT0STEP, SIGMASTEP, PT0MIN,
218  EXPPOWMIN, PROBATLOWB, BSTEP, BMAX, EXPMAX,
219  KCONVERGE, CONVERT2MB, ROOTMIN, ECMDEV, WTACCWARN,
220  SIGMAMBLIMIT;
221 
222  // Initialization data, read from Settings.
223  bool allowRescatter, allowDoubleRes, canVetoMPI, doPartonVertex, doVarEcm,
224  setAntiSame, setAntiSameNow, allowIDAswitch;
225  int pTmaxMatch, alphaSorder, alphaEMorder, alphaSnfmax, bProfile,
226  processLevel, bSelScale, rescatterMode, nQuarkIn, nSample,
227  enhanceScreening, pT0paramMode, reuseInit;
228  double alphaSvalue, Kfactor, pT0Ref, ecmRef, ecmPow, pTmin, coreRadius,
229  coreFraction, expPow, ySepResc, deltaYResc, sigmaPomP, mPomP, pPomP,
230  mMaxPertDiff, mMinPertDiff;
231  string initFile;
232 
233  // x-dependent matter profile:
234  // Constants.
235  static const int XDEP_BBIN;
236  static const double XDEP_A0, XDEP_A1, XDEP_BSTEP, XDEP_BSTEPINC,
237  XDEP_CUTOFF, XDEP_WARN, XDEP_SMB2FM;
238 
239  // Table of Int( dSigma/dX * O(b, X), dX ) in bins of b for fast
240  // integration. Only needed during initialisation.
241  vector <double> sigmaIntWgt, sigmaSumWgt;
242 
243  // a1: value of a1 constant, taken from settings database.
244  // a0now (a02now): tuned value of a0 (squared value).
245  // bstepNow: current size of bins in b.
246  // a2max: given an xMin, a maximal (squared) value of the
247  // width, to be used in overestimate Omax(b).
248  // enhanceBmax, retain enhanceB as enhancement factor for the hardest
249  // enhanceBnow: interaction. Use enhanceBmax as overestimate for fastPT2,
250  // and enhanceBnow to store the value for the current
251  // interaction.
252  double a1, a0now, a02now, bstepNow, a2max, b2now, enhanceBmax, enhanceBnow;
253 
254  // Storage for trial interactions.
255  int id1Save, id2Save;
256  double pT2Save, x1Save, x2Save, sHatSave, tHatSave, uHatSave,
257  alpSsave, alpEMsave, pT2FacSave, pT2RenSave, xPDF1nowSave,
258  xPDF2nowSave, scaleLimitPTsave;
259  SigmaProcessPtr dSigmaDtSelSave;
260 
261  // vsc1, vsc2: for minimum-bias events with trial interaction, store
262  // decision on whether hard interaction was valence or sea.
263  int vsc1, vsc2;
264 
265  // Other initialization data.
266  bool hasPomeronBeams, hasLowPow, globalRecoilFSR;
267  int iDiffSys, nMaxGlobalRecoilFSR, bSelHard;
268  double eCM, sCM, pT0, pT20, pT2min, pTmax, pT2max, pT20R, pT20minR,
269  pT20maxR, pT20min0maxR, pT2maxmin, sigmaND, pT4dSigmaMax,
270  pT4dProbMax, dSigmaApprox, sigmaInt, sudExpPT[101],
271  zeroIntCorr, normOverlap, nAvg, kNow, normPi, bAvg, bDiv,
272  probLowB, radius2B, radius2C, fracA, fracB, fracC, fracAhigh,
273  fracBhigh, fracChigh, fracABChigh, expRev, cDiv, cMax,
274  enhanceBavg;
275 
276  // Properties specific to current system.
277  bool bIsSet, bSetInFirst, isAtLowB, pickOtherSel;
278  int id1, id2, i1Sel, i2Sel, id1Sel, id2Sel, iPDFA, nPDFA;
279  vector<int> idAList;
280  double bNow, enhanceB, pT2, pT2shift, pT2Ren, pT2Fac, x1, x2, xT, xT2,
281  tau, y, sHat, tHat, uHat, alpS, alpEM, xPDF1now, xPDF2now,
282  dSigmaSum, x1Sel, x2Sel, sHatSel, tHatSel, uHatSel;
283 
284  // Local values for beam particle switch and mass interpolation.
285  int iPDFAsave, nStep, iStepFrom, iStepTo;
286  double eCMsave, eStepMin, eStepMax, eStepSize, eStepMix, eStepFrom, eStepTo;
287 
288  // Stored values for mass interpolation. First index projectile particle.
289  struct MPIInterpolationInfo {
290  int nStepSave;
291  double eStepMinSave, eStepMaxSave, eStepSizeSave;
292  vector<double> pT0Save, pT4dSigmaMaxSave, pT4dProbMaxSave,
293  sigmaIntSave, zeroIntCorrSave, normOverlapSave, kNowSave,
294  bAvgSave, bDivSave, probLowBSave,
295  fracAhighSave, fracBhighSave, fracChighSave,
296  fracABChighSave, cDivSave, cMaxSave;
297  vector<array<double, 101> > sudExpPTSave;
298 
299  void init(int nStepIn);
300  };
301 
302  vector<MPIInterpolationInfo> mpis;
303 
304  // Beam offset wrt. normal situation and other photon-related parameters.
305  int beamOffset;
306  double mGmGmMin, mGmGmMax;
307  bool hasGamma, isGammaGamma, isGammaHadron, isHadronGamma;
308 
309  // Pointer to assign space-time vertices during parton evolution.
310  PartonVertexPtr partonVertexPtr;
311 
312  // Collections of parton-level 2 -> 2 cross sections. Selected one.
313  SigmaMultiparton sigma2gg, sigma2qg, sigma2qqbarSame, sigma2qq;
314  SigmaMultiparton* sigma2Sel;
315  SigmaProcessPtr dSigmaDtSel;
316 
317  // Statistics on generated 2 -> 2 processes.
318  map<int, int> nGen;
319 
320  // alphaStrong and alphaEM calculations.
321  AlphaStrong alphaS;
322  AlphaEM alphaEM;
323 
324  // Scattered partons.
325  vector<int> scatteredA, scatteredB;
326 
327  // Determine constant in d(Prob)/d(pT2) < const / (pT2 + r * pT20)^2.
328  void upperEnvelope();
329 
330  // Integrate the parton-parton interaction cross section.
331  void jetCrossSection();
332 
333  // Read or write initialization data from/to file, to save startup time.
334  bool saveMPIdata();
335  bool loadMPIdata();
336 
337  // Evaluate "Sudakov form factor" for not having a harder interaction.
338  double sudakov(double pT2sud, double enhance = 1.);
339 
340  // Do a quick evolution towards the next smaller pT.
341  double fastPT2( double pT2beg);
342 
343  // Calculate the actual cross section, either for the first interaction
344  // (including at initialization) or for any subsequent in the sequence.
345  double sigmaPT2scatter(bool isFirst = false, bool doSymmetrize = false);
346 
347  // Find the partons that may rescatter.
348  void findScatteredPartons( Event& event);
349 
350  // Calculate the actual cross section for a rescattering.
351  double sigmaPT2rescatter( Event& event);
352 
353  // Calculate factor relating matter overlap and interaction rate.
354  void overlapInit();
355 
356  // Pick impact parameter and interaction rate enhancement,
357  // either before the first interaction (for nondiffractive) or after it.
358  void overlapFirst();
359  void overlapNext(Event& event, double pTscale, bool rehashB);
360 
361 };
362 
363 //==========================================================================
364 
365 } // end namespace Pythia8
366 
367 #endif // Pythia8_MultipartonInteractions_H
void prepare(Event &event, double pTscale=1000., bool rehashB=false)
Prepare system for evolution.
Definition: MultipartonInteractions.h:168
SigmaProcessPtr sigmaSel()
Return one subprocess, picked according to relative cross sections.
Definition: MultipartonInteractions.cc:243
Definition: StandardModel.h:23
Definition: PhysicsBase.h:27
double sigma(int id1, int id2, double x1, double x2, double sHat, double tHat, double uHat, double alpS, double alpEM, bool restore=false, bool pickOtherIn=false)
Calculate cross section summed over possibilities.
Definition: MultipartonInteractions.cc:182
Definition: Info.h:45
The Event class holds all info on the generated event.
Definition: Event.h:453
Definition: BeamParticle.h:133
bool init(int inState, int processLevel, Info *infoPtr, BeamParticle *beamAPtr, BeamParticle *beamBPtr)
Initialize list of processes.
Definition: MultipartonInteractions.cc:37
Definition: MultipartonInteractions.h:35
MultipartonInteractions()
Constructor.
Definition: MultipartonInteractions.h:105
void accumulate()
Definition: MultipartonInteractions.h:206
Definition: StandardModel.h:106
int nProc() const
Return code or name of a specified process, for statistics table.
Definition: MultipartonInteractions.h:67
void updateBeamIDs()
Switch to new beam particle identities.
Definition: MultipartonInteractions.h:50
Definition: Basics.h:386
void initSwitchID(const vector< int > &idAListIn)
Special setup to allow switching between beam PDFs for MPI handling.
Definition: MultipartonInteractions.h:142
bool pickedOther()
Return whether the other, rare processes were selected.
Definition: MultipartonInteractions.h:60
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
void setEmpty()
Set "empty" values to avoid query of undefined quantities.
Definition: MultipartonInteractions.h:178
Definition: MultipartonInteractions.h:100
int getBeamOffset() const
Definition: MultipartonInteractions.h:201
int getVSC1() const
Definition: MultipartonInteractions.h:196
double Q2Ren() const
Get some information on current interaction.
Definition: MultipartonInteractions.h:182
This class holds a map of all ParticleDataEntries.
Definition: ParticleData.h:422
void setBeamID(int iPDFAin)
Switch to new beam particle identities, and possibly PDFs.
Definition: MultipartonInteractions.h:147
SigmaMultiparton()
Constructor.
Definition: MultipartonInteractions.h:40