1 | // MultipartonInteractions.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 | // 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 "Basics.h" |
---|
14 | #include "BeamParticle.h" |
---|
15 | #include "Event.h" |
---|
16 | #include "Info.h" |
---|
17 | #include "PartonSystems.h" |
---|
18 | #include "PythiaStdlib.h" |
---|
19 | #include "Settings.h" |
---|
20 | #include "SigmaTotal.h" |
---|
21 | #include "SigmaProcess.h" |
---|
22 | #include "StandardModel.h" |
---|
23 | #include "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 | SigmaMultiparton() {} |
---|
39 | |
---|
40 | // Destructor. |
---|
41 | ~SigmaMultiparton() { |
---|
42 | for (int i = 0; i < int(sigmaT.size()); ++i) delete sigmaT[i]; |
---|
43 | 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 | bool pickedOther() {return pickOther;} |
---|
57 | |
---|
58 | // Return one subprocess, picked according to relative cross sections. |
---|
59 | SigmaProcess* sigmaSel(); |
---|
60 | bool swapTU() {return pickedU;} |
---|
61 | |
---|
62 | // Return code or name of a specified process, for statistics table. |
---|
63 | int nProc() const {return nChan;} |
---|
64 | int codeProc(int iProc) const {return sigmaT[iProc]->code();} |
---|
65 | 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 | class MultipartonInteractions { |
---|
96 | |
---|
97 | public: |
---|
98 | |
---|
99 | // Constructor. |
---|
100 | 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 minbias process. |
---|
114 | void pTfirst(); |
---|
115 | |
---|
116 | // Set up kinematics for first = hardest pT in minbias 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 | 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 | void setEmpty() {pT2Ren = alpS = alpEM = x1 = x2 = pT2Fac |
---|
134 | = 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 | double bMPI() const {return (bIsSet) ? bNow / bAvg : 0.;} |
---|
146 | 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 | int getVSC1() const {return vsc1;} |
---|
151 | int getVSC2() const {return vsc2;} |
---|
152 | |
---|
153 | // Update and print statistics on number of processes. |
---|
154 | // Note: currently only valid for MB systems, not diffraction?? |
---|
155 | void accumulate() { int iBeg = (infoPtr->isMinBias()) ? 0 : 1; |
---|
156 | for (int i = iBeg; i < infoPtr->nMPI(); ++i) |
---|
157 | ++nGen[ infoPtr->codeMPI(i) ];} |
---|
158 | void statistics(bool resetStat = false, ostream& os = cout); |
---|
159 | void resetStatistics() { for ( map<int, int>::iterator iter = nGen.begin(); |
---|
160 | 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, bProfile, processLevel, |
---|
173 | rescatterMode, nQuarkIn, nSample, enhanceScreening; |
---|
174 | double alphaSvalue, Kfactor, pT0Ref, ecmRef, ecmPow, pTmin, coreRadius, |
---|
175 | coreFraction, expPow, ySepResc, deltaYResc, sigmaPomP, mPomP, pPomP, |
---|
176 | mMaxPertDiff, mMinPertDiff; |
---|
177 | |
---|
178 | // x-dependent matter profile: |
---|
179 | // Constants. |
---|
180 | static const int XDEP_BBIN; |
---|
181 | static const double XDEP_A0, XDEP_A1, XDEP_BSTEP, XDEP_BSTEPINC, |
---|
182 | XDEP_CUTOFF, XDEP_WARN, XDEP_SMB2FM; |
---|
183 | |
---|
184 | // Table of Int( dSigma/dX * O(b, X), dX ) in bins of b for fast |
---|
185 | // integration. Only needed during initialisation. |
---|
186 | vector <double> sigmaIntWgt, sigmaSumWgt; |
---|
187 | |
---|
188 | // a1: value of a1 constant, taken from settings database. |
---|
189 | // a0now (a02now): tuned value of a0 (squared value). |
---|
190 | // bstepNow: current size of bins in b. |
---|
191 | // a2max: given an xMin, a maximal (squared) value of the |
---|
192 | // width, to be used in overestimate Omax(b). |
---|
193 | // enhanceBmax, retain enhanceB as enhancement factor for the hardest |
---|
194 | // enhanceBnow: interaction. Use enhanceBmax as overestimate for fastPT2, |
---|
195 | // and enhanceBnow to store the value for the current |
---|
196 | // interaction. |
---|
197 | double a1, a0now, a02now, bstepNow, a2max, b2now, enhanceBmax, enhanceBnow; |
---|
198 | |
---|
199 | // Storage for trial interactions. |
---|
200 | int id1Save, id2Save; |
---|
201 | double pT2Save, x1Save, x2Save, sHatSave, tHatSave, uHatSave, |
---|
202 | alpSsave, alpEMsave, pT2FacSave, pT2RenSave, xPDF1nowSave, |
---|
203 | xPDF2nowSave; |
---|
204 | SigmaProcess *dSigmaDtSelSave; |
---|
205 | |
---|
206 | // vsc1, vsc2: for minimum-bias events with trial interaction, store |
---|
207 | // decision on whether hard interaction was valence or sea. |
---|
208 | int vsc1, vsc2; |
---|
209 | |
---|
210 | // Other initialization data. |
---|
211 | bool hasBaryonBeams, hasLowPow, globalRecoilFSR; |
---|
212 | int iDiffSys, nMaxGlobalRecoilFSR; |
---|
213 | double eCM, sCM, pT0, pT20, pT2min, pTmax, pT2max, pT20R, pT20minR, |
---|
214 | pT20maxR, pT20min0maxR, pT2maxmin, sigmaND, pT4dSigmaMax, |
---|
215 | pT4dProbMax, dSigmaApprox, sigmaInt, sudExpPT[101], |
---|
216 | zeroIntCorr, normOverlap, nAvg, kNow, normPi, bAvg, bDiv, |
---|
217 | probLowB, radius2B, radius2C, fracA, fracB, fracC, fracAhigh, |
---|
218 | fracBhigh, fracChigh, fracABChigh, expRev, cDiv, cMax; |
---|
219 | |
---|
220 | // Properties specific to current system. |
---|
221 | bool bIsSet, bSetInFirst, isAtLowB, pickOtherSel; |
---|
222 | int id1, id2, i1Sel, i2Sel, id1Sel, id2Sel; |
---|
223 | double bNow, enhanceB, pT2, pT2shift, pT2Ren, pT2Fac, x1, x2, xT, xT2, |
---|
224 | tau, y, sHat, tHat, uHat, alpS, alpEM, xPDF1now, xPDF2now, |
---|
225 | dSigmaSum, x1Sel, x2Sel, sHatSel, tHatSel, uHatSel; |
---|
226 | |
---|
227 | // Stored values for mass interpolation for diffractive systems. |
---|
228 | int nStep, iStepFrom, iStepTo; |
---|
229 | double eCMsave, eStepSize, eStepSave, eStepFrom, eStepTo, pT0Save[5], |
---|
230 | pT4dSigmaMaxSave[5], pT4dProbMaxSave[5], sigmaIntSave[5], |
---|
231 | sudExpPTSave[5][101], zeroIntCorrSave[5], normOverlapSave[5], |
---|
232 | kNowSave[5], bAvgSave[5], bDivSave[5], probLowBSave[5], |
---|
233 | fracAhighSave[5], fracBhighSave[5], fracChighSave[5], |
---|
234 | fracABChighSave[5], cDivSave[5], cMaxSave[5]; |
---|
235 | |
---|
236 | // Pointer to various information on the generation. |
---|
237 | Info* infoPtr; |
---|
238 | |
---|
239 | // Pointer to the random number generator. |
---|
240 | Rndm* rndmPtr; |
---|
241 | |
---|
242 | // Pointers to the two incoming beams. |
---|
243 | BeamParticle* beamAPtr; |
---|
244 | BeamParticle* beamBPtr; |
---|
245 | |
---|
246 | // Pointers to Standard Model couplings. |
---|
247 | Couplings* couplingsPtr; |
---|
248 | |
---|
249 | // Pointer to information on subcollision parton locations. |
---|
250 | PartonSystems* partonSystemsPtr; |
---|
251 | |
---|
252 | // Pointer to total cross section parametrization. |
---|
253 | SigmaTotal* sigmaTotPtr; |
---|
254 | |
---|
255 | // Pointer to user hooks. |
---|
256 | UserHooks* userHooksPtr; |
---|
257 | |
---|
258 | // Collections of parton-level 2 -> 2 cross sections. Selected one. |
---|
259 | SigmaMultiparton sigma2gg, sigma2qg, sigma2qqbarSame, sigma2qq; |
---|
260 | SigmaMultiparton* sigma2Sel; |
---|
261 | SigmaProcess* dSigmaDtSel; |
---|
262 | |
---|
263 | // Statistics on generated 2 -> 2 processes. |
---|
264 | map<int, int> nGen; |
---|
265 | |
---|
266 | // alphaStrong and alphaEM calculations. |
---|
267 | AlphaStrong alphaS; |
---|
268 | AlphaEM alphaEM; |
---|
269 | |
---|
270 | // Scattered partons. |
---|
271 | vector<int> scatteredA, scatteredB; |
---|
272 | |
---|
273 | // Determine constant in d(Prob)/d(pT2) < const / (pT2 + r * pT20)^2. |
---|
274 | void upperEnvelope(); |
---|
275 | |
---|
276 | // Integrate the parton-parton interaction cross section. |
---|
277 | void jetCrossSection(); |
---|
278 | |
---|
279 | // Evaluate "Sudakov form factor" for not having a harder interaction. |
---|
280 | double sudakov(double pT2sud, double enhance = 1.); |
---|
281 | |
---|
282 | // Do a quick evolution towards the next smaller pT. |
---|
283 | double fastPT2( double pT2beg); |
---|
284 | |
---|
285 | // Calculate the actual cross section, either for the first interaction |
---|
286 | // (including at initialization) or for any subsequent in the sequence. |
---|
287 | double sigmaPT2scatter(bool isFirst = false); |
---|
288 | |
---|
289 | // Find the partons that may rescatter. |
---|
290 | void findScatteredPartons( Event& event); |
---|
291 | |
---|
292 | // Calculate the actual cross section for a rescattering. |
---|
293 | double sigmaPT2rescatter( Event& event); |
---|
294 | |
---|
295 | // Calculate factor relating matter overlap and interaction rate. |
---|
296 | void overlapInit(); |
---|
297 | |
---|
298 | // Pick impact parameter and interaction rate enhancement, |
---|
299 | // either before the first interaction (for minbias) or after it. |
---|
300 | void overlapFirst(); |
---|
301 | void overlapNext(Event& event, double pTscale); |
---|
302 | |
---|
303 | }; |
---|
304 | |
---|
305 | //========================================================================== |
---|
306 | |
---|
307 | } // end namespace Pythia8 |
---|
308 | |
---|
309 | #endif // Pythia8_MultipartonInteractions_H |
---|