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