Line data Source code
1 : // MultipartonInteractions.h is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 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 : // 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/PythiaStdlib.h"
19 : #include "Pythia8/Settings.h"
20 : #include "Pythia8/SigmaTotal.h"
21 : #include "Pythia8/SigmaProcess.h"
22 : #include "Pythia8/StandardModel.h"
23 : #include "Pythia8/UserHooks.h"
24 :
25 : namespace Pythia8 {
26 :
27 : //==========================================================================
28 :
29 : // SigmaMultiparton is a helper class to MultipartonInteractions.
30 : // It packs pointers to the allowed processes for different
31 : // flavour combinations and levels of ambition.
32 :
33 : class SigmaMultiparton {
34 :
35 : public:
36 :
37 : // Constructor.
38 0 : SigmaMultiparton() {}
39 :
40 : // Destructor.
41 0 : ~SigmaMultiparton() {
42 0 : for (int i = 0; i < int(sigmaT.size()); ++i) delete sigmaT[i];
43 0 : for (int i = 0; i < int(sigmaU.size()); ++i) delete sigmaU[i];}
44 :
45 : // Initialize list of processes.
46 : bool init(int inState, int processLevel, Info* infoPtr,
47 : Settings* settingsPtr, ParticleData* particleDataPtr, Rndm* rndmPtrIn,
48 : BeamParticle* beamAPtr, BeamParticle* beamBPtr, Couplings* couplingsPtr);
49 :
50 : // Calculate cross section summed over possibilities.
51 : double sigma( int id1, int id2, double x1, double x2, double sHat,
52 : double tHat, double uHat, double alpS, double alpEM,
53 : bool restore = false, bool pickOtherIn = false);
54 :
55 : // Return whether the other, rare processes were selected.
56 0 : bool pickedOther() {return pickOther;}
57 :
58 : // Return one subprocess, picked according to relative cross sections.
59 : SigmaProcess* sigmaSel();
60 0 : bool swapTU() {return pickedU;}
61 :
62 : // Return code or name of a specified process, for statistics table.
63 0 : int nProc() const {return nChan;}
64 0 : int codeProc(int iProc) const {return sigmaT[iProc]->code();}
65 0 : string nameProc(int iProc) const {return sigmaT[iProc]->name();}
66 :
67 : private:
68 :
69 : // Constants: could only be changed in the code itself.
70 : static const double MASSMARGIN, OTHERFRAC;
71 :
72 : // Number of processes. Some use massive matrix elements.
73 : int nChan;
74 : vector<bool> needMasses;
75 : vector<double> m3Fix, m4Fix, sHatMin;
76 :
77 : // Vector of process list, one for t-channel and one for u-channel.
78 : vector<SigmaProcess*> sigmaT, sigmaU;
79 :
80 : // Values of cross sections in process list above.
81 : vector<double> sigmaTval, sigmaUval;
82 : double sigmaTsum, sigmaUsum;
83 : bool pickOther, pickedU;
84 :
85 : // Pointer to the random number generator.
86 : Rndm* rndmPtr;
87 :
88 : };
89 :
90 : //==========================================================================
91 :
92 : // The MultipartonInteractions class contains the main methods for the
93 : // generation of multiparton parton-parton interactions in hadronic collisions.
94 :
95 0 : class MultipartonInteractions {
96 :
97 : public:
98 :
99 : // Constructor.
100 0 : MultipartonInteractions() : bIsSet(false) {}
101 :
102 : // Initialize the generation process for given beams.
103 : bool init( bool doMPIinit, int iDiffSysIn, Info* infoPtrIn,
104 : Settings& settings, ParticleData* particleDataPtr, Rndm* rndmPtrIn,
105 : BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
106 : Couplings* couplingsPtrIn, PartonSystems* partonSystemsPtrIn,
107 : SigmaTotal* sigmaTotPtrIn, UserHooks* userHooksPtrIn,
108 : ostream& os = cout);
109 :
110 : // Reset impact parameter choice and update the CM energy.
111 : void reset();
112 :
113 : // Select first = hardest pT in nondiffractive process.
114 : void pTfirst();
115 :
116 : // Set up kinematics for first = hardest pT in nondiffractive process.
117 : void setupFirstSys( Event& process);
118 :
119 : // Find whether to limit maximum scale of emissions.
120 : bool limitPTmax( Event& event);
121 :
122 : // Prepare system for evolution.
123 : void prepare(Event& event, double pTscale = 1000.) {
124 0 : if (!bSetInFirst) overlapNext(event, pTscale); }
125 :
126 : // Select next pT in downwards evolution.
127 : double pTnext( double pTbegAll, double pTendAll, Event& event);
128 :
129 : // Set up kinematics of acceptable interaction.
130 : bool scatter( Event& event);
131 :
132 : // Set "empty" values to avoid query of undefined quantities.
133 0 : void setEmpty() {pT2Ren = alpS = alpEM = x1 = x2 = pT2Fac
134 0 : = xPDF1now = xPDF2now = 0.; bIsSet = false;}
135 :
136 : // Get some information on current interaction.
137 : double Q2Ren() const {return pT2Ren;}
138 : double alphaSH() const {return alpS;}
139 : double alphaEMH() const {return alpEM;}
140 : double x1H() const {return x1;}
141 : double x2H() const {return x2;}
142 : double Q2Fac() const {return pT2Fac;}
143 : double pdf1() const {return xPDF1now;}
144 : double pdf2() const {return xPDF2now;}
145 0 : double bMPI() const {return (bIsSet) ? bNow / bAvg : 0.;}
146 0 : double enhanceMPI() const {return (bIsSet) ? enhanceB / zeroIntCorr : 1.;}
147 :
148 : // For x-dependent matter profile, return incoming valence/sea
149 : // decision from trial interactions.
150 0 : int getVSC1() const {return vsc1;}
151 0 : int getVSC2() const {return vsc2;}
152 :
153 : // Update and print statistics on number of processes.
154 : // Note: currently only valid for nondiffractive systems, not diffraction??
155 0 : void accumulate() { int iBeg = (infoPtr->isNonDiffractive()) ? 0 : 1;
156 0 : for (int i = iBeg; i < infoPtr->nMPI(); ++i)
157 0 : ++nGen[ infoPtr->codeMPI(i) ];}
158 : void statistics(bool resetStat = false, ostream& os = cout);
159 0 : void resetStatistics() { for ( map<int, int>::iterator iter = nGen.begin();
160 0 : iter != nGen.end(); ++iter) iter->second = 0; }
161 :
162 : private:
163 :
164 : // Constants: could only be changed in the code itself.
165 : static const bool SHIFTFACSCALE, PREPICKRESCATTER;
166 : static const double SIGMAFUDGE, RPT20, PT0STEP, SIGMASTEP, PT0MIN,
167 : EXPPOWMIN, PROBATLOWB, BSTEP, BMAX, EXPMAX,
168 : KCONVERGE, CONVERT2MB, ROOTMIN, ECMDEV, WTACCWARN;
169 :
170 : // Initialization data, read from Settings.
171 : bool allowRescatter, allowDoubleRes, canVetoMPI;
172 : int pTmaxMatch, alphaSorder, alphaEMorder, alphaSnfmax, bProfile,
173 : processLevel, bSelScale, rescatterMode, nQuarkIn, nSample,
174 : enhanceScreening;
175 : double alphaSvalue, Kfactor, pT0Ref, ecmRef, ecmPow, pTmin, coreRadius,
176 : coreFraction, expPow, ySepResc, deltaYResc, sigmaPomP, mPomP, pPomP,
177 : mMaxPertDiff, mMinPertDiff;
178 :
179 : // x-dependent matter profile:
180 : // Constants.
181 : static const int XDEP_BBIN;
182 : static const double XDEP_A0, XDEP_A1, XDEP_BSTEP, XDEP_BSTEPINC,
183 : XDEP_CUTOFF, XDEP_WARN, XDEP_SMB2FM;
184 :
185 : // Table of Int( dSigma/dX * O(b, X), dX ) in bins of b for fast
186 : // integration. Only needed during initialisation.
187 : vector <double> sigmaIntWgt, sigmaSumWgt;
188 :
189 : // a1: value of a1 constant, taken from settings database.
190 : // a0now (a02now): tuned value of a0 (squared value).
191 : // bstepNow: current size of bins in b.
192 : // a2max: given an xMin, a maximal (squared) value of the
193 : // width, to be used in overestimate Omax(b).
194 : // enhanceBmax, retain enhanceB as enhancement factor for the hardest
195 : // enhanceBnow: interaction. Use enhanceBmax as overestimate for fastPT2,
196 : // and enhanceBnow to store the value for the current
197 : // interaction.
198 : double a1, a0now, a02now, bstepNow, a2max, b2now, enhanceBmax, enhanceBnow;
199 :
200 : // Storage for trial interactions.
201 : int id1Save, id2Save;
202 : double pT2Save, x1Save, x2Save, sHatSave, tHatSave, uHatSave,
203 : alpSsave, alpEMsave, pT2FacSave, pT2RenSave, xPDF1nowSave,
204 : xPDF2nowSave;
205 : SigmaProcess *dSigmaDtSelSave;
206 :
207 : // vsc1, vsc2: for minimum-bias events with trial interaction, store
208 : // decision on whether hard interaction was valence or sea.
209 : int vsc1, vsc2;
210 :
211 : // Other initialization data.
212 : bool hasBaryonBeams, hasLowPow, globalRecoilFSR;
213 : int iDiffSys, nMaxGlobalRecoilFSR;
214 : double eCM, sCM, pT0, pT20, pT2min, pTmax, pT2max, pT20R, pT20minR,
215 : pT20maxR, pT20min0maxR, pT2maxmin, sigmaND, pT4dSigmaMax,
216 : pT4dProbMax, dSigmaApprox, sigmaInt, sudExpPT[101],
217 : zeroIntCorr, normOverlap, nAvg, kNow, normPi, bAvg, bDiv,
218 : probLowB, radius2B, radius2C, fracA, fracB, fracC, fracAhigh,
219 : fracBhigh, fracChigh, fracABChigh, expRev, cDiv, cMax;
220 :
221 : // Properties specific to current system.
222 : bool bIsSet, bSetInFirst, isAtLowB, pickOtherSel;
223 : int id1, id2, i1Sel, i2Sel, id1Sel, id2Sel;
224 : double bNow, enhanceB, pT2, pT2shift, pT2Ren, pT2Fac, x1, x2, xT, xT2,
225 : tau, y, sHat, tHat, uHat, alpS, alpEM, xPDF1now, xPDF2now,
226 : dSigmaSum, x1Sel, x2Sel, sHatSel, tHatSel, uHatSel;
227 :
228 : // Stored values for mass interpolation for diffractive systems.
229 : int nStep, iStepFrom, iStepTo;
230 : double eCMsave, eStepSize, eStepSave, eStepFrom, eStepTo, pT0Save[5],
231 : pT4dSigmaMaxSave[5], pT4dProbMaxSave[5], sigmaIntSave[5],
232 : sudExpPTSave[5][101], zeroIntCorrSave[5], normOverlapSave[5],
233 : kNowSave[5], bAvgSave[5], bDivSave[5], probLowBSave[5],
234 : fracAhighSave[5], fracBhighSave[5], fracChighSave[5],
235 : fracABChighSave[5], cDivSave[5], cMaxSave[5];
236 :
237 : // Pointer to various information on the generation.
238 : Info* infoPtr;
239 :
240 : // Pointer to the random number generator.
241 : Rndm* rndmPtr;
242 :
243 : // Pointers to the two incoming beams.
244 : BeamParticle* beamAPtr;
245 : BeamParticle* beamBPtr;
246 :
247 : // Pointers to Standard Model couplings.
248 : Couplings* couplingsPtr;
249 :
250 : // Pointer to information on subcollision parton locations.
251 : PartonSystems* partonSystemsPtr;
252 :
253 : // Pointer to total cross section parametrization.
254 : SigmaTotal* sigmaTotPtr;
255 :
256 : // Pointer to user hooks.
257 : UserHooks* userHooksPtr;
258 :
259 : // Collections of parton-level 2 -> 2 cross sections. Selected one.
260 : SigmaMultiparton sigma2gg, sigma2qg, sigma2qqbarSame, sigma2qq;
261 : SigmaMultiparton* sigma2Sel;
262 : SigmaProcess* dSigmaDtSel;
263 :
264 : // Statistics on generated 2 -> 2 processes.
265 : map<int, int> nGen;
266 :
267 : // alphaStrong and alphaEM calculations.
268 : AlphaStrong alphaS;
269 : AlphaEM alphaEM;
270 :
271 : // Scattered partons.
272 : vector<int> scatteredA, scatteredB;
273 :
274 : // Determine constant in d(Prob)/d(pT2) < const / (pT2 + r * pT20)^2.
275 : void upperEnvelope();
276 :
277 : // Integrate the parton-parton interaction cross section.
278 : void jetCrossSection();
279 :
280 : // Evaluate "Sudakov form factor" for not having a harder interaction.
281 : double sudakov(double pT2sud, double enhance = 1.);
282 :
283 : // Do a quick evolution towards the next smaller pT.
284 : double fastPT2( double pT2beg);
285 :
286 : // Calculate the actual cross section, either for the first interaction
287 : // (including at initialization) or for any subsequent in the sequence.
288 : double sigmaPT2scatter(bool isFirst = false);
289 :
290 : // Find the partons that may rescatter.
291 : void findScatteredPartons( Event& event);
292 :
293 : // Calculate the actual cross section for a rescattering.
294 : double sigmaPT2rescatter( Event& event);
295 :
296 : // Calculate factor relating matter overlap and interaction rate.
297 : void overlapInit();
298 :
299 : // Pick impact parameter and interaction rate enhancement,
300 : // either before the first interaction (for nondiffractive) or after it.
301 : void overlapFirst();
302 : void overlapNext(Event& event, double pTscale);
303 :
304 : };
305 :
306 : //==========================================================================
307 :
308 : } // end namespace Pythia8
309 :
310 : #endif // Pythia8_MultipartonInteractions_H
|
