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