PYTHIA8/pythia8145/src/ProcessLevel.cxx

   1 // ProcessLevel.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 ProcessLevel class.
   7
   8 #include "ProcessLevel.h"
   9
  10 namespace Pythia8 {
  11
  12 //==========================================================================
  13
  14 // The ProcessLevel class.
  15
  16 //--------------------------------------------------------------------------
  17
  18 // Constants: could be changed here if desired, but normally should not.
  19 // These are of technical nature, as described for each.
  20
  21 // Allow a few failures in final construction of events.
  22 const int ProcessLevel::MAXLOOP = 5;
  23
  24 //--------------------------------------------------------------------------
  25
  26 // Destructor.
  27
  28 ProcessLevel::~ProcessLevel() {
  29
  30   // Run through list of first hard processes and delete them.
  31   for (int i = 0; i < int(containerPtrs.size()); ++i)
  32     delete containerPtrs[i];
  33
  34   // Run through list of second hard processes and delete them.
  35   for (int i =0; i < int(container2Ptrs.size()); ++i)
  36     delete container2Ptrs[i];
  37
  38 }
  39
  40 //--------------------------------------------------------------------------
  41
  42 // Main routine to initialize generation process.
  43
  44 bool ProcessLevel::init( Info* infoPtrIn, Settings& settings,
  45   ParticleData* particleDataPtrIn, Rndm* rndmPtrIn,
  46   BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn, Couplings* couplingsPtrIn,
  47   SigmaTotal* sigmaTotPtrIn, bool doLHA, SusyLesHouches* slhaPtrIn,
  48   UserHooks* userHooksPtrIn, vector<SigmaProcess*>& sigmaPtrs, ostream& os) {
  49
  50   // Store input pointers for future use.
  51   infoPtr         = infoPtrIn;
  52   particleDataPtr = particleDataPtrIn;
  53   rndmPtr         = rndmPtrIn;
  54   beamAPtr        = beamAPtrIn;
  55   beamBPtr        = beamBPtrIn;
  56   couplingsPtr    = couplingsPtrIn;
  57   sigmaTotPtr     = sigmaTotPtrIn;
  58   userHooksPtr    = userHooksPtrIn;
  59   slhaPtr         = slhaPtrIn;
  60
  61   // Send on some input pointers.
  62   resonanceDecays.init( infoPtr, particleDataPtr, rndmPtr);
  63
  64   // Set up SigmaTotal. Store sigma_nondiffractive for future use.
  65   sigmaTotPtr->init( infoPtr, settings, particleDataPtr);
  66   int    idA = infoPtr->idA();
  67   int    idB = infoPtr->idB();
  68   double eCM = infoPtr->eCM();
  69   sigmaTotPtr->calc( idA, idB, eCM);
  70   sigmaND = sigmaTotPtr->sigmaND();
  71
  72   // Options to allow second hard interaction and resonance decays.
  73   doSecondHard  = settings.flag("SecondHard:generate");
  74   doSameCuts    = settings.flag("PhaseSpace:sameForSecond");
  75   doResDecays   = settings.flag("ProcessLevel:resonanceDecays");
  76   startColTag   = settings.mode("Event:startColTag");
  77
  78   // Mass and pT cuts for two hard processes.
  79   mHatMin1      = settings.parm("PhaseSpace:mHatMin");
  80   mHatMax1      = settings.parm("PhaseSpace:mHatMax");
  81   if (mHatMax1 < mHatMin1) mHatMax1 = eCM;
  82   pTHatMin1     = settings.parm("PhaseSpace:pTHatMin");
  83   pTHatMax1     = settings.parm("PhaseSpace:pTHatMax");
  84   if (pTHatMax1 < pTHatMin1) pTHatMax1 = eCM;
  85   mHatMin2      = settings.parm("PhaseSpace:mHatMinSecond");
  86   mHatMax2      = settings.parm("PhaseSpace:mHatMaxSecond");
  87   if (mHatMax2 < mHatMin2) mHatMax2 = eCM;
  88   pTHatMin2     = settings.parm("PhaseSpace:pTHatMinSecond");
  89   pTHatMax2     = settings.parm("PhaseSpace:pTHatMaxSecond");
  90   if (pTHatMax2 < pTHatMin2) pTHatMax2 = eCM;
  91
  92   // Check whether mass and pT ranges agree or overlap.
  93   cutsAgree     = doSameCuts;
  94   if (mHatMin2 == mHatMin1 && mHatMax2 == mHatMax1 && pTHatMin2 == pTHatMin1
  95       && pTHatMax2 == pTHatMax1) cutsAgree = true;
  96   cutsOverlap   = cutsAgree;
  97   if (!cutsAgree) {
  98     bool mHatOverlap = (max( mHatMin1, mHatMin2)
  99                       < min( mHatMax1, mHatMax2));
 100     bool pTHatOverlap = (max( pTHatMin1, pTHatMin2)
 101                        < min( pTHatMax1, pTHatMax2));
 102     if (mHatOverlap && pTHatOverlap) cutsOverlap = true;
 103   }
 104
 105   // Set up containers for all the internal hard processes.
 106   SetupContainers setupContainers;
 107   setupContainers.init(containerPtrs, settings, particleDataPtr, couplingsPtr);
 108
 109   // Append containers for external hard processes, if any.
 110   if (sigmaPtrs.size() > 0) {
 111     for (int iSig = 0; iSig < int(sigmaPtrs.size()); ++iSig)
 112       containerPtrs.push_back( new ProcessContainer(sigmaPtrs[iSig],
 113       true) );
 114   }
 115
 116   // Append single container for Les Houches processes, if any.
 117   if (doLHA) {
 118     SigmaProcess* sigmaPtr = new SigmaLHAProcess();
 119     containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
 120
 121     // Store location of this container, and send in LHA pointer.
 122     iLHACont = containerPtrs.size() - 1;
 123     containerPtrs[iLHACont]->setLHAPtr(lhaUpPtr);
 124   }
 125
 126   // If no processes found then refuse to do anything.
 127   if ( int(containerPtrs.size()) == 0) {
 128     infoPtr->errorMsg("Error in ProcessLevel::init: "
 129       "no process switched on");
 130     return false;
 131   }
 132
 133   // Check that SUSY couplings were indeed initialized where necessary.
 134   bool hasSUSY = false;
 135   for (int i = 0; i < int(containerPtrs.size()); ++i)
 136     if (containerPtrs[i]->isSUSY()) hasSUSY = true;
 137
 138   //If SUSY processes requested but no SUSY couplings present
 139   if(hasSUSY && !couplingsPtr->isSUSY) {
 140     infoPtr->errorMsg("Error in ProcessLevel::init: "
 141       "SUSY process switched on but no SUSY couplings found");
 142     return false;
 143   }
 144
 145   // Initialize each process.
 146   int numberOn = 0;
 147   for (int i = 0; i < int(containerPtrs.size()); ++i)
 148     if (containerPtrs[i]->init(true, infoPtr, settings, particleDataPtr,
 149       rndmPtr, beamAPtr, beamBPtr, couplingsPtr, sigmaTotPtr,
 150       &resonanceDecays, slhaPtr, userHooksPtr)) ++numberOn;
 151
 152   // Sum maxima for Monte Carlo choice.
 153   sigmaMaxSum = 0.;
 154   for (int i = 0; i < int(containerPtrs.size()); ++i)
 155     sigmaMaxSum += containerPtrs[i]->sigmaMax();
 156
 157   // Option to pick a second hard interaction: repeat as above.
 158   int number2On = 0;
 159   if (doSecondHard) {
 160     setupContainers.init2(container2Ptrs, settings);
 161     if ( int(container2Ptrs.size()) == 0) {
 162       infoPtr->errorMsg("Error in ProcessLevel::init: "
 163         "no second hard process switched on");
 164       return false;
 165     }
 166     for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
 167       if (container2Ptrs[i2]->init(false, infoPtr, settings, particleDataPtr,
 168         rndmPtr, beamAPtr, beamBPtr, couplingsPtr, sigmaTotPtr,
 169         &resonanceDecays, slhaPtr, userHooksPtr)) ++number2On;
 170     sigma2MaxSum = 0.;
 171     for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
 172       sigma2MaxSum += container2Ptrs[i2]->sigmaMax();
 173   }
 174
 175   // Construct string with incoming beams and for cm energy.
 176   string collision = "We collide " + particleDataPtr->name(idA)
 177     + " with " + particleDataPtr->name(idB) + " at a CM energy of ";
 178   string pad( 51 - collision.length(), ' ');
 179
 180   // Print initialization information: header.
 181   os << "\n *-------  PYTHIA Process Initialization  ---------"
 182      << "-----------------*\n"
 183      << " |                                                   "
 184      << "               |\n"
 185      << " | " << collision << scientific << setprecision(3)<< setw(9) << eCM
 186      << " GeV" << pad << " |\n"
 187      << " |                                                   "
 188      << "               |\n"
 189      << " |---------------------------------------------------"
 190      << "---------------|\n"
 191      << " |                                                   "
 192      << " |             |\n"
 193      << " | Subprocess                                    Code"
 194      << " |   Estimated |\n"
 195      << " |                                                   "
 196      << " |    max (mb) |\n"
 197      << " |                                                   "
 198      << " |             |\n"
 199      << " |---------------------------------------------------"
 200      << "---------------|\n"
 201      << " |                                                   "
 202      << " |             |\n";
 203
 204
 205   // Loop over existing processes: print individual process info.
 206   for (int i = 0; i < int(containerPtrs.size()); ++i)
 207   os << " | " << left << setw(45) << containerPtrs[i]->name()
 208      << right << setw(5) << containerPtrs[i]->code() << " | "
 209      << scientific << setprecision(3) << setw(11)
 210      << containerPtrs[i]->sigmaMax() << " |\n";
 211
 212   // Loop over second hard processes, if any, and repeat as above.
 213   if (doSecondHard) {
 214     os << " |                                                   "
 215        << " |             |\n"
 216        << " |---------------------------------------------------"
 217        <<"---------------|\n"
 218        << " |                                                   "
 219        <<" |             |\n";
 220     for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
 221     os << " | " << left << setw(45) << container2Ptrs[i2]->name()
 222        << right << setw(5) << container2Ptrs[i2]->code() << " | "
 223        << scientific << setprecision(3) << setw(11)
 224        << container2Ptrs[i2]->sigmaMax() << " |\n";
 225   }
 226
 227   // Listing finished.
 228   os << " |                                                     "
 229      << "             |\n"
 230      << " *-------  End PYTHIA Process Initialization ----------"
 231      <<"-------------*" << endl;
 232
 233   // If sum of maxima vanishes then refuse to do anything.
 234   if ( numberOn == 0  || sigmaMaxSum <= 0.) {
 235     infoPtr->errorMsg("Error in ProcessLevel::init: "
 236       "all processes have vanishing cross sections");
 237     return false;
 238   }
 239   if ( doSecondHard && (number2On == 0  || sigma2MaxSum <= 0.) ) {
 240     infoPtr->errorMsg("Error in ProcessLevel::init: "
 241       "all second hard processes have vanishing cross sections");
 242     return false;
 243   }
 244
 245   // If two hard processes then check whether some (but not all) agree.
 246   allHardSame  = true;
 247   noneHardSame = true;
 248   if (doSecondHard) {
 249     bool foundMatch = false;
 250
 251     // Check for each first process if matched in second.
 252     for (int i = 0; i < int(containerPtrs.size()); ++i) {
 253       foundMatch = false;
 254       if (cutsOverlap)
 255       for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
 256         if (container2Ptrs[i2]->code() == containerPtrs[i]->code())
 257           foundMatch = true;
 258       containerPtrs[i]->isSame( foundMatch );
 259       if (!foundMatch)  allHardSame = false;
 260       if ( foundMatch) noneHardSame = false;
 261     }
 262
 263     // Check for each second process if matched in first.
 264     for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2) {
 265       foundMatch = false;
 266       if (cutsOverlap)
 267       for (int i = 0; i < int(containerPtrs.size()); ++i)
 268         if (containerPtrs[i]->code() == container2Ptrs[i2]->code())
 269           foundMatch = true;
 270       container2Ptrs[i2]->isSame( foundMatch );
 271       if (!foundMatch)  allHardSame = false;
 272       if ( foundMatch) noneHardSame = false;
 273     }
 274   }
 275
 276   // Concluding classification, including cuts.
 277   if (!cutsAgree) allHardSame = false;
 278   someHardSame = !allHardSame && !noneHardSame;
 279
 280   // Reset counters for average impact-parameter enhancement.
 281   nImpact       = 0;
 282   sumImpactFac  = 0.;
 283   sum2ImpactFac = 0.;
 284
 285   // Statistics for LHA events.
 286   codeLHA.resize(0);
 287   nEvtLHA.resize(0);
 288
 289   // Done.
 290   return true;
 291 }
 292
 293 //--------------------------------------------------------------------------
 294
 295 // Main routine to generate the hard process.
 296
 297 bool ProcessLevel::next( Event& process) {
 298
 299   // Generate the next event with two or one hard interactions.
 300   bool physical = (doSecondHard) ? nextTwo( process) : nextOne( process);
 301
 302   // Check that colour assignments make sense.
 303   if (physical) physical = checkColours( process);
 304
 305   // Done.
 306   return physical;
 307 }
 308
 309 //--------------------------------------------------------------------------
 310
 311 // Accumulate and update statistics (after possible user veto).
 312
 313 void ProcessLevel::accumulate() {
 314
 315   // Increase number of accepted events.
 316   containerPtrs[iContainer]->accumulate();
 317
 318   // Provide current generated cross section estimate.
 319   long   nTrySum    = 0;
 320   long   nSelSum    = 0;
 321   long   nAccSum    = 0;
 322   double sigmaSum   = 0.;
 323   double delta2Sum  = 0.;
 324   double sigSelSum  = 0.;
 325   for (int i = 0; i < int(containerPtrs.size()); ++i)
 326   if (containerPtrs[i]->sigmaMax() != 0.) {
 327     nTrySum        += containerPtrs[i]->nTried();
 328     nSelSum        += containerPtrs[i]->nSelected();
 329     nAccSum        += containerPtrs[i]->nAccepted();
 330     sigmaSum       += containerPtrs[i]->sigmaMC();
 331     delta2Sum      += pow2(containerPtrs[i]->deltaMC());
 332     sigSelSum      += containerPtrs[i]->sigmaSelMC();
 333   }
 334
 335   // For Les Houches events find subprocess type and update counter.
 336   if (infoPtr->isLHA()) {
 337     int codeLHANow = infoPtr->codeSub();
 338     int iFill = -1;
 339     for (int i = 0; i < int(codeLHA.size()); ++i)
 340       if (codeLHANow == codeLHA[i]) iFill = i;
 341     if (iFill >= 0) ++nEvtLHA[iFill];
 342
 343     // Add new process when new code and then arrange in order.
 344     else {
 345       codeLHA.push_back(codeLHANow);
 346       nEvtLHA.push_back(1);
 347       for (int i = int(codeLHA.size()) - 1; i > 0; --i) {
 348         if (codeLHA[i] < codeLHA[i - 1]) {
 349           swap(codeLHA[i], codeLHA[i - 1]);
 350           swap(nEvtLHA[i], nEvtLHA[i - 1]);
 351         }
 352         else break;
 353       }
 354     }
 355   }
 356
 357   // Normally only one hard interaction. Then store info and done.
 358   if (!doSecondHard) {
 359     double deltaSum = sqrtpos(delta2Sum);
 360     infoPtr->setSigma( nTrySum, nSelSum, nAccSum, sigmaSum, deltaSum);
 361     return;
 362   }
 363
 364   // Increase counter for a second hard interaction.
 365   container2Ptrs[i2Container]->accumulate();
 366
 367   // Update statistics on average impact factor.
 368   ++nImpact;
 369   sumImpactFac     += infoPtr->enhanceMI();
 370   sum2ImpactFac    += pow2(infoPtr->enhanceMI());
 371
 372   // Cross section estimate for second hard process.
 373   double sigma2Sum  = 0.;
 374   double sig2SelSum = 0.;
 375   for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
 376   if (container2Ptrs[i2]->sigmaMax() != 0.) {
 377     nTrySum        += container2Ptrs[i2]->nTried();
 378     sigma2Sum      += container2Ptrs[i2]->sigmaMC();
 379     sig2SelSum     += container2Ptrs[i2]->sigmaSelMC();
 380   }
 381
 382   // Average impact-parameter factor and error.
 383   double invN       = 1. / max(1, nImpact);
 384   double impactFac  = max( 1., sumImpactFac * invN);
 385   double impactErr2 = ( sum2ImpactFac * invN / pow2(impactFac) - 1.) * invN;
 386
 387   // Cross section estimate for combination of first and second process.
 388   // Combine two possible ways and take average.
 389   double sigmaComb  = 0.5 * (sigmaSum * sig2SelSum + sigSelSum * sigma2Sum);
 390   sigmaComb        *= impactFac / sigmaND;
 391   if (allHardSame) sigmaComb *= 0.5;
 392   double deltaComb  = sqrtpos(2. / nAccSum + impactErr2) * sigmaComb;
 393
 394   // Store info and done.
 395   infoPtr->setSigma( nTrySum, nSelSum, nAccSum, sigmaComb, deltaComb);
 396
 397 }
 398
 399 //--------------------------------------------------------------------------
 400
 401 // Print statistics on cross sections and number of events.
 402
 403 void ProcessLevel::statistics(bool reset, ostream& os) {
 404
 405   // Special processing if two hard interactions selected.
 406   if (doSecondHard) {
 407     statistics2(reset, os);
 408     return;
 409   }
 410
 411   // Header.
 412   os << "\n *-------  PYTHIA Event and Cross Section Statistics  ------"
 413      << "-------------------------------------------------------*\n"
 414      << " |                                                            "
 415      << "                                                     |\n"
 416      << " | Subprocess                                    Code |       "
 417      << "     Number of events       |      sigma +- delta    |\n"
 418      << " |                                                    |       "
 419      << "Tried   Selected   Accepted |     (estimated) (mb)   |\n"
 420      << " |                                                    |       "
 421      << "                            |                        |\n"
 422      << " |------------------------------------------------------------"
 423      << "-----------------------------------------------------|\n"
 424      << " |                                                    |       "
 425      << "                            |                        |\n";
 426
 427   // Reset sum counters.
 428   long   nTrySum   = 0;
 429   long   nSelSum   = 0;
 430   long   nAccSum   = 0;
 431   double sigmaSum  = 0.;
 432   double delta2Sum = 0.;
 433
 434   // Loop over existing processes.
 435   for (int i = 0; i < int(containerPtrs.size()); ++i)
 436   if (containerPtrs[i]->sigmaMax() != 0.) {
 437
 438     // Read info for process. Sum counters.
 439     long   nTry    = containerPtrs[i]->nTried();
 440     long   nSel    = containerPtrs[i]->nSelected();
 441     long   nAcc    = containerPtrs[i]->nAccepted();
 442     double sigma   = containerPtrs[i]->sigmaMC();
 443     double delta   = containerPtrs[i]->deltaMC();
 444     nTrySum       += nTry;
 445     nSelSum       += nSel;
 446     nAccSum       += nAcc;
 447     sigmaSum      += sigma;
 448     delta2Sum     += pow2(delta);
 449
 450     // Print individual process info.
 451     os << " | " << left << setw(45) << containerPtrs[i]->name()
 452        << right << setw(5) << containerPtrs[i]->code() << " | "
 453        << setw(11) << nTry << " " << setw(10) << nSel << " "
 454        << setw(10) << nAcc << " | " << scientific << setprecision(3)
 455        << setw(11) << sigma << setw(11) << delta << " |\n";
 456
 457     // Print subdivision by user code for Les Houches process.
 458     if (containerPtrs[i]->code() == 9999)
 459     for (int j = 0; j < int(codeLHA.size()); ++j)
 460       os << " |    ... whereof user classification code " << setw(10)
 461          << codeLHA[j] << " |                        " << setw(10)
 462          << nEvtLHA[j] << " |                        | \n";
 463   }
 464
 465   // Print summed process info.
 466   os << " |                                                    |       "
 467      << "                            |                        |\n"
 468      << " | " << left << setw(50) << "sum" << right << " | " << setw(11)
 469      << nTrySum << " " << setw(10) << nSelSum << " " << setw(10)
 470      << nAccSum << " | " << scientific << setprecision(3) << setw(11)
 471      << sigmaSum << setw(11) << sqrtpos(delta2Sum) << " |\n";
 472
 473   // Listing finished.
 474   os << " |                                                            "
 475      << "                                                     |\n"
 476      << " *-------  End PYTHIA Event and Cross Section Statistics -----"
 477      << "-----------------------------------------------------*" << endl;
 478
 479   // Optionally reset statistics contants.
 480   if (reset) for (int i = 0; i < int(containerPtrs.size()); ++i)
 481     containerPtrs[i]->reset();
 482
 483 }
 484
 485 //--------------------------------------------------------------------------
 486
 487 // Generate the next event with one interaction.
 488
 489 bool ProcessLevel::nextOne( Event& process) {
 490
 491   // Update CM energy for phase space selection.
 492   double eCM = infoPtr->eCM();
 493   for (int i = 0; i < int(containerPtrs.size()); ++i)
 494     containerPtrs[i]->newECM(eCM);
 495
 496   // Outer loop in case of rare failures.
 497   bool physical = true;
 498   for (int loop = 0; loop < MAXLOOP; ++loop) {
 499     if (!physical) process.clear();
 500     physical = true;
 501
 502     // Loop over tries until trial event succeeds.
 503     for ( ; ; ) {
 504
 505       // Pick one of the subprocesses.
 506       double sigmaMaxNow = sigmaMaxSum * rndmPtr->flat();
 507       int iMax = containerPtrs.size() - 1;
 508       iContainer = -1;
 509       do sigmaMaxNow -= containerPtrs[++iContainer]->sigmaMax();
 510       while (sigmaMaxNow > 0. && iContainer < iMax);
 511
 512       // Do a trial event of this subprocess; accept or not.
 513       if (containerPtrs[iContainer]->trialProcess()) break;
 514
 515       // Check for end-of-file condition for Les Houches events.
 516       if (infoPtr->atEndOfFile()) return false;
 517     }
 518
 519     // Update sum of maxima if current maximum violated.
 520     if (containerPtrs[iContainer]->newSigmaMax()) {
 521       sigmaMaxSum = 0.;
 522       for (int i = 0; i < int(containerPtrs.size()); ++i)
 523         sigmaMaxSum += containerPtrs[i]->sigmaMax();
 524     }
 525
 526     // Construct kinematics of acceptable process.
 527     if ( !containerPtrs[iContainer]->constructProcess( process) )
 528       physical = false;
 529
 530     // Do all resonance decays.
 531     if ( physical && doResDecays
 532       && !containerPtrs[iContainer]->decayResonances( process) )
 533       physical = false;
 534
 535     // Add any junctions to the process event record list.
 536     if (physical) findJunctions( process);
 537
 538     // Outer loop should normally work first time around.
 539     if (physical) break;
 540   }
 541
 542   // Done.
 543   return physical;
 544 }
 545
 546 //--------------------------------------------------------------------------
 547
 548 // Generate the next event with two hard interactions.
 549
 550 bool ProcessLevel::nextTwo( Event& process) {
 551
 552   // Update CM energy for phase space selection.
 553   double eCM = infoPtr->eCM();
 554   for (int i = 0; i < int(containerPtrs.size()); ++i)
 555     containerPtrs[i]->newECM(eCM);
 556   for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
 557     container2Ptrs[i2]->newECM(eCM);
 558
 559   // Outer loop in case of rare failures.
 560   bool physical = true;
 561   for (int loop = 0; loop < MAXLOOP; ++loop) {
 562     if (!physical) process.clear();
 563     physical = true;
 564
 565     // Loop over both hard processes to find consistent common kinematics.
 566     for ( ; ; ) {
 567
 568       // Loop internally over tries for hardest process until succeeds.
 569       for ( ; ; ) {
 570
 571         // Pick one of the subprocesses.
 572         double sigmaMaxNow = sigmaMaxSum * rndmPtr->flat();
 573         int iMax = containerPtrs.size() - 1;
 574         iContainer = -1;
 575         do sigmaMaxNow -= containerPtrs[++iContainer]->sigmaMax();
 576         while (sigmaMaxNow > 0. && iContainer < iMax);
 577
 578         // Do a trial event of this subprocess; accept or not.
 579         if (containerPtrs[iContainer]->trialProcess()) break;
 580
 581         // Check for end-of-file condition for Les Houches events.
 582         if (infoPtr->atEndOfFile()) return false;
 583       }
 584
 585       // Update sum of maxima if current maximum violated.
 586       if (containerPtrs[iContainer]->newSigmaMax()) {
 587         sigmaMaxSum = 0.;
 588         for (int i = 0; i < int(containerPtrs.size()); ++i)
 589           sigmaMaxSum += containerPtrs[i]->sigmaMax();
 590       }
 591
 592       // Loop internally over tries for second hardest process until succeeds.
 593       for ( ; ; ) {
 594
 595         // Pick one of the subprocesses.
 596         double sigma2MaxNow = sigma2MaxSum * rndmPtr->flat();
 597         int i2Max = container2Ptrs.size() - 1;
 598         i2Container = -1;
 599         do sigma2MaxNow -= container2Ptrs[++i2Container]->sigmaMax();
 600         while (sigma2MaxNow > 0. && i2Container < i2Max);
 601
 602         // Do a trial event of this subprocess; accept or not.
 603         if (container2Ptrs[i2Container]->trialProcess()) break;
 604       }
 605
 606       // Update sum of maxima if current maximum violated.
 607       if (container2Ptrs[i2Container]->newSigmaMax()) {
 608         sigma2MaxSum = 0.;
 609         for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
 610           sigma2MaxSum += container2Ptrs[i2]->sigmaMax();
 611       }
 612
 613       // Check whether common set of x values is kinematically possible.
 614       double xA1      = containerPtrs[iContainer]->x1();
 615       double xB1      = containerPtrs[iContainer]->x2();
 616       double xA2      = container2Ptrs[i2Container]->x1();
 617       double xB2      = container2Ptrs[i2Container]->x2();
 618       if (xA1 + xA2 >= 1. || xB1 + xB2 >= 1.) continue;
 619
 620       // Reset beam contents. Naive parton densities for second interaction.
 621       // (Subsequent procedure could be symmetrized, but would be overkill.)
 622       beamAPtr->clear();
 623       beamBPtr->clear();
 624       int    idA1     = containerPtrs[iContainer]->id1();
 625       int    idB1     = containerPtrs[iContainer]->id2();
 626       int    idA2     = container2Ptrs[i2Container]->id1();
 627       int    idB2     = container2Ptrs[i2Container]->id2();
 628       double Q2Fac1   = containerPtrs[iContainer]->Q2Fac();
 629       double Q2Fac2   = container2Ptrs[i2Container]->Q2Fac();
 630       double pdfA2Raw = beamAPtr->xf( idA2, xA2,Q2Fac2);
 631       double pdfB2Raw = beamBPtr->xf( idB2, xB2,Q2Fac2);
 632
 633       // Remove partons in first interaction from beams.
 634       beamAPtr->append( 3, idA1, xA1);
 635       beamAPtr->xfISR( 0, idA1, xA1, Q2Fac1);
 636       beamAPtr->pickValSeaComp();
 637       beamBPtr->append( 4, idB1, xB1);
 638       beamBPtr->xfISR( 0, idB1, xB1, Q2Fac1);
 639       beamBPtr->pickValSeaComp();
 640
 641       // Reevaluate pdf's for second interaction and weight by reduction.
 642       double pdfA2Mod = beamAPtr->xfMI( idA2, xA2,Q2Fac2);
 643       double pdfB2Mod = beamBPtr->xfMI( idB2, xB2,Q2Fac2);
 644       double wtPdfMod = (pdfA2Mod * pdfB2Mod) / (pdfA2Raw * pdfB2Raw);
 645       if (wtPdfMod < rndmPtr->flat()) continue;
 646
 647       // Reduce by a factor of 2 for identical processes when others not,
 648       // and when in same phase space region.
 649       bool toLoop = false;
 650       if ( someHardSame && containerPtrs[iContainer]->isSame()
 651         && container2Ptrs[i2Container]->isSame()) {
 652         if (cutsAgree) {
 653           if (rndmPtr->flat() > 0.5) toLoop = true;
 654         } else {
 655         double mHat1 = containerPtrs[iContainer]->mHat();
 656         double pTHat1 = containerPtrs[iContainer]->pTHat();
 657         double mHat2 = container2Ptrs[i2Container]->mHat();
 658         double pTHat2 = container2Ptrs[i2Container]->pTHat();
 659         if (mHat1 > mHatMin2 && mHat1 < mHatMax2
 660            && pTHat1 > pTHatMin2 && pTHat1 < pTHatMax2
 661            && mHat2 > mHatMin1 && mHat2 < mHatMax1
 662            && pTHat2 > pTHatMin1 && pTHat2 < pTHatMax1
 663            && rndmPtr->flat() > 0.5) toLoop = true;
 664         }
 665       }
 666       if (toLoop) continue;
 667
 668       // If come this far then acceptable event.
 669       break;
 670     }
 671
 672     // Construct kinematics of acceptable processes.
 673     Event process2;
 674     process2.init( "(second hard)", particleDataPtr, startColTag);
 675     process2.initColTag();
 676     if ( !containerPtrs[iContainer]->constructProcess( process) )
 677       physical = false;
 678     if (physical && !container2Ptrs[i2Container]->constructProcess( process2,
 679       false) ) physical = false;
 680
 681     // Do all resonance decays.
 682     if ( physical && doResDecays
 683       &&  !containerPtrs[iContainer]->decayResonances( process) )
 684       physical = false;
 685     if ( physical && doResDecays
 686       &&  !container2Ptrs[i2Container]->decayResonances( process2) )
 687       physical = false;
 688
 689     // Append second hard interaction to normal process object.
 690     if (physical) combineProcessRecords( process, process2);
 691
 692     // Add any junctions to the process event record list.
 693     if (physical) findJunctions( process);
 694
 695     // Outer loop should normally work first time around.
 696     if (physical) break;
 697   }
 698
 699   // Done.
 700   return physical;
 701 }
 702
 703 //--------------------------------------------------------------------------
 704
 705 // Append second hard interaction to normal process object.
 706 // Complication: all resonance decay chains must be put at the end.
 707
 708 void ProcessLevel::combineProcessRecords( Event& process, Event& process2) {
 709
 710   // Find first event record size, excluding resonances.
 711   int nSize = process.size();
 712   int nHard = 5;
 713   while (nHard < nSize && process[nHard].mother1() == 3) ++nHard;
 714
 715   // Save resonance products temporarily elsewhere.
 716   vector<Particle> resProd;
 717   if (nSize > nHard) {
 718     for (int i = nHard; i < nSize; ++i) resProd.push_back( process[i] );
 719     process.popBack(nSize - nHard);
 720   }
 721
 722   // Find second event record size, excluding resonances.
 723   int nSize2 = process2.size();
 724   int nHard2 = 5;
 725   while (nHard2 < nSize2 && process2[nHard2].mother1() == 3) ++nHard2;
 726
 727   // Find amount of necessary position and colour offset for second process.
 728   int addPos  = nHard  - 3;
 729   int addCol  = process.lastColTag() - startColTag;
 730
 731   // Loop over all particles (except beams) from second process.
 732   for (int i = 3; i < nSize2; ++i) {
 733
 734     // Offset mother and daughter pointers and colour tags of particle.
 735     process2[i].offsetHistory( 2, addPos, 2, addPos);
 736     process2[i].offsetCol( addCol);
 737
 738     // Append hard-process particles from process2 to process.
 739     if (i < nHard2) process.append( process2[i] );
 740   }
 741
 742   // Reinsert resonance decay chains of first hard process.
 743   int addPos2 = nHard2 - 3;
 744   if (nHard < nSize) {
 745
 746     // Offset daughter pointers of unmoved mothers.
 747     for (int i = 5; i < nHard; ++i)
 748       process[i].offsetHistory( 0, 0, nHard - 1, addPos2);
 749
 750     // Modify history of resonance products when restoring.
 751     for (int i = 0; i < int(resProd.size()); ++i) {
 752       resProd[i].offsetHistory( nHard - 1, addPos2, nHard - 1, addPos2);
 753       process.append( resProd[i] );
 754     }
 755   }
 756
 757   // Insert resonance decay chains of second hard process.
 758   if (nHard2 < nSize2) {
 759     int nHard3  = nHard + nHard2 - 3;
 760     int addPos3 = nSize - nHard;
 761
 762     // Offset daughter pointers of second-process mothers.
 763     for (int i = nHard + 2; i < nHard3; ++i)
 764       process[i].offsetHistory( 0, 0, nHard3 - 1, addPos3);
 765
 766     // Modify history of second-process resonance products and insert.
 767     for (int i = nHard2; i < nSize2; ++i) {
 768       process2[i].offsetHistory( nHard3 - 1, addPos3, nHard3 - 1, addPos3);
 769       process.append( process2[i] );
 770     }
 771   }
 772
 773   // Store PDF scale for second interaction.
 774   process.scaleSecond( process2.scale() );
 775
 776 }
 777
 778 //--------------------------------------------------------------------------
 779
 780 // Add any junctions to the process event record list.
 781 // Also check that do not doublebook if called repeatedly.
 782
 783 void ProcessLevel::findJunctions( Event& junEvent) {
 784
 785   // Find the total number of unmatched colours & anticolours in the event
 786   map<int,int> colMap, acolMap;
 787   for (int i = 0; i < junEvent.size(); ++i) {
 788     // Consider all external lines (including initial incoming ones)
 789     int col  = junEvent[i].col();
 790     int acol = junEvent[i].acol();
 791     if (junEvent[i].isFinal() || abs(junEvent[i].status()) == 21) {
 792       // Cross initial-state colors
 793       if (junEvent[i].status() == -21) {
 794         col  = junEvent[i].acol();
 795         acol = junEvent[i].col();
 796       }
 797       if (col > 0) {
 798         // If unmatched (so far), add end. Else erase matching parton.
 799         if (acolMap.find(col) == acolMap.end() ) colMap[col] = i;
 800         else acolMap.erase(col);
 801       }
 802       if (acol > 0) {
 803         // If unmatched (so far), add end. Else erase matching parton.
 804         if (colMap.find(acol) == colMap.end()) acolMap[acol] = i;
 805         else colMap.erase(acol);
 806       }
 807     }
 808     // For internal particles, check for internal junction lines
 809     // by counting total number of triplets at each vertex
 810     else if (abs(junEvent[i].status()) == 22) {
 811       int cSum = 0;
 812       int idRes = junEvent[i].id();
 813       // Only consider triplets for the time being, since these are the ones
 814       // most likely to cause trouble.
 815       // (E.g., internal stop in resonant stop production)
 816       if (junEvent[i].colType() != 1) continue;
 817       // Decaying particle crossed => triplet <-> antitriplet
 818       if (idRes < 0) cSum += (junEvent[i].colType() % 2) ;
 819       else cSum -= (junEvent[i].colType() % 2) ;
 820       int iDau1 = junEvent[i].daughter1();
 821       int iDau2 = max(iDau1,junEvent[i].daughter2());
 822       for (int iDau=iDau1; iDau<=iDau2; iDau++) {
 823         if (junEvent[iDau].id() > 0) cSum += (junEvent[iDau].colType() % 2) ;
 824         else cSum -= (junEvent[iDau].colType() % 2) ;
 825       }
 826       // Add dangling internal tag if not already present
 827       if (cSum != 0) {
 828         if (idRes > 0 && acolMap.find(col) == acolMap.end()) acolMap[col]=i;
 829         if (idRes < 0 && colMap.find(acol) == colMap.end()) colMap[acol]=i;
 830       }
 831     }
 832   }
 833
 834   // Number of expected junctions and antijunctions
 835   // (Add 2 to account for max 2 junction-antijunction connections)
 836   int nJunc     = (colMap.size()+2)/3;
 837   int nAntiJunc = (acolMap.size()+2)/3;
 838
 839   //cout<<" nJunc = "<<nJunc<<" nAntiJunc = "<<nAntiJunc<<endl;
 840
 841   // If no junctions expected, return without doing anything more
 842   if (max(nJunc,nAntiJunc) == 0) return;
 843
 844   // Junctions expected. Check (and skip) any that have already been added
 845   for (int iJun = 0; iJun < junEvent.sizeJunction(); ++iJun) {
 846     // Remove the tags corresponding to each of the 3 existing junction legs
 847     for (int j = 0; j < 3; ++j) {
 848       int colNow = junEvent.colJunction(iJun, j);
 849       if (junEvent.kindJunction(iJun) % 2 == 1) colMap.erase(colNow);
 850       else acolMap.erase(colNow);
 851     }
 852   }
 853
 854   // Update Number of expected junctions and antijunctions
 855   nJunc     = (colMap.size()+2)/3;
 856   nAntiJunc = (acolMap.size()+2)/3;
 857
 858   // cout<<" AR nJunc = "<<nJunc<<" nAntiJunc = "<<nAntiJunc<<endl;
 859
 860   // If no (more) junctions expected, return without doing anything more
 861   if (max(nJunc,nAntiJunc) == 0) return;
 862
 863   // Loop: first check colors, then anticolors
 864   for (int mAcol = 0; mAcol <= 1; mAcol++) {
 865     map<int,int>& cMap = colMap;
 866     if (mAcol == 1) cMap = acolMap;
 867
 868     // Trace colours to find junction vertices and missing lines from
 869     // internal junction-antijunction connections (such connections *must*
 870     // be represented by internal lines connecting the junction vertices)
 871     for (map<int,int>::iterator it = cMap.begin(); it != cMap.end(); it++) {
 872       int col1 = it->first;
 873       int i1   = it->second;
 874       int i2   = i1;
 875       bool goForward = false;
 876       // Reverse direction for particles in initial state
 877       if ( !junEvent[i1].isFinal() ) goForward = true;
 878       // Skip colour tags that have already been associated with a junction
 879       bool skipCol = false;
 880       for (int iJun = 0; iJun < junEvent.sizeJunction(); ++iJun) {
 881         // Only look at other (anti)junctions
 882         if (junEvent.kindJunction(iJun) % 2 == mAcol) continue;
 883         for (int j = 0; j < 3; ++j) {
 884           int colNow = junEvent.colJunction(iJun, j);
 885           if (colNow == col1) skipCol = true;
 886         }
 887       }
 888       if (skipCol) {
 889         //cout << "Skipping (a)col : "<<col1<<" at "<<i1<<endl;
 890         continue;
 891       }
 892       //cout << "Tracing (a)col : "<<col1<<" from "<<i1<<endl;
 893       bool search = true;
 894       while ( search ) {
 895         // Step to mother (or daughter, if doing forward search)
 896         if (goForward) {
 897           i1     = i2;
 898           // For now, assume daughter1 always is the one that leads to the
 899           // hard interaction. Otherwise, may need to search among daughters.
 900           i2     = junEvent[i1].daughter1();
 901         } else {
 902           i2     = junEvent[i1].mother2();
 903           i1     = junEvent[i1].mother1();
 904         }
 905         int iDau1  = junEvent[i1].daughter1();
 906         int iDau2  = max(iDau1,junEvent[i1].daughter2());
 907         int barSum = 0;
 908         // If we went all the way back to the beam without finding a violation,
 909         // skip and try next
 910         if (i1 <= 0 || abs(junEvent[i1].status()) <= 20) search = false;
 911         // Check baryon number of vertex, counting incoming
 912         // partons with a minus sign
 913         else {
 914           //cout << "    at : "<<i1<<" colType = "<<junEvent[i1].colType()<<endl;
 915           if (abs(junEvent[i1].colType()) == 1) barSum -= junEvent[i1].colType();
 916           if (i2 != 0 && abs(junEvent[i2].colType()) == 1 )
 917             barSum -= junEvent[i2].colType();
 918           for (int iDau = iDau1; iDau <= iDau2; iDau++)
 919             if ( abs(junEvent[iDau].colType()) == 1 )
 920               barSum += junEvent[iDau].colType();
 921         }
 922         // If still searching, and if B is violated at this vertex => junction
 923         // Still needs additional refinement for complicated cases ???
 924         if (search && barSum != 0) {
 925           //cout<<" (anti)junction found at : "<<i1<<endl;
 926           // Find the 3 (anti)colour tags for this junction
 927           // If going forward, now change i2 from daughter to other mother
 928           if (goForward && i2 == iDau1) {
 929             if (junEvent[i2].mother1() == i1) i2 = junEvent[i2].mother2();
 930             else i2 = junEvent[i2].mother1();
 931           }
 932           //cout <<" Checking vertex: i1 = "<<i1<<" i2 = "<<i2<<endl;
 933           // Determine which (anti)color tags go to junction
 934           vector<int> colJu;
 935           int c1(0), c2(0), cMot1(0), cMot2(0);
 936           // First determine number of incoming lines
 937           if (mAcol == 0) {
 938             // For junctions: look at incoming anticolors
 939             c1 = junEvent[i1].acol();
 940             c2 = 0;
 941             // Also store incoming colors, to check for flow-through later
 942             cMot1  = junEvent[i1].col();
 943             if ( i2 != 0 ) {
 944               c2     = junEvent[i2].acol();
 945               cMot2  = junEvent[i1].col();
 946               // Remove annihilation-type flows
 947               if (c1 != 0 && junEvent[i2].col() == c1) c1 = 0;
 948               if (c2 != 0 && junEvent[i1].col() == c2) c2 = 0;
 949               // Remove anticolors flowing through diagram
 950               for (int iDau = iDau1; iDau <= iDau2; iDau++) {
 951                 int acolDau = junEvent[iDau].acol();
 952                 if (c1 != 0 && acolDau == c1) c1 = 0;
 953                 if (c2 != 0 && acolDau == c2) c2 = 0;
 954               }
 955             }
 956           } else {
 957             // For antijunctions: look at incoming colors
 958             c1 = junEvent[i1].col();
 959             c2 = 0;
 960             // Also store incoming anticolors, to check for flow-through later
 961             cMot1  = junEvent[i1].acol();
 962             if ( i2 != 0 ) {
 963               c2    = junEvent[i2].col();
 964               cMot2 = junEvent[i2].acol();
 965               // Remove annihilation-type flows
 966               if (c1 != 0 && junEvent[i2].acol() == c1) c1 = 0;
 967               if (c2 != 0 && junEvent[i1].acol() == c2) c2 = 0;
 968               // Remove colors flowing through diagram
 969               for (int iDau = iDau1; iDau <= iDau2; iDau++) {
 970                 int colDau = junEvent[iDau].col();
 971                 if (c1 != 0 && colDau == c1) c1 = 0;
 972                 if (c2 != 0 && colDau == c2) c2 = 0;
 973               }
 974             }
 975           }
 976           // Add non-annihilated colors
 977           if (c1 != 0) colJu.push_back(c1);
 978           if (c2 != 0) colJu.push_back(c2);
 979           // Determine junction kind from number of incoming lines
 980           int kindJun = 1 + mAcol + 2*colJu.size();
 981           //cout<<" Junction type : "<<kindJun<<endl;
 982           // Now add daughter color tags to junction
 983           for (int iDau = iDau1; iDau <= iDau2; iDau++) {
 984             c1 = 0;
 985             // For junctions: look at final-state colors, and vice-versa
 986             if (mAcol == 0) c1 = junEvent[iDau].col();
 987             else c1 = junEvent[iDau].acol();
 988             // Remove non-existent and (anti)colors flowing through diagram
 989             if (c1 == 0 || c1 == cMot1 || c1 == cMot2) continue;
 990             // Remove final-state matching color-anticolor pairs
 991             bool matchedPair = false;
 992             for (int jDau = iDau1; jDau <= iDau2; jDau++)
 993               if (mAcol==0 && c1 == junEvent[jDau].acol()) matchedPair = true;
 994               else if (mAcol==1 && c1 == junEvent[jDau].col()) matchedPair = true;
 995             if (matchedPair) continue;
 996             // Else add anticolor to junction ends
 997             else colJu.push_back(c1);
 998           }
 999           if (colJu.size() != 3) {
1000             cout<< " (ProcessLevel::findJunctions:) problem identifying junction "<<endl;
1001             search = false;
1002           } else {
1003             //cout<< " adding junction type "
1004             //     <<kindJun<<" "<<colJu[0]<<" "<<colJu[1]<<" "<<colJu[2]<<endl;
1005             junEvent.appendJunction( kindJun, colJu[0], colJu[1], colJu[2]);
1006             search = false;
1007           }
1008         }
1009       }
1010     }
1011   }
1012
1013   // Done.
1014 }
1015
1016 //--------------------------------------------------------------------------
1017
1018 // Check that colours match up.
1019
1020 bool ProcessLevel::checkColours( Event& process) {
1021
1022   // Variables and arrays for common usage.
1023   bool physical = true;
1024   bool match;
1025   int colType, col, acol, iPos, iNow, iNowA;
1026   vector<int> colTags, colPos, acolPos;
1027
1028   // Check that each particle has the kind of colours expected of it.
1029   for (int i = 0; i < process.size(); ++i) {
1030     colType = process[i].colType();
1031     col     = process[i].col();
1032     acol    = process[i].acol();
1033     if      (colType ==  0 && (col != 0 || acol != 0)) physical = false;
1034     else if (colType ==  1 && (col <= 0 || acol != 0)) physical = false;
1035     else if (colType == -1 && (col != 0 || acol <= 0)) physical = false;
1036     else if (colType ==  2 && (col <= 0 || acol <= 0)) physical = false;
1037     else if (colType < -1 || colType > 2)              physical = false;
1038
1039     // Add to the list of colour tags.
1040     if (col > 0) {
1041       match = false;
1042       for (int ic = 0; ic < int(colTags.size()) ; ++ic)
1043         if (col == colTags[ic]) match = true;
1044       if (!match) colTags.push_back(col);
1045     } else if (acol > 0) {
1046       match = false;
1047       for (int ic = 0; ic < int(colTags.size()) ; ++ic)
1048         if (acol == colTags[ic]) match = true;
1049       if (!match) colTags.push_back(acol);
1050     }
1051   }
1052
1053   // Warn and give up if particles did not have the expected colours.
1054   if (!physical) {
1055     infoPtr->errorMsg("Error in ProcessLevel::checkColours: "
1056       "incorrect colour assignment");
1057     return false;
1058   }
1059
1060   // Remove (anti)colours coming from an (anti)junction.
1061
1062   //cout<<" Removing junction tags in:"<<endl;
1063   //process.list();
1064   //cout<<"   Size of junctions "<<process.sizeJunction()<<endl;
1065
1066   for (int iJun = 0; iJun < process.sizeJunction(); ++iJun) {
1067     for (int j = 0; j < 3; ++j) {
1068       int colJun = process.colJunction(iJun, j);
1069       for (int ic = 0; ic < int(colTags.size()) ; ++ic)
1070         if (colJun == colTags[ic]) {
1071           colTags[ic] = colTags[colTags.size() - 1];
1072           colTags.pop_back();
1073           break;
1074         }
1075     }
1076   }
1077
1078   // Loop through all colour tags and find their positions (by sign).
1079   for (int ic = 0; ic < int(colTags.size()); ++ic) {
1080     col = colTags[ic];
1081     colPos.resize(0);
1082     acolPos.resize(0);
1083     for (int i = 0; i < process.size(); ++i) {
1084       if (process[i].col() == col) colPos.push_back(i);
1085       if (process[i].acol() == col) acolPos.push_back(i);
1086     }
1087
1088     // Trace colours back through decays; remove daughters.
1089     while (colPos.size() > 1) {
1090       iPos = colPos.size() - 1;
1091       iNow = colPos[iPos];
1092       if ( process[iNow].mother1() == colPos[iPos - 1]
1093         && process[iNow].mother2() == 0) colPos.pop_back();
1094       else break;
1095     }
1096     while (acolPos.size() > 1) {
1097       iPos = acolPos.size() - 1;
1098       iNow = acolPos[iPos];
1099       if ( process[iNow].mother1() == acolPos[iPos - 1]
1100         && process[iNow].mother2() == 0) acolPos.pop_back();
1101       else break;
1102     }
1103
1104     // Now colour should exist in only 2 copies.
1105     if (colPos.size() + acolPos.size() != 2) physical = false;
1106
1107     // If both colours or both anticolours then one mother of the other.
1108     else if (colPos.size() == 2) {
1109       iNow = colPos[1];
1110       if ( process[iNow].mother1() != colPos[0]
1111         && process[iNow].mother2() != colPos[0] ) physical = false;
1112     }
1113     else if (acolPos.size() == 2) {
1114       iNowA = acolPos[1];
1115       if ( process[iNowA].mother1() != acolPos[0]
1116         && process[iNowA].mother2() != acolPos[0] ) physical = false;
1117     }
1118
1119     // If one of each then should have same mother(s), or point to beams.
1120     else {
1121       iNow  = colPos[0];
1122       iNowA = acolPos[0];
1123       if ( process[iNow].status() == -21 &&  process[iNowA].status() == -21 );
1124       else if ( (process[iNow].mother1() != process[iNowA].mother1())
1125                 || (process[iNow].mother2() != process[iNowA].mother2()) ) {
1126         cout<<" Argh!"<<endl;
1127         physical = false;
1128       }
1129     }
1130
1131   }
1132
1133   // Error message if problem found. Done.
1134   if (!physical) infoPtr->errorMsg("Error in ProcessLevel::checkColours: "
1135     "unphysical colour flow");
1136   return physical;
1137
1138 }
1139
1140 //--------------------------------------------------------------------------
1141
1142 // Print statistics when two hard processes allowed.
1143
1144 void ProcessLevel::statistics2(bool reset, ostream& os) {
1145
1146   // Average impact-parameter factor and error.
1147   double invN          = 1. / max(1, nImpact);
1148   double impactFac     = max( 1., sumImpactFac * invN);
1149   double impactErr2    = ( sum2ImpactFac * invN / pow2(impactFac) - 1.) * invN;
1150
1151   // Derive scaling factor to be applied to first set of processes.
1152   double sigma2SelSum  = 0.;
1153   int    n2SelSum      = 0;
1154   for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2) {
1155     sigma2SelSum      += container2Ptrs[i2]->sigmaSelMC();
1156     n2SelSum          += container2Ptrs[i2]->nSelected();
1157   }
1158   double factor1       = impactFac * sigma2SelSum / sigmaND;
1159   double rel1Err       = sqrt(1. / max(1, n2SelSum) + impactErr2);
1160   if (allHardSame) factor1 *= 0.5;
1161
1162   // Derive scaling factor to be applied to second set of processes.
1163   double sigma1SelSum  = 0.;
1164   int    n1SelSum      = 0;
1165   for (int i = 0; i < int(containerPtrs.size()); ++i) {
1166     sigma1SelSum      += containerPtrs[i]->sigmaSelMC();
1167     n1SelSum          += containerPtrs[i]->nSelected();
1168   }
1169   double factor2       = impactFac * sigma1SelSum / sigmaND;
1170   if (allHardSame) factor2 *= 0.5;
1171   double rel2Err       = sqrt(1. / max(1, n1SelSum) + impactErr2);
1172
1173   // Header.
1174   os << "\n *-------  PYTHIA Event and Cross Section Statistics  ------"
1175      << "--------------------------------------------------*\n"
1176      << " |                                                            "
1177      << "                                                |\n"
1178      << " | Subprocess                               Code |            "
1179      << "Number of events       |      sigma +- delta    |\n"
1180      << " |                                               |       Tried"
1181      << "   Selected   Accepted |     (estimated) (mb)   |\n"
1182      << " |                                               |            "
1183      << "                       |                        |\n"
1184      << " |------------------------------------------------------------"
1185      << "------------------------------------------------|\n"
1186      << " |                                               |            "
1187      << "                       |                        |\n"
1188      << " | First hard process:                           |            "
1189      << "                       |                        |\n"
1190      << " |                                               |            "
1191      << "                       |                        |\n";
1192
1193   // Reset sum counters.
1194   long   nTrySum   = 0;
1195   long   nSelSum   = 0;
1196   long   nAccSum   = 0;
1197   double sigmaSum  = 0.;
1198   double delta2Sum = 0.;
1199
1200   // Loop over existing first processes.
1201   for (int i = 0; i < int(containerPtrs.size()); ++i)
1202   if (containerPtrs[i]->sigmaMax() != 0.) {
1203
1204     // Read info for process. Sum counters.
1205     long   nTry    = containerPtrs[i]->nTried();
1206     long   nSel    = containerPtrs[i]->nSelected();
1207     long   nAcc    = containerPtrs[i]->nAccepted();
1208     double sigma   = containerPtrs[i]->sigmaMC() * factor1;
1209     double delta2  = pow2( containerPtrs[i]->deltaMC() * factor1 );
1210     nTrySum       += nTry;
1211     nSelSum       += nSel;
1212     nAccSum       += nAcc;
1213     sigmaSum      += sigma;
1214     delta2Sum     += delta2;
1215     delta2        += pow2( sigma * rel1Err );
1216
1217     // Print individual process info.
1218     os << " | " << left << setw(40) << containerPtrs[i]->name()
1219        << right << setw(5) << containerPtrs[i]->code() << " | "
1220        << setw(11) << nTry << " " << setw(10) << nSel << " "
1221        << setw(10) << nAcc << " | " << scientific << setprecision(3)
1222        << setw(11) << sigma << setw(11) << sqrtpos(delta2) << " |\n";
1223   }
1224
1225   // Print summed info for first processes.
1226   delta2Sum       += pow2( sigmaSum * rel1Err );
1227   os << " |                                               |            "
1228      << "                       |                        |\n"
1229      << " | " << left << setw(45) << "sum" << right << " | " << setw(11)
1230      << nTrySum << " " << setw(10) << nSelSum << " " << setw(10)
1231      << nAccSum << " | " << scientific << setprecision(3) << setw(11)
1232      << sigmaSum << setw(11) << sqrtpos(delta2Sum) << " |\n";
1233
1234
1235   // Separation lines to second hard processes.
1236   os << " |                                               |            "
1237      << "                       |                        |\n"
1238      << " |------------------------------------------------------------"
1239      << "------------------------------------------------|\n"
1240      << " |                                               |            "
1241      << "                       |                        |\n"
1242      << " | Second hard process:                          |            "
1243      << "                       |                        |\n"
1244      << " |                                               |            "
1245      << "                       |                        |\n";
1246
1247   // Reset sum counters.
1248   nTrySum   = 0;
1249   nSelSum   = 0;
1250   nAccSum   = 0;
1251   sigmaSum  = 0.;
1252   delta2Sum = 0.;
1253
1254   // Loop over existing second processes.
1255   for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
1256   if (container2Ptrs[i2]->sigmaMax() != 0.) {
1257
1258     // Read info for process. Sum counters.
1259     long   nTry    = container2Ptrs[i2]->nTried();
1260     long   nSel    = container2Ptrs[i2]->nSelected();
1261     long   nAcc    = container2Ptrs[i2]->nAccepted();
1262     double sigma   = container2Ptrs[i2]->sigmaMC() * factor2;
1263     double delta2  = pow2( container2Ptrs[i2]->deltaMC() * factor2 );
1264     nTrySum       += nTry;
1265     nSelSum       += nSel;
1266     nAccSum       += nAcc;
1267     sigmaSum      += sigma;
1268     delta2Sum     += delta2;
1269     delta2        += pow2( sigma * rel2Err );
1270
1271     // Print individual process info.
1272     os << " | " << left << setw(40) << container2Ptrs[i2]->name()
1273        << right << setw(5) << container2Ptrs[i2]->code() << " | "
1274        << setw(11) << nTry << " " << setw(10) << nSel << " "
1275        << setw(10) << nAcc << " | " << scientific << setprecision(3)
1276        << setw(11) << sigma << setw(11) << sqrtpos(delta2) << " |\n";
1277   }
1278
1279   // Print summed info for second processes.
1280   delta2Sum       += pow2( sigmaSum * rel2Err );
1281   os << " |                                               |            "
1282      << "                       |                        |\n"
1283      << " | " << left << setw(45) << "sum" << right << " | " << setw(11)
1284      << nTrySum << " " << setw(10) << nSelSum << " " << setw(10)
1285      << nAccSum << " | " << scientific << setprecision(3) << setw(11)
1286      << sigmaSum << setw(11) << sqrtpos(delta2Sum) << " |\n";
1287
1288   // Print information on how the two processes were combined.
1289   os << " |                                               |            "
1290      << "                       |                        |\n"
1291      << " |------------------------------------------------------------"
1292      << "------------------------------------------------|\n"
1293      << " |                                                            "
1294      << "                                                |\n"
1295      << " | Uncombined cross sections for the two event sets were "
1296      << setw(10) << sigma1SelSum << " and " << sigma2SelSum << " mb, "
1297      << "respectively, combined  |\n"
1298      << " | using a sigma(nonDiffractive) of " << setw(10) << sigmaND
1299      << " mb and an impact-parameter enhancement factor of "
1300      << setw(10) << impactFac << ".   |\n";
1301
1302   // Listing finished.
1303   os << " |                                                            "
1304      << "                                                |\n"
1305      << " *-------  End PYTHIA Event and Cross Section Statistics -----"
1306      << "------------------------------------------------*" << endl;
1307
1308   // Optionally reset statistics contants.
1309   if (reset) {
1310     for (int i = 0; i < int(containerPtrs.size()); ++i)
1311       containerPtrs[i]->reset();
1312     for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
1313       container2Ptrs[i2]->reset();
1314   }
1315
1316 }
1317
1318 //==========================================================================
1319
1320 } // end namespace Pythia8