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