PYTHIA8/pythia8130/src/ProcessLevel.cxx

   1 // ProcessLevel.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2008 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 // Destructor.
  19
  20 ProcessLevel::~ProcessLevel() {
  21
  22   // Run through list of first hard processes and delete them.
  23   for (int i = 0; i < int(containerPtrs.size()); ++i)
  24     delete containerPtrs[i];
  25
  26   // Run through list of second hard processes and delete them.
  27   for (int i =0; i < int(container2Ptrs.size()); ++i)
  28     delete container2Ptrs[i];
  29
  30 }
  31
  32 //*********
  33
  34 // Main routine to initialize generation process.
  35
  36 bool ProcessLevel::init( Info* infoPtrIn, BeamParticle* beamAPtrIn,
  37   BeamParticle* beamBPtrIn, SigmaTotal* sigmaTotPtrIn, bool doLHA,
  38   SusyLesHouches* slhaPtrIn, UserHooks* userHooksPtrIn,
  39   vector<SigmaProcess*>& sigmaPtrs, ostream& os) {
  40
  41   // Store input pointers for future use.
  42   infoPtr       = infoPtrIn;
  43   beamAPtr      = beamAPtrIn;
  44   beamBPtr      = beamBPtrIn;
  45   sigmaTotPtr   = sigmaTotPtrIn;
  46   userHooksPtr  = userHooksPtrIn;
  47   slhaPtr       = slhaPtrIn;
  48
  49   // Send on some input pointers.
  50   resonanceDecays.init( infoPtr);
  51
  52   // Options to allow second hard interaction and resonance decays.
  53   doSecondHard  = Settings::flag("SecondHard:generate");
  54   doResDecays   = Settings::flag("ProcessLevel:resonanceDecays");
  55
  56   // Set up SigmaTotal. Store sigma_nondiffractive for future use.
  57   sigmaTotPtr->init( infoPtr);
  58   int    idA = infoPtr->idA();
  59   int    idB = infoPtr->idB();
  60   double eCM = infoPtr->eCM();
  61   sigmaTotPtr->calc( idA, idB, eCM);
  62   sigmaND = sigmaTotPtr->sigmaND();
  63
  64   // Initialize SUSY Les Houches Accord data
  65   if (!initSLHA()) return false;
  66
  67   // Set up containers for all the internal hard processes.
  68   SetupContainers setupContainers;
  69   setupContainers.init(containerPtrs);
  70
  71   // Append containers for external hard processes, if any.
  72   if (sigmaPtrs.size() > 0) {
  73     for (int iSig = 0; iSig < int(sigmaPtrs.size()); ++iSig)
  74       containerPtrs.push_back( new ProcessContainer(sigmaPtrs[iSig]) );
  75   }
  76
  77   // Append single container for Les Houches processes, if any.
  78   if (doLHA) {
  79     SigmaProcess* sigmaPtr = new SigmaLHAProcess();
  80     containerPtrs.push_back( new ProcessContainer(sigmaPtr) );
  81
  82     // Store location of this container, and send in LHA pointer.
  83     iLHACont = containerPtrs.size() - 1;
  84     containerPtrs[iLHACont]->setLHAPtr(lhaUpPtr);
  85   }
  86
  87   // If no processes found then refuse to do anything.
  88   if ( int(containerPtrs.size()) == 0) {
  89     infoPtr->errorMsg("Error in ProcessLevel::init: "
  90       "no process switched on");
  91     return false;
  92   }
  93
  94   // Initialize alphaStrong generation for SigmaProcess.
  95   double alphaSvalue  = Settings::parm("SigmaProcess:alphaSvalue");
  96   int    alphaSorder  = Settings::mode("SigmaProcess:alphaSorder");
  97   alphaS.init( alphaSvalue, alphaSorder);
  98
  99   // Initialize alphaEM generation for SigmaProcess.
 100   int    alphaEMorder = Settings::mode("SigmaProcess:alphaEMorder");
 101   alphaEM.init( alphaEMorder);
 102
 103   // Initialize each process.
 104   int numberOn = 0;
 105   for (int i = 0; i < int(containerPtrs.size()); ++i)
 106     if (containerPtrs[i]->init(infoPtr, beamAPtr, beamBPtr, &alphaS,
 107       &alphaEM, sigmaTotPtr, &resonanceDecays, slhaPtr, userHooksPtr))
 108       ++numberOn;
 109
 110   // Sum maxima for Monte Carlo choice.
 111   sigmaMaxSum = 0.;
 112   for (int i = 0; i < int(containerPtrs.size()); ++i)
 113     sigmaMaxSum += containerPtrs[i]->sigmaMax();
 114
 115   // Option to pick a second hard interaction: repeat as above.
 116   int number2On = 0;
 117   if (doSecondHard) {
 118     setupContainers.init2(container2Ptrs);
 119     if ( int(container2Ptrs.size()) == 0) {
 120       infoPtr->errorMsg("Error in ProcessLevel::init: "
 121         "no second hard process switched on");
 122       return false;
 123     }
 124     for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
 125       if (container2Ptrs[i2]->init(infoPtr, beamAPtr, beamBPtr, &alphaS,
 126         &alphaEM, sigmaTotPtr, &resonanceDecays, slhaPtr, userHooksPtr))
 127         ++number2On;
 128     sigma2MaxSum = 0.;
 129     for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
 130       sigma2MaxSum += container2Ptrs[i2]->sigmaMax();
 131   }
 132
 133   // Construct string with incoming beams and for cm energy.
 134   string collision = "We collide " + ParticleDataTable::name(idA)
 135     + " with " + ParticleDataTable::name(idB) + " at a CM energy of ";
 136   string pad( 51 - collision.length(), ' ');
 137
 138   // Print initialization information: header.
 139   os << "\n *-------  PYTHIA Process Initialization  ---------"
 140      << "-----------------*\n"
 141      << " |                                                   "
 142      << "               |\n"
 143      << " | " << collision << scientific << setprecision(3)<< setw(9) << eCM
 144      << " GeV" << pad << " |\n"
 145      << " |                                                   "
 146      << "               |\n"
 147      << " |---------------------------------------------------"
 148      << "---------------|\n"
 149      << " |                                                   "
 150      << " |             |\n"
 151      << " | Subprocess                                    Code"
 152      << " |   Estimated |\n"
 153      << " |                                                   "
 154      << " |    max (mb) |\n"
 155      << " |                                                   "
 156      << " |             |\n"
 157      << " |---------------------------------------------------"
 158      << "---------------|\n"
 159      << " |                                                   "
 160      << " |             |\n";
 161
 162
 163   // Loop over existing processes: print individual process info.
 164   for (int i = 0; i < int(containerPtrs.size()); ++i)
 165   os << " | " << left << setw(45) << containerPtrs[i]->name()
 166      << right << setw(5) << containerPtrs[i]->code() << " | "
 167      << scientific << setprecision(3) << setw(11)
 168      << containerPtrs[i]->sigmaMax() << " |\n";
 169
 170   // Loop over second hard processes, if any, and repeat as above.
 171   if (doSecondHard) {
 172     os << " |                                                   "
 173        << " |             |\n"
 174        << " |---------------------------------------------------"
 175        <<"---------------|\n"
 176        << " |                                                   "
 177        <<" |             |\n";
 178     for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
 179     os << " | " << left << setw(45) << container2Ptrs[i2]->name()
 180        << right << setw(5) << container2Ptrs[i2]->code() << " | "
 181        << scientific << setprecision(3) << setw(11)
 182        << container2Ptrs[i2]->sigmaMax() << " |\n";
 183   }
 184
 185   // Listing finished.
 186   os << " |                                                     "
 187      << "             |\n"
 188      << " *-------  End PYTHIA Process Initialization ----------"
 189      <<"-------------*" << endl;
 190
 191   // If sum of maxima vanishes then refuse to do anything.
 192   if ( numberOn == 0  || sigmaMaxSum <= 0.) {
 193     infoPtr->errorMsg("Error in ProcessLevel::init: "
 194       "all processes have vanishing cross sections");
 195     return false;
 196   }
 197   if ( doSecondHard && (number2On == 0  || sigma2MaxSum <= 0.) ) {
 198     infoPtr->errorMsg("Error in ProcessLevel::init: "
 199       "all second hard processes have vanishing cross sections");
 200     return false;
 201   }
 202
 203   // If two hard processes then check whether some (but not all) agree.
 204   allHardSame  = true;
 205   noneHardSame = true;
 206   if (doSecondHard) {
 207     bool foundMatch = false;
 208
 209     // Check for each first process if matched in second.
 210     for (int i = 0; i < int(containerPtrs.size()); ++i) {
 211       foundMatch = false;
 212       for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
 213         if (container2Ptrs[i2]->code() == containerPtrs[i]->code())
 214           foundMatch = true;
 215       containerPtrs[i]->isSame( foundMatch );
 216       if (!foundMatch)  allHardSame = false;
 217       if ( foundMatch) noneHardSame = false;
 218     }
 219
 220     // Check for each second process if matched in first.
 221     for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2) {
 222       foundMatch = false;
 223       for (int i = 0; i < int(containerPtrs.size()); ++i)
 224         if (containerPtrs[i]->code() == container2Ptrs[i2]->code())
 225           foundMatch = true;
 226       container2Ptrs[i2]->isSame( foundMatch );
 227       if (!foundMatch)  allHardSame = false;
 228       if ( foundMatch) noneHardSame = false;
 229     }
 230   }
 231   someHardSame = !allHardSame && !noneHardSame;
 232
 233   // Reset counters for average impact-parameter enhancement.
 234   nImpact       = 0;
 235   sumImpactFac  = 0.;
 236   sum2ImpactFac = 0.;
 237
 238   // Statistics for LHA events.
 239   codeLHA.resize(0);
 240   nEvtLHA.resize(0);
 241
 242   // Done.
 243   return true;
 244 }
 245
 246 //*********
 247
 248 // Main routine to generate the hard process.
 249
 250 bool ProcessLevel::next( Event& process) {
 251
 252   // Generate the next event with two or one hard interactions.
 253   bool physical = (doSecondHard) ? nextTwo( process) : nextOne( process);
 254
 255   // Check that colour assignments make sense.
 256   if (physical) physical = checkColours( process);
 257
 258   // Done.
 259   return physical;
 260 }
 261
 262 //*********
 263
 264 // Accumulate and update statistics (after possible user veto).
 265
 266 void ProcessLevel::accumulate() {
 267
 268   // Increase number of accepted events.
 269   containerPtrs[iContainer]->accumulate();
 270
 271   // Provide current generated cross section estimate.
 272   long   nTrySum    = 0;
 273   long   nSelSum    = 0;
 274   long   nAccSum    = 0;
 275   double sigmaSum   = 0.;
 276   double delta2Sum  = 0.;
 277   double sigSelSum  = 0.;
 278   for (int i = 0; i < int(containerPtrs.size()); ++i)
 279   if (containerPtrs[i]->sigmaMax() != 0.) {
 280     nTrySum        += containerPtrs[i]->nTried();
 281     nSelSum        += containerPtrs[i]->nSelected();
 282     nAccSum        += containerPtrs[i]->nAccepted();
 283     sigmaSum       += containerPtrs[i]->sigmaMC();
 284     delta2Sum      += pow2(containerPtrs[i]->deltaMC());
 285     sigSelSum      += containerPtrs[i]->sigmaSelMC();
 286   }
 287
 288   // For Les Houches events find subprocess type and update counter.
 289   if (infoPtr->isLHA()) {
 290     int codeLHANow = infoPtr->codeSub();
 291     int iFill = -1;
 292     for (int i = 0; i < int(codeLHA.size()); ++i)
 293       if (codeLHANow == codeLHA[i]) iFill = i;
 294     if (iFill >= 0) ++nEvtLHA[iFill];
 295
 296     // Add new process when new code and then arrange in order.
 297     else {
 298       codeLHA.push_back(codeLHANow);
 299       nEvtLHA.push_back(1);
 300       for (int i = int(codeLHA.size()) - 1; i > 0; --i) {
 301         if (codeLHA[i] < codeLHA[i - 1]) {
 302           swap(codeLHA[i], codeLHA[i - 1]);
 303           swap(nEvtLHA[i], nEvtLHA[i - 1]);
 304         }
 305         else break;
 306       }
 307     }
 308   }
 309
 310   // Normally only one hard interaction. Then store info and done.
 311   if (!doSecondHard) {
 312     double deltaSum = sqrtpos(delta2Sum);
 313     infoPtr->setSigma( nTrySum, nSelSum, nAccSum, sigmaSum, deltaSum);
 314     return;
 315   }
 316
 317   // Increase counter for a second hard interaction.
 318   container2Ptrs[i2Container]->accumulate();
 319
 320   // Update statistics on average impact factor.
 321   ++nImpact;
 322   sumImpactFac     += infoPtr->enhanceMI();
 323   sum2ImpactFac    += pow2(infoPtr->enhanceMI());
 324
 325   // Cross section estimate for second hard process.
 326   double sigma2Sum  = 0.;
 327   double sig2SelSum = 0.;
 328   for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
 329   if (container2Ptrs[i2]->sigmaMax() != 0.) {
 330     nTrySum        += container2Ptrs[i2]->nTried();
 331     sigma2Sum      += container2Ptrs[i2]->sigmaMC();
 332     sig2SelSum     += container2Ptrs[i2]->sigmaSelMC();
 333   }
 334
 335   // Average impact-parameter factor and error.
 336   double invN       = 1. / max(1, nImpact);
 337   double impactFac  = max( 1., sumImpactFac * invN);
 338   double impactErr2 = ( sum2ImpactFac * invN / pow2(impactFac) - 1.) * invN;
 339
 340   // Cross section estimate for combination of first and second process.
 341   // Combine two possible ways and take average.
 342   double sigmaComb  = 0.5 * (sigmaSum * sig2SelSum + sigSelSum * sigma2Sum);
 343   sigmaComb        *= impactFac / sigmaND;
 344   if (allHardSame) sigmaComb *= 0.5;
 345   double deltaComb  = sqrtpos(2. / nAccSum + impactErr2) * sigmaComb;
 346
 347   // Store info and done.
 348   infoPtr->setSigma( nTrySum, nSelSum, nAccSum, sigmaComb, deltaComb);
 349
 350 }
 351
 352 //*********
 353
 354 // Print statistics on cross sections and number of events.
 355
 356 void ProcessLevel::statistics(bool reset, ostream& os) {
 357
 358   // Special processing if two hard interactions selected.
 359   if (doSecondHard) {
 360     statistics2(reset, os);
 361     return;
 362   }
 363
 364   // Header.
 365   os << "\n *-------  PYTHIA Event and Cross Section Statistics  ------"
 366      << "-------------------------------------------------------*\n"
 367      << " |                                                            "
 368      << "                                                     |\n"
 369      << " | Subprocess                                    Code |       "
 370      << "     Number of events       |      sigma +- delta    |\n"
 371      << " |                                                    |       "
 372      << "Tried   Selected   Accepted |     (estimated) (mb)   |\n"
 373      << " |                                                    |       "
 374      << "                            |                        |\n"
 375      << " |------------------------------------------------------------"
 376      << "-----------------------------------------------------|\n"
 377      << " |                                                    |       "
 378      << "                            |                        |\n";
 379
 380   // Reset sum counters.
 381   long   nTrySum   = 0;
 382   long   nSelSum   = 0;
 383   long   nAccSum   = 0;
 384   double sigmaSum  = 0.;
 385   double delta2Sum = 0.;
 386
 387   // Loop over existing processes.
 388   for (int i = 0; i < int(containerPtrs.size()); ++i)
 389   if (containerPtrs[i]->sigmaMax() != 0.) {
 390
 391     // Read info for process. Sum counters.
 392     long   nTry    = containerPtrs[i]->nTried();
 393     long   nSel    = containerPtrs[i]->nSelected();
 394     long   nAcc    = containerPtrs[i]->nAccepted();
 395     double sigma   = containerPtrs[i]->sigmaMC();
 396     double delta   = containerPtrs[i]->deltaMC();
 397     nTrySum       += nTry;
 398     nSelSum       += nSel;
 399     nAccSum       += nAcc;
 400     sigmaSum      += sigma;
 401     delta2Sum     += pow2(delta);
 402
 403     // Print individual process info.
 404     os << " | " << left << setw(45) << containerPtrs[i]->name()
 405        << right << setw(5) << containerPtrs[i]->code() << " | "
 406        << setw(11) << nTry << " " << setw(10) << nSel << " "
 407        << setw(10) << nAcc << " | " << scientific << setprecision(3)
 408        << setw(11) << sigma << setw(11) << delta << " |\n";
 409
 410     // Print subdivision by user code for Les Houches process.
 411     if (containerPtrs[i]->code() == 9999)
 412     for (int j = 0; j < int(codeLHA.size()); ++j)
 413       os << " |    ... whereof user classification code " << setw(10)
 414          << codeLHA[j] << " |                        " << setw(10)
 415          << nEvtLHA[j] << " |                        | \n";
 416   }
 417
 418   // Print summed process info.
 419   os << " |                                                    |       "
 420      << "                            |                        |\n"
 421      << " | " << left << setw(50) << "sum" << right << " | " << setw(11)
 422      << nTrySum << " " << setw(10) << nSelSum << " " << setw(10)
 423      << nAccSum << " | " << scientific << setprecision(3) << setw(11)
 424      << sigmaSum << setw(11) << sqrtpos(delta2Sum) << " |\n";
 425
 426   // Listing finished.
 427   os << " |                                                            "
 428      << "                                                     |\n"
 429      << " *-------  End PYTHIA Event and Cross Section Statistics -----"
 430      << "-----------------------------------------------------*" << endl;
 431
 432   // Optionally reset statistics contants.
 433   if (reset) for (int i = 0; i < int(containerPtrs.size()); ++i)
 434     containerPtrs[i]->reset();
 435
 436 }
 437
 438 //*********
 439
 440 // Initialize SUSY Les Houches Accord data.
 441
 442 bool ProcessLevel::initSLHA() {
 443
 444   // Check whether SUSY is on.
 445   if ( !Settings::flag("SUSY") ) return true;
 446
 447   // Read SUSY Les Houches Accord File.
 448   string slhaFile = Settings::word("SUSY:SusyLesHouchesFile");
 449   int ifail = slhaPtr->readFile(slhaFile);
 450
 451   // In case of problems, print error and fail init.
 452   if (ifail != 0) {
 453     infoPtr->errorMsg("Error from Pythia::initSLHA: "
 454       "problem reading SLHA file", slhaFile);
 455     return false;
 456   };
 457
 458   // Update particle data.
 459   int id = slhaPtr->mass.first();
 460   for (int i = 1; i <= slhaPtr->mass.size() ; i++) {
 461     double mass = abs(slhaPtr->mass(id));
 462     ParticleDataTable::m0(id,mass);
 463     id = slhaPtr->mass.next();
 464   };
 465
 466   // Check spectrum for consistency. Switch off SUSY if necessary.
 467   ifail = slhaPtr->checkSpectrum();
 468   if (ifail != 0) {
 469     infoPtr->errorMsg("Warning from Pythia::initSLHA: "
 470       "Problem with SLHA spectrum.",
 471       "\n Only using masses and switching off SUSY.");
 472     Settings::flag("SUSY", false);
 473     slhaPtr->printSpectrum();
 474     return true;
 475   };
 476
 477   // Print spectrum. Done.
 478   slhaPtr->printSpectrum();
 479   return true;
 480
 481 }
 482
 483 //*********
 484
 485 // Generate the next event with one interaction.
 486
 487 bool ProcessLevel::nextOne( Event& process) {
 488
 489   // Update CM energy for phase space selection.
 490   double eCM = infoPtr->eCM();
 491   for (int i = 0; i < int(containerPtrs.size()); ++i)
 492     containerPtrs[i]->newECM(eCM);
 493
 494   // Loop over tries until trial event succeeds.
 495   for ( ; ; ) {
 496
 497     // Pick one of the subprocesses.
 498     double sigmaMaxNow = sigmaMaxSum * Rndm::flat();
 499     int iMax = containerPtrs.size() - 1;
 500     iContainer = -1;
 501     do sigmaMaxNow -= containerPtrs[++iContainer]->sigmaMax();
 502     while (sigmaMaxNow > 0. && iContainer < iMax);
 503
 504     // Do a trial event of this subprocess; accept or not.
 505     if (containerPtrs[iContainer]->trialProcess()) break;
 506
 507     // Check for end-of-file condition for Les Houches events.
 508     if (infoPtr->atEndOfFile()) return false;
 509   }
 510
 511   // Update sum of maxima if current maximum violated.
 512   if (containerPtrs[iContainer]->newSigmaMax()) {
 513     sigmaMaxSum = 0.;
 514     for (int i = 0; i < int(containerPtrs.size()); ++i)
 515       sigmaMaxSum += containerPtrs[i]->sigmaMax();
 516   }
 517
 518   // Construct kinematics of acceptable process.
 519   if ( !containerPtrs[iContainer]->constructProcess( process) ) return false;
 520
 521   // Do all resonance decays.
 522   if ( doResDecays && !containerPtrs[iContainer]->decayResonances(
 523     process) ) return false;
 524
 525   // Add any junctions to the process event record list.
 526   findJunctions( process);
 527
 528   // Done.
 529   return true;
 530 }
 531
 532 //*********
 533
 534 // Generate the next event with two hard interactions.
 535
 536 bool ProcessLevel::nextTwo( Event& process) {
 537
 538   // Update CM energy for phase space selection.
 539   double eCM = infoPtr->eCM();
 540   for (int i = 0; i < int(containerPtrs.size()); ++i)
 541     containerPtrs[i]->newECM(eCM);
 542   for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
 543     container2Ptrs[i2]->newECM(eCM);
 544
 545   // Loop over both hard processes to find consistent common kinematics.
 546   for ( ; ; ) {
 547
 548     // Loop internally over tries for hardest process until succeeds.
 549     for ( ; ; ) {
 550
 551       // Pick one of the subprocesses.
 552       double sigmaMaxNow = sigmaMaxSum * Rndm::flat();
 553       int iMax = containerPtrs.size() - 1;
 554       iContainer = -1;
 555       do sigmaMaxNow -= containerPtrs[++iContainer]->sigmaMax();
 556       while (sigmaMaxNow > 0. && iContainer < iMax);
 557
 558       // Do a trial event of this subprocess; accept or not.
 559       if (containerPtrs[iContainer]->trialProcess()) break;
 560
 561       // Check for end-of-file condition for Les Houches events.
 562       if (infoPtr->atEndOfFile()) return false;
 563     }
 564
 565     // Update sum of maxima if current maximum violated.
 566     if (containerPtrs[iContainer]->newSigmaMax()) {
 567       sigmaMaxSum = 0.;
 568       for (int i = 0; i < int(containerPtrs.size()); ++i)
 569         sigmaMaxSum += containerPtrs[i]->sigmaMax();
 570     }
 571
 572     // Loop internally over tries for second hardest process until succeeds.
 573     for ( ; ; ) {
 574
 575       // Pick one of the subprocesses.
 576       double sigma2MaxNow = sigma2MaxSum * Rndm::flat();
 577       int i2Max = container2Ptrs.size() - 1;
 578       i2Container = -1;
 579       do sigma2MaxNow -= container2Ptrs[++i2Container]->sigmaMax();
 580       while (sigma2MaxNow > 0. && i2Container < i2Max);
 581
 582       // Do a trial event of this subprocess; accept or not.
 583       if (container2Ptrs[i2Container]->trialProcess()) break;
 584     }
 585
 586     // Update sum of maxima if current maximum violated.
 587     if (container2Ptrs[i2Container]->newSigmaMax()) {
 588       sigma2MaxSum = 0.;
 589       for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
 590         sigma2MaxSum += container2Ptrs[i2]->sigmaMax();
 591     }
 592
 593     // Check whether common set of x values is kinematically possible.
 594     double xA1      = containerPtrs[iContainer]->x1();
 595     double xB1      = containerPtrs[iContainer]->x2();
 596     double xA2      = container2Ptrs[i2Container]->x1();
 597     double xB2      = container2Ptrs[i2Container]->x2();
 598     if (xA1 + xA2 >= 1. || xB1 + xB2 >= 1.) continue;
 599
 600     // Reset beam contents. Naive parton densities for second interaction.
 601     // (Subsequent procedure could be symmetrized, but would be overkill.)
 602     beamAPtr->clear();
 603     beamBPtr->clear();
 604     int    idA1     = containerPtrs[iContainer]->id1();
 605     int    idB1     = containerPtrs[iContainer]->id2();
 606     int    idA2     = container2Ptrs[i2Container]->id1();
 607     int    idB2     = container2Ptrs[i2Container]->id2();
 608     double Q2Fac1   = containerPtrs[iContainer]->Q2Fac();
 609     double Q2Fac2   = container2Ptrs[i2Container]->Q2Fac();
 610     double pdfA2Raw = beamAPtr->xf( idA2, xA2,Q2Fac2);
 611     double pdfB2Raw = beamBPtr->xf( idB2, xB2,Q2Fac2);
 612
 613     // Remove partons in first interaction from beams.
 614     beamAPtr->append( 3, idA1, xA1);
 615     beamAPtr->xfISR( 0, idA1, xA1, Q2Fac1);
 616     beamAPtr->pickValSeaComp();
 617     beamBPtr->append( 4, idB1, xB1);
 618     beamBPtr->xfISR( 0, idB1, xB1, Q2Fac1);
 619     beamBPtr->pickValSeaComp();
 620
 621     // Reevaluate pdf's for second interaction and weight by reduction.
 622     double pdfA2Mod = beamAPtr->xfMI( idA2, xA2,Q2Fac2);
 623     double pdfB2Mod = beamBPtr->xfMI( idB2, xB2,Q2Fac2);
 624     double wtPdfMod = (pdfA2Mod * pdfB2Mod) / (pdfA2Raw * pdfB2Raw);
 625     if (wtPdfMod < Rndm::flat()) continue;
 626
 627     // Reduce by a factor of 2 for identical processes when others not.
 628     if ( someHardSame && containerPtrs[iContainer]->isSame()
 629       && container2Ptrs[i2Container]->isSame() && Rndm::flat() > 0.5)
 630       continue;
 631
 632     // If come this far then acceptable event.
 633     break;
 634   }
 635
 636   // Construct kinematics of acceptable processes.
 637   Event process2;
 638   process2.initColTag();
 639   startColTag2 = process2.lastColTag();
 640   if ( !containerPtrs[iContainer]->constructProcess( process) ) return false;
 641   if ( !container2Ptrs[i2Container]->constructProcess( process2, false) )
 642     return false;
 643
 644   // Do all resonance decays.
 645   if ( doResDecays &&  !containerPtrs[iContainer]->decayResonances(
 646     process) ) return false;
 647   if ( doResDecays &&  !container2Ptrs[i2Container]->decayResonances(
 648     process2) ) return false;
 649
 650   // Append second hard interaction to normal process object.
 651   combineProcessRecords( process, process2);
 652
 653   // Add any junctions to the process event record list.
 654   findJunctions( process);
 655
 656   // Done.
 657   return true;
 658 }
 659
 660 //*********
 661
 662 // Append second hard interaction to normal process object.
 663 // Complication: all resonance decay chains must be put at the end.
 664
 665 void ProcessLevel::combineProcessRecords( Event& process, Event& process2) {
 666
 667   // Find first event record size, excluding resonances.
 668   int nSize = process.size();
 669   int nHard = 5;
 670   while (nHard < nSize && process[nHard].mother1() == 3) ++nHard;
 671
 672   // Save resonance products temporarily elsewhere.
 673   vector<Particle> resProd;
 674   if (nSize > nHard) {
 675     for (int i = nHard; i < nSize; ++i) resProd.push_back( process[i] );
 676     process.popBack(nSize - nHard);
 677   }
 678
 679   // Find second event record size, excluding resonances.
 680   int nSize2 = process2.size();
 681   int nHard2 = 5;
 682   while (nHard2 < nSize2 && process2[nHard2].mother1() == 3) ++nHard2;
 683
 684   // Find amount of necessary position and colour offset for second process.
 685   int addPos  = nHard  - 3;
 686   int addCol  = process.lastColTag() - startColTag2;
 687
 688   // Loop over all particles (except beams) from second process.
 689   for (int i = 3; i < nSize2; ++i) {
 690
 691     // Offset mother and daughter pointers and colour tags of particle.
 692     process2[i].offsetHistory( 2, addPos, 2, addPos);
 693     process2[i].offsetCol( addCol);
 694
 695     // Append hard-process particles from process2 to process.
 696     if (i < nHard2) process.append( process2[i] );
 697   }
 698
 699   // Reinsert resonance decay chains of first hard process.
 700   int addPos2 = nHard2 - 3;
 701   if (nHard < nSize) {
 702
 703     // Offset daughter pointers of unmoved mothers.
 704     for (int i = 5; i < nHard; ++i)
 705       process[i].offsetHistory( 0, 0, nHard - 1, addPos2);
 706
 707     // Modify history of resonance products when restoring.
 708     for (int i = 0; i < int(resProd.size()); ++i) {
 709       resProd[i].offsetHistory( nHard - 1, addPos2, nHard - 1, addPos2);
 710       process.append( resProd[i] );
 711     }
 712   }
 713
 714   // Insert resonance decay chains of second hard process.
 715   if (nHard2 < nSize2) {
 716     int nHard3  = nHard + nHard2 - 3;
 717     int addPos3 = nSize - nHard;
 718
 719     // Offset daughter pointers of second-process mothers.
 720     for (int i = nHard + 2; i < nHard3; ++i)
 721       process[i].offsetHistory( 0, 0, nHard3 - 1, addPos3);
 722
 723     // Modify history of second-process resonance products and insert.
 724     for (int i = nHard2; i < nSize2; ++i) {
 725       process2[i].offsetHistory( nHard3 - 1, addPos3, nHard3 - 1, addPos3);
 726       process.append( process2[i] );
 727     }
 728   }
 729
 730   // Store PDF scale for second interaction.
 731   process.scaleSecond( process2.scale() );
 732
 733 }
 734
 735 //*********
 736
 737 // Add any junctions to the process event record list.
 738 // First try, so still incomplete. ??
 739 // Also check that do not doublebook if called repeatedly.
 740
 741 void ProcessLevel::findJunctions( Event& junEvent) {
 742
 743   // Loop though event; isolate all uncoloured particles.
 744   for (int i = 0; i < junEvent.size(); ++i)
 745   if ( junEvent[i].col() == 0 && junEvent[i].acol() == 0) {
 746
 747     // Find all daughters and store daughter colours and anticolours.
 748     vector<int> daughters = junEvent.daughterList(i);
 749     vector<int> cols, acols;
 750     for (int j = 0; j < int(daughters.size()); ++j) {
 751       int colDau  = junEvent[ daughters[j] ].col();
 752       int acolDau = junEvent[ daughters[j] ].acol();
 753       if (colDau > 0)  cols.push_back( colDau);
 754       if (acolDau > 0) acols.push_back( acolDau);
 755     }
 756
 757     // Remove all matching colour-anticolour pairs.
 758     bool foundPair = true;
 759     while (foundPair && cols.size() > 0 && acols.size() > 0) {
 760       foundPair = false;
 761       for (int j = 0; j < int(cols.size()); ++j) {
 762         for (int k = 0; k < int(acols.size()); ++k) {
 763           if (acols[k] == cols[j]) {
 764             cols[j]  = cols.back();
 765             cols.pop_back();
 766             acols[k] = acols.back();
 767             acols.pop_back();
 768             foundPair = true;
 769             break;
 770           }
 771         } if (foundPair) break;
 772       }
 773     }
 774
 775     // Store an (anti)junction when three (anti)coloured daughters.
 776     if (cols.size() == 3 && acols.size() == 0)
 777       junEvent.appendJunction( 1, cols[0], cols[1], cols[2]);
 778     if (acols.size() == 3 && cols.size() == 0)
 779       junEvent.appendJunction( 2, acols[0], acols[1], acols[2]);
 780   }
 781
 782   // Done.
 783 }
 784
 785 //*********
 786
 787 // Check that colours match up.
 788
 789 bool ProcessLevel::checkColours( Event& process) {
 790
 791   // Variables and arrays for common usage.
 792   bool physical = true;
 793   bool match;
 794   int colType, col, acol, iPos, iNow, iNowA;
 795   vector<int> colTags, colPos, acolPos;
 796
 797   // Check that each particle has the kind of colours expected of it.
 798   for (int i = 0; i < process.size(); ++i) {
 799     colType = process[i].colType();
 800     col     = process[i].col();
 801     acol    = process[i].acol();
 802     if      (colType ==  0 && (col != 0 || acol != 0)) physical = false;
 803     else if (colType ==  1 && (col <= 0 || acol != 0)) physical = false;
 804     else if (colType == -1 && (col != 0 || acol <= 0)) physical = false;
 805     else if (colType ==  2 && (col <= 0 || acol <= 0)) physical = false;
 806     else if (colType < -1 || colType > 2)              physical = false;
 807
 808     // Add to the list of colour tags.
 809     if (col > 0) {
 810       match = false;
 811       for (int ic = 0; ic < int(colTags.size()) ; ++ic)
 812         if (col == colTags[ic]) match = true;
 813       if (!match) colTags.push_back(col);
 814     } else if (acol > 0) {
 815       match = false;
 816       for (int ic = 0; ic < int(colTags.size()) ; ++ic)
 817         if (acol == colTags[ic]) match = true;
 818       if (!match) colTags.push_back(acol);
 819     }
 820   }
 821
 822   // Warn and give up if particles did not have the expected colours.
 823   if (!physical) {
 824     infoPtr->errorMsg("Error in ProcessLevel::checkColours: "
 825       "incorrect colour assignment");
 826     return false;
 827   }
 828
 829   // Loop through all colour tags and find their positions (by sign).
 830   for (int ic = 0; ic < int(colTags.size()); ++ic) {
 831     col = colTags[ic];
 832     colPos.resize(0);
 833     acolPos.resize(0);
 834     for (int i = 0; i < process.size(); ++i) {
 835       if (process[i].col() == col) colPos.push_back(i);
 836       if (process[i].acol() == col) acolPos.push_back(i);
 837     }
 838
 839     // Trace colours back through decays; remove daughters.
 840     while (colPos.size() > 1) {
 841       iPos = colPos.size() - 1;
 842       iNow = colPos[iPos];
 843       if ( process[iNow].mother1() == colPos[iPos - 1]
 844         && process[iNow].mother2() == 0) colPos.pop_back();
 845       else break;
 846     }
 847     while (acolPos.size() > 1) {
 848       iPos = acolPos.size() - 1;
 849       iNow = acolPos[iPos];
 850       if ( process[iNow].mother1() == acolPos[iPos - 1]
 851         && process[iNow].mother2() == 0) acolPos.pop_back();
 852       else break;
 853     }
 854
 855     // Now colour should exist in only 2 copies.
 856     if (colPos.size() + acolPos.size() != 2) physical = false;
 857
 858     // If both colours or both anticolours then one mother of the other.
 859     else if (colPos.size() == 2) {
 860       iNow = colPos[1];
 861       if ( process[iNow].mother1() != colPos[0]
 862         && process[iNow].mother2() != colPos[0] ) physical = false;
 863     }
 864     else if (acolPos.size() == 2) {
 865       iNowA = acolPos[1];
 866       if ( process[iNowA].mother1() != acolPos[0]
 867         && process[iNowA].mother2() != acolPos[0] ) physical = false;
 868     }
 869
 870     // If one of each then should have same mother(s), or point to beams.
 871     else {
 872       iNow  = colPos[0];
 873       iNowA = acolPos[0];
 874       if ( process[iNow].status() == -21 &&  process[iNowA].status() == -21 );
 875       else if ( (process[iNow].mother1() != process[iNowA].mother1())
 876              || (process[iNow].mother2() != process[iNowA].mother2()) )
 877              physical = false;
 878     }
 879
 880   }
 881
 882   // Error message if problem found. Done.
 883   if (!physical) infoPtr->errorMsg("Error in ProcessLevel::checkColours: "
 884     "unphysical colour flow");
 885   return physical;
 886
 887 }
 888
 889 //*********
 890
 891 // Print statistics when two hard processes allowed.
 892
 893 void ProcessLevel::statistics2(bool reset, ostream& os) {
 894
 895   // Average impact-parameter factor and error.
 896   double invN          = 1. / max(1, nImpact);
 897   double impactFac     = max( 1., sumImpactFac * invN);
 898   double impactErr2    = ( sum2ImpactFac * invN / pow2(impactFac) - 1.) * invN;
 899
 900   // Derive scaling factor to be applied to first set of processes.
 901   double sigma2SelSum  = 0.;
 902   int    n2SelSum      = 0;
 903   for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2) {
 904     sigma2SelSum      += container2Ptrs[i2]->sigmaSelMC();
 905     n2SelSum          += container2Ptrs[i2]->nSelected();
 906   }
 907   double factor1       = impactFac * sigma2SelSum / sigmaND;
 908   double rel1Err       = sqrt(1. / max(1, n2SelSum) + impactErr2);
 909   if (allHardSame) factor1 *= 0.5;
 910
 911   // Derive scaling factor to be applied to second set of processes.
 912   double sigma1SelSum  = 0.;
 913   int    n1SelSum      = 0;
 914   for (int i = 0; i < int(containerPtrs.size()); ++i) {
 915     sigma1SelSum      += containerPtrs[i]->sigmaSelMC();
 916     n1SelSum          += containerPtrs[i]->nSelected();
 917   }
 918   double factor2       = impactFac * sigma1SelSum / sigmaND;
 919   if (allHardSame) factor2 *= 0.5;
 920   double rel2Err       = sqrt(1. / max(1, n1SelSum) + impactErr2);
 921
 922   // Header.
 923   os << "\n *-------  PYTHIA Event and Cross Section Statistics  ------"
 924      << "--------------------------------------------------*\n"
 925      << " |                                                            "
 926      << "                                                |\n"
 927      << " | Subprocess                               Code |            "
 928      << "Number of events       |      sigma +- delta    |\n"
 929      << " |                                               |       Tried"
 930      << "   Selected   Accepted |     (estimated) (mb)   |\n"
 931      << " |                                               |            "
 932      << "                       |                        |\n"
 933      << " |------------------------------------------------------------"
 934      << "------------------------------------------------|\n"
 935      << " |                                               |            "
 936      << "                       |                        |\n"
 937      << " | First hard process:                           |            "
 938      << "                       |                        |\n"
 939      << " |                                               |            "
 940      << "                       |                        |\n";
 941
 942   // Reset sum counters.
 943   long   nTrySum   = 0;
 944   long   nSelSum   = 0;
 945   long   nAccSum   = 0;
 946   double sigmaSum  = 0.;
 947   double delta2Sum = 0.;
 948
 949   // Loop over existing first processes.
 950   for (int i = 0; i < int(containerPtrs.size()); ++i)
 951   if (containerPtrs[i]->sigmaMax() != 0.) {
 952
 953     // Read info for process. Sum counters.
 954     long   nTry    = containerPtrs[i]->nTried();
 955     long   nSel    = containerPtrs[i]->nSelected();
 956     long   nAcc    = containerPtrs[i]->nAccepted();
 957     double sigma   = containerPtrs[i]->sigmaMC() * factor1;
 958     double delta2  = pow2( containerPtrs[i]->deltaMC() * factor1 );
 959     nTrySum       += nTry;
 960     nSelSum       += nSel;
 961     nAccSum       += nAcc;
 962     sigmaSum      += sigma;
 963     delta2Sum     += delta2;
 964     delta2        += pow2( sigma * rel1Err );
 965
 966     // Print individual process info.
 967     os << " | " << left << setw(40) << containerPtrs[i]->name()
 968        << right << setw(5) << containerPtrs[i]->code() << " | "
 969        << setw(11) << nTry << " " << setw(10) << nSel << " "
 970        << setw(10) << nAcc << " | " << scientific << setprecision(3)
 971        << setw(11) << sigma << setw(11) << sqrtpos(delta2) << " |\n";
 972   }
 973
 974   // Print summed info for first processes.
 975   delta2Sum       += pow2( sigmaSum * rel1Err );
 976   os << " |                                               |            "
 977      << "                       |                        |\n"
 978      << " | " << left << setw(45) << "sum" << right << " | " << setw(11)
 979      << nTrySum << " " << setw(10) << nSelSum << " " << setw(10)
 980      << nAccSum << " | " << scientific << setprecision(3) << setw(11)
 981      << sigmaSum << setw(11) << sqrtpos(delta2Sum) << " |\n";
 982
 983
 984   // Separation lines to second hard processes.
 985   os << " |                                               |            "
 986      << "                       |                        |\n"
 987      << " |------------------------------------------------------------"
 988      << "------------------------------------------------|\n"
 989      << " |                                               |            "
 990      << "                       |                        |\n"
 991      << " | Second hard process:                          |            "
 992      << "                       |                        |\n"
 993      << " |                                               |            "
 994      << "                       |                        |\n";
 995
 996   // Reset sum counters.
 997   nTrySum   = 0;
 998   nSelSum   = 0;
 999   nAccSum   = 0;
