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