]>
Commit | Line | Data |
---|---|---|
c6b60c38 | 1 | // MultipartonInteractions.cc is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2013 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 | // Function definitions (not found in the header) for the | |
7 | // SigmaMultiparton and MultipartonInteractions classes. | |
8 | ||
9 | #include "MultipartonInteractions.h" | |
10 | ||
11 | // Internal headers for special processes. | |
12 | #include "SigmaQCD.h" | |
13 | #include "SigmaEW.h" | |
14 | #include "SigmaOnia.h" | |
15 | ||
16 | namespace Pythia8 { | |
17 | ||
18 | //========================================================================== | |
19 | ||
20 | // The SigmaMultiparton class. | |
21 | ||
22 | //-------------------------------------------------------------------------- | |
23 | ||
24 | // Constants: could be changed here if desired, but normally should not. | |
25 | // These are of technical nature, as described for each. | |
26 | ||
27 | // The sum of outgoing masses must not be too close to the cm energy. | |
28 | const double SigmaMultiparton::MASSMARGIN = 0.1; | |
29 | ||
30 | // Fraction of time not the dominant "gluon t-channel" process is picked. | |
31 | const double SigmaMultiparton::OTHERFRAC = 0.2; | |
32 | ||
33 | //-------------------------------------------------------------------------- | |
34 | ||
35 | // Initialize the generation process for given beams. | |
36 | ||
37 | bool SigmaMultiparton::init(int inState, int processLevel, Info* infoPtr, | |
38 | Settings* settingsPtr, ParticleData* particleDataPtr, Rndm* rndmPtrIn, | |
39 | BeamParticle* beamAPtr, BeamParticle* beamBPtr, Couplings* couplingsPtr) { | |
40 | ||
41 | // Store input pointer for future use. | |
42 | rndmPtr = rndmPtrIn; | |
43 | ||
44 | // Reset vector sizes (necessary in case of re-initialization). | |
45 | if (sigmaT.size() > 0) { | |
46 | for (int i = 0; i < int(sigmaT.size()); ++i) delete sigmaT[i]; | |
47 | sigmaT.resize(0); | |
48 | } | |
49 | if (sigmaU.size() > 0) { | |
50 | for (int i = 0; i < int(sigmaU.size()); ++i) delete sigmaU[i]; | |
51 | sigmaU.resize(0); | |
52 | } | |
53 | ||
54 | // Always store mimimal set of processes:QCD 2 -> 2 t-channel. | |
55 | ||
56 | // Gluon-gluon instate. | |
57 | if (inState == 0) { | |
58 | sigmaT.push_back( new Sigma2gg2gg() ); | |
59 | sigmaU.push_back( new Sigma2gg2gg() ); | |
60 | ||
61 | // Quark-gluon instate. | |
62 | } else if (inState == 1) { | |
63 | sigmaT.push_back( new Sigma2qg2qg() ); | |
64 | sigmaU.push_back( new Sigma2qg2qg() ); | |
65 | ||
66 | // Quark-(anti)quark instate. | |
67 | } else { | |
68 | sigmaT.push_back( new Sigma2qq2qq() ); | |
69 | sigmaU.push_back( new Sigma2qq2qq() ); | |
70 | } | |
71 | ||
72 | // Normally store QCD processes to new flavour. | |
73 | if (processLevel > 0) { | |
74 | if (inState == 0) { | |
75 | sigmaT.push_back( new Sigma2gg2qqbar() ); | |
76 | sigmaU.push_back( new Sigma2gg2qqbar() ); | |
77 | sigmaT.push_back( new Sigma2gg2QQbar(4, 121) ); | |
78 | sigmaU.push_back( new Sigma2gg2QQbar(4, 121) ); | |
79 | sigmaT.push_back( new Sigma2gg2QQbar(5, 123) ); | |
80 | sigmaU.push_back( new Sigma2gg2QQbar(5, 123) ); | |
81 | } else if (inState == 2) { | |
82 | sigmaT.push_back( new Sigma2qqbar2gg() ); | |
83 | sigmaU.push_back( new Sigma2qqbar2gg() ); | |
84 | sigmaT.push_back( new Sigma2qqbar2qqbarNew() ); | |
85 | sigmaU.push_back( new Sigma2qqbar2qqbarNew() ); | |
86 | sigmaT.push_back( new Sigma2qqbar2QQbar(4, 122) ); | |
87 | sigmaU.push_back( new Sigma2qqbar2QQbar(4, 122) ); | |
88 | sigmaT.push_back( new Sigma2qqbar2QQbar(5, 124) ); | |
89 | sigmaU.push_back( new Sigma2qqbar2QQbar(5, 124) ); | |
90 | } | |
91 | } | |
92 | ||
93 | // Optionally store electroweak processes, mainly photon production. | |
94 | if (processLevel > 1) { | |
95 | if (inState == 0) { | |
96 | sigmaT.push_back( new Sigma2gg2ggamma() ); | |
97 | sigmaU.push_back( new Sigma2gg2ggamma() ); | |
98 | sigmaT.push_back( new Sigma2gg2gammagamma() ); | |
99 | sigmaU.push_back( new Sigma2gg2gammagamma() ); | |
100 | } else if (inState == 1) { | |
101 | sigmaT.push_back( new Sigma2qg2qgamma() ); | |
102 | sigmaU.push_back( new Sigma2qg2qgamma() ); | |
103 | } else if (inState == 2) { | |
104 | sigmaT.push_back( new Sigma2qqbar2ggamma() ); | |
105 | sigmaU.push_back( new Sigma2qqbar2ggamma() ); | |
106 | sigmaT.push_back( new Sigma2ffbar2gammagamma() ); | |
107 | sigmaU.push_back( new Sigma2ffbar2gammagamma() ); | |
108 | sigmaT.push_back( new Sigma2ffbar2ffbarsgm() ); | |
109 | sigmaU.push_back( new Sigma2ffbar2ffbarsgm() ); | |
110 | } | |
111 | if (inState >= 2) { | |
112 | sigmaT.push_back( new Sigma2ff2fftgmZ() ); | |
113 | sigmaU.push_back( new Sigma2ff2fftgmZ() ); | |
114 | sigmaT.push_back( new Sigma2ff2fftW() ); | |
115 | sigmaU.push_back( new Sigma2ff2fftW() ); | |
116 | } | |
117 | } | |
118 | ||
119 | // Optionally store charmonium and bottomonium production. | |
120 | if (processLevel > 2) { | |
121 | if (inState == 0) { | |
122 | sigmaT.push_back( new Sigma2gg2QQbar3S11g(4, 401) ); | |
123 | sigmaU.push_back( new Sigma2gg2QQbar3S11g(4, 401) ); | |
124 | sigmaT.push_back( new Sigma2gg2QQbar3PJ1g(4, 0, 402) ); | |
125 | sigmaU.push_back( new Sigma2gg2QQbar3PJ1g(4, 0, 402) ); | |
126 | sigmaT.push_back( new Sigma2gg2QQbar3PJ1g(4, 1, 403) ); | |
127 | sigmaU.push_back( new Sigma2gg2QQbar3PJ1g(4, 1, 403) ); | |
128 | sigmaT.push_back( new Sigma2gg2QQbar3PJ1g(4, 2, 404) ); | |
129 | sigmaU.push_back( new Sigma2gg2QQbar3PJ1g(4, 2, 404) ); | |
130 | sigmaT.push_back( new Sigma2gg2QQbarX8g(4, 0, 411) ); | |
131 | sigmaU.push_back( new Sigma2gg2QQbarX8g(4, 0, 411) ); | |
132 | sigmaT.push_back( new Sigma2gg2QQbarX8g(4, 1, 412) ); | |
133 | sigmaU.push_back( new Sigma2gg2QQbarX8g(4, 1, 412) ); | |
134 | sigmaT.push_back( new Sigma2gg2QQbarX8g(4, 2, 413) ); | |
135 | sigmaU.push_back( new Sigma2gg2QQbarX8g(4, 2, 413) ); | |
136 | sigmaT.push_back( new Sigma2gg2QQbar3S11g(5, 501) ); | |
137 | sigmaU.push_back( new Sigma2gg2QQbar3S11g(5, 501) ); | |
138 | sigmaT.push_back( new Sigma2gg2QQbar3PJ1g(5, 0, 502) ); | |
139 | sigmaU.push_back( new Sigma2gg2QQbar3PJ1g(5, 0, 502) ); | |
140 | sigmaT.push_back( new Sigma2gg2QQbar3PJ1g(5, 1, 503) ); | |
141 | sigmaU.push_back( new Sigma2gg2QQbar3PJ1g(5, 1, 503) ); | |
142 | sigmaT.push_back( new Sigma2gg2QQbar3PJ1g(5, 2, 504) ); | |
143 | sigmaU.push_back( new Sigma2gg2QQbar3PJ1g(5, 2, 504) ); | |
144 | sigmaT.push_back( new Sigma2gg2QQbarX8g(5, 0, 511) ); | |
145 | sigmaU.push_back( new Sigma2gg2QQbarX8g(5, 0, 511) ); | |
146 | sigmaT.push_back( new Sigma2gg2QQbarX8g(5, 1, 512) ); | |
147 | sigmaU.push_back( new Sigma2gg2QQbarX8g(5, 1, 512) ); | |
148 | sigmaT.push_back( new Sigma2gg2QQbarX8g(5, 2, 513) ); | |
149 | sigmaU.push_back( new Sigma2gg2QQbarX8g(5, 2, 513) ); | |
150 | } else if (inState == 1) { | |
151 | sigmaT.push_back( new Sigma2qg2QQbar3PJ1q(4, 0, 405) ); | |
152 | sigmaU.push_back( new Sigma2qg2QQbar3PJ1q(4, 0, 405) ); | |
153 | sigmaT.push_back( new Sigma2qg2QQbar3PJ1q(4, 1, 406) ); | |
154 | sigmaU.push_back( new Sigma2qg2QQbar3PJ1q(4, 1, 406) ); | |
155 | sigmaT.push_back( new Sigma2qg2QQbar3PJ1q(4, 2, 407) ); | |
156 | sigmaU.push_back( new Sigma2qg2QQbar3PJ1q(4, 2, 407) ); | |
157 | sigmaT.push_back( new Sigma2qg2QQbarX8q(4, 0, 414) ); | |
158 | sigmaU.push_back( new Sigma2qg2QQbarX8q(4, 0, 414) ); | |
159 | sigmaT.push_back( new Sigma2qg2QQbarX8q(4, 1, 415) ); | |
160 | sigmaU.push_back( new Sigma2qg2QQbarX8q(4, 1, 415) ); | |
161 | sigmaT.push_back( new Sigma2qg2QQbarX8q(4, 2, 416) ); | |
162 | sigmaU.push_back( new Sigma2qg2QQbarX8q(4, 2, 416) ); | |
163 | sigmaT.push_back( new Sigma2qg2QQbar3PJ1q(5, 0, 505) ); | |
164 | sigmaU.push_back( new Sigma2qg2QQbar3PJ1q(5, 0, 505) ); | |
165 | sigmaT.push_back( new Sigma2qg2QQbar3PJ1q(5, 1, 506) ); | |
166 | sigmaU.push_back( new Sigma2qg2QQbar3PJ1q(5, 1, 506) ); | |
167 | sigmaT.push_back( new Sigma2qg2QQbar3PJ1q(5, 2, 507) ); | |
168 | sigmaU.push_back( new Sigma2qg2QQbar3PJ1q(5, 2, 507) ); | |
169 | sigmaT.push_back( new Sigma2qg2QQbarX8q(5, 0, 514) ); | |
170 | sigmaU.push_back( new Sigma2qg2QQbarX8q(5, 0, 514) ); | |
171 | sigmaT.push_back( new Sigma2qg2QQbarX8q(5, 1, 515) ); | |
172 | sigmaU.push_back( new Sigma2qg2QQbarX8q(5, 1, 515) ); | |
173 | sigmaT.push_back( new Sigma2qg2QQbarX8q(5, 2, 516) ); | |
174 | sigmaU.push_back( new Sigma2qg2QQbarX8q(5, 2, 516) ); | |
175 | } else if (inState == 2) { | |
176 | sigmaT.push_back( new Sigma2qqbar2QQbar3PJ1g(4, 0, 408) ); | |
177 | sigmaU.push_back( new Sigma2qqbar2QQbar3PJ1g(4, 0, 408) ); | |
178 | sigmaT.push_back( new Sigma2qqbar2QQbar3PJ1g(4, 1, 409) ); | |
179 | sigmaU.push_back( new Sigma2qqbar2QQbar3PJ1g(4, 1, 409) ); | |
180 | sigmaT.push_back( new Sigma2qqbar2QQbar3PJ1g(4, 2, 410) ); | |
181 | sigmaU.push_back( new Sigma2qqbar2QQbar3PJ1g(4, 2, 410) ); | |
182 | sigmaT.push_back( new Sigma2qqbar2QQbarX8g(4, 0, 417) ); | |
183 | sigmaU.push_back( new Sigma2qqbar2QQbarX8g(4, 0, 417) ); | |
184 | sigmaT.push_back( new Sigma2qqbar2QQbarX8g(4, 1, 418) ); | |
185 | sigmaU.push_back( new Sigma2qqbar2QQbarX8g(4, 1, 418) ); | |
186 | sigmaT.push_back( new Sigma2qqbar2QQbarX8g(4, 2, 419) ); | |
187 | sigmaU.push_back( new Sigma2qqbar2QQbarX8g(4, 2, 419) ); | |
188 | sigmaT.push_back( new Sigma2qqbar2QQbar3PJ1g(5, 0, 508) ); | |
189 | sigmaU.push_back( new Sigma2qqbar2QQbar3PJ1g(5, 0, 508) ); | |
190 | sigmaT.push_back( new Sigma2qqbar2QQbar3PJ1g(5, 1, 509) ); | |
191 | sigmaU.push_back( new Sigma2qqbar2QQbar3PJ1g(5, 1, 509) ); | |
192 | sigmaT.push_back( new Sigma2qqbar2QQbar3PJ1g(5, 2, 510) ); | |
193 | sigmaU.push_back( new Sigma2qqbar2QQbar3PJ1g(5, 2, 510) ); | |
194 | sigmaT.push_back( new Sigma2qqbar2QQbarX8g(5, 0, 517) ); | |
195 | sigmaU.push_back( new Sigma2qqbar2QQbarX8g(5, 0, 517) ); | |
196 | sigmaT.push_back( new Sigma2qqbar2QQbarX8g(5, 1, 518) ); | |
197 | sigmaU.push_back( new Sigma2qqbar2QQbarX8g(5, 1, 518) ); | |
198 | sigmaT.push_back( new Sigma2qqbar2QQbarX8g(5, 2, 519) ); | |
199 | sigmaU.push_back( new Sigma2qqbar2QQbarX8g(5, 2, 519) ); | |
200 | } | |
201 | } | |
202 | ||
203 | // Resize arrays to match sizes above. | |
204 | nChan = sigmaT.size(); | |
205 | needMasses.resize(nChan); | |
206 | m3Fix.resize(nChan); | |
207 | m4Fix.resize(nChan); | |
208 | sHatMin.resize(nChan); | |
209 | sigmaTval.resize(nChan); | |
210 | sigmaUval.resize(nChan); | |
211 | ||
212 | // Initialize the processes. | |
213 | for (int i = 0; i < nChan; ++i) { | |
214 | sigmaT[i]->init( infoPtr, settingsPtr, particleDataPtr, rndmPtr, | |
215 | beamAPtr, beamBPtr, couplingsPtr); | |
216 | sigmaT[i]->initProc(); | |
217 | sigmaU[i]->init( infoPtr, settingsPtr, particleDataPtr, rndmPtr, | |
218 | beamAPtr, beamBPtr, couplingsPtr); | |
219 | sigmaU[i]->initProc(); | |
220 | ||
221 | // Prepare for massive kinematics (but fixed masses!) where required. | |
222 | needMasses[i] = false; | |
223 | int id3Mass = sigmaT[i]->id3Mass(); | |
224 | int id4Mass = sigmaT[i]->id4Mass(); | |
225 | m3Fix[i] = 0.; | |
226 | m4Fix[i] = 0.; | |
227 | if (id3Mass > 0 || id4Mass > 0) { | |
228 | needMasses[i] = true; | |
229 | m3Fix[i] = particleDataPtr->m0(id3Mass); | |
230 | m4Fix[i] = particleDataPtr->m0(id4Mass); | |
231 | } | |
232 | sHatMin[i] = pow2( m3Fix[i] + m4Fix[i] + MASSMARGIN); | |
233 | } | |
234 | ||
235 | // Done. | |
236 | return true; | |
237 | ||
238 | } | |
239 | ||
240 | //-------------------------------------------------------------------------- | |
241 | ||
242 | // Calculate cross section summed over possibilities. | |
243 | ||
244 | double SigmaMultiparton::sigma( int id1, int id2, double x1, double x2, | |
245 | double sHat, double tHat, double uHat, double alpS, double alpEM, | |
246 | bool restore, bool pickOtherIn) { | |
247 | ||
248 | // Choose either the dominant process (in slot 0) or the rest of them. | |
249 | if (restore) pickOther = pickOtherIn; | |
250 | else pickOther = (rndmPtr->flat() < OTHERFRAC); | |
251 | ||
252 | // Iterate over all subprocesses. | |
253 | sigmaTsum = 0.; | |
254 | sigmaUsum = 0.; | |
255 | for (int i = 0; i < nChan; ++i) { | |
256 | sigmaTval[i] = 0.; | |
257 | sigmaUval[i] = 0.; | |
258 | ||
259 | // Skip the not chosen processes. | |
260 | if (i == 0 && pickOther) continue; | |
261 | if (i > 0 && !pickOther) continue; | |
262 | ||
263 | // t-channel-sampling contribution. | |
264 | if (sHat > sHatMin[i]) { | |
265 | sigmaT[i]->set2KinMPI( x1, x2, sHat, tHat, uHat, | |
266 | alpS, alpEM, needMasses[i], m3Fix[i], m4Fix[i]); | |
267 | sigmaTval[i] = sigmaT[i]->sigmaHatWrap(id1, id2); | |
268 | sigmaT[i]->pickInState(id1, id2); | |
269 | // Correction factor for tHat rescaling in massive kinematics. | |
270 | if (needMasses[i]) sigmaTval[i] *= sigmaT[i]->sHBetaMPI() / sHat; | |
271 | sigmaTsum += sigmaTval[i]; | |
272 | } | |
273 | ||
274 | // u-channel-sampling contribution. | |
275 | if (sHat > sHatMin[i]) { | |
276 | sigmaU[i]->set2KinMPI( x1, x2, sHat, uHat, tHat, | |
277 | alpS, alpEM, needMasses[i], m3Fix[i], m4Fix[i]); | |
278 | sigmaUval[i] = sigmaU[i]->sigmaHatWrap( id1, id2); | |
279 | sigmaU[i]->pickInState(id1, id2); | |
280 | // Correction factor for tHat rescaling in massive kinematics. | |
281 | if (needMasses[i]) sigmaUval[i] *= sigmaU[i]->sHBetaMPI() / sHat; | |
282 | sigmaUsum += sigmaUval[i]; | |
283 | } | |
284 | ||
285 | // Average of t- and u-channel sampling; corrected for not selected channels. | |
286 | } | |
287 | double sigmaAvg = 0.5 * (sigmaTsum + sigmaUsum); | |
288 | if (pickOther) sigmaAvg /= OTHERFRAC; | |
289 | if (!pickOther) sigmaAvg /= (1. - OTHERFRAC); | |
290 | return sigmaAvg; | |
291 | ||
292 | } | |
293 | ||
294 | //-------------------------------------------------------------------------- | |
295 | ||
296 | // Return one subprocess, picked according to relative cross sections. | |
297 | ||
298 | SigmaProcess* SigmaMultiparton::sigmaSel() { | |
299 | ||
300 | // Decide between t- and u-channel-sampled kinematics. | |
301 | pickedU = (rndmPtr->flat() * (sigmaTsum + sigmaUsum) < sigmaUsum); | |
302 | ||
303 | // Pick one of t-channel-sampled processes. | |
304 | if (!pickedU) { | |
305 | double sigmaRndm = sigmaTsum * rndmPtr->flat(); | |
306 | int iPick = -1; | |
307 | do sigmaRndm -= sigmaTval[++iPick]; | |
308 | while (sigmaRndm > 0.); | |
309 | return sigmaT[iPick]; | |
310 | ||
311 | // Pick one of u-channel-sampled processes. | |
312 | } else { | |
313 | double sigmaRndm = sigmaUsum * rndmPtr->flat(); | |
314 | int iPick = -1; | |
315 | do sigmaRndm -= sigmaUval[++iPick]; | |
316 | while (sigmaRndm > 0.); | |
317 | return sigmaU[iPick]; | |
318 | } | |
319 | ||
320 | } | |
321 | ||
322 | //========================================================================== | |
323 | ||
324 | // The MultipartonInteractions class. | |
325 | ||
326 | //-------------------------------------------------------------------------- | |
327 | ||
328 | // Constants: could be changed here if desired, but normally should not. | |
329 | // These are of technical nature, as described for each. | |
330 | ||
331 | // Factorization scale pT2 by default, but could be shifted to pT2 + pT02. | |
332 | // (A priori possible, but problems for flavour threshold interpretation.) | |
333 | const bool MultipartonInteractions::SHIFTFACSCALE = false; | |
334 | ||
335 | // Pick one parton to represent rescattering cross section to speed up. | |
336 | const bool MultipartonInteractions::PREPICKRESCATTER = true; | |
337 | ||
338 | // Naive upper estimate of cross section too pessimistic, so reduce by this. | |
339 | const double MultipartonInteractions::SIGMAFUDGE = 0.8; | |
340 | ||
341 | // The r value above, picked to allow a flatter correct/trial cross section. | |
342 | const double MultipartonInteractions::RPT20 = 0.25; | |
343 | ||
344 | // Reduce pT0 by factor pT0STEP if sigmaInt < SIGMASTEP * sigmaND. | |
345 | const double MultipartonInteractions::PT0STEP = 0.9; | |
346 | const double MultipartonInteractions::SIGMASTEP = 1.1; | |
347 | ||
348 | // Stop if pT0 or pTmin fall below this, or alpha_s blows up. | |
349 | const double MultipartonInteractions::PT0MIN = 0.2; | |
350 | ||
351 | // Refuse too low expPow in impact parameter profile. | |
352 | const double MultipartonInteractions::EXPPOWMIN = 0.4; | |
353 | ||
354 | // Define low-b region by interaction rate above given number. | |
355 | const double MultipartonInteractions::PROBATLOWB = 0.6; | |
356 | ||
357 | // Basic step size for b integration; sometimes modified. | |
358 | const double MultipartonInteractions::BSTEP = 0.01; | |
359 | ||
360 | // Stop b integration when integrand dropped enough. | |
361 | const double MultipartonInteractions::BMAX = 1e-8; | |
362 | ||
363 | // Do not allow too large argument to exp function. | |
364 | const double MultipartonInteractions::EXPMAX = 50.; | |
365 | ||
366 | // Convergence criterion for k iteration. | |
367 | const double MultipartonInteractions::KCONVERGE = 1e-7; | |
368 | ||
369 | // Conversion of GeV^{-2} to mb for cross section. | |
370 | const double MultipartonInteractions::CONVERT2MB = 0.389380; | |
371 | ||
372 | // Stay away from division by zero in Jacobian for tHat -> pT2. | |
373 | const double MultipartonInteractions::ROOTMIN = 0.01; | |
374 | ||
375 | // No need to reinitialize parameters if energy close to previous. | |
376 | const double MultipartonInteractions::ECMDEV = 0.01; | |
377 | ||
378 | // Settings for x-dependent matter profile: | |
379 | // Number of bins in b (with these settings, no bStep increase and | |
380 | // reintegration needed with a1 ~ 0.20 up to ECM ~ 40TeV). | |
381 | const int MultipartonInteractions::XDEP_BBIN = 500; | |
382 | // Initial value of a0. | |
383 | const double MultipartonInteractions::XDEP_A0 = 1.0; | |
384 | // Width of form ( XDEP_A1 + a1 * log(1 / x) ). | |
385 | const double MultipartonInteractions::XDEP_A1 = 1.0; | |
386 | // Initial step size in b and increment. | |
387 | const double MultipartonInteractions::XDEP_BSTEP = 0.02; | |
388 | const double MultipartonInteractions::XDEP_BSTEPINC = 0.01; | |
389 | // Overlap-weighted cross section in last bin of b must be beneath | |
390 | // XDEP_CUTOFF * sigmaInt. | |
391 | const double MultipartonInteractions::XDEP_CUTOFF = 1e-4; | |
392 | // a0 is calculated in units of sqrt(mb), so convert to fermi. | |
393 | const double MultipartonInteractions::XDEP_SMB2FM = sqrt(0.1); | |
394 | ||
395 | // Only write warning when weight clearly above unity. | |
396 | const double MultipartonInteractions::WTACCWARN = 1.1; | |
397 | ||
398 | //-------------------------------------------------------------------------- | |
399 | ||
400 | // Initialize the generation process for given beams. | |
401 | ||
402 | bool MultipartonInteractions::init( bool doMPIinit, int iDiffSysIn, | |
403 | Info* infoPtrIn, Settings& settings, ParticleData* particleDataPtr, | |
404 | Rndm* rndmPtrIn, BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn, | |
405 | Couplings* couplingsPtrIn, PartonSystems* partonSystemsPtrIn, | |
406 | SigmaTotal* sigmaTotPtrIn, UserHooks* userHooksPtrIn, ostream& os) { | |
407 | ||
408 | // Store input pointers for future use. Done if no initialization. | |
409 | iDiffSys = iDiffSysIn; | |
410 | infoPtr = infoPtrIn; | |
411 | rndmPtr = rndmPtrIn; | |
412 | beamAPtr = beamAPtrIn; | |
413 | beamBPtr = beamBPtrIn; | |
414 | couplingsPtr = couplingsPtrIn; | |
415 | partonSystemsPtr = partonSystemsPtrIn; | |
416 | sigmaTotPtr = sigmaTotPtrIn; | |
417 | userHooksPtr = userHooksPtrIn; | |
418 | if (!doMPIinit) return false; | |
419 | ||
420 | // If both beams are baryons then softer PDF's than for mesons/Pomerons. | |
421 | hasBaryonBeams = ( beamAPtr->isBaryon() && beamBPtr->isBaryon() ); | |
422 | ||
423 | // Matching in pT of hard interaction to further interactions. | |
424 | pTmaxMatch = settings.mode("MultipartonInteractions:pTmaxMatch"); | |
425 | ||
426 | // Parameters of alphaStrong generation. | |
427 | alphaSvalue = settings.parm("MultipartonInteractions:alphaSvalue"); | |
428 | alphaSorder = settings.mode("MultipartonInteractions:alphaSorder"); | |
429 | ||
430 | // Parameters of alphaEM generation. | |
431 | alphaEMorder = settings.mode("MultipartonInteractions:alphaEMorder"); | |
432 | ||
433 | // Parameters of cross section generation. | |
434 | Kfactor = settings.parm("MultipartonInteractions:Kfactor"); | |
435 | ||
436 | // Regularization of QCD evolution for pT -> 0. | |
437 | pT0Ref = settings.parm("MultipartonInteractions:pT0Ref"); | |
438 | ecmRef = settings.parm("MultipartonInteractions:ecmRef"); | |
439 | ecmPow = settings.parm("MultipartonInteractions:ecmPow"); | |
440 | pTmin = settings.parm("MultipartonInteractions:pTmin"); | |
441 | ||
442 | // Impact parameter profile: nondiffractive topologies. | |
443 | if (iDiffSys == 0) { | |
444 | bProfile = settings.mode("MultipartonInteractions:bProfile"); | |
445 | coreRadius = settings.parm("MultipartonInteractions:coreRadius"); | |
446 | coreFraction = settings.parm("MultipartonInteractions:coreFraction"); | |
447 | expPow = settings.parm("MultipartonInteractions:expPow"); | |
448 | expPow = max(EXPPOWMIN, expPow); | |
449 | // x-dependent impact parameter profile. | |
450 | a1 = settings.parm("MultipartonInteractions:a1"); | |
451 | ||
452 | // Impact parameter profile: diffractive topologies. | |
453 | } else { | |
454 | bProfile = settings.mode("Diffraction:bProfile"); | |
455 | coreRadius = settings.parm("Diffraction:coreRadius"); | |
456 | coreFraction = settings.parm("Diffraction:coreFraction"); | |
457 | expPow = settings.parm("Diffraction:expPow"); | |
458 | expPow = max(EXPPOWMIN, expPow); | |
459 | } | |
460 | ||
461 | // Common choice of "pT" scale for determining impact parameter. | |
462 | bSelScale = settings.mode("MultipartonInteractions:bSelScale"); | |
463 | ||
464 | // Process sets to include in machinery. | |
465 | processLevel = settings.mode("MultipartonInteractions:processLevel"); | |
466 | ||
467 | // Parameters of rescattering description. | |
468 | allowRescatter = settings.flag("MultipartonInteractions:allowRescatter"); | |
469 | allowDoubleRes | |
470 | = settings.flag("MultipartonInteractions:allowDoubleRescatter"); | |
471 | rescatterMode = settings.mode("MultipartonInteractions:rescatterMode"); | |
472 | ySepResc = settings.parm("MultipartonInteractions:ySepRescatter"); | |
473 | deltaYResc = settings.parm("MultipartonInteractions:deltaYRescatter"); | |
474 | ||
475 | // Rescattering not yet implemented for x-dependent impact profile. | |
476 | if (bProfile == 4) allowRescatter = false; | |
477 | ||
478 | // A global recoil FSR stategy restricts rescattering. | |
479 | globalRecoilFSR = settings.flag("TimeShower:globalRecoil"); | |
480 | nMaxGlobalRecoilFSR = settings.mode("TimeShower:nMaxGlobalRecoil"); | |
481 | ||
482 | // Various other parameters. | |
483 | nQuarkIn = settings.mode("MultipartonInteractions:nQuarkIn"); | |
484 | nSample = settings.mode("MultipartonInteractions:nSample"); | |
485 | ||
486 | // Optional dampening at small pT's when large multiplicities. | |
487 | enhanceScreening = settings.mode("MultipartonInteractions:enhanceScreening"); | |
488 | ||
489 | // Parameters for diffractive systems. | |
490 | sigmaPomP = settings.parm("Diffraction:sigmaRefPomP"); | |
491 | mPomP = settings.parm("Diffraction:mRefPomP"); | |
492 | pPomP = settings.parm("Diffraction:mPowPomP"); | |
493 | mMinPertDiff = settings.parm("Diffraction:mMinPert"); | |
494 | ||
495 | // Possibility to allow user veto of MPI | |
496 | canVetoMPI = (userHooksPtr != 0) ? userHooksPtr->canVetoMPIEmission() | |
497 | : false; | |
498 | ||
499 | // Some common combinations for double Gaussian, as shorthand. | |
500 | if (bProfile == 2) { | |
501 | fracA = pow2(1. - coreFraction); | |
502 | fracB = 2. * coreFraction * (1. - coreFraction); | |
503 | fracC = pow2(coreFraction); | |
504 | radius2B = 0.5 * (1. + pow2(coreRadius)); | |
505 | radius2C = pow2(coreRadius); | |
506 | ||
507 | // Some common combinations for exp(b^pow), as shorthand. | |
508 | } else if (bProfile == 3) { | |
509 | hasLowPow = (expPow < 2.); | |
510 | expRev = 2. / expPow - 1.; | |
511 | } | |
512 | ||
513 | // Initialize alpha_strong generation. | |
514 | alphaS.init( alphaSvalue, alphaSorder); | |
515 | double Lambda3 = alphaS.Lambda3(); | |
516 | ||
517 | // Initialize alphaEM generation. | |
518 | alphaEM.init( alphaEMorder, &settings); | |
519 | ||
520 | // Attach matrix-element calculation objects. | |
521 | sigma2gg.init( 0, processLevel, infoPtr, &settings, particleDataPtr, | |
522 | rndmPtr, beamAPtr, beamBPtr, couplingsPtr); | |
523 | sigma2qg.init( 1, processLevel, infoPtr, &settings, particleDataPtr, | |
524 | rndmPtr, beamAPtr, beamBPtr, couplingsPtr); | |
525 | sigma2qqbarSame.init( 2, processLevel, infoPtr, &settings, particleDataPtr, | |
526 | rndmPtr, beamAPtr, beamBPtr, couplingsPtr); | |
527 | sigma2qq.init( 3, processLevel, infoPtr, &settings, particleDataPtr, | |
528 | rndmPtr, beamAPtr, beamBPtr, couplingsPtr); | |
529 | ||
530 | // Calculate invariant mass of system. | |
531 | eCM = infoPtr->eCM(); | |
532 | sCM = eCM * eCM; | |
533 | mMaxPertDiff = eCM; | |
534 | eCMsave = eCM; | |
535 | ||
536 | // Get the total inelastic and nondiffractive cross section. | |
537 | if (!sigmaTotPtr->hasSigmaTot()) return false; | |
538 | bool isNonDiff = (iDiffSys == 0); | |
539 | sigmaND = sigmaTotPtr->sigmaND(); | |
540 | double sigmaMaxViol = 0.; | |
541 | ||
542 | // Output initialization info - first part. | |
543 | bool showMPI = settings.flag("Init:showMultipartonInteractions"); | |
544 | if (showMPI) { | |
545 | os << "\n *------- PYTHIA Multiparton Interactions Initialization " | |
546 | << "---------* \n" | |
547 | << " | " | |
548 | << " | \n"; | |
549 | if (isNonDiff) | |
550 | os << " | minbias, sigmaNonDiffractive = " << fixed | |
551 | << setprecision(2) << setw(7) << sigmaND << " mb | \n"; | |
552 | else if (iDiffSys == 1) | |
553 | os << " | diffraction XB " | |
554 | << " | \n"; | |
555 | else if (iDiffSys == 2) | |
556 | os << " | diffraction AX " | |
557 | << " | \n"; | |
558 | else if (iDiffSys == 3) | |
559 | os << " | diffraction AXB " | |
560 | << " | \n"; | |
561 | os << " | " | |
562 | << " | \n"; | |
563 | } | |
564 | ||
565 | // For diffraction need to cover range of diffractive masses. | |
566 | nStep = (iDiffSys == 0) ? 1 : 5; | |
567 | eStepSize = (nStep < 2) ? 1. | |
568 | : log(mMaxPertDiff / mMinPertDiff) / (nStep - 1.); | |
569 | for (int iStep = 0; iStep < nStep; ++iStep) { | |
570 | ||
571 | // Update and output current diffractive mass and | |
572 | // fictitious Pomeron-proton cross section for normalization. | |
573 | if (nStep > 1) { | |
574 | eCM = mMinPertDiff * pow( mMaxPertDiff / mMinPertDiff, | |
575 | iStep / (nStep - 1.) ); | |
576 | sCM = eCM * eCM; | |
577 | sigmaND = sigmaPomP * pow( eCM / mPomP, pPomP); | |
578 | if (showMPI) os << " | diffractive mass = " << scientific | |
579 | << setprecision(3) << setw(9) << eCM << " GeV and sigmaNorm = " | |
580 | << fixed << setw(6) << sigmaND << " mb | \n"; | |
581 | } | |
582 | ||
583 | // Set current pT0 scale. | |
584 | pT0 = pT0Ref * pow(eCM / ecmRef, ecmPow); | |
585 | ||
586 | // The pT0 value may need to be decreased, if sigmaInt < sigmaND. | |
587 | double pT4dSigmaMaxBeg = 0.; | |
588 | for ( ; ; ) { | |
589 | ||
590 | // Derived pT kinematics combinations. | |
591 | pT20 = pT0*pT0; | |
592 | pT2min = pTmin*pTmin; | |
593 | pTmax = 0.5*eCM; | |
594 | pT2max = pTmax*pTmax; | |
595 | pT20R = RPT20 * pT20; | |
596 | pT20minR = pT2min + pT20R; | |
597 | pT20maxR = pT2max + pT20R; | |
598 | pT20min0maxR = pT20minR * pT20maxR; | |
599 | pT2maxmin = pT2max - pT2min; | |
600 | ||
601 | // Provide upper estimate of interaction rate d(Prob)/d(pT2). | |
602 | upperEnvelope(); | |
603 | ||
604 | // Setup binning in b for x-dependent matter profile. | |
605 | if (bProfile == 4) { | |
606 | sigmaIntWgt.resize(XDEP_BBIN); | |
607 | sigmaSumWgt.resize(XDEP_BBIN); | |
608 | bstepNow = XDEP_BSTEP; | |
609 | } | |
610 | ||
611 | // Integrate the parton-parton interaction cross section. | |
612 | pT4dSigmaMaxBeg = pT4dSigmaMax; | |
613 | jetCrossSection(); | |
614 | ||
615 | // If the overlap-weighted cross section has not fallen below | |
616 | // cutoff, then increase bin size in b and reintegrate. | |
617 | while (bProfile == 4 | |
618 | && sigmaIntWgt[XDEP_BBIN - 1] > XDEP_CUTOFF * sigmaInt) { | |
619 | bstepNow += XDEP_BSTEPINC; | |
620 | jetCrossSection(); | |
621 | } | |
622 | ||
623 | // Sufficiently big SigmaInt or reduce pT0; maybe also pTmin. | |
624 | if (sigmaInt > SIGMASTEP * sigmaND) break; | |
625 | if (showMPI) os << fixed << setprecision(2) << " | pT0 = " | |
626 | << setw(5) << pT0 << " gives sigmaInteraction = " << setw(7) | |
627 | << sigmaInt << " mb: rejected | \n"; | |
628 | if (pTmin > pT0) pTmin *= PT0STEP; | |
629 | pT0 *= PT0STEP; | |
630 | ||
631 | // Give up if pT0 and pTmin fall too low. | |
632 | if ( max(pT0, pTmin) < max(PT0MIN, Lambda3) ) { | |
633 | infoPtr->errorMsg("Error in MultipartonInteractions::init:" | |
634 | " failed to find acceptable pT0 and pTmin"); | |
635 | infoPtr->setTooLowPTmin(true); | |
636 | return false; | |
637 | } | |
638 | } | |
639 | ||
640 | // Output for accepted pT0. | |
641 | if (showMPI) os << fixed << setprecision(2) << " | pT0 = " | |
642 | << setw(5) << pT0 << " gives sigmaInteraction = "<< setw(7) | |
643 | << sigmaInt << " mb: accepted | \n"; | |
644 | ||
645 | // Calculate factor relating matter overlap and interaction rate. | |
646 | overlapInit(); | |
647 | ||
648 | // Maximum violation relative to first estimate. | |
649 | sigmaMaxViol = max( sigmaMaxViol, pT4dSigmaMax / pT4dSigmaMaxBeg); | |
650 | ||
651 | // Save values calculated. | |
652 | if (nStep > 1) { | |
653 | pT0Save[iStep] = pT0; | |
654 | pT4dSigmaMaxSave[iStep] = pT4dSigmaMax; | |
655 | pT4dProbMaxSave[iStep] = pT4dProbMax; | |
656 | sigmaIntSave[iStep] = sigmaInt; | |
657 | for (int j = 0; j <= 100; ++j) sudExpPTSave[iStep][j] = sudExpPT[j]; | |
658 | zeroIntCorrSave[iStep] = zeroIntCorr; | |
659 | normOverlapSave[iStep] = normOverlap; | |
660 | kNowSave[iStep] = kNow; | |
661 | bAvgSave[iStep] = bAvg; | |
662 | bDivSave[iStep] = bDiv; | |
663 | probLowBSave[iStep] = probLowB; | |
664 | fracAhighSave[iStep] = fracAhigh; | |
665 | fracBhighSave[iStep] = fracBhigh; | |
666 | fracChighSave[iStep] = fracBhigh; | |
667 | fracABChighSave[iStep] = fracABChigh; | |
668 | cDivSave[iStep] = cDiv; | |
669 | cMaxSave[iStep] = cMax; | |
670 | } | |
671 | ||
672 | // End of loop over diffractive masses. | |
673 | } | |
674 | ||
675 | // Output details for x-dependent matter profile. | |
676 | if (bProfile == 4 && showMPI) | |
677 | os << " | " | |
678 | << " | \n" | |
679 | << fixed << setprecision(2) | |
680 | << " | x-dependent matter profile: a1 = " << a1 << ", " | |
681 | << "a0 = " << a0now * XDEP_SMB2FM << ", bStep = " | |
682 | << bstepNow << " | \n"; | |
683 | ||
684 | // End initialization printout. | |
685 | if (showMPI) os << " | " | |
686 | << " | \n" | |
687 | << " *------- End PYTHIA Multiparton Interactions Initialization" | |
688 | << " -----* " << endl; | |
689 | ||
690 | // Amount of violation from upperEnvelope to jetCrossSection. | |
691 | if (sigmaMaxViol > 1.) { | |
692 | ostringstream osWarn; | |
693 | osWarn << "by factor " << fixed << setprecision(3) << sigmaMaxViol; | |
694 | infoPtr->errorMsg("Warning in MultipartonInteractions::init:" | |
695 | " maximum increased", osWarn.str()); | |
696 | } | |
697 | ||
698 | // Reset statistics. | |
699 | SigmaMultiparton* dSigma; | |
700 | for (int i = 0; i < 4; ++i) { | |
701 | if (i == 0) dSigma = &sigma2gg; | |
702 | else if (i == 1) dSigma = &sigma2qg; | |
703 | else if (i == 2) dSigma = &sigma2qqbarSame; | |
704 | else dSigma = &sigma2qq; | |
705 | int nProc = dSigma->nProc(); | |
706 | for (int iProc = 0; iProc < nProc; ++iProc) | |
707 | nGen[ dSigma->codeProc(iProc) ] = 0; | |
708 | } | |
709 | ||
710 | // Additional setup for x-dependent matter profile. | |
711 | if (bProfile == 4) { | |
712 | sigmaIntWgt.clear(); | |
713 | sigmaSumWgt.clear(); | |
714 | } | |
715 | // No preselection of sea/valence content and initialise a0. | |
716 | vsc1 = 0; | |
717 | vsc2 = 0; | |
718 | ||
719 | // Done. | |
720 | return true; | |
721 | } | |
722 | ||
723 | //-------------------------------------------------------------------------- | |
724 | ||
725 | // Reset impact parameter choice and update the CM energy. | |
726 | // For diffraction also interpolate parameters to current CM energy. | |
727 | ||
728 | void MultipartonInteractions::reset( ) { | |
729 | ||
730 | // Reset impact parameter choice. | |
731 | bIsSet = false; | |
732 | bSetInFirst = false; | |
733 | ||
734 | // Update CM energy. Done if not diffraction and not new energy. | |
735 | eCM = infoPtr->eCM(); | |
736 | sCM = eCM * eCM; | |
737 | if (nStep == 1 || abs( eCM / eCMsave - 1.) < ECMDEV) return; | |
738 | ||
739 | // Set fictitious Pomeron-proton cross section for diffractive system. | |
740 | sigmaND = sigmaPomP * pow( eCM / mPomP, pPomP); | |
741 | ||
742 | // Current interpolation point. | |
743 | eCMsave = eCM; | |
744 | eStepSave = log(eCM / mMinPertDiff) / eStepSize; | |
745 | iStepFrom = max( 0, min( nStep - 2, int( eStepSave) ) ); | |
746 | iStepTo = iStepFrom + 1; | |
747 | eStepTo = max( 0., min( 1., eStepSave - iStepFrom) ); | |
748 | eStepFrom = 1. - eStepTo; | |
749 | ||
750 | // Update pT0 and combinations derived from it. | |
751 | pT0 = eStepFrom * pT0Save[iStepFrom] | |
752 | + eStepTo * pT0Save[iStepTo]; | |
753 | pT20 = pT0*pT0; | |
754 | pT2min = pTmin*pTmin; | |
755 | pTmax = 0.5*eCM; | |
756 | pT2max = pTmax*pTmax; | |
757 | pT20R = RPT20 * pT20; | |
758 | pT20minR = pT2min + pT20R; | |
759 | pT20maxR = pT2max + pT20R; | |
760 | pT20min0maxR = pT20minR * pT20maxR; | |
761 | pT2maxmin = pT2max - pT2min; | |
762 | ||
763 | // Update other parameters used in pT choice. | |
764 | pT4dSigmaMax = eStepFrom * pT4dSigmaMaxSave[iStepFrom] | |
765 | + eStepTo * pT4dSigmaMaxSave[iStepTo]; | |
766 | pT4dProbMax = eStepFrom * pT4dProbMaxSave[iStepFrom] | |
767 | + eStepTo * pT4dProbMaxSave[iStepTo]; | |
768 | sigmaInt = eStepFrom * sigmaIntSave[iStepFrom] | |
769 | + eStepTo * sigmaIntSave[iStepTo]; | |
770 | for (int j = 0; j <= 100; ++j) | |
771 | sudExpPT[j] = eStepFrom * sudExpPTSave[iStepFrom][j] | |
772 | + eStepTo * sudExpPTSave[iStepTo][j]; | |
773 | ||
774 | // Update parameters related to the impact-parameter picture. | |
775 | zeroIntCorr = eStepFrom * zeroIntCorrSave[iStepFrom] | |
776 | + eStepTo * zeroIntCorrSave[iStepTo]; | |
777 | normOverlap = eStepFrom * normOverlapSave[iStepFrom] | |
778 | + eStepTo * normOverlapSave[iStepTo]; | |
779 | kNow = eStepFrom * kNowSave[iStepFrom] | |
780 | + eStepTo * kNowSave[iStepTo]; | |
781 | bAvg = eStepFrom * bAvgSave[iStepFrom] | |
782 | + eStepTo * bAvgSave[iStepTo]; | |
783 | bDiv = eStepFrom * bDivSave[iStepFrom] | |
784 | + eStepTo * bDivSave[iStepTo]; | |
785 | probLowB = eStepFrom * probLowBSave[iStepFrom] | |
786 | + eStepTo * probLowBSave[iStepTo]; | |
787 | fracAhigh = eStepFrom * fracAhighSave[iStepFrom] | |
788 | + eStepTo * fracAhighSave[iStepTo]; | |
789 | fracBhigh = eStepFrom * fracBhighSave[iStepFrom] | |
790 | + eStepTo * fracBhighSave[iStepTo]; | |
791 | fracChigh = eStepFrom * fracChighSave[iStepFrom] | |
792 | + eStepTo * fracChighSave[iStepTo]; | |
793 | fracABChigh = eStepFrom * fracABChighSave[iStepFrom] | |
794 | + eStepTo * fracABChighSave[iStepTo]; | |
795 | cDiv = eStepFrom * cDivSave[iStepFrom] | |
796 | + eStepTo * cDivSave[iStepTo]; | |
797 | cMax = eStepFrom * cMaxSave[iStepFrom] | |
798 | + eStepTo * cMaxSave[iStepTo]; | |
799 | ||
800 | } | |
801 | ||
802 | //-------------------------------------------------------------------------- | |
803 | ||
804 | // Select first = hardest pT in minbias process. | |
805 | // Requires separate treatment at low and high b values. | |
806 | ||
807 | void MultipartonInteractions::pTfirst() { | |
808 | ||
809 | // Pick impact parameter and thereby interaction rate enhancement. | |
810 | // This is not used for the x-dependent matter profile, which | |
811 | // instead uses trial interactions. | |
812 | if (bProfile == 4) isAtLowB = false; | |
813 | else overlapFirst(); | |
814 | bSetInFirst = true; | |
815 | double WTacc; | |
816 | ||
817 | // At low b values evolve downwards with Sudakov. | |
818 | if (isAtLowB) { | |
819 | pT2 = pT2max; | |
820 | do { | |
821 | ||
822 | // Pick a pT using a quick-and-dirty cross section estimate. | |
823 | pT2 = fastPT2(pT2); | |
824 | ||
825 | // If fallen below lower cutoff then need to restart at top. | |
826 | if (pT2 < pT2min) { | |
827 | pT2 = pT2max; | |
828 | WTacc = 0.; | |
829 | ||
830 | // Else pick complete kinematics and evaluate cross-section correction. | |
831 | } else { | |
832 | WTacc = sigmaPT2scatter(true) / dSigmaApprox; | |
833 | if (WTacc > WTACCWARN) infoPtr->errorMsg("Warning in " | |
834 | "MultipartonInteractions::pTfirst: weight above unity"); | |
835 | } | |
836 | ||
837 | // Loop until acceptable pT and acceptable kinematics. | |
838 | } while (WTacc < rndmPtr->flat() || !dSigmaDtSel->final2KinMPI()); | |
839 | ||
840 | // At high b values make preliminary pT choice without Sudakov factor. | |
841 | } else { | |
842 | ||
843 | while (true) { | |
844 | do { | |
845 | pT2 = pT20min0maxR / (pT20minR + rndmPtr->flat() * pT2maxmin) - pT20R; | |
846 | ||
847 | // Evaluate upper estimate of cross section for this pT2 choice. | |
848 | dSigmaApprox = pT4dSigmaMax / pow2(pT2 + pT20R); | |
849 | ||
850 | // Pick complete kinematics and evaluate cross-section correction. | |
851 | WTacc = sigmaPT2scatter(true) / dSigmaApprox; | |
852 | ||
853 | // Evaluate and include Sudakov factor. | |
854 | if (bProfile != 4) WTacc *= sudakov( pT2, enhanceB); | |
855 | ||
856 | // Warn for weight above unity | |
857 | if (WTacc > WTACCWARN) infoPtr->errorMsg("Warning in " | |
858 | "MultipartonInteractions::pTfirst: weight above unity"); | |
859 | ||
860 | // Loop until acceptable pT and acceptable kinematics. | |
861 | } while (WTacc < rndmPtr->flat() || !dSigmaDtSel->final2KinMPI()); | |
862 | ||
863 | // For x-dependent matter profile, use trial interactions to | |
864 | // generate Sudakov, otherwise done. | |
865 | if (bProfile != 4) break; | |
866 | else { | |
867 | // Save details of the original hard interaction. | |
868 | pT2Save = pT2; | |
869 | id1Save = id1; | |
870 | id2Save = id2; | |
871 | x1Save = x1; | |
872 | x2Save = x2; | |
873 | sHatSave = sHat; | |
874 | tHatSave = tHat; | |
875 | uHatSave = uHat; | |
876 | alpSsave = alpS; | |
877 | alpEMsave = alpEM; | |
878 | pT2FacSave = pT2Fac; | |
879 | pT2RenSave = pT2Ren; | |
880 | xPDF1nowSave = xPDF1now; | |
881 | xPDF2nowSave = xPDF2now; | |
882 | // Save accepted kinematics and pointer to SigmaProcess. | |
883 | dSigmaDtSel->saveKin(); | |
884 | dSigmaDtSelSave = dSigmaDtSel; | |
885 | ||
886 | // Put x1, x2 information into beam pointers to get correct | |
887 | // PDF rescaling in trial interaction (for hard process, can | |
888 | // be sea or valence, not companion). | |
889 | beamAPtr->append( 0, id1, x1); | |
890 | beamAPtr->xfISR( 0, id1, x1, pT2Fac * pT2Fac); | |
891 | vsc1 = beamAPtr->pickValSeaComp(); | |
892 | beamBPtr->append( 0, id2, x2); | |
893 | beamBPtr->xfISR( 0, id2, x2, pT2Fac * pT2Fac); | |
894 | vsc2 = beamBPtr->pickValSeaComp(); | |
895 | ||
896 | // Pick b according to O(b, x1, x2). | |
897 | double w1 = XDEP_A1 + a1 * log(1. / x1); | |
898 | double w2 = XDEP_A1 + a1 * log(1. / x2); | |
899 | double fac = a02now * (w1 * w1 + w2 * w2); | |
900 | double expb2 = rndmPtr->flat(); | |
901 | b2now = - fac * log(expb2); | |
902 | bNow = sqrt(b2now); | |
903 | ||
904 | // Enhancement factor for the hard process and overestimate | |
905 | // for fastPT2. Note that existing framework has a (1. / sigmaND) | |
906 | // present. | |
907 | enhanceB = sigmaND / M_PI / fac * expb2; | |
908 | enhanceBmax = sigmaND / 2. / M_PI / a02now * | |
909 | exp( -b2now / 2. / a2max ); | |
910 | ||
911 | // Trial interaction with dummy event. | |
912 | Event evDummy; | |
913 | double pTtrial = pTnext(pTmax, pTmin, evDummy); | |
914 | ||
915 | // Restore beams. | |
916 | beamAPtr->clear(); | |
917 | beamBPtr->clear(); | |
918 | ||
919 | // Accept if fallen beneath factorisation scale. | |
920 | if (pTtrial < sqrt(pT2FacSave)) { | |
921 | // Restore previous values and original sigma. | |
922 | swap(pT2, pT2Save); | |
923 | swap(pT2Fac, pT2FacSave); | |
924 | swap(pT2Ren, pT2RenSave); | |
925 | swap(id1, id1Save); | |
926 | swap(id2, id2Save); | |
927 | swap(x1, x1Save); | |
928 | swap(x2, x2Save); | |
929 | swap(sHat, sHatSave); | |
930 | swap(tHat, tHatSave); | |
931 | swap(uHat, uHatSave); | |
932 | swap(alpS, alpSsave); | |
933 | swap(alpEM, alpEMsave); | |
934 | swap(xPDF1now, xPDF1nowSave); | |
935 | swap(xPDF2now, xPDF2nowSave); | |
936 | if (dSigmaDtSel == dSigmaDtSelSave) dSigmaDtSel->swapKin(); | |
937 | else swap(dSigmaDtSel, dSigmaDtSelSave); | |
938 | ||
939 | // Accept. | |
940 | bIsSet = true; | |
941 | break; | |
942 | } | |
943 | } // if (bProfile == 4) | |
944 | } // while (true) | |
945 | ||
946 | // End handling for high b. | |
947 | } | |
948 | ||
949 | } | |
950 | ||
951 | //-------------------------------------------------------------------------- | |
952 | ||
953 | // Set up kinematics for first = hardest pT in minbias process. | |
954 | ||
955 | void MultipartonInteractions::setupFirstSys( Event& process) { | |
956 | ||
957 | // Last beam-status particles. Offset relative to normal beam locations. | |
958 | int sizeProc = process.size(); | |
959 | int nBeams = 3; | |
960 | for (int i = 3; i < sizeProc; ++i) | |
961 | if (process[i].statusAbs() < 20) nBeams = i + 1; | |
962 | int nOffset = nBeams - 3; | |
963 | ||
964 | // Remove any partons of previous failed interactions. | |
965 | if (sizeProc > nBeams) { | |
966 | process.popBack( sizeProc - nBeams); | |
967 | process.initColTag(); | |
968 | } | |
969 | ||
970 | // Entries 3 and 4, now to be added, come from 1 and 2. | |
971 | process[1 + nOffset].daughter1(3 + nOffset); | |
972 | process[2 + nOffset].daughter1(4 + nOffset); | |
973 | ||
974 | // Negate beam status, if not already done. (Case with offset beams.) | |
975 | process[1 + nOffset].statusNeg(); | |
976 | process[2 + nOffset].statusNeg(); | |
977 | ||
978 | // Loop over four partons and offset info relative to subprocess itself. | |
979 | int colOffset = process.lastColTag(); | |
980 | for (int i = 1; i <= 4; ++i) { | |
981 | Particle parton = dSigmaDtSel->getParton(i); | |
982 | if (i <= 2 ) parton.mothers( i + nOffset, 0); | |
983 | else parton.mothers( 3 + nOffset, 4 + nOffset); | |
984 | if (i <= 2 ) parton.daughters( 5 + nOffset, 6 + nOffset); | |
985 | else parton.daughters( 0, 0); | |
986 | int col = parton.col(); | |
987 | if (col > 0) parton.col( col + colOffset); | |
988 | int acol = parton.acol(); | |
989 | if (acol > 0) parton.acol( acol + colOffset); | |
990 | ||
991 | // Put the partons into the event record. | |
992 | process.append(parton); | |
993 | } | |
994 | ||
995 | // Set scale from which to begin evolution. | |
996 | process.scale( sqrt(pT2Fac) ); | |
997 | ||
998 | // Info on subprocess - specific to mimimum-bias events. | |
999 | string nameSub = dSigmaDtSel->name(); | |
1000 | int codeSub = dSigmaDtSel->code(); | |
1001 | int nFinalSub = dSigmaDtSel->nFinal(); | |
1002 | double pTMPI = dSigmaDtSel->pTMPIFin(); | |
1003 | infoPtr->setSubType( iDiffSys, nameSub, codeSub, nFinalSub); | |
1004 | if (iDiffSys == 0) infoPtr->setTypeMPI( codeSub, pTMPI, 0, 0, | |
1005 | enhanceB / zeroIntCorr); | |
1006 | ||
1007 | // Further standard info on process. | |
1008 | infoPtr->setPDFalpha( iDiffSys, id1, id2, x1, x2, xPDF1now, xPDF2now, | |
1009 | pT2Fac, alpEM, alpS, pT2Ren); | |
1010 | double m3 = dSigmaDtSel->m(3); | |
1011 | double m4 = dSigmaDtSel->m(4); | |
1012 | double theta = dSigmaDtSel->thetaMPI(); | |
1013 | double phi = dSigmaDtSel->phiMPI(); | |
1014 | infoPtr->setKin( iDiffSys, id1, id2, x1, x2, sHat, tHat, uHat, sqrt(pT2), | |
1015 | m3, m4, theta, phi); | |
1016 | ||
1017 | } | |
1018 | ||
1019 | //-------------------------------------------------------------------------- | |
1020 | ||
1021 | // Find whether to limit maximum scale of emissions. | |
1022 | ||
1023 | bool MultipartonInteractions::limitPTmax( Event& event) { | |
1024 | ||
1025 | // User-set cases. | |
1026 | if (pTmaxMatch == 1) return true; | |
1027 | if (pTmaxMatch == 2) return false; | |
1028 | ||
1029 | // Look if only quarks (u, d, s, c, b), gluons and photons in final state. | |
1030 | bool onlyQGP1 = true; | |
1031 | bool onlyQGP2 = true; | |
1032 | int n21 = 0; | |
1033 | for (int i = 5; i < event.size(); ++i) { | |
1034 | if (event[i].status() == -21) ++n21; | |
1035 | else if (n21 == 0) { | |
1036 | int idAbs = event[i].idAbs(); | |
1037 | if (idAbs > 5 && idAbs != 21 && idAbs != 22) onlyQGP1 = false; | |
1038 | } else if (n21 == 2) { | |
1039 | int idAbs = event[i].idAbs(); | |
1040 | if (idAbs > 5 && idAbs != 21 && idAbs != 22) onlyQGP2 = false; | |
1041 | } | |
1042 | } | |
1043 | ||
1044 | // If two hard interactions then limit if one only contains q/g/gamma. | |
1045 | bool onlyQGP = (n21 == 2) ? (onlyQGP1 || onlyQGP2) : onlyQGP1; | |
1046 | return (onlyQGP); | |
1047 | ||
1048 | } | |
1049 | ||
1050 | //-------------------------------------------------------------------------- | |
1051 | ||
1052 | // Select next pT in downwards evolution. | |
1053 | ||
1054 | double MultipartonInteractions::pTnext( double pTbegAll, double pTendAll, | |
1055 | Event& event) { | |
1056 | ||
1057 | // Initial values. | |
1058 | bool pickRescatter, acceptKin; | |
1059 | double dSigmaScatter, dSigmaRescatter, WTacc; | |
1060 | double pT2end = pow2( max(pTmin, pTendAll) ); | |
1061 | ||
1062 | // With the x-dependent matter profile, it is possible to reuse | |
1063 | // values that have been stored during trial interactions for a | |
1064 | // slight speedup. bIsSet is false during trial interactions, | |
1065 | // counter 21 in case partonLevel is retried and counter 22 for | |
1066 | // the first pass through partonLevel. | |
1067 | if (bProfile == 4 && bIsSet && infoPtr->getCounter(21) == 1 | |
1068 | && infoPtr->getCounter(22) == 1) { | |
1069 | ||
1070 | // Minimum bias. | |
1071 | if (bSetInFirst) { | |
1072 | if (pT2Save < pT2end) return 0.; | |
1073 | pT2 = pT2Save; | |
1074 | pT2Fac = pT2FacSave; | |
1075 | pT2Ren = pT2RenSave; | |
1076 | id1 = id1Save; | |
1077 | id2 = id2Save; | |
1078 | x1 = x1Save; | |
1079 | x2 = x2Save; | |
1080 | sHat = sHatSave; | |
1081 | tHat = tHatSave; | |
1082 | uHat = uHatSave; | |
1083 | alpS = alpSsave; | |
1084 | alpEM = alpEMsave; | |
1085 | xPDF1now = xPDF1nowSave; | |
1086 | xPDF2now = xPDF2nowSave; | |
1087 | if (dSigmaDtSel == dSigmaDtSelSave) dSigmaDtSel->swapKin(); | |
1088 | else dSigmaDtSel = dSigmaDtSelSave; | |
1089 | return sqrt(pT2); | |
1090 | ||
1091 | // Hard process. | |
1092 | } else { | |
1093 | return (pT2 < pT2end) ? 0. : sqrt(pT2); | |
1094 | } | |
1095 | } | |
1096 | ||
1097 | // Do not allow rescattering while sill FSR with global recoil. | |
1098 | bool allowRescatterNow = allowRescatter; | |
1099 | if (globalRecoilFSR && partonSystemsPtr->sizeOut(0) <= nMaxGlobalRecoilFSR) | |
1100 | allowRescatterNow = false; | |
1101 | ||
1102 | // Initial pT2 value. | |
1103 | pT2 = pow2(pTbegAll); | |
1104 | ||
1105 | // Find the set of already scattered partons on the two sides. | |
1106 | if (allowRescatterNow) findScatteredPartons( event); | |
1107 | ||
1108 | // Pick a pT2 using a quick-and-dirty cross section estimate. | |
1109 | do { | |
1110 | do { | |
1111 | pT2 = fastPT2(pT2); | |
1112 | if (pT2 < pT2end) return 0.; | |
1113 | ||
1114 | // Initial values: no rescattering. | |
1115 | i1Sel = 0; | |
1116 | i2Sel = 0; | |
1117 | dSigmaSum = 0.; | |
1118 | pickRescatter = false; | |
1119 | ||
1120 | // Pick complete kinematics and evaluate interaction cross-section. | |
1121 | dSigmaScatter = sigmaPT2scatter(false); | |
1122 | ||
1123 | // Also cross section from rescattering if allowed. | |
1124 | dSigmaRescatter = (allowRescatterNow) ? sigmaPT2rescatter( event) : 0.; | |
1125 | ||
1126 | // Normalize to dSigmaApprox, which was set in fastPT2 above. | |
1127 | WTacc = (dSigmaScatter + dSigmaRescatter) / dSigmaApprox; | |
1128 | if (WTacc > WTACCWARN) infoPtr->errorMsg("Warning in " | |
1129 | "MultipartonInteractions::pTnext: weight above unity"); | |
1130 | ||
1131 | // Idea suggested by Gosta Gustafson: increased screening in events | |
1132 | // with large activity can be simulated by pT0_eff = sqrt(n) * pT0. | |
1133 | if (enhanceScreening > 0) { | |
1134 | int nSysNow = infoPtr->nMPI() + 1; | |
1135 | if (enhanceScreening == 2) nSysNow += infoPtr->nISR(); | |
1136 | double WTscreen = pow2( (pT2 + pT20) / (pT2 + nSysNow * pT20) ); | |
1137 | WTacc *= WTscreen; | |
1138 | } | |
1139 | ||
1140 | // x-dependent matter profile overlap weighting. | |
1141 | if (bProfile == 4) { | |
1142 | double w1 = XDEP_A1 + a1 * log(1. / x1); | |
1143 | double w2 = XDEP_A1 + a1 * log(1. / x2); | |
1144 | double fac = a02now * (w1 * w1 + w2 * w2); | |
1145 | // Correct enhancement factor and weight | |
1146 | enhanceBnow = sigmaND / M_PI / fac * exp( - b2now / fac); | |
1147 | double oWgt = enhanceBnow / enhanceBmax; | |
1148 | if (oWgt > 1.0000000001) infoPtr->errorMsg("Warning in Multiparton" | |
1149 | "Interactions::pTnext: overlap weight above unity"); | |
1150 | WTacc *= oWgt; | |
1151 | } | |
1152 | ||
1153 | // Decide whether to keep the event based on weight. | |
1154 | } while (WTacc < rndmPtr->flat()); | |
1155 | ||
1156 | // When rescattering possible: new interaction or rescattering? | |
1157 | if (allowRescatterNow) { | |
1158 | pickRescatter = (i1Sel > 0 || i2Sel > 0); | |
1159 | ||
1160 | // Restore kinematics for selected scattering/rescattering. | |
1161 | id1 = id1Sel; | |
1162 | id2 = id2Sel; | |
1163 | x1 = x1Sel; | |
1164 | x2 = x2Sel; | |
1165 | sHat = sHatSel; | |
1166 | tHat = tHatSel; | |
1167 | uHat = uHatSel; | |
1168 | sigma2Sel->sigma( id1, id2, x1, x2, sHat, tHat, uHat, alpS, alpEM, | |
1169 | true, pickOtherSel); | |
1170 | } | |
1171 | ||
1172 | // Pick one of the possible channels summed above. | |
1173 | dSigmaDtSel = sigma2Sel->sigmaSel(); | |
1174 | if (sigma2Sel->swapTU()) swap( tHat, uHat); | |
1175 | ||
1176 | // Decide to keep event based on kinematics (normally always OK). | |
1177 | // Rescattering: need to provide incoming four-vectors and masses. | |
1178 | if (pickRescatter) { | |
1179 | Vec4 p1Res = (i1Sel == 0) ? 0.5 * eCM * x1Sel * Vec4( 0., 0., 1., 1.) | |
1180 | : event[i1Sel].p(); | |
1181 | Vec4 p2Res = (i2Sel == 0) ? 0.5 * eCM * x2Sel * Vec4( 0., 0., -1., 1.) | |
1182 | : event[i2Sel].p(); | |
1183 | double m1Res = (i1Sel == 0) ? 0. : event[i1Sel].m(); | |
1184 | double m2Res = (i2Sel == 0) ? 0. : event[i2Sel].m(); | |
1185 | acceptKin = dSigmaDtSel->final2KinMPI( i1Sel, i2Sel, p1Res, p2Res, | |
1186 | m1Res, m2Res); | |
1187 | // New interaction: already stored values enough. | |
1188 | } else acceptKin = dSigmaDtSel->final2KinMPI(); | |
1189 | } while (!acceptKin); | |
1190 | ||
1191 | // Done. | |
1192 | return sqrt(pT2); | |
1193 | ||
1194 | } | |
1195 | ||
1196 | //-------------------------------------------------------------------------- | |
1197 | ||
1198 | // Set up the kinematics of the 2 -> 2 scattering process, | |
1199 | // and store the scattering in the event record. | |
1200 | ||
1201 | bool MultipartonInteractions::scatter( Event& event) { | |
1202 | ||
1203 | // Last beam-status particles. Offset relative to normal beam locations. | |
1204 | int sizeProc = event.size(); | |
1205 | int nBeams = 3; | |
1206 | for (int i = 3; i < sizeProc; ++i) | |
1207 | if (event[i].statusAbs() < 20) nBeams = i + 1; | |
1208 | int nOffset = nBeams - 3; | |
1209 | ||
1210 | // Loop over four partons and offset info relative to subprocess itself. | |
1211 | int motherOffset = event.size(); | |
1212 | int colOffset = event.lastColTag(); | |
1213 | for (int i = 1; i <= 4; ++i) { | |
1214 | Particle parton = dSigmaDtSel->getParton(i); | |
1215 | if (i <= 2 ) parton.mothers( i + nOffset, 0); | |
1216 | else parton.mothers( motherOffset, motherOffset + 1); | |
1217 | if (i <= 2 ) parton.daughters( motherOffset + 2, motherOffset + 3); | |
1218 | else parton.daughters( 0, 0); | |
1219 | int col = parton.col(); | |
1220 | if (col > 0) parton.col( col + colOffset); | |
1221 | int acol = parton.acol(); | |
1222 | if (acol > 0) parton.acol( acol + colOffset); | |
1223 | ||
1224 | // Put the partons into the event record. | |
1225 | event.append(parton); | |
1226 | } | |
1227 | ||
1228 | // Allow veto of MPI. If so restore event record to before scatter. | |
1229 | if (canVetoMPI && userHooksPtr->doVetoMPIEmission(sizeProc, event)) { | |
1230 | event.popBack(event.size() - sizeProc); | |
1231 | return false; | |
1232 | } | |
1233 | ||
1234 | // Store participating partons as a new set in list of all systems. | |
1235 | int iSys = partonSystemsPtr->addSys(); | |
1236 | partonSystemsPtr->setInA(iSys, motherOffset); | |
1237 | partonSystemsPtr->setInB(iSys, motherOffset + 1); | |
1238 | partonSystemsPtr->addOut(iSys, motherOffset + 2); | |
1239 | partonSystemsPtr->addOut(iSys, motherOffset + 3); | |
1240 | partonSystemsPtr->setSHat(iSys, sHat); | |
1241 | ||
1242 | // Tag double rescattering graphs that annihilate one initial colour. | |
1243 | bool annihil1 = false; | |
1244 | bool annihil2 = false; | |
1245 | if (i1Sel > 0 && i2Sel > 0) { | |
1246 | if (event[motherOffset].col() == event[motherOffset + 1].acol() | |
1247 | && event[motherOffset].col() > 0) annihil1 = true; | |
1248 | if (event[motherOffset].acol() == event[motherOffset + 1].col() | |
1249 | && event[motherOffset].acol() > 0) annihil2 = true; | |
1250 | } | |
1251 | ||
1252 | // Beam remnant A: add scattered partons to list. | |
1253 | BeamParticle& beamA = *beamAPtr; | |
1254 | int iA = beamA.append( motherOffset, id1, x1); | |
1255 | ||
1256 | // Find whether incoming partons are valence or sea, so prepared for ISR. | |
1257 | if (i1Sel == 0) { | |
1258 | beamA.xfISR( iA, id1, x1, pT2); | |
1259 | beamA.pickValSeaComp(); | |
1260 | ||
1261 | // Remove rescattered parton from final state and change history. | |
1262 | // Propagate existing colour labels throught graph. | |
1263 | } else { | |
1264 | beamA[iA].companion(-10); | |
1265 | event[i1Sel].statusNeg(); | |
1266 | event[i1Sel].daughters( motherOffset, motherOffset); | |
1267 | event[motherOffset].mothers( i1Sel, i1Sel); | |
1268 | int colOld = event[i1Sel].col(); | |
1269 | if (colOld > 0) { | |
1270 | int colNew = event[motherOffset].col(); | |
1271 | for (int i = motherOffset; i < motherOffset + 4; ++i) { | |
1272 | if (event[i].col() == colNew) event[i].col( colOld); | |
1273 | if (event[i].acol() == colNew) event[i].acol( colOld); | |
1274 | } | |
1275 | } | |
1276 | int acolOld = event[i1Sel].acol(); | |
1277 | if (acolOld > 0) { | |
1278 | int acolNew = event[motherOffset].acol(); | |
1279 | for (int i = motherOffset; i < motherOffset + 4; ++i) { | |
1280 | if (event[i].col() == acolNew) event[i].col( acolOld); | |
1281 | if (event[i].acol() == acolNew) event[i].acol( acolOld); | |
1282 | } | |
1283 | } | |
1284 | } | |
1285 | ||
1286 | // Beam remnant B: add scattered partons to list. | |
1287 | BeamParticle& beamB = *beamBPtr; | |
1288 | int iB = beamB.append( motherOffset + 1, id2, x2); | |
1289 | ||
1290 | // Find whether incoming partons are valence or sea, so prepared for ISR. | |
1291 | if (i2Sel == 0) { | |
1292 | beamB.xfISR( iB, id2, x2, pT2); | |
1293 | beamB.pickValSeaComp(); | |
1294 | ||
1295 | // Remove rescattered parton from final state and change history. | |
1296 | // Propagate existing colour labels throught graph. | |
1297 | } else { | |
1298 | beamB[iB].companion(-10); | |
1299 | event[i2Sel].statusNeg(); | |
1300 | event[i2Sel].daughters( motherOffset + 1, motherOffset + 1); | |
1301 | event[motherOffset + 1].mothers( i2Sel, i2Sel); | |
1302 | int colOld = event[i2Sel].col(); | |
1303 | if (colOld > 0) { | |
1304 | int colNew = event[motherOffset + 1].col(); | |
1305 | for (int i = motherOffset; i < motherOffset + 4; ++i) { | |
1306 | if (event[i].col() == colNew) event[i].col( colOld); | |
1307 | if (event[i].acol() == colNew) event[i].acol( colOld); | |
1308 | } | |
1309 | } | |
1310 | int acolOld = event[i2Sel].acol(); | |
1311 | if (acolOld > 0) { | |
1312 | int acolNew = event[motherOffset + 1].acol(); | |
1313 | for (int i = motherOffset; i < motherOffset + 4; ++i) { | |
1314 | if (event[i].col() == acolNew) event[i].col( acolOld); | |
1315 | if (event[i].acol() == acolNew) event[i].acol( acolOld); | |
1316 | } | |
1317 | } | |
1318 | } | |
1319 | ||
1320 | // Annihilating colour in double rescattering requires relabelling | |
1321 | // of one colour into the other in the whole preceding event. | |
1322 | if (annihil1 || annihil2) { | |
1323 | int colLeft = (annihil1) ? event[motherOffset].col() | |
1324 | : event[motherOffset].acol(); | |
1325 | int mother1 = event[motherOffset].mother1(); | |
1326 | int mother2 = event[motherOffset + 1].mother1(); | |
1327 | int colLost = (annihil1) | |
1328 | ? event[mother1].col() + event[mother2].acol() - colLeft | |
1329 | : event[mother1].acol() + event[mother2].col() - colLeft; | |
1330 | for (int i = 0; i < motherOffset; ++i) { | |
1331 | if (event[i].col() == colLost) event[i].col( colLeft ); | |
1332 | if (event[i].acol() == colLost) event[i].acol( colLeft ); | |
1333 | } | |
1334 | } | |
1335 | ||
1336 | // Store info on subprocess code and rescattered partons. | |
1337 | int codeMPI = dSigmaDtSel->code(); | |
1338 | double pTMPI = dSigmaDtSel->pTMPIFin(); | |
1339 | if (iDiffSys == 0) infoPtr->setTypeMPI( codeMPI, pTMPI, i1Sel, i2Sel, | |
1340 | enhanceBnow / zeroIntCorr); | |
1341 | partonSystemsPtr->setPTHat(iSys, pTMPI); | |
1342 | ||
1343 | // Done. | |
1344 | return true; | |
1345 | } | |
1346 | ||
1347 | //-------------------------------------------------------------------------- | |
1348 | ||
1349 | // Determine constant in d(Prob)/d(pT2) < const / (pT2 + r * pT20)^2. | |
1350 | ||
1351 | void MultipartonInteractions::upperEnvelope() { | |
1352 | ||
1353 | // Initially determine constant in jet cross section upper estimate | |
1354 | // d(sigma_approx)/d(pT2) < const / (pT2 + r * pT20)^2. | |
1355 | pT4dSigmaMax = 0.; | |
1356 | ||
1357 | // Loop thorough allowed pT range logarithmically evenly. | |
1358 | for (int iPT = 0; iPT < 100; ++iPT) { | |
1359 | double pT = pTmin * pow( pTmax/pTmin, 0.01 * (iPT + 0.5) ); | |
1360 | pT2 = pT*pT; | |
1361 | pT2shift = pT2 + pT20; | |
1362 | pT2Ren = pT2shift; | |
1363 | pT2Fac = (SHIFTFACSCALE) ? pT2shift : pT2; | |
1364 | xT = 2. * pT / eCM; | |
1365 | ||
1366 | // Evaluate parton density sums at x1 = x2 = xT. | |
1367 | double xPDF1sumMax = (9./4.) * beamAPtr->xf(21, xT, pT2Fac); | |
1368 | for (int id = 1; id <= nQuarkIn; ++id) | |
1369 | xPDF1sumMax += beamAPtr->xf( id, xT, pT2Fac) | |
1370 | + beamAPtr->xf(-id, xT, pT2Fac); | |
1371 | double xPDF2sumMax = (9./4.) * beamBPtr->xf(21, xT, pT2Fac); | |
1372 | for (int id = 1; id <= nQuarkIn; ++id) | |
1373 | xPDF2sumMax += beamBPtr->xf( id, xT, pT2Fac) | |
1374 | + beamBPtr->xf(-id, xT, pT2Fac); | |
1375 | ||
1376 | // Evaluate alpha_strong and _EM, matrix element and phase space volume. | |
1377 | alpS = alphaS.alphaS(pT2Ren); | |
1378 | alpEM = alphaEM.alphaEM(pT2Ren); | |
1379 | double dSigmaPartonApprox = CONVERT2MB * Kfactor * 0.5 * M_PI | |
1380 | * pow2(alpS / pT2shift); | |
1381 | double yMax = log(1./xT + sqrt(1./(xT*xT) - 1.)); | |
1382 | double volumePhSp = pow2(2. * yMax); | |
1383 | ||
1384 | // Final comparison to determine upper estimate. | |
1385 | double dSigmaApproxNow = SIGMAFUDGE * xPDF1sumMax * xPDF2sumMax | |
1386 | * dSigmaPartonApprox * volumePhSp; | |
1387 | double pT4dSigmaNow = pow2(pT2 + pT20R) * dSigmaApproxNow; | |
1388 | if ( pT4dSigmaNow > pT4dSigmaMax) pT4dSigmaMax = pT4dSigmaNow; | |
1389 | } | |
1390 | ||
1391 | // Get wanted constant by dividing by the nondiffractive cross section. | |
1392 | pT4dProbMax = pT4dSigmaMax / sigmaND; | |
1393 | ||
1394 | } | |
1395 | ||
1396 | //-------------------------------------------------------------------------- | |
1397 | ||
1398 | // Integrate the parton-parton interaction cross section, | |
1399 | // using stratified Monte Carlo sampling. | |
1400 | // Store result in pT bins for use as Sudakov form factors. | |
1401 | ||
1402 | void MultipartonInteractions::jetCrossSection() { | |
1403 | ||
1404 | // Common factor from bin size in dpT2 / (pT2 + r * pT20)^2 and statistics. | |
1405 | double sigmaFactor = (1. / pT20minR - 1. / pT20maxR) / (100. * nSample); | |
1406 | ||
1407 | // Reset overlap-weighted cross section for x-dependent matter profile. | |
1408 | if (bProfile == 4) for (int bBin = 0; bBin < XDEP_BBIN; bBin++) | |
1409 | sigmaIntWgt[bBin] = 0.; | |
1410 | ||
1411 | // Loop through allowed pT range evenly in dpT2 / (pT2 + r * pT20)^2. | |
1412 | sigmaInt = 0.; | |
1413 | double dSigmaMax = 0.; | |
1414 | sudExpPT[100] = 0.; | |
1415 | ||
1416 | for (int iPT = 99; iPT >= 0; --iPT) { | |
1417 | double sigmaSum = 0.; | |
1418 | ||
1419 | // Reset pT-binned overlap-weigted integration. | |
1420 | if (bProfile == 4) for (int bBin = 0; bBin < XDEP_BBIN; bBin++) | |
1421 | sigmaSumWgt[bBin] = 0.; | |
1422 | ||
1423 | // In each pT bin sample a number of random pT values. | |
1424 | for (int iSample = 0; iSample < nSample; ++iSample) { | |
1425 | double mappedPT2 = 1. - 0.01 * (iPT + rndmPtr->flat()); | |
1426 | pT2 = pT20min0maxR / (pT20minR + mappedPT2 * pT2maxmin) - pT20R; | |
1427 | ||
1428 | // Evaluate cross section dSigma/dpT2 in phase space point. | |
1429 | double dSigma = sigmaPT2scatter(true); | |
1430 | ||
1431 | // Multiply by (pT2 + r * pT20)^2 to compensate for pT sampling. Sum. | |
1432 | dSigma *= pow2(pT2 + pT20R); | |
1433 | sigmaSum += dSigma; | |
1434 | if (dSigma > dSigmaMax) dSigmaMax = dSigma; | |
1435 | ||
1436 | // Overlap-weighted cross section for x-dependent matter profile. | |
1437 | // Note that dSigma can be 0. when points are rejected. | |
1438 | if (bProfile == 4 && dSigma > 0.) { | |
1439 | double w1 = XDEP_A1 + a1 * log(1. / x1); | |
1440 | double w2 = XDEP_A1 + a1 * log(1. / x2); | |
1441 | double fac = XDEP_A0 * XDEP_A0 * (w1 * w1 + w2 * w2); | |
1442 | double b = 0.5 * bstepNow; | |
1443 | for (int bBin = 0; bBin < XDEP_BBIN; bBin++) { | |
1444 | double wgt = exp( - b * b / fac ) / fac / M_PI; | |
1445 | sigmaSumWgt[bBin] += dSigma * wgt; | |
1446 | b += bstepNow; | |
1447 | } | |
1448 | } | |
1449 | } | |
1450 | ||
1451 | // Store total cross section and exponent of Sudakov. | |
1452 | sigmaSum *= sigmaFactor; | |
1453 | sigmaInt += sigmaSum; | |
1454 | sudExpPT[iPT] = sudExpPT[iPT + 1] + sigmaSum / sigmaND; | |
1455 | ||
1456 | // Sum overlap-weighted cross section | |
1457 | if (bProfile == 4) for (int bBin = 0; bBin < XDEP_BBIN; bBin++) { | |
1458 | sigmaSumWgt[bBin] *= sigmaFactor; | |
1459 | sigmaIntWgt[bBin] += sigmaSumWgt[bBin]; | |
1460 | } | |
1461 | ||
1462 | // End of loop over pT values. | |
1463 | } | |
1464 | ||
1465 | // Update upper estimate of differential cross section. Done. | |
1466 | if (dSigmaMax > pT4dSigmaMax) { | |
1467 | pT4dSigmaMax = dSigmaMax; | |
1468 | pT4dProbMax = dSigmaMax / sigmaND; | |
1469 | } | |
1470 | ||
1471 | } | |
1472 | ||
1473 | //-------------------------------------------------------------------------- | |
1474 | ||
1475 | // Evaluate "Sudakov form factor" for not having a harder interaction | |
1476 | // at the selected b value, given the pT scale of the event. | |
1477 | ||
1478 | double MultipartonInteractions::sudakov(double pT2sud, double enhance) { | |
1479 | ||
1480 | // Find bin the pT2 scale falls in. | |
1481 | double xBin = (pT2sud - pT2min) * pT20maxR | |
1482 | / (pT2maxmin * (pT2sud + pT20R)); | |
1483 | xBin = max(1e-6, min(100. - 1e-6, 100. * xBin) ); | |
1484 | int iBin = int(xBin); | |
1485 | ||
1486 | // Interpolate inside bin. Optionally include enhancement factor. | |
1487 | double sudExp = sudExpPT[iBin] + (xBin - iBin) | |
1488 | * (sudExpPT[iBin + 1] - sudExpPT[iBin]); | |
1489 | return exp( -enhance * sudExp); | |
1490 | ||
1491 | } | |
1492 | ||
1493 | //-------------------------------------------------------------------------- | |
1494 | ||
1495 | // Pick a trial next pT, based on a simple upper estimate of the | |
1496 | // d(sigma)/d(pT2) spectrum. | |
1497 | ||
1498 | double MultipartonInteractions::fastPT2( double pT2beg) { | |
1499 | ||
1500 | // Use d(Prob)/d(pT2) < pT4dProbMax / (pT2 + r * pT20)^2. | |
1501 | double pT20begR = pT2beg + pT20R; | |
1502 | double pT4dProbMaxNow = pT4dProbMax * enhanceBmax; | |
1503 | double pT2try = pT4dProbMaxNow * pT20begR | |
1504 | / (pT4dProbMaxNow - pT20begR * log(rndmPtr->flat())) - pT20R; | |
1505 | ||
1506 | // Save cross section associated with ansatz above. Done. | |
1507 | dSigmaApprox = pT4dSigmaMax / pow2(pT2try + pT20R); | |
1508 | return pT2try; | |
1509 | ||
1510 | } | |
1511 | ||
1512 | //-------------------------------------------------------------------------- | |
1513 | ||
1514 | // Calculate the actual cross section to decide whether fast choice is OK. | |
1515 | // Select flavours and kinematics for interaction at given pT. | |
1516 | // Slightly different treatment for first interaction and subsequent ones. | |
1517 | ||
1518 | double MultipartonInteractions::sigmaPT2scatter(bool isFirst) { | |
1519 | ||
1520 | // Derive recormalization and factorization scales, amd alpha_strong/em. | |
1521 | pT2shift = pT2 + pT20; | |
1522 | pT2Ren = pT2shift; | |
1523 | pT2Fac = (SHIFTFACSCALE) ? pT2shift : pT2; | |
1524 | alpS = alphaS.alphaS(pT2Ren); | |
1525 | alpEM = alphaEM.alphaEM(pT2Ren); | |
1526 | ||
1527 | // Derive rapidity limits from chosen pT2. | |
1528 | xT = 2. * sqrt(pT2) / eCM; | |
1529 | if (xT >= 1.) return 0.; | |
1530 | xT2 = xT*xT; | |
1531 | double yMax = log(1./xT + sqrt(1./xT2 - 1.)); | |
1532 | ||
1533 | // Select rapidities y3 and y4 of the two produced partons. | |
1534 | double y3 = yMax * (2. * rndmPtr->flat() - 1.); | |
1535 | double y4 = yMax * (2. * rndmPtr->flat() - 1.); | |
1536 | y = 0.5 * (y3 + y4); | |
1537 | ||
1538 | // Reject some events at large rapidities to improve efficiency. | |
1539 | // (Works for baryons, not pions or Pomerons if they have hard PDF's.) | |
1540 | double WTy = (hasBaryonBeams) | |
1541 | ? (1. - pow2(y3/yMax)) * (1. - pow2(y4/yMax)) : 1.; | |
1542 | if (WTy < rndmPtr->flat()) return 0.; | |
1543 | ||
1544 | // Failure if x1 or x2 exceed what is left in respective beam. | |
1545 | x1 = 0.5 * xT * (exp(y3) + exp(y4)); | |
1546 | x2 = 0.5 * xT * (exp(-y3) + exp(-y4)); | |
1547 | if (isFirst) { | |
1548 | if (x1 > 1. || x2 > 1.) return 0.; | |
1549 | } else { | |
1550 | if (x1 > beamAPtr->xMax() || x2 > beamBPtr->xMax()) return 0.; | |
1551 | } | |
1552 | tau = x1 * x2; | |
1553 | ||
1554 | // Begin evaluate parton densities at actual x1 and x2. | |
07530286 | 1555 | double xPDF1[21] = {0.}; |
c6b60c38 | 1556 | double xPDF1sum = 0.; |
07530286 | 1557 | double xPDF2[21] = {0.}; |
c6b60c38 | 1558 | double xPDF2sum = 0.; |
1559 | ||
1560 | // For first interaction use normal densities. | |
1561 | if (isFirst) { | |
1562 | for (int id = -nQuarkIn; id <= nQuarkIn; ++id) { | |
1563 | if (id == 0) xPDF1[10] = (9./4.) * beamAPtr->xf(21, x1, pT2Fac); | |
1564 | else xPDF1[id+10] = beamAPtr->xf(id, x1, pT2Fac); | |
1565 | xPDF1sum += xPDF1[id+10]; | |
1566 | } | |
1567 | for (int id = -nQuarkIn; id <= nQuarkIn; ++id) { | |
1568 | if (id == 0) xPDF2[10] = (9./4.) * beamBPtr->xf(21, x2, pT2Fac); | |
1569 | else xPDF2[id+10] = beamBPtr->xf(id, x2, pT2Fac); | |
1570 | xPDF2sum += xPDF2[id+10]; | |
1571 | } | |
1572 | ||
1573 | // For subsequent interactions use rescaled densities. | |
1574 | } else { | |
1575 | for (int id = -nQuarkIn; id <= nQuarkIn; ++id) { | |
1576 | if (id == 0) xPDF1[10] = (9./4.) * beamAPtr->xfMPI(21, x1, pT2Fac); | |
1577 | else xPDF1[id+10] = beamAPtr->xfMPI(id, x1, pT2Fac); | |
1578 | xPDF1sum += xPDF1[id+10]; | |
1579 | } | |
1580 | for (int id = -nQuarkIn; id <= nQuarkIn; ++id) { | |
1581 | if (id == 0) xPDF2[10] = (9./4.) * beamBPtr->xfMPI(21, x2, pT2Fac); | |
1582 | else xPDF2[id+10] = beamBPtr->xfMPI(id, x2, pT2Fac); | |
1583 | xPDF2sum += xPDF2[id+10]; | |
1584 | } | |
1585 | } | |
1586 | ||
1587 | // Select incoming flavours according to actual PDF's. | |
1588 | id1 = -nQuarkIn - 1; | |
1589 | double temp = xPDF1sum * rndmPtr->flat(); | |
1590 | do { xPDF1now = xPDF1[(++id1) + 10]; temp -= xPDF1now; } | |
1591 | while (temp > 0. && id1 < nQuarkIn); | |
1592 | if (id1 == 0) id1 = 21; | |
1593 | id2 = -nQuarkIn-1; | |
1594 | temp = xPDF2sum * rndmPtr->flat(); | |
1595 | do { xPDF2now = xPDF2[(++id2) + 10]; temp -= xPDF2now;} | |
1596 | while (temp > 0. && id2 < nQuarkIn); | |
1597 | if (id2 == 0) id2 = 21; | |
1598 | ||
1599 | // Assign pointers to processes relevant for incoming flavour choice: | |
1600 | // g + g, q + g, q + qbar (same flavour), q + q(bar) (the rest). | |
1601 | // Factor 4./9. per incoming gluon to compensate for preweighting. | |
1602 | SigmaMultiparton* sigma2Tmp; | |
1603 | double gluFac = 1.; | |
1604 | if (id1 == 21 && id2 == 21) { | |
1605 | sigma2Tmp = &sigma2gg; | |
1606 | gluFac = 16. / 81.; | |
1607 | } else if (id1 == 21 || id2 == 21) { | |
1608 | sigma2Tmp = &sigma2qg; | |
1609 | gluFac = 4. / 9.; | |
1610 | } else if (id1 == -id2) sigma2Tmp = &sigma2qqbarSame; | |
1611 | else sigma2Tmp = &sigma2qq; | |
1612 | ||
1613 | // Prepare to generate differential cross sections. | |
1614 | sHat = tau * sCM; | |
1615 | double root = sqrtpos(1. - xT2 / tau); | |
1616 | tHat = -0.5 * sHat * (1. - root); | |
1617 | uHat = -0.5 * sHat * (1. + root); | |
1618 | ||
1619 | // Evaluate cross sections, include possibility of K factor. | |
1620 | // Use kinematics for two symmetrical configurations (tHat <-> uHat). | |
1621 | // (Not necessary, but makes for better MC efficiency.) | |
1622 | double dSigmaPartonCorr = Kfactor * gluFac | |
1623 | * sigma2Tmp->sigma( id1, id2, x1, x2, sHat, tHat, uHat, alpS, alpEM); | |
1624 | ||
1625 | // Combine cross section, pdf's and phase space integral. | |
1626 | double volumePhSp = pow2(2. * yMax) / WTy; | |
1627 | double dSigmaScat = dSigmaPartonCorr * xPDF1sum * xPDF2sum * volumePhSp; | |
1628 | ||
1629 | // Dampen cross section at small pT values; part of formalism. | |
1630 | // Use pT2 corrected for massive kinematics at this step.?? | |
1631 | // double pT2massive = dSigmaDtSel->pT2MPI(); | |
1632 | double pT2massive = pT2; | |
1633 | dSigmaScat *= pow2( pT2massive / (pT2massive + pT20) ); | |
1634 | ||
1635 | // Sum up total contribution for all scatterings and rescatterings. | |
1636 | dSigmaSum += dSigmaScat; | |
1637 | ||
1638 | // Save values for comparison with rescattering processes. | |
1639 | i1Sel = 0; | |
1640 | i2Sel = 0; | |
1641 | id1Sel = id1; | |
1642 | id2Sel = id2; | |
1643 | x1Sel = x1; | |
1644 | x2Sel = x2; | |
1645 | sHatSel = sHat; | |
1646 | tHatSel = tHat; | |
1647 | uHatSel = uHat; | |
1648 | sigma2Sel = sigma2Tmp; | |
1649 | pickOtherSel = sigma2Tmp->pickedOther(); | |
1650 | ||
1651 | // For first interaction: pick one of the possible channels summed above. | |
1652 | if (isFirst) { | |
1653 | dSigmaDtSel = sigma2Tmp->sigmaSel(); | |
1654 | if (sigma2Tmp->swapTU()) swap( tHat, uHat); | |
1655 | } | |
1656 | ||
1657 | // Done. | |
1658 | return dSigmaScat; | |
1659 | } | |
1660 | ||
1661 | //-------------------------------------------------------------------------- | |
1662 | ||
1663 | // Find the partons that are allowed to rescatter. | |
1664 | ||
1665 | void MultipartonInteractions::findScatteredPartons( Event& event) { | |
1666 | ||
1667 | // Reset arrays. | |
1668 | scatteredA.resize(0); | |
1669 | scatteredB.resize(0); | |
1670 | double yTmp, probA, probB; | |
1671 | ||
1672 | // Loop though the event record and catch "final" partons. | |
1673 | for (int i = 0; i < event.size(); ++i) | |
1674 | if (event[i].isFinal() && (event[i].idAbs() <= nQuarkIn | |
1675 | || event[i].id() == 21)) { | |
1676 | yTmp = event[i].y(); | |
1677 | ||
1678 | // Different strategies to determine which partons may rescatter. | |
1679 | switch(rescatterMode) { | |
1680 | ||
1681 | // Case 0: step function at origin | |
1682 | case 0: | |
1683 | if ( yTmp > 0.) scatteredA.push_back( i); | |
1684 | if (-yTmp > 0.) scatteredB.push_back( i); | |
1685 | break; | |
1686 | ||
1687 | // Case 1: step function as ySepResc. | |
1688 | case 1: | |
1689 | if ( yTmp > ySepResc) scatteredA.push_back( i); | |
1690 | if (-yTmp > ySepResc) scatteredB.push_back( i); | |
1691 | break; | |
1692 | ||
1693 | // Case 2: linear rise from ySep - deltaY to ySep + deltaY. | |
1694 | case 2: | |
1695 | probA = 0.5 * (1. + ( yTmp - ySepResc) / deltaYResc); | |
1696 | if (probA > rndmPtr->flat()) scatteredA.push_back( i); | |
1697 | probB = 0.5 * (1. + (-yTmp - ySepResc) / deltaYResc); | |
1698 | if (probB > rndmPtr->flat()) scatteredB.push_back( i); | |
1699 | break; | |
1700 | ||
1701 | // Case 3: rise like (1/2) * ( 1 + tanh((y - ySep) / deltaY) ). | |
1702 | // Use that (1/2) (1 + tanh(x)) = 1 / (1 + exp(-2x)). | |
1703 | case 3: | |
1704 | probA = 1. / (1. + exp(-2. * ( yTmp - ySepResc) / deltaYResc)); | |
1705 | if (probA > rndmPtr->flat()) scatteredA.push_back( i); | |
1706 | probB = 1. / (1. + exp(-2. * (-yTmp - ySepResc) / deltaYResc)); | |
1707 | if (probB > rndmPtr->flat()) scatteredB.push_back( i); | |
1708 | break; | |
1709 | ||
1710 | // Case 4 and undefined values: all partons can rescatter. | |
1711 | default: | |
1712 | scatteredA.push_back( i); | |
1713 | scatteredB.push_back( i); | |
1714 | break; | |
1715 | ||
1716 | // End of loop over partons. Done. | |
1717 | } | |
1718 | } | |
1719 | ||
1720 | } | |
1721 | ||
1722 | //-------------------------------------------------------------------------- | |
1723 | ||
1724 | // Rescattering contribution summed over all already scattered partons. | |
1725 | // Calculate the actual cross section to decide whether fast choice is OK. | |
1726 | // Select flavours and kinematics for interaction at given pT. | |
1727 | ||
1728 | double MultipartonInteractions::sigmaPT2rescatter( Event& event) { | |
1729 | ||
1730 | // Derive xT scale from chosen pT2. | |
1731 | xT = 2. * sqrt(pT2) / eCM; | |
1732 | if (xT >= 1.) return 0.; | |
1733 | xT2 = xT*xT; | |
1734 | ||
1735 | // Pointer to cross section and total rescatter contribution. | |
1736 | SigmaMultiparton* sigma2Tmp; | |
1737 | double dSigmaResc = 0.; | |
1738 | ||
1739 | // Normally save time by picking one random scattered parton from side A. | |
1740 | int nScatA = scatteredA.size(); | |
1741 | int iScatA = -1; | |
1742 | if (PREPICKRESCATTER && nScatA > 0) iScatA = min( nScatA - 1, | |
1743 | int( rndmPtr->flat() * double(nScatA) ) ); | |
1744 | ||
1745 | // Loop over all already scattered partons from side A. | |
1746 | for (int iScat = 0; iScat < nScatA; ++iScat) { | |
1747 | if (PREPICKRESCATTER && iScat != iScatA) continue; | |
1748 | int iA = scatteredA[iScat]; | |
1749 | int id1Tmp = event[iA].id(); | |
1750 | ||
1751 | // Calculate maximum possible sHat and check whether big enough. | |
1752 | double x1Tmp = event[iA].pPos() / eCM; | |
1753 | double sHatMax = x1Tmp * beamBPtr->xMax() * sCM; | |
1754 | if (sHatMax < 4. * pT2) continue; | |
1755 | ||
1756 | // Select rapidity separation between two produced partons. | |
1757 | double dyMax = log( sqrt(0.25 * sHatMax / pT2) | |
1758 | * (1. + sqrtpos(1. - 4. * pT2 / sHatMax)) ); | |
1759 | double dy = dyMax * (2. * rndmPtr->flat() - 1.); | |
1760 | ||
1761 | // Reconstruct kinematical variables, especially x2. | |
1762 | // Incoming c/b masses handled better if tau != x1 * x2. | |
1763 | double m2Tmp = event[iA].m2(); | |
1764 | double sHatTmp = m2Tmp + 4. * pT2 * pow2(cosh(dy)); | |
1765 | double x2Tmp = (sHatTmp - m2Tmp) /(x1Tmp * sCM); | |
1766 | double tauTmp = sHatTmp / sCM; | |
1767 | double root = sqrtpos(1. - xT2 / tauTmp); | |
1768 | double tHatTmp = -0.5 * sHatTmp * (1. - root); | |
1769 | double uHatTmp = -0.5 * sHatTmp * (1. + root); | |
1770 | ||
1771 | // Begin evaluate parton densities at actual x2. | |
e74d4fa4 | 1772 | double xPDF2[21] = {0.}; |
c6b60c38 | 1773 | double xPDF2sum = 0.; |
1774 | ||
1775 | // Use rescaled densities, with preweighting 9/4 for gluons. | |
1776 | for (int id = -nQuarkIn; id <= nQuarkIn; ++id) { | |
1777 | if (id == 0) xPDF2[10] = (9./4.) * beamBPtr->xfMPI(21, x2Tmp, pT2Fac); | |
1778 | else xPDF2[id+10] = beamBPtr->xfMPI(id, x2Tmp, pT2Fac); | |
1779 | xPDF2sum += xPDF2[id+10]; | |
1780 | } | |
1781 | ||
1782 | // Select incoming flavour according to actual PDF's. | |
1783 | int id2Tmp = -nQuarkIn - 1; | |
1784 | double temp = xPDF2sum * rndmPtr->flat(); | |
1785 | do { xPDF2now = xPDF2[(++id2Tmp) + 10]; temp -= xPDF2now;} | |
1786 | while (temp > 0. && id2Tmp < nQuarkIn); | |
1787 | if (id2Tmp == 0) id2Tmp = 21; | |
1788 | ||
1789 | // Assign pointers to processes relevant for incoming flavour choice: | |
1790 | // g + g, q + g, q + qbar (same flavour), q + q(bar) (the rest). | |
1791 | // Factor 4./9. for incoming gluon 2 to compensate for preweighting. | |
1792 | if (id1Tmp == 21 && id2Tmp == 21) sigma2Tmp = &sigma2gg; | |
1793 | else if (id1Tmp == 21 || id2Tmp == 21) sigma2Tmp = &sigma2qg; | |
1794 | else if (id1Tmp == -id2Tmp) sigma2Tmp = &sigma2qqbarSame; | |
1795 | else sigma2Tmp = &sigma2qq; | |
1796 | double gluFac = (id2Tmp == 21) ? 4. / 9. : 1.; | |
1797 | ||
1798 | // Evaluate cross sections, include possibility of K factor. | |
1799 | // Use kinematics for two symmetrical configurations (tHat <-> uHat). | |
1800 | // (Not necessary, but makes for better MC efficiency.) | |
1801 | double dSigmaPartonCorr = Kfactor * gluFac | |
1802 | * sigma2Tmp->sigma( id1Tmp, id2Tmp, x1Tmp, x2Tmp, sHatTmp, tHatTmp, | |
1803 | uHatTmp, alpS, alpEM); | |
1804 | ||
1805 | // Combine cross section, pdf's and phase space integral. | |
1806 | double volumePhSp = 4. *dyMax; | |
1807 | double dSigmaCorr = dSigmaPartonCorr * xPDF2sum * volumePhSp; | |
1808 | ||
1809 | // Dampen cross section at small pT values; part of formalism. | |
1810 | // Use pT2 corrected for massive kinematics at this step. | |
1811 | //?? double pT2massive = dSigmaDtSel->pT2MPI(); | |
1812 | double pT2massive = pT2; | |
1813 | dSigmaCorr *= pow2( pT2massive / (pT2massive + pT20) ); | |
1814 | ||
1815 | // Compensate for prepicked rescattering if appropriate. | |
1816 | if (PREPICKRESCATTER) dSigmaCorr *= nScatA; | |
1817 | ||
1818 | // Sum up total contribution for all scatterings or rescattering only. | |
1819 | dSigmaSum += dSigmaCorr; | |
1820 | dSigmaResc += dSigmaCorr; | |
1821 | ||
1822 | // Determine whether current rescattering should be selected. | |
1823 | if (dSigmaCorr > rndmPtr->flat() * dSigmaSum) { | |
1824 | i1Sel = iA; | |
1825 | i2Sel = 0; | |
1826 | id1Sel = id1Tmp; | |
1827 | id2Sel = id2Tmp; | |
1828 | x1Sel = x1Tmp; | |
1829 | x2Sel = x2Tmp; | |
1830 | sHatSel = sHatTmp; | |
1831 | tHatSel = tHatTmp; | |
1832 | uHatSel = uHatTmp; | |
1833 | sigma2Sel = sigma2Tmp; | |
1834 | pickOtherSel = sigma2Tmp->pickedOther(); | |
1835 | } | |
1836 | } | |
1837 | ||
1838 | // Normally save time by picking one random scattered parton from side B. | |
1839 | int nScatB = scatteredB.size(); | |
1840 | int iScatB = -1; | |
1841 | if (PREPICKRESCATTER && nScatB > 0) iScatB = min( nScatB - 1, | |
1842 | int( rndmPtr->flat() * double(nScatB) ) ); | |
1843 | ||
1844 | // Loop over all already scattered partons from side B. | |
1845 | for (int iScat = 0; iScat < nScatB; ++iScat) { | |
1846 | if (PREPICKRESCATTER && iScat != iScatB) continue; | |
1847 | int iB = scatteredB[iScat]; | |
1848 | int id2Tmp = event[iB].id(); | |
1849 | ||
1850 | // Calculate maximum possible sHat and check whether big enough. | |
1851 | double x2Tmp = event[iB].pNeg() / eCM; | |
1852 | double sHatMax = beamAPtr->xMax() * x2Tmp * sCM; | |
1853 | if (sHatMax < 4. * pT2) continue; | |
1854 | ||
1855 | // Select rapidity separation between two produced partons. | |
1856 | double dyMax = log( sqrt(0.25 * sHatMax / pT2) | |
1857 | * (1. + sqrtpos(1. - 4. * pT2 / sHatMax)) ); | |
1858 | double dy = dyMax * (2. * rndmPtr->flat() - 1.); | |
1859 | ||
1860 | // Reconstruct kinematical variables, especially x1. | |
1861 | // Incoming c/b masses handled better if tau != x1 * x2. | |
1862 | double m2Tmp = event[iB].m2(); | |
1863 | double sHatTmp = m2Tmp + 4. * pT2 * pow2(cosh(dy)); | |
1864 | double x1Tmp = (sHatTmp - m2Tmp) /(x2Tmp * sCM); | |
1865 | double tauTmp = sHatTmp / sCM; | |
1866 | double root = sqrtpos(1. - xT2 / tauTmp); | |
1867 | double tHatTmp = -0.5 * sHatTmp * (1. - root); | |
1868 | double uHatTmp = -0.5 * sHatTmp * (1. + root); | |
1869 | ||
1870 | // Begin evaluate parton densities at actual x1. | |
d099bed0 | 1871 | double xPDF1[21] = {0}; |
c6b60c38 | 1872 | double xPDF1sum = 0.; |
1873 | ||
1874 | // Use rescaled densities, with preweighting 9/4 for gluons. | |
1875 | for (int id = -nQuarkIn; id <= nQuarkIn; ++id) { | |
1876 | if (id == 0) xPDF1[10] = (9./4.) * beamAPtr->xfMPI(21, x1Tmp, pT2Fac); | |
1877 | else xPDF1[id+10] = beamAPtr->xfMPI(id, x1Tmp, pT2Fac); | |
1878 | xPDF1sum += xPDF1[id+10]; | |
1879 | } | |
1880 | ||
1881 | // Select incoming flavour according to actual PDF's. | |
1882 | int id1Tmp = -nQuarkIn - 1; | |
1883 | double temp = xPDF1sum * rndmPtr->flat(); | |
1884 | do { xPDF1now = xPDF1[(++id1Tmp) + 10]; temp -= xPDF1now;} | |
1885 | while (temp > 0. && id1Tmp < nQuarkIn); | |
1886 | if (id1Tmp == 0) id1Tmp = 21; | |
1887 | ||
1888 | // Assign pointers to processes relevant for incoming flavour choice: | |
1889 | // g + g, q + g, q + qbar (same flavour), q + q(bar) (the rest). | |
1890 | // Factor 4./9. for incoming gluon 2 to compensate for preweighting. | |
1891 | if (id1Tmp == 21 && id2Tmp == 21) sigma2Tmp = &sigma2gg; | |
1892 | else if (id1Tmp == 21 || id2Tmp == 21) sigma2Tmp = &sigma2qg; | |
1893 | else if (id1Tmp == -id2Tmp) sigma2Tmp = &sigma2qqbarSame; | |
1894 | else sigma2Tmp = &sigma2qq; | |
1895 | double gluFac = (id1Tmp == 21) ? 4. / 9. : 1.; | |
1896 | ||
1897 | // Evaluate cross sections, include possibility of K factor. | |
1898 | // Use kinematics for two symmetrical configurations (tHat <-> uHat). | |
1899 | // (Not necessary, but makes for better MC efficiency.) | |
1900 | double dSigmaPartonCorr = Kfactor * gluFac | |
1901 | * sigma2Tmp->sigma( id1Tmp, id2Tmp, x1Tmp, x2Tmp, sHatTmp, tHatTmp, | |
1902 | uHatTmp, alpS, alpEM); | |
1903 | ||
1904 | // Combine cross section, pdf's and phase space integral. | |
1905 | double volumePhSp = 4. *dyMax; | |
1906 | double dSigmaCorr = dSigmaPartonCorr * xPDF1sum * volumePhSp; | |
1907 | ||
1908 | // Dampen cross section at small pT values; part of formalism. | |
1909 | // Use pT2 corrected for massive kinematics at this step. | |
1910 | //?? double pT2massive = dSigmaDtSel->pT2MPI(); | |
1911 | double pT2massive = pT2; | |
1912 | dSigmaCorr *= pow2( pT2massive / (pT2massive + pT20) ); | |
1913 | ||
1914 | // Compensate for prepicked rescattering if appropriate. | |
1915 | if (PREPICKRESCATTER) dSigmaCorr *= nScatB; | |
1916 | ||
1917 | // Sum up total contribution for all scatterings or rescattering only. | |
1918 | dSigmaSum += dSigmaCorr; | |
1919 | dSigmaResc += dSigmaCorr; | |
1920 | ||
1921 | // Determine whether current rescattering should be selected. | |
1922 | if (dSigmaCorr > rndmPtr->flat() * dSigmaSum) { | |
1923 | i1Sel = 0; | |
1924 | i2Sel = iB; | |
1925 | id1Sel = id1Tmp; | |
1926 | id2Sel = id2Tmp; | |
1927 | x1Sel = x1Tmp; | |
1928 | x2Sel = x2Tmp; | |
1929 | sHatSel = sHatTmp; | |
1930 | tHatSel = tHatTmp; | |
1931 | uHatSel = uHatTmp; | |
1932 | sigma2Sel = sigma2Tmp; | |
1933 | pickOtherSel = sigma2Tmp->pickedOther(); | |
1934 | } | |
1935 | } | |
1936 | ||
1937 | // Double rescatter: loop over already scattered partons from both sides. | |
1938 | // As before, allow to use only one parton per side to speed up. | |
1939 | if (allowDoubleRes) { | |
1940 | for (int iScat1 = 0; iScat1 < nScatA; ++iScat1) { | |
1941 | if (PREPICKRESCATTER && iScat1 != iScatA) continue; | |
1942 | int iA = scatteredA[iScat1]; | |
1943 | int id1Tmp = event[iA].id(); | |
1944 | for (int iScat2 = 0; iScat2 < nScatB; ++iScat2) { | |
1945 | if (PREPICKRESCATTER && iScat2 != iScatB) continue; | |
1946 | int iB = scatteredB[iScat2]; | |
1947 | int id2Tmp = event[iB].id(); | |
1948 | ||
1949 | // Calculate current sHat and check whether big enough. | |
1950 | double sHatTmp = (event[iA].p() + event[iB].p()).m2Calc(); | |
1951 | if (sHatTmp < 4. * pT2) continue; | |
1952 | ||
1953 | // Reconstruct other kinematical variables. (x values not needed.) | |
1954 | double x1Tmp = event[iA].pPos() / eCM; | |
1955 | double x2Tmp = event[iB].pNeg() / eCM; | |
1956 | double tauTmp = sHatTmp / sCM; | |
1957 | double root = sqrtpos(1. - xT2 / tauTmp); | |
1958 | double tHatTmp = -0.5 * sHatTmp * (1. - root); | |
1959 | double uHatTmp = -0.5 * sHatTmp * (1. + root); | |
1960 | ||
1961 | // Assign pointers to processes relevant for incoming flavour choice: | |
1962 | // g + g, q + g, q + qbar (same flavour), q + q(bar) (the rest). | |
1963 | if (id1Tmp == 21 && id2Tmp == 21) sigma2Tmp = &sigma2gg; | |
1964 | else if (id1Tmp == 21 || id2Tmp == 21) sigma2Tmp = &sigma2qg; | |
1965 | else if (id1Tmp == -id2Tmp) sigma2Tmp = &sigma2qqbarSame; | |
1966 | else sigma2Tmp = &sigma2qq; | |
1967 | ||
1968 | // Evaluate cross sections, include possibility of K factor. | |
1969 | // Use kinematics for two symmetrical configurations (tHat <-> uHat). | |
1970 | // (Not necessary, but makes for better MC efficiency.) | |
1971 | double dSigmaPartonCorr = Kfactor | |
1972 | * sigma2Tmp->sigma( id1Tmp, id2Tmp, x1Tmp, x2Tmp, sHatTmp, tHatTmp, | |
1973 | uHatTmp, alpS, alpEM); | |
1974 | ||
1975 | // Combine cross section and Jacobian tHat -> pT2 | |
1976 | // (with safety minvalue). | |
1977 | double dSigmaCorr = dSigmaPartonCorr / max(ROOTMIN, root); | |
1978 | ||
1979 | // Dampen cross section at small pT values; part of formalism. | |
1980 | // Use pT2 corrected for massive kinematics at this step. | |
1981 | //?? double pT2massive = dSigmaDtSel->pT2MPI(); | |
1982 | double pT2massive = pT2; | |
1983 | dSigmaCorr *= pow2( pT2massive / (pT2massive + pT20) ); | |
1984 | ||
1985 | // Compensate for prepicked rescattering if appropriate. | |
1986 | if (PREPICKRESCATTER) dSigmaCorr *= nScatA * nScatB; | |
1987 | ||
1988 | // Sum up total contribution for all scatterings or rescattering only. | |
1989 | dSigmaSum += dSigmaCorr; | |
1990 | dSigmaResc += dSigmaCorr; | |
1991 | ||
1992 | // Determine whether current rescattering should be selected. | |
1993 | if (dSigmaCorr > rndmPtr->flat() * dSigmaSum) { | |
1994 | i1Sel = iA; | |
1995 | i2Sel = iB; | |
1996 | id1Sel = id1Tmp; | |
1997 | id2Sel = id2Tmp; | |
1998 | x1Sel = x1Tmp; | |
1999 | x2Sel = x2Tmp; | |
2000 | sHatSel = sHatTmp; | |
2001 | tHatSel = tHatTmp; | |
2002 | uHatSel = uHatTmp; | |
2003 | sigma2Sel = sigma2Tmp; | |
2004 | pickOtherSel = sigma2Tmp->pickedOther(); | |
2005 | } | |
2006 | } | |
2007 | } | |
2008 | } | |
2009 | ||
2010 | // Done. | |
2011 | return dSigmaResc; | |
2012 | } | |
2013 | ||
2014 | ||
2015 | //-------------------------------------------------------------------------- | |
2016 | ||
2017 | // Calculate factor relating matter overlap and interaction rate, | |
2018 | // i.e. k in <n_interaction(b)> = k * overlap(b) (neglecting | |
2019 | // n_int = 0 corrections and energy-momentum conservation effects). | |
2020 | ||
2021 | void MultipartonInteractions::overlapInit() { | |
2022 | ||
2023 | // Initial values for iteration. Step size of b integration. | |
2024 | nAvg = sigmaInt / sigmaND; | |
2025 | kNow = 0.5; | |
2026 | int stepDir = 1; | |
2027 | double deltaB = BSTEP; | |
2028 | if (bProfile == 2) deltaB *= min( 0.5, 2.5 * coreRadius); | |
2029 | if (bProfile == 3) deltaB *= max(1., pow(2. / expPow, 1. / expPow)); | |
2030 | ||
2031 | // Further variables, with dummy initial values. | |
2032 | double nNow = 0.; | |
2033 | double kLow = 0.; | |
2034 | double nLow = 0.; | |
2035 | double kHigh = 0.; | |
2036 | double nHigh = 0.; | |
2037 | double overlapNow = 0.; | |
2038 | double probNow = 0.; | |
2039 | double overlapInt = 0.5; | |
2040 | double probInt = 0.; | |
2041 | double probOverlapInt = 0.; | |
2042 | double bProbInt = 0.; | |
2043 | normPi = 1. / (2. * M_PI); | |
2044 | ||
2045 | // Subdivision into low-b and high-b region by interaction rate. | |
2046 | bool pastBDiv = false; | |
2047 | double overlapHighB = 0.; | |
2048 | ||
2049 | // For x-dependent matter profile, try to find a0 rather than k. | |
2050 | // Existing framework is used, but with substitutions: | |
2051 | // a0 tuned according to Int( Pint(b), d^2b ) = sigmaND, | |
2052 | // nAvg = sigmaND, kNow = a0now, kLow = a0low, kHigh = a0high, | |
2053 | // nNow = probInt, nLow = probIntLow, nHigh = probIntHigh. | |
2054 | double rescale2 = 1.; | |
2055 | if (bProfile == 4) { | |
2056 | nAvg = sigmaND; | |
2057 | kNow = XDEP_A0 / 2.0; | |
2058 | } | |
2059 | ||
2060 | // First close k into an interval by binary steps, | |
2061 | // then find k by successive interpolation. | |
2062 | do { | |
2063 | if (stepDir == 1) kNow *= 2.; | |
2064 | else if (stepDir == -1) kNow *= 0.5; | |
2065 | else kNow = kLow + (nAvg - nLow) * (kHigh - kLow) / (nHigh - nLow); | |
2066 | ||
2067 | // Overlap trivial if no impact parameter dependence. | |
2068 | if (bProfile <= 0 || bProfile > 4) { | |
2069 | probInt = 0.5 * M_PI * (1. - exp(-kNow)); | |
2070 | probOverlapInt = probInt / M_PI; | |
2071 | bProbInt = probInt; | |
2072 | ||
2073 | // Ratio of b-integrated k * overlap / (1 - exp( - k * overlap)). | |
2074 | nNow = M_PI * kNow * overlapInt / probInt; | |
2075 | ||
2076 | // Else integrate overlap over impact parameter. | |
2077 | } else if (bProfile < 4) { | |
2078 | ||
2079 | // Reset integrals. | |
2080 | overlapInt = (bProfile == 3) ? 0. : 0.5; | |
2081 | probInt = 0.; | |
2082 | probOverlapInt = 0.; | |
2083 | bProbInt = 0.; | |
2084 | pastBDiv = false; | |
2085 | overlapHighB = 0.; | |
2086 | ||
2087 | // Step in b space. | |
2088 | double b = -0.5 * deltaB; | |
2089 | double bArea = 0.; | |
2090 | do { | |
2091 | b += deltaB; | |
2092 | bArea = 2. * M_PI * b * deltaB; | |
2093 | ||
2094 | // Evaluate overlap at current b value. | |
2095 | if (bProfile == 1) { | |
2096 | overlapNow = normPi * exp( -b*b); | |
2097 | } else if (bProfile == 2) { | |
2098 | overlapNow = normPi * ( fracA * exp( -min(EXPMAX, b*b)) | |
2099 | + fracB * exp( -min(EXPMAX, b*b / radius2B)) / radius2B | |
2100 | + fracC * exp( -min(EXPMAX, b*b / radius2C)) / radius2C ); | |
2101 | } else { | |
2102 | overlapNow = normPi * exp( -pow( b, expPow)); | |
2103 | overlapInt += bArea * overlapNow; | |
2104 | } | |
2105 | if (pastBDiv) overlapHighB += bArea * overlapNow; | |
2106 | ||
2107 | // Calculate interaction probability and integrate. | |
2108 | probNow = 1. - exp( -min(EXPMAX, M_PI * kNow * overlapNow)); | |
2109 | probInt += bArea * probNow; | |
2110 | probOverlapInt += bArea * overlapNow * probNow; | |
2111 | bProbInt += b * bArea * probNow; | |
2112 | ||
2113 | // Check when interaction probability has dropped sufficiently. | |
2114 | if (!pastBDiv && probNow < PROBATLOWB) { | |
2115 | bDiv = b + 0.5 * deltaB; | |
2116 | pastBDiv = true; | |
2117 | } | |
2118 | ||
2119 | // Continue out in b until overlap too small. | |
2120 | } while (b < 1. || b * probNow > BMAX); | |
2121 | ||
2122 | // Ratio of b-integrated k * overlap / (1 - exp( - k * overlap)). | |
2123 | nNow = M_PI * kNow * overlapInt / probInt; | |
2124 | ||
2125 | // x-dependent matter profile. | |
2126 | } else if (bProfile == 4) { | |
2127 | rescale2 = pow2(kNow / XDEP_A0); | |
2128 | probInt = 0.; | |
2129 | double b = 0.5 * bstepNow; | |
2130 | for (int bBin = 0; bBin < XDEP_BBIN; bBin++) { | |
2131 | double bArea = 2. * M_PI * b * bstepNow; | |
2132 | double pIntNow = 1 - exp( -min(EXPMAX, sigmaIntWgt[bBin] / rescale2) ); | |
2133 | probInt += bArea * rescale2 * pIntNow; | |
2134 | b += bstepNow; | |
2135 | } | |
2136 | nNow = probInt; | |
2137 | } | |
2138 | ||
2139 | // Replace lower or upper limit of k. | |
2140 | if (nNow < nAvg) { | |
2141 | kLow = kNow; | |
2142 | nLow = nNow; | |
2143 | if (stepDir == -1) stepDir = 0; | |
2144 | } else { | |
2145 | kHigh = kNow; | |
2146 | nHigh = nNow; | |
2147 | if (stepDir == 1) stepDir = -1; | |
2148 | } | |
2149 | ||
2150 | // Continue iteration until convergence. | |
2151 | } while (abs(nNow - nAvg) > KCONVERGE * nAvg); | |
2152 | ||
2153 | // Save relevant final numbers for overlap values. | |
2154 | if (bProfile >= 0 && bProfile < 4) { | |
2155 | double avgOverlap = probOverlapInt / probInt; | |
2156 | zeroIntCorr = probOverlapInt / overlapInt; | |
2157 | normOverlap = normPi * zeroIntCorr / avgOverlap; | |
2158 | bAvg = bProbInt / probInt; | |
2159 | ||
2160 | // Values for x-dependent matter profile. | |
2161 | } else if (bProfile == 4) { | |
2162 | // bAvg = Int(b * Pint(b), d2b) / sigmaND. | |
2163 | // zeroIntCorr = Int(<n(b)> * Pint(b), d2b) / sigmaInt. | |
2164 | bAvg = 0.; | |
2165 | zeroIntCorr = 0.; | |
2166 | double b = 0.5 * bstepNow; | |
2167 | for (int bBin = 0; bBin < XDEP_BBIN; bBin++) { | |
2168 | double bArea = 2. * M_PI * b * bstepNow; | |
2169 | double pIntNow = 1. - exp( -min(EXPMAX, sigmaIntWgt[bBin] / rescale2) ); | |
2170 | bAvg += sqrt(rescale2) * b * bArea * rescale2 * pIntNow; | |
2171 | zeroIntCorr += bArea * sigmaIntWgt[bBin] * pIntNow; | |
2172 | b += bstepNow; | |
2173 | } | |
2174 | bAvg /= nNow; | |
2175 | zeroIntCorr /= sigmaInt; | |
2176 | ||
2177 | // Other required values. | |
2178 | a0now = kNow; | |
2179 | infoPtr->seta0MPI(a0now * XDEP_SMB2FM); | |
2180 | a02now = a0now * a0now; | |
2181 | double xMin = 2. * pTmin / infoPtr->eCM(); | |
2182 | a2max = a0now * (XDEP_A1 + a1 * log(1. / xMin)); | |
2183 | a2max *= a2max; | |
2184 | } | |
2185 | ||
2186 | // Relative rates for preselection of low-b and high-b region. | |
2187 | // Other useful combinations for subsequent selection. | |
2188 | if (bProfile > 0 && bProfile <= 3) { | |
2189 | probLowB = M_PI * bDiv*bDiv; | |
2190 | double probHighB = M_PI * kNow * overlapHighB; | |
2191 | if (bProfile == 1) probHighB = M_PI * kNow * 0.5 * exp( -bDiv*bDiv); | |
2192 | else if (bProfile == 2) { | |
2193 | fracAhigh = fracA * exp( -bDiv*bDiv); | |
2194 | fracBhigh = fracB * exp( -bDiv*bDiv / radius2B); | |
2195 | fracChigh = fracC * exp( -bDiv*bDiv / radius2C); | |
2196 | fracABChigh = fracAhigh + fracBhigh + fracChigh; | |
2197 | probHighB = M_PI * kNow * 0.5 * fracABChigh; | |
2198 | } else { | |
2199 | cDiv = pow( bDiv, expPow); | |
2200 | cMax = max(2. * expRev, cDiv); | |
2201 | } | |
2202 | probLowB /= (probLowB + probHighB); | |
2203 | } | |
2204 | ||
2205 | } | |
2206 | ||
2207 | //-------------------------------------------------------------------------- | |
2208 | ||
2209 | // Pick impact parameter and interaction rate enhancement beforehand, | |
2210 | // i.e. before even the hardest interaction for minimum-bias events. | |
2211 | ||
2212 | void MultipartonInteractions::overlapFirst() { | |
2213 | ||
2214 | // Trivial values if no impact parameter dependence. | |
2215 | if (bProfile <= 0 || bProfile > 4) { | |
2216 | bNow = bAvg; | |
2217 | enhanceB = zeroIntCorr; | |
2218 | bIsSet = true; | |
2219 | isAtLowB = true; | |
2220 | return; | |
2221 | } | |
2222 | ||
2223 | // Preliminary choice between and inside low-b and high-b regions. | |
2224 | double overlapNow = 0.; | |
2225 | double probAccept = 0.; | |
2226 | do { | |
2227 | ||
2228 | // Treatment in low-b region: pick b flat in area. | |
2229 | if (rndmPtr->flat() < probLowB) { | |
2230 | isAtLowB = true; | |
2231 | bNow = bDiv * sqrt(rndmPtr->flat()); | |
2232 | ||
2233 | // Evaluate overlap and from that acceptance probability. | |
2234 | if (bProfile == 1) overlapNow = normPi * exp( -bNow*bNow); | |
2235 | else if (bProfile == 2) overlapNow = normPi * | |
2236 | ( fracA * exp( -bNow*bNow) | |
2237 | + fracB * exp( -bNow*bNow / radius2B) / radius2B | |
2238 | + fracC * exp( -bNow*bNow / radius2C) / radius2C ); | |
2239 | else overlapNow = normPi * exp( -pow( bNow, expPow)); | |
2240 | probAccept = 1. - exp( -min(EXPMAX, M_PI * kNow * overlapNow)); | |
2241 | ||
2242 | // Treatment in high-b region: pick b according to overlap. | |
2243 | } else { | |
2244 | isAtLowB = false; | |
2245 | ||
2246 | // For simple and double Gaussian pick b according to exp(-b^2 / r^2). | |
2247 | if (bProfile == 1) { | |
2248 | bNow = sqrt(bDiv*bDiv - log(rndmPtr->flat())); | |
2249 | overlapNow = normPi * exp( -min(EXPMAX, bNow*bNow)); | |
2250 | } else if (bProfile == 2) { | |
2251 | double pickFrac = rndmPtr->flat() * fracABChigh; | |
2252 | if (pickFrac < fracAhigh) | |
2253 | bNow = sqrt(bDiv*bDiv - log(rndmPtr->flat())); | |
2254 | else if (pickFrac < fracAhigh + fracBhigh) | |
2255 | bNow = sqrt(bDiv*bDiv - radius2B * log(rndmPtr->flat())); | |
2256 | else bNow = sqrt(bDiv*bDiv - radius2C * log(rndmPtr->flat())); | |
2257 | overlapNow = normPi * ( fracA * exp( -min(EXPMAX, bNow*bNow)) | |
2258 | + fracB * exp( -min(EXPMAX, bNow*bNow / radius2B)) / radius2B | |
2259 | + fracC * exp( -min(EXPMAX, bNow*bNow / radius2C)) / radius2C ); | |
2260 | ||
2261 | // For exp( - b^expPow) transform to variable c = b^expPow so that | |
2262 | // f(b) = b * exp( - b^expPow) -> f(c) = c^r * exp(-c) with r = expRev. | |
2263 | // case hasLowPow: expPow < 2 <=> r > 0: preselect according to | |
2264 | // f(c) < N exp(-c/2) and then accept with N' * c^r * exp(-c/2). | |
2265 | } else if (hasLowPow) { | |
2266 | double cNow, acceptC; | |
2267 | do { | |
2268 | cNow = cDiv - 2. * log(rndmPtr->flat()); | |
2269 | acceptC = pow(cNow / cMax, expRev) * exp( -0.5 * (cNow - cMax)); | |
2270 | } while (acceptC < rndmPtr->flat()); | |
2271 | bNow = pow( cNow, 1. / expPow); | |
2272 | overlapNow = normPi * exp( -cNow); | |
2273 | ||
2274 | // case !hasLowPow: expPow >= 2 <=> - 1 < r < 0: preselect according to | |
2275 | // f(c) < N exp(-c) and then accept with N' * c^r. | |
2276 | } else { | |
2277 | double cNow, acceptC; | |
2278 | do { | |
2279 | cNow = cDiv - log(rndmPtr->flat()); | |
2280 | acceptC = pow(cNow / cDiv, expRev); | |
2281 | } while (acceptC < rndmPtr->flat()); | |
2282 | bNow = pow( cNow, 1. / expPow); | |
2283 | overlapNow = normPi * exp( -cNow); | |
2284 | } | |
2285 | double temp = M_PI * kNow *overlapNow; | |
2286 | probAccept = (1. - exp( -min(EXPMAX, temp))) / temp; | |
2287 | } | |
2288 | ||
2289 | // Confirm choice of b value. Derive enhancement factor. | |
2290 | } while (probAccept < rndmPtr->flat()); | |
2291 | ||
2292 | // Same enhancement for hardest process and all subsequent MPI | |
2293 | enhanceB = enhanceBmax = enhanceBnow = (normOverlap / normPi) * overlapNow; | |
2294 | ||
2295 | // Done. | |
2296 | bIsSet = true; | |
2297 | ||
2298 | } | |
2299 | ||
2300 | //-------------------------------------------------------------------------- | |
2301 | ||
2302 | // Pick impact parameter and interaction rate enhancement afterwards, | |
2303 | // i.e. after a hard interaction is known but before rest of MPI treatment. | |
2304 | ||
2305 | void MultipartonInteractions::overlapNext(Event& event, double pTscale) { | |
2306 | ||
2307 | // Default, valid for bProfile = 0. Also initial Sudakov. | |
2308 | enhanceB = zeroIntCorr; | |
2309 | if (bProfile <= 0 || bProfile > 4) return; | |
2310 | ||
2311 | // Alternative choices of event scale for Sudakov in (pT, b) space. | |
2312 | if (bSelScale == 1) { | |
2313 | vector<double> mmT; | |
2314 | for (int i = 5; i < event.size(); ++i) if (event[i].isFinal()) { | |
2315 | mmT.push_back( event[i].m() + event[i].mT() ); | |
2316 | for (int j = int(mmT.size()) - 1; j > 0; --j) | |
2317 | if (mmT[j] > mmT[j - 1]) swap( mmT[j], mmT[j - 1] ); | |
2318 | } | |
2319 | pTscale = 0.5 * mmT[0]; | |
2320 | for (int j = 1; j < int(mmT.size()); ++j) pTscale += mmT[j] / (j + 1.); | |
2321 | } else if (bSelScale == 2) pTscale = event.scale(); | |
2322 | double pT2scale = pTscale*pTscale; | |
2323 | ||
2324 | // Use trial interaction for x-dependent matter profile. | |
2325 | if (bProfile == 4) { | |
2326 | double pTtrial = 0.; | |
2327 | do { | |
2328 | // Pick b according to wanted O(b, x1, x2). | |
2329 | double expb2 = rndmPtr->flat(); | |
2330 | double w1 = XDEP_A1 + a1 * log(1. / infoPtr->x1()); | |
2331 | double w2 = XDEP_A1 + a1 * log(1. / infoPtr->x2()); | |
2332 | double fac = a02now * (w1 * w1 + w2 * w2); | |
2333 | b2now = - fac * log(expb2); | |
2334 | bNow = sqrt(b2now); | |
2335 | ||
2336 | // Enhancement factor for the hard process and overestimate | |
2337 | // for fastPT2. Note that existing framework has a (1. / sigmaND) | |
2338 | // present. | |
2339 | enhanceB = sigmaND / M_PI / fac * expb2; | |
2340 | enhanceBmax = sigmaND / 2. / M_PI / a02now | |
2341 | * exp( -b2now / 2. / a2max ); | |
2342 | ||
2343 | // Trial interaction. Keep going until pTtrial < pTscale. | |
2344 | pTtrial = pTnext(pTmax, pTmin, event); | |
2345 | } while (pTtrial > pTscale); | |
2346 | bIsSet = true; | |
2347 | return; | |
2348 | } | |
2349 | ||
2350 | // Begin loop over pT-dependent rejection of b value. | |
2351 | do { | |
2352 | ||
2353 | // Flat enhancement distribution for simple Gaussian. | |
2354 | if (bProfile == 1) { | |
2355 | double expb2 = rndmPtr->flat(); | |
2356 | // Same enhancement for hardest process and all subsequent MPI. | |
2357 | enhanceB = enhanceBmax = enhanceBnow = normOverlap * expb2; | |
2358 | bNow = sqrt( -log(expb2)); | |
2359 | ||
2360 | // For double Gaussian go the way via b, according to each Gaussian. | |
2361 | } else if (bProfile == 2) { | |
2362 | double bType = rndmPtr->flat(); | |
2363 | double b2 = -log( rndmPtr->flat() ); | |
2364 | if (bType < fracA) ; | |
2365 | else if (bType < fracA + fracB) b2 *= radius2B; | |
2366 | else b2 *= radius2C; | |
2367 | // Same enhancement for hardest process and all subsequent MPI. | |
2368 | enhanceB = enhanceBmax = enhanceBnow = normOverlap * | |
2369 | ( fracA * exp( -min(EXPMAX, b2)) | |
2370 | + fracB * exp( -min(EXPMAX, b2 / radius2B)) / radius2B | |
2371 | + fracC * exp( -min(EXPMAX, b2 / radius2C)) / radius2C ); | |
2372 | bNow = sqrt(b2); | |
2373 | ||
2374 | // For exp( - b^expPow) transform to variable c = b^expPow so that | |
2375 | // f(b) = b * exp( - b^expPow) -> f(c) = c^r * exp(-c) with r = expRev. | |
2376 | // case hasLowPow: expPow < 2 <=> r > 0: | |
2377 | // f(c) < r^r exp(-r) for c < 2r, < (2r)^r exp(-r-c/2) for c > 2r. | |
2378 | } else if (bProfile == 3 && hasLowPow) { | |
2379 | double cNow, acceptC; | |
2380 | double probLowC = expRev / (expRev + pow(2., expRev) * exp( - expRev)); | |
2381 | do { | |
2382 | if (rndmPtr->flat() < probLowC) { | |
2383 | cNow = 2. * expRev * rndmPtr->flat(); | |
2384 | acceptC = pow( cNow / expRev, expRev) * exp(expRev - cNow); | |
2385 | } else { | |
2386 | cNow = 2. * (expRev - log( rndmPtr->flat() )); | |
2387 | acceptC = pow(0.5 * cNow / expRev, expRev) | |
2388 | * exp(expRev - 0.5 * cNow); | |
2389 | } | |
2390 | } while (acceptC < rndmPtr->flat()); | |
2391 | // Same enhancement for hardest process and all subsequent MPI. | |
2392 | enhanceB = enhanceBmax = enhanceBnow = normOverlap *exp(-cNow); | |
2393 | bNow = pow( cNow, 1. / expPow); | |
2394 | ||
2395 | // case !hasLowPow: expPow >= 2 <=> - 1 < r < 0: | |
2396 | // f(c) < c^r for c < 1, f(c) < exp(-c) for c > 1. | |
2397 | } else if (bProfile == 3 && !hasLowPow) { | |
2398 | double cNow, acceptC; | |
2399 | double probLowC = expPow / (2. * exp(-1.) + expPow); | |
2400 | do { | |
2401 | if (rndmPtr->flat() < probLowC) { | |
2402 | cNow = pow( rndmPtr->flat(), 0.5 * expPow); | |
2403 | acceptC = exp(-cNow); | |
2404 | } else { | |
2405 | cNow = 1. - log( rndmPtr->flat() ); | |
2406 | acceptC = pow( cNow, expRev); | |
2407 | } | |
2408 | } while (acceptC < rndmPtr->flat()); | |
2409 | // Same enhancement for hardest process and all subsequent MPI. | |
2410 | enhanceB = enhanceBmax = enhanceBnow = normOverlap * exp(-cNow); | |
2411 | bNow = pow( cNow, 1. / expPow); | |
2412 | } | |
2413 | ||
2414 | // Evaluate "Sudakov form factor" for not having a harder interaction. | |
2415 | } while (sudakov(pT2scale, enhanceB) < rndmPtr->flat()); | |
2416 | ||
2417 | // Done. | |
2418 | bIsSet = true; | |
2419 | } | |
2420 | ||
2421 | //-------------------------------------------------------------------------- | |
2422 | ||
2423 | // Printe statistics on number of multiparton-interactions processes. | |
2424 | ||
2425 | void MultipartonInteractions::statistics(bool resetStat, ostream& os) { | |
2426 | ||
2427 | // Header. | |
2428 | os << "\n *------- PYTHIA Multiparton Interactions Statistics -----" | |
2429 | << "---*\n" | |
2430 | << " | " | |
2431 | << " |\n" | |
2432 | << " | Note: excludes hardest subprocess if already listed above " | |
2433 | << " |\n" | |
2434 | << " | " | |
2435 | << " |\n" | |
2436 | << " | Subprocess Code | Times" | |
2437 | << " |\n" | |
2438 | << " | | " | |
2439 | << " |\n" | |
2440 | << " |------------------------------------------------------------" | |
2441 | << "-|\n" | |
2442 | << " | | " | |
2443 | << " |\n"; | |
2444 | ||
2445 | // Loop over existing processes. Sum of all subprocesses. | |
2446 | int numberSum = 0; | |
2447 | for ( map<int, int>::iterator iter = nGen.begin(); iter != nGen.end(); | |
2448 | ++iter) { | |
2449 | int code = iter -> first; | |
2450 | int number = iter->second; | |
2451 | numberSum += number; | |
2452 | ||
2453 | // Find process name that matches code. | |
2454 | string name = " "; | |
2455 | bool foundName = false; | |
2456 | SigmaMultiparton* dSigma; | |
2457 | for (int i = 0; i < 4; ++i) { | |
2458 | if (i == 0) dSigma = &sigma2gg; | |
2459 | else if (i == 1) dSigma = &sigma2qg; | |
2460 | else if (i == 2) dSigma = &sigma2qqbarSame; | |
2461 | else dSigma = &sigma2qq; | |
2462 | int nProc = dSigma->nProc(); | |
2463 | for (int iProc = 0; iProc < nProc; ++iProc) | |
2464 | if (dSigma->codeProc(iProc) == code) { | |
2465 | name = dSigma->nameProc(iProc); | |
2466 | foundName = true; | |
2467 | } | |
2468 | if (foundName) break; | |
2469 | } | |
2470 | ||
2471 | // Print individual process info. | |
2472 | os << " | " << left << setw(40) << name << right << setw(5) << code | |
2473 | << " | " << setw(11) << number << " |\n"; | |
2474 | } | |
2475 | ||
2476 | // Print summed process info. | |
2477 | os << " | " | |
2478 | << " |\n" | |
2479 | << " | " << left << setw(45) << "sum" << right << " | " << setw(11) | |
2480 | << numberSum << " |\n"; | |
2481 | ||
2482 | // Listing finished. | |
2483 | os << " | | " | |
2484 | << " |\n" | |
2485 | << " *------- End PYTHIA Multiparton Interactions Statistics ----" | |
2486 | << "-*" << endl; | |
2487 | ||
2488 | // Optionally reset statistics contents. | |
2489 | if (resetStat) resetStatistics(); | |
2490 | ||
2491 | } | |
2492 | ||
2493 | //========================================================================== | |
2494 | ||
2495 | } // end namespace Pythia8 |