1000   sigmaSum  = 0.;
1001   delta2Sum = 0.;
1002
1003   // Loop over existing second processes.
1004   for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
1005   if (container2Ptrs[i2]->sigmaMax() != 0.) {
1006
1007     // Read info for process. Sum counters.
1008     long   nTry    = container2Ptrs[i2]->nTried();
1009     long   nSel    = container2Ptrs[i2]->nSelected();
1010     long   nAcc    = container2Ptrs[i2]->nAccepted();
1011     double sigma   = container2Ptrs[i2]->sigmaMC() * factor2;
1012     double delta2  = pow2( container2Ptrs[i2]->deltaMC() * factor2 );
1013     nTrySum       += nTry;
1014     nSelSum       += nSel;
1015     nAccSum       += nAcc;
1016     sigmaSum      += sigma;
1017     delta2Sum     += delta2;
1018     delta2        += pow2( sigma * rel2Err );
1019
1020     // Print individual process info.
1021     os << " | " << left << setw(40) << container2Ptrs[i2]->name()
1022        << right << setw(5) << container2Ptrs[i2]->code() << " | "
1023        << setw(11) << nTry << " " << setw(10) << nSel << " "
1024        << setw(10) << nAcc << " | " << scientific << setprecision(3)
1025        << setw(11) << sigma << setw(11) << sqrtpos(delta2) << " |\n";
1026   }
1027
1028   // Print summed info for second processes.
1029   delta2Sum       += pow2( sigmaSum * rel2Err );
1030   os << " |                                               |            "
1031      << "                       |                        |\n"
1032      << " | " << left << setw(45) << "sum" << right << " | " << setw(11)
1033      << nTrySum << " " << setw(10) << nSelSum << " " << setw(10)
1034      << nAccSum << " | " << scientific << setprecision(3) << setw(11)
1035      << sigmaSum << setw(11) << sqrtpos(delta2Sum) << " |\n";
1036
1037   // Print information on how the two processes were combined.
1038   os << " |                                               |            "
1039      << "                       |                        |\n"
1040      << " |------------------------------------------------------------"
1041      << "------------------------------------------------|\n"
1042      << " |                                                            "
1043      << "                                                |\n"
1044      << " | Uncombined cross sections for the two event sets were "
1045      << setw(10) << sigma1SelSum << " and " << sigma2SelSum << " mb, "
1046      << "respectively, combined  |\n"
1047      << " | using a sigma(nonDiffractive) of " << setw(10) << sigmaND
1048      << " mb and an impact-parameter enhancement factor of "
1049      << setw(10) << impactFac << ".   |\n";
1050
1051   // Listing finished.
1052   os << " |                                                            "
1053      << "                                                |\n"
1054      << " *-------  End PYTHIA Event and Cross Section Statistics -----"
1055      << "------------------------------------------------*" << endl;
1056
1057   // Optionally reset statistics contants.
1058   if (reset) {
1059     for (int i = 0; i < int(containerPtrs.size()); ++i)
1060       containerPtrs[i]->reset();
1061     for (int i2 = 0; i2 < int(container2Ptrs.size()); ++i2)
1062       container2Ptrs[i2]->reset();
1063   }
1064
1065 }
1066
1067 //**************************************************************************
1068
1069 } // end namespace Pythia8