PYTHIA8/pythia8145/src/MultipleInteractions.cxx

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