PYTHIA8/pythia8170/src/ProcessLevel.cxx

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