StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MultipartonInteractions.h
1 // MultipartonInteractions.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 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 
35 class SigmaMultiparton {
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  // Destructor.
46  ~SigmaMultiparton() {
47  for (int i = 0; i < int(sigmaT.size()); ++i) delete sigmaT[i];
48  for (int i = 0; i < int(sigmaU.size()); ++i) delete sigmaU[i];}
49 
50  // Initialize list of processes.
51  bool init(int inState, int processLevel, Info* infoPtr,
52  BeamParticle* beamAPtr, BeamParticle* beamBPtr);
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  SigmaProcess* 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<SigmaProcess*> 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 
100 class MultipartonInteractions : public PhysicsBase {
101 
102 public:
103 
104  // Constructor.
105  MultipartonInteractions() : allowRescatter(), allowDoubleRes(), canVetoMPI(),
106  doPartonVertex(), doVarEcm(), pTmaxMatch(), alphaSorder(), alphaEMorder(),
107  alphaSnfmax(), bProfile(), processLevel(), bSelScale(), rescatterMode(),
108  nQuarkIn(), nSample(), enhanceScreening(), pT0paramMode(), alphaSvalue(),
109  Kfactor(), pT0Ref(), ecmRef(), ecmPow(), pTmin(), coreRadius(),
110  coreFraction(), expPow(), ySepResc(), deltaYResc(), sigmaPomP(), mPomP(),
111  pPomP(), mMaxPertDiff(), mMinPertDiff(), a1(), a0now(), a02now(),
112  bstepNow(), a2max(), b2now(), enhanceBmax(), enhanceBnow(), id1Save(),
113  id2Save(), pT2Save(), x1Save(), x2Save(), sHatSave(), tHatSave(),
114  uHatSave(), alpSsave(), alpEMsave(), pT2FacSave(), pT2RenSave(),
115  xPDF1nowSave(), xPDF2nowSave(), scaleLimitPTsave(), dSigmaDtSelSave(),
116  vsc1(), vsc2(), hasBaryonBeams(), hasPomeronBeams(), hasLowPow(),
117  globalRecoilFSR(), iDiffSys(), nMaxGlobalRecoilFSR(), bSelHard(), eCM(),
118  sCM(), pT0(), pT20(), pT2min(), pTmax(), pT2max(), pT20R(), pT20minR(),
119  pT20maxR(), pT20min0maxR(), pT2maxmin(), sigmaND(), pT4dSigmaMax(),
120  pT4dProbMax(), dSigmaApprox(), sigmaInt(), sudExpPT(), zeroIntCorr(),
121  normOverlap(), nAvg(), kNow(), normPi(), bAvg(), bDiv(), probLowB(),
122  radius2B(), radius2C(), fracA(), fracB(), fracC(), fracAhigh(),
123  fracBhigh(), fracChigh(), fracABChigh(), expRev(), cDiv(), cMax(),
124  enhanceBavg(), bIsSet(false), bSetInFirst(), isAtLowB(), pickOtherSel(),
125  id1(), id2(), i1Sel(), i2Sel(), id1Sel(), id2Sel(), bNow(), enhanceB(),
126  pT2(), pT2shift(), pT2Ren(), pT2Fac(), x1(), x2(), xT(), xT2(), tau(),
127  y(), sHat(), tHat(), uHat(), alpS(), alpEM(), xPDF1now(), xPDF2now(),
128  dSigmaSum(), x1Sel(), x2Sel(), sHatSel(), tHatSel(), uHatSel(), nStep(),
129  iStepFrom(), iStepTo(), eCMsave(), eStepMin(), eStepMax(), eStepSize(),
130  eStepSave(), eStepFrom(), eStepTo(), pT0Save(), pT4dSigmaMaxSave(),
131  pT4dProbMaxSave(), sigmaIntSave(), sudExpPTSave(), zeroIntCorrSave(),
132  normOverlapSave(), kNowSave(), bAvgSave(), bDivSave(), probLowBSave(),
133  fracAhighSave(), fracBhighSave(), fracChighSave(), fracABChighSave(),
134  cDivSave(), cMaxSave(), beamOffset(), mGmGmMin(), mGmGmMax(), hasGamma(),
135  isGammaGamma(), isGammaHadron(), isHadronGamma(),
136  partonVertexPtr(), sigma2Sel(), dSigmaDtSel() {}
137 
138  // Initialize the generation process for given beams.
139  bool init( bool doMPIinit, int iDiffSysIn,
140  BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
141  PartonVertexPtr partonVertexPtrIn, bool hasGammaIn = false);
142 
143  // Reset impact parameter choice and update the CM energy.
144  void reset();
145 
146  // Select first = hardest pT in nondiffractive process.
147  void pTfirst();
148 
149  // Set up kinematics for first = hardest pT in nondiffractive process.
150  void setupFirstSys( Event& process);
151 
152  // Find whether to limit maximum scale of emissions.
153  // Provide sum pT / 2 as potential limit where relevant.
154  bool limitPTmax( Event& event);
155  double scaleLimitPT() const {return scaleLimitPTsave;}
156 
157  // Prepare system for evolution.
158  void prepare(Event& event, double pTscale = 1000., bool rehashB = false) {
159  if (!bSetInFirst) overlapNext(event, pTscale, rehashB); }
160 
161  // Select next pT in downwards evolution.
162  double pTnext( double pTbegAll, double pTendAll, Event& event);
163 
164  // Set up kinematics of acceptable interaction.
165  bool scatter( Event& event);
166 
167  // Set "empty" values to avoid query of undefined quantities.
168  void setEmpty() {pT2Ren = alpS = alpEM = x1 = x2 = pT2Fac
169  = xPDF1now = xPDF2now = 0.; bIsSet = false;}
170 
171  // Get some information on current interaction.
172  double Q2Ren() const {return pT2Ren;}
173  double alphaSH() const {return alpS;}
174  double alphaEMH() const {return alpEM;}
175  double x1H() const {return x1;}
176  double x2H() const {return x2;}
177  double Q2Fac() const {return pT2Fac;}
178  double pdf1() const {return (id1 == 21 ? 4./9. : 1.) * xPDF1now;}
179  double pdf2() const {return (id2 == 21 ? 4./9. : 1.) * xPDF2now;}
180  double bMPI() const {return (bIsSet) ? bNow : 0.;}
181  double enhanceMPI() const {return (bIsSet) ? enhanceB / zeroIntCorr : 1.;}
182  double enhanceMPIavg() const {return enhanceBavg;}
183 
184  // For x-dependent matter profile, return incoming valence/sea
185  // decision from trial interactions.
186  int getVSC1() const {return vsc1;}
187  int getVSC2() const {return vsc2;}
188 
189  // Set the offset wrt. to normal beam particle positions for hard diffraction
190  // and for photon beam from lepton.
191  int getBeamOffset() const {return beamOffset;}
192  void setBeamOffset(int offsetIn) {beamOffset = offsetIn;}
193 
194  // Update and print statistics on number of processes.
195  // Note: currently only valid for nondiffractive systems, not diffraction??
196  void accumulate() { int iBeg = (infoPtr->isNonDiffractive()) ? 0 : 1;
197  for (int i = iBeg; i < infoPtr->nMPI(); ++i)
198  ++nGen[ infoPtr->codeMPI(i) ];}
199  void statistics(bool resetStat = false);
200  void resetStatistics() { for ( map<int, int>::iterator iter = nGen.begin();
201  iter != nGen.end(); ++iter) iter->second = 0; }
202 
203 private:
204 
205  // Constants: could only be changed in the code itself.
206  static const bool SHIFTFACSCALE, PREPICKRESCATTER;
207  static const double SIGMAFUDGE, RPT20, PT0STEP, SIGMASTEP, PT0MIN,
208  EXPPOWMIN, PROBATLOWB, BSTEP, BMAX, EXPMAX,
209  KCONVERGE, CONVERT2MB, ROOTMIN, ECMDEV, WTACCWARN,
210  SIGMAMBLIMIT;
211 
212  // Initialization data, read from Settings.
213  bool allowRescatter, allowDoubleRes, canVetoMPI, doPartonVertex, doVarEcm;
214  int pTmaxMatch, alphaSorder, alphaEMorder, alphaSnfmax, bProfile,
215  processLevel, bSelScale, rescatterMode, nQuarkIn, nSample,
216  enhanceScreening, pT0paramMode;
217  double alphaSvalue, Kfactor, pT0Ref, ecmRef, ecmPow, pTmin, coreRadius,
218  coreFraction, expPow, ySepResc, deltaYResc, sigmaPomP, mPomP, pPomP,
219  mMaxPertDiff, mMinPertDiff;
220 
221  // x-dependent matter profile:
222  // Constants.
223  static const int XDEP_BBIN;
224  static const double XDEP_A0, XDEP_A1, XDEP_BSTEP, XDEP_BSTEPINC,
225  XDEP_CUTOFF, XDEP_WARN, XDEP_SMB2FM;
226 
227  // Table of Int( dSigma/dX * O(b, X), dX ) in bins of b for fast
228  // integration. Only needed during initialisation.
229  vector <double> sigmaIntWgt, sigmaSumWgt;
230 
231  // a1: value of a1 constant, taken from settings database.
232  // a0now (a02now): tuned value of a0 (squared value).
233  // bstepNow: current size of bins in b.
234  // a2max: given an xMin, a maximal (squared) value of the
235  // width, to be used in overestimate Omax(b).
236  // enhanceBmax, retain enhanceB as enhancement factor for the hardest
237  // enhanceBnow: interaction. Use enhanceBmax as overestimate for fastPT2,
238  // and enhanceBnow to store the value for the current
239  // interaction.
240  double a1, a0now, a02now, bstepNow, a2max, b2now, enhanceBmax, enhanceBnow;
241 
242  // Storage for trial interactions.
243  int id1Save, id2Save;
244  double pT2Save, x1Save, x2Save, sHatSave, tHatSave, uHatSave,
245  alpSsave, alpEMsave, pT2FacSave, pT2RenSave, xPDF1nowSave,
246  xPDF2nowSave, scaleLimitPTsave;
247  SigmaProcess *dSigmaDtSelSave;
248 
249  // vsc1, vsc2: for minimum-bias events with trial interaction, store
250  // decision on whether hard interaction was valence or sea.
251  int vsc1, vsc2;
252 
253  // Other initialization data.
254  bool hasBaryonBeams, hasPomeronBeams, hasLowPow, globalRecoilFSR;
255  int iDiffSys, nMaxGlobalRecoilFSR, bSelHard;
256  double eCM, sCM, pT0, pT20, pT2min, pTmax, pT2max, pT20R, pT20minR,
257  pT20maxR, pT20min0maxR, pT2maxmin, sigmaND, pT4dSigmaMax,
258  pT4dProbMax, dSigmaApprox, sigmaInt, sudExpPT[101],
259  zeroIntCorr, normOverlap, nAvg, kNow, normPi, bAvg, bDiv,
260  probLowB, radius2B, radius2C, fracA, fracB, fracC, fracAhigh,
261  fracBhigh, fracChigh, fracABChigh, expRev, cDiv, cMax,
262  enhanceBavg;
263 
264  // Properties specific to current system.
265  bool bIsSet, bSetInFirst, isAtLowB, pickOtherSel;
266  int id1, id2, i1Sel, i2Sel, id1Sel, id2Sel;
267  double bNow, enhanceB, pT2, pT2shift, pT2Ren, pT2Fac, x1, x2, xT, xT2,
268  tau, y, sHat, tHat, uHat, alpS, alpEM, xPDF1now, xPDF2now,
269  dSigmaSum, x1Sel, x2Sel, sHatSel, tHatSel, uHatSel;
270 
271  // Stored values for mass interpolation for diffractive systems.
272  int nStep, iStepFrom, iStepTo;
273  double eCMsave, eStepMin, eStepMax, eStepSize, eStepSave, eStepFrom, eStepTo,
274  pT0Save[20], pT4dSigmaMaxSave[20], pT4dProbMaxSave[20],
275  sigmaIntSave[20], sudExpPTSave[20][101], zeroIntCorrSave[20],
276  normOverlapSave[20], kNowSave[20], bAvgSave[20], bDivSave[20],
277  probLowBSave[20], fracAhighSave[20], fracBhighSave[20],
278  fracChighSave[20], fracABChighSave[20], cDivSave[20], cMaxSave[20];
279 
280  // Beam offset wrt. normal situation and other photon-related parameters.
281  int beamOffset;
282  double mGmGmMin, mGmGmMax;
283  bool hasGamma, isGammaGamma, isGammaHadron, isHadronGamma;
284 
285  // Pointer to assign space-time vertices during parton evolution.
286  PartonVertexPtr partonVertexPtr;
287 
288  // Collections of parton-level 2 -> 2 cross sections. Selected one.
289  SigmaMultiparton sigma2gg, sigma2qg, sigma2qqbarSame, sigma2qq;
290  SigmaMultiparton* sigma2Sel;
291  SigmaProcess* dSigmaDtSel;
292 
293  // Statistics on generated 2 -> 2 processes.
294  map<int, int> nGen;
295 
296  // alphaStrong and alphaEM calculations.
297  AlphaStrong alphaS;
298  AlphaEM alphaEM;
299 
300  // Scattered partons.
301  vector<int> scatteredA, scatteredB;
302 
303  // Determine constant in d(Prob)/d(pT2) < const / (pT2 + r * pT20)^2.
304  void upperEnvelope();
305 
306  // Integrate the parton-parton interaction cross section.
307  void jetCrossSection();
308 
309  // Evaluate "Sudakov form factor" for not having a harder interaction.
310  double sudakov(double pT2sud, double enhance = 1.);
311 
312  // Do a quick evolution towards the next smaller pT.
313  double fastPT2( double pT2beg);
314 
315  // Calculate the actual cross section, either for the first interaction
316  // (including at initialization) or for any subsequent in the sequence.
317  double sigmaPT2scatter(bool isFirst = false);
318 
319  // Find the partons that may rescatter.
320  void findScatteredPartons( Event& event);
321 
322  // Calculate the actual cross section for a rescattering.
323  double sigmaPT2rescatter( Event& event);
324 
325  // Calculate factor relating matter overlap and interaction rate.
326  void overlapInit();
327 
328  // Pick impact parameter and interaction rate enhancement,
329  // either before the first interaction (for nondiffractive) or after it.
330  void overlapFirst();
331  void overlapNext(Event& event, double pTscale, bool rehashB);
332 
333 };
334 
335 //==========================================================================
336 
337 } // end namespace Pythia8
338 
339 #endif // Pythia8_MultipartonInteractions_H
Definition: AgUStep.h:26