PYTHIA8/pythia8170/src/Pythia.cxx

   1 // Pythia.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 Pythia class.
   7
   8 #include "Pythia.h"
   9
  10 // Access time information.
  11 #include <ctime>
  12
  13 // Allow string and character manipulation.
  14 #include <cctype>
  15
  16 namespace Pythia8 {
  17
  18 //==========================================================================
  19
  20 // The Pythia class.
  21
  22 //--------------------------------------------------------------------------
  23
  24 // The current Pythia (sub)version number, to agree with XML version.
  25 const double Pythia::VERSIONNUMBERCODE = 8.170;
  26
  27 //--------------------------------------------------------------------------
  28
  29 // Constants: could be changed here if desired, but normally should not.
  30 // These are of technical nature, as described for each.
  31
  32 // Maximum number of tries to produce parton level from given input.
  33 const int Pythia::NTRY          = 10;
  34
  35 // Negative integer to denote that no subrun has been set.
  36 const int Pythia::SUBRUNDEFAULT = -999;
  37
  38 //--------------------------------------------------------------------------
  39
  40 // Constructor.
  41
  42 Pythia::Pythia(string xmlDir) {
  43
  44   // Initial values for pointers to PDF's.
  45   useNewPdfA      = false;
  46   useNewPdfB      = false;
  47   useNewPdfHard   = false;
  48   useNewPdfPomA   = false;
  49   useNewPdfPomB   = false;
  50   pdfAPtr         = 0;
  51   pdfBPtr         = 0;
  52   pdfHardAPtr     = 0;
  53   pdfHardBPtr     = 0;
  54   pdfPomAPtr      = 0;
  55   pdfPomBPtr      = 0;
  56
  57   // Initial values for pointers to Les Houches Event objects.
  58   doLHA           = false;
  59   useNewLHA       = false;
  60   lhaUpPtr        = 0;
  61
  62   //Initial value for couplings pointer
  63   couplingsPtr    = &couplings;
  64
  65   // Initial value for pointer to external decay handler.
  66   decayHandlePtr  = 0;
  67
  68   // Initial value for pointer to user hooks.
  69   userHooksPtr    = 0;
  70
  71   // Initial value for pointer to merging hooks.
  72   doMerging          = false;
  73   hasMergingHooks    = false;
  74   hasOwnMergingHooks = false;
  75   mergingHooksPtr    = 0;
  76
  77   // Initial value for pointer to beam shape.
  78   useNewBeamShape = false;
  79   beamShapePtr    = 0;
  80
  81   // Initial values for pointers to timelike and spacelike showers.
  82   useNewTimes     = false;
  83   useNewSpace     = false;
  84   timesDecPtr     = 0;
  85   timesPtr        = 0;
  86   spacePtr        = 0;
  87
  88   // Find path to data files, i.e. xmldoc directory location.
  89   // Environment variable takes precedence, else use constructor input.
  90   xmlPath     = "";
  91   const char* PYTHIA8DATA = "PYTHIA8DATA";
  92   char* envPath = getenv(PYTHIA8DATA);
  93   if (envPath != 0 && *envPath != '\0') {
  94     int i = 0;
  95     while (*(envPath+i) != '\0') xmlPath += *(envPath+(i++));
  96   }
  97   else xmlPath = xmlDir;
  98   if (xmlPath[ xmlPath.length() - 1 ] != '/') xmlPath += "/";
  99
 100   // Read in files with all flags, modes, parms and words.
 101   settings.initPtr( &info);
 102   string initFile = xmlPath + "Index.xml";
 103   isConstructed = settings.init( initFile);
 104   if (!isConstructed) {
 105     info.errorMsg("Abort from Pythia::Pythia: settings unavailable");
 106     return;
 107   }
 108
 109   // Check that XML version number matches code version number.
 110   double versionNumberXML = parm("Pythia:versionNumber");
 111   isConstructed = (abs(versionNumberXML - VERSIONNUMBERCODE) < 0.0005);
 112   if (!isConstructed) {
 113     ostringstream errCode;
 114     errCode << fixed << setprecision(3) << ": in code " << VERSIONNUMBERCODE
 115             << " but in XML " << versionNumberXML;
 116     info.errorMsg("Abort from Pythia::Pythia: unmatched version numbers",
 117       errCode.str());
 118     return;
 119   }
 120
 121   // Read in files with all particle data.
 122   particleData.initPtr( &info, &settings, &rndm, couplingsPtr);
 123   string dataFile = xmlPath + "ParticleData.xml";
 124   isConstructed = particleData.init( dataFile);
 125   if (!isConstructed) {
 126     info.errorMsg("Abort from Pythia::Pythia: particle data unavailable");
 127     return;
 128   }
 129
 130   // Write the Pythia banner to output.
 131   banner();
 132
 133   // Not initialized until at the end of the init() call.
 134   isInit = false;
 135   info.addCounter(0);
 136
 137 }
 138
 139 //--------------------------------------------------------------------------
 140
 141 // Destructor.
 142
 143 Pythia::~Pythia() {
 144
 145   // Delete the PDF's created with new.
 146   if (useNewPdfHard && pdfHardAPtr != pdfAPtr) delete pdfHardAPtr;
 147   if (useNewPdfHard && pdfHardBPtr != pdfBPtr) delete pdfHardBPtr;
 148   if (useNewPdfA) delete pdfAPtr;
 149   if (useNewPdfB) delete pdfBPtr;
 150   if (useNewPdfPomA) delete pdfPomAPtr;
 151   if (useNewPdfPomB) delete pdfPomBPtr;
 152
 153   // Delete the Les Houches object created with new.
 154   if (useNewLHA) delete lhaUpPtr;
 155
 156   // Delete the MergingHooks object created with new.
 157   if (hasOwnMergingHooks) delete mergingHooksPtr;
 158
 159   // Delete the BeamShape object created with new.
 160   if (useNewBeamShape) delete beamShapePtr;
 161
 162   // Delete the timelike and spacelike showers created with new.
 163   if (useNewTimes) delete timesPtr;
 164   if (useNewSpace) delete spacePtr;
 165
 166 }
 167
 168 //--------------------------------------------------------------------------
 169
 170 // Read in one update for a setting or particle data from a single line.
 171
 172 bool Pythia::readString(string line, bool warn) {
 173
 174   // Check that constructor worked.
 175   if (!isConstructed) return false;
 176
 177   // If empty line then done.
 178   if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos) return true;
 179
 180   // If first character is not a letter/digit, then taken to be a comment.
 181   int firstChar = line.find_first_not_of(" \n\t\v\b\r\f\a");
 182   if (!isalnum(line[firstChar])) return true;
 183
 184   // Send on particle data to the ParticleData database.
 185   if (isdigit(line[firstChar]))
 186     return particleData.readString(line, warn);
 187
 188   // Everything else sent on to Settings.
 189   return settings.readString(line, warn);
 190
 191 }
 192
 193 //--------------------------------------------------------------------------
 194
 195 // Read in updates for settings or particle data from user-defined file.
 196
 197 bool Pythia::readFile(string fileName, bool warn, int subrun) {
 198
 199   // Check that constructor worked.
 200   if (!isConstructed) return false;
 201
 202   // Open file for reading.
 203   const char* cstring = fileName.c_str();
 204   ifstream is(cstring);
 205   if (!is.good()) {
 206     info.errorMsg("Error in Pythia::readFile: did not find file", fileName);
 207     return false;
 208   }
 209
 210   // Hand over real work to next method.
 211   return readFile( is, warn, subrun);
 212
 213 }
 214
 215 //--------------------------------------------------------------------------
 216
 217 // Read in updates for settings or particle data
 218 // from user-defined stream (or file).
 219
 220 bool Pythia::readFile(istream& is, bool warn, int subrun) {
 221
 222   // Check that constructor worked.
 223   if (!isConstructed) return false;
 224
 225   // Read in one line at a time.
 226   string line;
 227   bool accepted = true;
 228   int subrunNow = SUBRUNDEFAULT;
 229   while ( getline(is, line) ) {
 230
 231     // Check whether entered new subrun.
 232     int subrunLine = readSubrun( line, warn);
 233     if (subrunLine >= 0) subrunNow = subrunLine;
 234
 235     // Process the line if in correct subrun.
 236     if ( (subrunNow == subrun || subrunNow == SUBRUNDEFAULT)
 237        && !readString( line, warn) ) accepted = false;
 238
 239   // Reached end of input file.
 240   };
 241   return accepted;
 242
 243 }
 244
 245 //--------------------------------------------------------------------------
 246
 247 // Routine to pass in pointers to PDF's. Usage optional.
 248
 249 bool Pythia::setPDFPtr( PDF* pdfAPtrIn, PDF* pdfBPtrIn, PDF* pdfHardAPtrIn,
 250   PDF* pdfHardBPtrIn, PDF* pdfPomAPtrIn, PDF* pdfPomBPtrIn) {
 251
 252   // Delete any PDF's created in a previous init call.
 253   if (useNewPdfHard && pdfHardAPtr != pdfAPtr) delete pdfHardAPtr;
 254   if (useNewPdfHard && pdfHardBPtr != pdfBPtr) delete pdfHardBPtr;
 255   if (useNewPdfA) delete pdfAPtr;
 256   if (useNewPdfB) delete pdfBPtr;
 257   if (useNewPdfPomA) delete pdfPomAPtr;
 258   if (useNewPdfPomB) delete pdfPomBPtr;
 259
 260   // Reset pointers to be empty.
 261   useNewPdfA    = false;
 262   useNewPdfB    = false;
 263   useNewPdfHard = false;
 264   useNewPdfPomA = false;
 265   useNewPdfPomB = false;
 266   pdfAPtr       = 0;
 267   pdfBPtr       = 0;
 268   pdfHardAPtr   = 0;
 269   pdfHardBPtr   = 0;
 270   pdfPomAPtr    = 0;
 271   pdfPomBPtr    = 0;
 272
 273   // Switch off external PDF's by zero as input.
 274   if (pdfAPtrIn == 0 && pdfBPtrIn == 0) return true;
 275
 276   // The two PDF objects cannot be one and the same.
 277   if (pdfAPtrIn == pdfBPtrIn) return false;
 278
 279   // Save pointers.
 280   pdfAPtr       = pdfAPtrIn;
 281   pdfBPtr       = pdfBPtrIn;
 282
 283   // By default same pointers for hard-process PDF's.
 284   pdfHardAPtr   = pdfAPtrIn;
 285   pdfHardBPtr   = pdfBPtrIn;
 286
 287   // Optionally allow separate pointers for hard process.
 288   if (pdfHardAPtrIn != 0 && pdfHardBPtrIn != 0) {
 289     if (pdfHardAPtrIn == pdfHardBPtrIn) return false;
 290     pdfHardAPtr = pdfHardAPtrIn;
 291     pdfHardBPtr = pdfHardBPtrIn;
 292   }
 293
 294   // Optionally allow pointers for Pomerons in the proton.
 295   if (pdfPomAPtrIn != 0 && pdfPomBPtrIn != 0) {
 296     if (pdfPomAPtrIn == pdfPomBPtrIn) return false;
 297     pdfPomAPtr  = pdfPomAPtrIn;
 298     pdfPomBPtr  = pdfPomBPtrIn;
 299   }
 300
 301   // Done.
 302   return true;
 303 }
 304
 305 //--------------------------------------------------------------------------
 306
 307 // Routine to initialize with the variable values of the Beams kind.
 308
 309 bool Pythia::init() {
 310
 311   // Check that constructor worked.
 312   isInit = false;
 313   if (!isConstructed) {
 314     info.errorMsg("Abort from Pythia::init: constructur "
 315       "initialization failed");
 316     return false;
 317   }
 318
 319   // Begin initialization. Find which frame type to use.
 320   info.addCounter(1);
 321   frameType = mode("Beams:frameType");
 322
 323   // Initialization with internal processes: read in and set values.
 324   if (frameType < 4 ) {
 325     doLHA     = false;
 326     boostType = frameType;
 327     idA       = mode("Beams:idA");
 328     idB       = mode("Beams:idB");
 329     eCM       = parm("Beams:eCM");
 330     eA        = parm("Beams:eA");
 331     eB        = parm("Beams:eB");
 332     pxA       = parm("Beams:pxA");
 333     pyA       = parm("Beams:pyA");
 334     pzA       = parm("Beams:pzA");
 335     pxB       = parm("Beams:pxB");
 336     pyB       = parm("Beams:pyB");
 337     pzB       = parm("Beams:pzB");
 338
 339    // Initialization with a Les Houches Event File or an LHAup object.
 340   } else {
 341     doLHA     = true;
 342     boostType = 2;
 343     string lhef        = word("Beams:LHEF");
 344     string lhefHeader  = word("Beams:LHEFheader");
 345     bool   readHeaders = flag("Beams:readLHEFheaders");
 346     bool   skipInit    = flag("Beams:newLHEFsameInit");
 347     int    nSkipAtInit = mode("Beams:nSkipLHEFatInit");
 348
 349     // For file input: renew file stream or (re)new Les Houches object.
 350     if (frameType == 4) {
 351       const char* cstring1 = lhef.c_str();
 352       if (useNewLHA && skipInit) lhaUpPtr->newEventFile(cstring1);
 353       else {
 354         if (useNewLHA) delete lhaUpPtr;
 355         // Header is optional, so use NULL pointer to indicate no value.
 356         const char* cstring2 = (lhefHeader == "void") ?
 357                                NULL : lhefHeader.c_str();
 358         lhaUpPtr   = new LHAupLHEF(cstring1, cstring2, readHeaders);
 359         useNewLHA  = true;
 360       }
 361
 362       // Check that file was properly opened.
 363       if (!lhaUpPtr->fileFound()) {
 364         info.errorMsg("Abort from Pythia::init: "
 365           "Les Houches Event File not found");
 366         return false;
 367       }
 368
 369     // For object input: at least check that not null pointer.
 370     } else {
 371       if (lhaUpPtr == 0) {
 372         info.errorMsg("Abort from Pythia::init: "
 373           "LHAup object not found");
 374         return false;
 375       }
 376
 377       // LHAup object generic abort using fileFound() routine
 378       if (!lhaUpPtr->fileFound()) {
 379         info.errorMsg("Abort from Pythia::init: "
 380           "LHAup initialisation error");
 381         return false;
 382       }
 383     }
 384
 385     // Send in pointer to info. Store or replace LHA pointer in other classes.
 386     lhaUpPtr->setPtr( &info);
 387     processLevel.setLHAPtr( lhaUpPtr);
 388
 389     // If second time around, only with new file, then simplify.
 390     // Optionally skip ahead a number of events at beginning of file.
 391     if (skipInit) {
 392       isInit = true;
 393       info.addCounter(2);
 394       if (nSkipAtInit > 0) lhaUpPtr->skipEvent(nSkipAtInit);
 395       return true;
 396     }
 397
 398     // Set up values related to merging hooks.
 399     if (frameType == 4) {
 400       // Set up values related to CKKW-L merging.
 401       doUserMerging     = settings.flag("Merging:doUserMerging");
 402       doMGMerging       = settings.flag("Merging:doMGMerging");
 403       doKTMerging       = settings.flag("Merging:doKTMerging");
 404       doPTLundMerging   = settings.flag("Merging:doPTLundMerging");
 405       doCutBasedMerging = settings.flag("Merging:doCutBasedMerging");
 406       doMerging         = doUserMerging || doMGMerging || doKTMerging
 407                        || doPTLundMerging || doCutBasedMerging;
 408       // Set up MergingHooks object
 409       if (!doUserMerging){
 410         if (hasOwnMergingHooks && mergingHooksPtr) delete mergingHooksPtr;
 411         mergingHooksPtr = new MergingHooks();
 412         hasOwnMergingHooks = true;
 413       }
 414       hasMergingHooks  = (mergingHooksPtr != 0);
 415       // Merging hooks required for merging. If no merging hooks pointer is
 416       // available, exit.
 417       if (!hasMergingHooks) {
 418         info.errorMsg("Abort from Pythia::init: "
 419           "no merging hooks object has been provided");
 420         return false;
 421       } else mergingHooksPtr->setLHEInputFile( lhef);
 422       // Initialise counting of Les Houches Events significantly above the
 423       // merging scale.
 424       info.setCounter(41,0);
 425     }
 426
 427     // Set LHAinit information (in some external program).
 428     if (settings.flag("ProcessLevel:all") && !lhaUpPtr->setInit()) {
 429       info.errorMsg("Abort from Pythia::init: "
 430         "Les Houches initialization failed");
 431       return false;
 432     }
 433
 434     // Extract beams from values set in an LHAinit object.
 435     idA = lhaUpPtr->idBeamA();
 436     idB = lhaUpPtr->idBeamB();
 437     eA  = lhaUpPtr->eBeamA();
 438     eB  = lhaUpPtr->eBeamB();
 439     // Optionally skip ahead a number of events at beginning of file.
 440     if (nSkipAtInit > 0) lhaUpPtr->skipEvent(nSkipAtInit);
 441   }
 442
 443   // Set up values related to user hooks.
 444   hasUserHooks  = (userHooksPtr != 0);
 445   doVetoProcess = false;
 446   doVetoPartons = false;
 447   if (hasUserHooks) {
 448     userHooksPtr->initPtr( &info, &settings, &particleData, &rndm,
 449         &beamA, &beamB, &beamPomA, &beamPomB, couplingsPtr, &partonSystems,
 450         &sigmaTot);
 451     if (!userHooksPtr->initAfterBeams()) {
 452       info.errorMsg("Abort from Pythia::init: could not initialise UserHooks");
 453       return false;
 454     }
 455     doVetoProcess = userHooksPtr->canVetoProcessLevel();
 456     doVetoPartons = userHooksPtr->canVetoPartonLevel();
 457   }
 458
 459   // Back to common initialization. Reset error counters.
 460   nErrEvent = 0;
 461   info.errorReset();
 462   info.setTooLowPTmin(false);
 463   info.sigmaReset();
 464
 465   // Initialize data members extracted from database.
 466   doProcessLevel   = settings.flag("ProcessLevel:all");
 467   doPartonLevel    = settings.flag("PartonLevel:all");
 468   doHadronLevel    = settings.flag("HadronLevel:all");
 469   doDiffraction    = settings.flag("SoftQCD:all")
 470                   || settings.flag("SoftQCD:singleDiffractive")
 471                   || settings.flag("SoftQCD:doubleDiffractive")
 472                   || settings.flag("SoftQCD:centralDiffractive");
 473   doResDec         = settings.flag("Standalone:allowResDec");
 474   doFSRinRes       = doPartonLevel && settings.flag("PartonLevel:FSR")
 475                   && settings.flag("PartonLevel:FSRinResonances");
 476   decayRHadrons    = settings.flag("RHadrons:allowDecay");
 477   doMomentumSpread = settings.flag("Beams:allowMomentumSpread");
 478   doVertexSpread   = settings.flag("Beams:allowVertexSpread");
 479   abortIfVeto      = settings.flag("Check:abortIfVeto");
 480   checkEvent       = settings.flag("Check:event");
 481   checkHistory     = settings.flag("Check:history");
 482   nErrList         = settings.mode("Check:nErrList");
 483   epTolErr         = settings.parm("Check:epTolErr");
 484   epTolWarn        = settings.parm("Check:epTolWarn");
 485
 486
 487   // Initialise merging hooks.
 488   if (hasMergingHooks && doMerging )
 489     mergingHooksPtr->init( settings, &info, &particleData );
 490
 491   // Initialize the random number generator.
 492   if ( settings.flag("Random:setSeed") )
 493     rndm.init( settings.mode("Random:seed") );
 494
 495   // Check that combinations of settings are allowed; change if not.
 496   checkSettings();
 497
 498   // Initialize couplings (needed to initialize resonances).
 499   // Check if SUSY couplings need to be read in
 500   if( !initSLHA()) info.errorMsg("Error in Pythia::init: "
 501     "Could not read SLHA file");
 502   if (couplingsPtr->isSUSY) {
 503     // Initialize the SM and SUSY.
 504     coupSUSY.init( settings, &rndm);
 505     coupSUSY.initSUSY(&slha, &settings, &particleData);
 506     couplingsPtr = (Couplings *) &coupSUSY;
 507   } else {
 508     // Initialize the SM couplings.
 509     couplingsPtr->init( settings, &rndm);
 510   }
 511
 512   // Reset couplingsPtr to the correct place.
 513   particleData.initPtr( &info, &settings, &rndm, couplingsPtr);
 514
 515   // Set headers to distinguish the two event listing kinds.
 516   int startColTag = settings.mode("Event:startColTag");
 517   process.init("(hard process)", &particleData, startColTag);
 518   event.init("(complete event)", &particleData, startColTag);
 519
 520   // Final setup stage of particle data, notably resonance widths.
 521   particleData.initWidths( resonancePtrs);
 522
 523   // Set up R-hadrons particle data, where relevant.
 524   rHadrons.init( &info, settings, &particleData, &rndm);
 525
 526   // Set up objects for timelike and spacelike showers.
 527   if (timesDecPtr == 0 || timesPtr == 0) {
 528     TimeShower* timesNow = new TimeShower();
 529     if (timesDecPtr == 0) timesDecPtr = timesNow;
 530     if (timesPtr == 0) timesPtr = timesNow;
 531     useNewTimes = true;
 532   }
 533   if (spacePtr == 0) {
 534     spacePtr    = new SpaceShower();
 535     useNewSpace = true;
 536   }
 537
 538   // Initialize showers, especially for simple showers in decays.
 539   timesPtr->initPtr( &info, &settings, &particleData, &rndm, couplingsPtr,
 540     &partonSystems, userHooksPtr);
 541   timesDecPtr->initPtr( &info, &settings, &particleData, &rndm, couplingsPtr,
 542     &partonSystems, userHooksPtr);
 543   spacePtr->initPtr( &info, &settings, &particleData, &rndm,
 544     &partonSystems, userHooksPtr);
 545   timesDecPtr->init( 0, 0);
 546
 547   // Set up values related to beam shape.
 548   if (beamShapePtr == 0) {
 549     beamShapePtr    = new BeamShape();
 550     useNewBeamShape = true;
 551   }
 552   beamShapePtr->init( settings, &rndm);
 553
 554   // Check that beams and beam combination can be handled.
 555   if (!checkBeams()) {
 556     info.errorMsg("Abort from Pythia::init: "
 557       "checkBeams initialization failed");
 558     return false;
 559   }
 560
 561   // Do not set up beam kinematics when no process level.
 562   if (!doProcessLevel) boostType = 1;
 563   else {
 564
 565     // Set up beam kinematics.
 566     if (!initKinematics()) {
 567       info.errorMsg("Abort from Pythia::init: "
 568         "kinematics initialization failed");
 569       return false;
 570     }
 571
 572     // Set up pointers to PDFs.
 573     if (!initPDFs()) {
 574       info.errorMsg("Abort from Pythia::init: PDF initialization failed");
 575       return false;
 576     }
 577
 578     // Set up the two beams and the common remnant system.
 579     StringFlav* flavSelPtr = hadronLevel.getStringFlavPtr();
 580     beamA.init( idA, pzAcm, eA, mA, &info, settings, &particleData, &rndm,
 581       pdfAPtr, pdfHardAPtr, isUnresolvedA, flavSelPtr);
 582     beamB.init( idB, pzBcm, eB, mB, &info, settings, &particleData, &rndm,
 583       pdfBPtr, pdfHardBPtr, isUnresolvedB, flavSelPtr);
 584
 585     // Optionally set up new alternative beams for these Pomerons.
 586     if ( doDiffraction) {
 587       beamPomA.init( 990,  0.5 * eCM, 0.5 * eCM, 0., &info, settings,
 588         &particleData, &rndm, pdfPomAPtr, pdfPomAPtr, false, flavSelPtr);
 589       beamPomB.init( 990, -0.5 * eCM, 0.5 * eCM, 0., &info, settings,
 590         &particleData, &rndm, pdfPomBPtr, pdfPomBPtr, false, flavSelPtr);
 591     }
 592   }
 593
 594   // Send info/pointers to process level for initialization.
 595   if ( doProcessLevel && !processLevel.init( &info, settings, &particleData,
 596     &rndm, &beamA, &beamB, couplingsPtr, &sigmaTot, doLHA, &slha, userHooksPtr,
 597     sigmaPtrs) ) {
 598     info.errorMsg("Abort from Pythia::init: "
 599       "processLevel initialization failed");
 600     return false;
 601   }
 602
 603   // Optionally only initialize resonance decays.
 604   if ( !doProcessLevel && doResDec) processLevel.initDecays( &info,
 605     &particleData, &rndm, lhaUpPtr);
 606
 607   // Send info/pointers to parton level for initialization.
 608   if ( doPartonLevel && doProcessLevel && !partonLevel.init( &info, settings,
 609     &particleData, &rndm, &beamA, &beamB, &beamPomA, &beamPomB, couplingsPtr,
 610     &partonSystems, &sigmaTot, timesDecPtr, timesPtr, spacePtr, &rHadrons,
 611     userHooksPtr, mergingHooksPtr, false) ) {
 612     info.errorMsg("Abort from Pythia::init: "
 613       "partonLevel initialization failed" );
 614     return false;
 615   }
 616
 617   // Optionally only initialize final-state showers in resonance decays.
 618   if ( (!doProcessLevel || !doPartonLevel) && doFSRinRes) partonLevel.init(
 619     &info, settings, &particleData, &rndm, 0, 0, 0, 0, couplingsPtr,
 620     &partonSystems, 0, timesDecPtr, 0, 0, &rHadrons, 0, 0, false);
 621
 622   // Send info/pointers to parton level for trial shower initialization.
 623   if ( doMerging
 624     && !trialPartonLevel.init( &info, settings, &particleData,
 625       &rndm, &beamA, &beamB, &beamPomA, &beamPomB, couplingsPtr,
 626       &partonSystems, &sigmaTot, timesDecPtr, timesPtr, spacePtr, &rHadrons,
 627       NULL, mergingHooksPtr, true) ) {
 628     info.errorMsg("Abort from Pythia::init: "
 629       "trialPartonLevel initialization failed");
 630     return false;
 631   }
 632
 633   // Send info/pointers to hadron level for initialization.
 634   // Note: forceHadronLevel() can come, so we must always initialize.
 635   if ( !hadronLevel.init( &info, settings, &particleData, &rndm,
 636     couplingsPtr, timesDecPtr, &rHadrons, decayHandlePtr,
 637     handledParticles) ) {
 638     info.errorMsg("Abort from Pythia::init: "
 639       "hadronLevel initialization failed");
 640     return false;
 641   }
 642
 643   // Optionally check particle data table for inconsistencies.
 644   if ( settings.flag("Check:particleData") )
 645     particleData.checkTable( settings.mode("Check:levelParticleData") );
 646
 647   // Optionally show settings and particle data, changed or all.
 648   bool showCS  = settings.flag("Init:showChangedSettings");
 649   bool showAS  = settings.flag("Init:showAllSettings");
 650   bool showCPD = settings.flag("Init:showChangedParticleData");
 651   bool showCRD = settings.flag("Init:showChangedResonanceData");
 652   bool showAPD = settings.flag("Init:showAllParticleData");
 653   int  show1PD = settings.mode("Init:showOneParticleData");
 654   bool showPro = settings.flag("Init:showProcesses");
 655   if (showCS)      settings.listChanged();
 656   if (showAS)      settings.listAll();
 657   if (show1PD > 0) particleData.list(show1PD);
 658   if (showCPD)     particleData.listChanged(showCRD);
 659   if (showAPD)     particleData.listAll();
 660
 661   // Listing options for the next() routine.
 662   nCount       = settings.mode("Next:numberCount");
 663   nShowLHA     = settings.mode("Next:numberShowLHA");
 664   nShowInfo    = settings.mode("Next:numberShowInfo");
 665   nShowProc    = settings.mode("Next:numberShowProcess");
 666   nShowEvt     = settings.mode("Next:numberShowEvent");
 667   showSaV      = settings.flag("Next:showScaleAndVertex");
 668   showMaD      = settings.flag("Next:showMothersAndDaughters");
 669
 670   // Succeeded.
 671   isInit = true;
 672   info.addCounter(2);
 673   if (useNewLHA && showPro) lhaUpPtr->listInit();
 674   return true;
 675
 676 }
 677
 678 //--------------------------------------------------------------------------
 679
 680 // Routine to initialize with CM energy.
 681
 682 bool Pythia::init( int idAin, int idBin, double eCMin) {
 683
 684   // Set arguments in Settings database.
 685   settings.mode("Beams:idA",  idAin);
 686   settings.mode("Beams:idB",  idBin);
 687   settings.mode("Beams:frameType",  1);
 688   settings.parm("Beams:eCM", eCMin);
 689
 690   // Send on to common initialization.
 691   return init();
 692
 693 }
 694
 695 //--------------------------------------------------------------------------
 696
 697 // Routine to initialize with two collinear beams,  energies specified.
 698
 699 bool Pythia::init( int idAin, int idBin, double eAin, double eBin) {
 700
 701   // Set arguments in Settings database.
 702   settings.mode("Beams:idA",  idAin);
 703   settings.mode("Beams:idB",  idBin);
 704   settings.mode("Beams:frameType",  2);
 705   settings.parm("Beams:eA",      eAin);
 706   settings.parm("Beams:eB",      eBin);
 707
 708   // Send on to common initialization.
 709   return init();
 710
 711 }
 712
 713 //--------------------------------------------------------------------------
 714
 715 // Routine to initialize with two beams specified by three-momenta.
 716
 717 bool Pythia::init( int idAin, int idBin, double pxAin, double pyAin,
 718   double pzAin, double pxBin, double pyBin, double pzBin) {
 719
 720   // Set arguments in Settings database.
 721   settings.mode("Beams:idA",  idAin);
 722   settings.mode("Beams:idB",  idBin);
 723   settings.mode("Beams:frameType",  3);
 724   settings.parm("Beams:pxA",    pxAin);
 725   settings.parm("Beams:pyA",    pyAin);
 726   settings.parm("Beams:pzA",    pzAin);
 727   settings.parm("Beams:pxB",    pxBin);
 728   settings.parm("Beams:pyB",    pyBin);
 729   settings.parm("Beams:pzB",    pzBin);
 730
 731   // Send on to common initialization.
 732   return init();
 733
 734 }
 735
 736 //--------------------------------------------------------------------------
 737
 738 // Routine to initialize when all info is given in a Les Houches Event File.
 739
 740 bool Pythia::init( string LesHouchesEventFile, bool skipInit) {
 741
 742   // Set arguments in Settings database.
 743   settings.mode("Beams:frameType",              4);
 744   settings.word("Beams:LHEF", LesHouchesEventFile);
 745   settings.flag("Beams:newLHEFsameInit", skipInit);
 746
 747   // Send on to common initialization.
 748   return init();
 749
 750 }
 751
 752 //--------------------------------------------------------------------------
 753
 754 // Routine to initialize when beam info is given in an LHAup object.
 755
 756 bool Pythia::init( LHAup* lhaUpPtrIn) {
 757
 758   // Store LHAup object and set arguments in Settings database.
 759   setLHAupPtr( lhaUpPtrIn);
 760   settings.mode("Beams:frameType", 5);
 761
 762   // Send on to common initialization.
 763   return init();
 764
 765 }
 766
 767 //--------------------------------------------------------------------------
 768
 769 // Check that combinations of settings are allowed; change if not.
 770
 771 void Pythia::checkSettings() {
 772
 773   // Double rescattering not allowed if ISR or FSR.
 774   if ((settings.flag("PartonLevel:ISR") || settings.flag("PartonLevel:FSR"))
 775     && settings.flag("MultipartonInteractions:allowDoubleRescatter")) {
 776     info.errorMsg("Warning in Pythia::checkSettings: "
 777         "double rescattering switched off since showering is on");
 778     settings.flag("MultipartonInteractions:allowDoubleRescatter", false);
 779   }
 780
 781 }
 782
 783 //--------------------------------------------------------------------------
 784
 785 // Check that beams and beam combination can be handled. Set up unresolved.
 786
 787 bool Pythia::checkBeams() {
 788
 789   // Absolute flavours. If not to do process level then no check needed.
 790   int idAabs = abs(idA);
 791   int idBabs = abs(idB);
 792   if (!doProcessLevel) return true;
 793
 794   // Neutrino beams always unresolved, charged lepton ones conditionally.
 795   bool isLeptonA  = (idAabs > 10 && idAabs < 17);
 796   bool isLeptonB  = (idBabs > 10 && idBabs < 17);
 797   bool isUnresLep = !settings.flag("PDF:lepton");
 798   isUnresolvedA   = isLeptonA && (idAabs%2 == 0 || isUnresLep);
 799   isUnresolvedB   = isLeptonB && (idBabs%2 == 0 || isUnresLep);
 800
 801   // Lepton-lepton collisions OK (including neutrinos) if both (un)resolved.
 802   if (isLeptonA && isLeptonB && isUnresolvedA == isUnresolvedB) return true;
 803
 804   // MBR model only implemented for pp/ppbar/pbarp collisions.
 805   int PomFlux     = settings.mode("Diffraction:PomFlux");
 806   if (PomFlux == 5) {
 807     bool ispp       = (idAabs == 2212 && idBabs == 2212);
 808     bool ispbarpbar = (idA == -2212 && idB == -2212);
 809     if (ispp && !ispbarpbar) return true;
 810     info.errorMsg("Error in Pythia::init: cannot handle this beam combination"
 811       " with PomFlux == 5");
 812     return false;
 813   }
 814
 815   // Hadron-hadron collisions OK, with Pomeron counted as hadron.
 816   bool isHadronA = (idAabs == 2212) || (idA == 111) || (idAabs == 211)
 817                 || (idA == 990);
 818   bool isHadronB = (idBabs == 2212) || (idA == 111)|| (idBabs == 211)
 819                 || (idB == 990);
 820   if (isHadronA && isHadronB) return true;
 821
 822   // If no case above then failed.
 823   info.errorMsg("Error in Pythia::init: cannot handle this beam combination");
 824   return false;
 825
 826 }
 827
 828 //--------------------------------------------------------------------------
 829
 830 // Calculate kinematics at initialization. Store beam four-momenta.
 831
 832 bool Pythia::initKinematics() {
 833
 834   // Find masses. Initial guess that we are in CM frame.
 835   mA       = particleData.m0(idA);
 836   mB       = particleData.m0(idB);
 837   betaZ    = 0.;
 838   gammaZ   = 1.;
 839
 840   // Collinear beams not in CM frame: find CM energy.
 841   if (boostType == 2) {
 842     eA     = max(eA, mA);
 843     eB     = max(eB, mB);
 844     pzA    = sqrt(eA*eA - mA*mA);
 845     pzB    = -sqrt(eB*eB - mB*mB);
 846     pAinit = Vec4( 0., 0., pzA, eA);
 847     pBinit = Vec4( 0., 0., pzB, eB);
 848     eCM    = sqrt( pow2(eA + eB) - pow2(pzA + pzB) );
 849
 850     // Find boost to rest frame.
 851     betaZ  = (pzA + pzB) / (eA + eB);
 852     gammaZ = (eA + eB) / eCM;
 853     if (abs(betaZ) < 1e-10) boostType = 1;
 854   }
 855
 856   // Completely general beam directions: find CM energy.
 857   else if (boostType == 3) {
 858     eA     = sqrt( pxA*pxA + pyA*pyA + pzA*pzA + mA*mA);
 859     eB     = sqrt( pxB*pxB + pyB*pyB + pzB*pzB + mB*mB);
 860     pAinit = Vec4( pxA, pyA, pzA, eA);
 861     pBinit = Vec4( pxB, pyB, pzB, eB);
 862     eCM = (pAinit + pBinit).mCalc();
 863
 864     // Find boost+rotation needed to move from/to CM frame.
 865     MfromCM.reset();
 866     MfromCM.fromCMframe( pAinit, pBinit);
 867     MtoCM = MfromCM;
 868     MtoCM.invert();
 869   }
 870
 871   // Fail if CM energy below beam masses.
 872   if (eCM < mA + mB) {
 873     info.errorMsg("Error in Pythia::initKinematics: too low energy");
 874     return false;
 875   }
 876
 877   // Set up CM-frame kinematics with beams along +-z axis.
 878   pzAcm    = 0.5 * sqrtpos( (eCM + mA + mB) * (eCM - mA - mB)
 879            * (eCM - mA + mB) * (eCM + mA - mB) ) / eCM;
 880   pzBcm    = -pzAcm;
 881   eA       = sqrt(mA*mA + pzAcm*pzAcm);
 882   eB       = sqrt(mB*mB + pzBcm*pzBcm);
 883
 884   // If in CM frame then store beam four-vectors (else already done above).
 885   if (boostType != 2 && boostType != 3) {
 886     pAinit = Vec4( 0., 0., pzAcm, eA);
 887     pBinit = Vec4( 0., 0., pzBcm, eB);
 888   }
 889
 890   // Store main info for access in process generation.
 891   info.setBeamA( idA, pzAcm, eA, mA);
 892   info.setBeamB( idB, pzBcm, eB, mB);
 893   info.setECM( eCM);
 894
 895   // Must allow for generic boost+rotation when beam momentum spread.
 896   if (doMomentumSpread) boostType = 3;
 897
 898   // Done.
 899   return true;
 900
 901 }
 902
 903 //--------------------------------------------------------------------------
 904
 905 // Set up pointers to PDFs.
 906
 907 bool Pythia::initPDFs() {
 908
 909   // Delete any PDF's created in a previous init call.
 910   if (useNewPdfHard) {
 911     if (pdfHardAPtr != pdfAPtr) {
 912       delete pdfHardAPtr;
 913       pdfHardAPtr = 0;
 914     }
 915     if (pdfHardBPtr != pdfBPtr) {
 916       delete pdfHardBPtr;
 917       pdfHardBPtr = 0;
 918     }
 919     useNewPdfHard = false;
 920   }
 921   if (useNewPdfA) {
 922     delete pdfAPtr;
 923     useNewPdfA    = false;
 924     pdfAPtr       = 0;
 925   }
 926   if (useNewPdfB) {
 927     delete pdfBPtr;
 928     useNewPdfB    = false;
 929     pdfBPtr       = 0;
 930   }
 931   if (useNewPdfPomA) {
 932     delete pdfPomAPtr;
 933     useNewPdfPomA = false;
 934     pdfPomAPtr    = 0;
 935   }
 936   if (useNewPdfPomB) {
 937     delete pdfPomBPtr;
 938     useNewPdfPomB = false;
 939     pdfPomBPtr    = 0;
 940   }
 941
 942   // Set up the PDF's, if not already done.
 943   if (pdfAPtr == 0) {
 944     pdfAPtr     = getPDFPtr(idA);
 945     if (pdfAPtr == 0 || !pdfAPtr->isSetup()) {
 946       info.errorMsg("Error in Pythia::init: "
 947         "could not set up PDF for beam A");
 948       return false;
 949     }
 950     pdfHardAPtr = pdfAPtr;
 951     useNewPdfA  = true;
 952   }
 953   if (pdfBPtr == 0) {
 954     pdfBPtr     = getPDFPtr(idB);
 955     if (pdfBPtr == 0 || !pdfBPtr->isSetup()) {
 956       info.errorMsg("Error in Pythia::init: "
 957         "could not set up PDF for beam B");
 958       return false;
 959     }
 960     pdfHardBPtr = pdfBPtr;
 961     useNewPdfB  = true;
 962   }
 963
 964   // Optionally set up separate PDF's for hard process.
 965   if (settings.flag("PDF:useHard") && useNewPdfA && useNewPdfB) {
 966     pdfHardAPtr = getPDFPtr(idA, 2);
 967     if (!pdfHardAPtr->isSetup()) return false;
 968     pdfHardBPtr = getPDFPtr(idB, 2);
 969     if (!pdfHardBPtr->isSetup()) return false;
 970     useNewPdfHard = true;
 971   }
 972
 973   // Optionally set up Pomeron PDF's for diffractive physics.
 974   if ( doDiffraction) {
 975     if (pdfPomAPtr == 0) {
 976       pdfPomAPtr    = getPDFPtr(990);
 977       useNewPdfPomA = true;
 978     }
 979     if (pdfPomBPtr == 0) {
 980       pdfPomBPtr    = getPDFPtr(990);
 981       useNewPdfPomB = true;
 982     }
 983   }
 984
 985   // Done.
 986   return true;
 987
 988 }
 989
 990 //--------------------------------------------------------------------------
 991
 992 // Main routine to generate the next event, using internal machinery.
 993
 994 bool Pythia::next() {
 995
 996   // Regularly print how many events have been generated.
 997   int nPrevious = info.getCounter(3);
 998   if (nCount > 0 && nPrevious > 0 && nPrevious%nCount == 0)
 999     cout << "\n Pythia::next(): " << nPrevious
1000          << " events have been generated " << endl;
1001
1002   // Set/reset info counters specific to each event.
1003   info.addCounter(3);
1004   for (int i = 10; i < 13; ++i) info.setCounter(i);
1005
1006   // Simpler option when no hard process, i.e. mainly hadron level.
1007   if (!doProcessLevel) {
1008
1009     // Optionally fetch in resonance decays from LHA interface.
1010     if (doLHA && !processLevel.nextLHAdec( event)) {
1011       if (info.atEndOfFile()) info.errorMsg("Abort from Pythia::next:"
1012         " reached end of Les Houches Events File");
1013       return false;
1014     }
1015
1016     // Reset info array (while event record contains data).
1017     info.clear();
1018
1019     // Set correct energy for system.
1020     Vec4 pSum = 0.;
1021     for (int i = 1; i < event.size(); ++i)
1022       if (event[i].isFinal()) pSum += event[i].p();
1023     event[0].p( pSum );
1024     event[0].m( pSum.mCalc() );
1025
1026     // Generate hadronization and decays.
1027     bool status = forceHadronLevel();
1028     if (status) info.addCounter(4);
1029     if (status && nPrevious < nShowEvt)  event.list(showSaV, showMaD);
1030     return status;
1031   }
1032
1033   // Reset arrays.
1034   info.clear();
1035   process.clear();
1036   event.clear();
1037   partonSystems.clear();
1038   beamA.clear();
1039   beamB.clear();
1040   beamPomA.clear();
1041   beamPomB.clear();
1042
1043   // Can only generate event if initialization worked.
1044   if (!isInit) {
1045     info.errorMsg("Abort from Pythia::next: "
1046       "not properly initialized so cannot generate events");
1047     return false;
1048   }
1049
1050   // Pick beam momentum spread and beam vertex.
1051   if (doMomentumSpread || doVertexSpread) beamShapePtr->pick();
1052
1053   // Recalculate kinematics when beam momentum spread.
1054   if (doMomentumSpread) nextKinematics();
1055
1056   // Outer loop over hard processes; only relevant for user-set vetoes.
1057   for ( ; ; ) {
1058     info.addCounter(10);
1059     bool hasVetoed = false;
1060
1061     // Provide the hard process that starts it off. Only one try.
1062     info.clear();
1063     process.clear();
1064
1065     if ( !processLevel.next( process) ) {
1066       if (doLHA && info.atEndOfFile()) info.errorMsg("Abort from "
1067         "Pythia::next: reached end of Les Houches Events File");
1068       else info.errorMsg("Abort from Pythia::next: "
1069         "processLevel failed; giving up");
1070       return false;
1071     }
1072
1073     info.addCounter(11);
1074
1075     // Possibility for a user veto of the process-level event.
1076     if (doVetoProcess) {
1077       hasVetoed = userHooksPtr->doVetoProcessLevel( process);
1078       if (hasVetoed) {
1079         if (abortIfVeto) return false;
1080         continue;
1081       }
1082     }
1083
1084     // Possibility to perform CKKW-L merging on this event.
1085     if(doMerging) {
1086       hasVetoed = mergeProcess();
1087       if (hasVetoed) {
1088         if (abortIfVeto) return false;
1089         continue;
1090       }
1091     }
1092
1093     // Possibility to stop the generation at this stage.
1094     if (!doPartonLevel) {
1095       boostAndVertex( true, true);
1096       processLevel.accumulate();
1097       info.addCounter(4);
1098       if (doLHA && nPrevious < nShowLHA) lhaUpPtr->listEvent();
1099       if (nPrevious < nShowInfo) info.list();
1100       if (nPrevious < nShowProc) process.list(showSaV, showMaD);
1101       return true;
1102     }
1103
1104     // Save spare copy of process record in case of problems.
1105     Event processSave = process;
1106     int sizeMPI       = info.sizeMPIarrays();
1107     info.addCounter(12);
1108     for (int i = 14; i < 19; ++i) info.setCounter(i);
1109
1110     // Allow up to ten tries for parton- and hadron-level processing.
1111     bool physical   = true;
1112     for (int iTry = 0; iTry < NTRY; ++ iTry) {
1113       info.addCounter(14);
1114       physical = true;
1115
1116       // Restore original process record if problems.
1117       if (iTry > 0) process = processSave;
1118       if (iTry > 0) info.resizeMPIarrays( sizeMPI);
1119
1120       // Reset event record and (extracted partons from) beam remnants.
1121       event.clear();
1122       beamA.clear();
1123       beamB.clear();
1124       beamPomA.clear();
1125       beamPomB.clear();
1126       partonSystems.clear();
1127
1128       // Parton-level evolution: ISR, FSR, MPI.
1129       if ( !partonLevel.next( process, event) ) {
1130         // Skip to next hard process for failure owing to deliberate veto.
1131         hasVetoed = partonLevel.hasVetoed();
1132         if (hasVetoed) {
1133           if (abortIfVeto) return false;
1134           break;
1135         }
1136         // Else make a new try for other failures.
1137         info.errorMsg("Error in Pythia::next: "
1138           "partonLevel failed; try again");
1139         physical = false;
1140         continue;
1141       }
1142       info.addCounter(15);
1143
1144       // Possibility for a user veto of the parton-level event.
1145       if (doVetoPartons) {
1146         hasVetoed = userHooksPtr->doVetoPartonLevel( event);
1147         if (hasVetoed) {
1148           if (abortIfVeto) return false;
1149           break;
1150         }
1151       }
1152
1153       // Boost to lab frame (before decays, for vertices).
1154       boostAndVertex( true, true);
1155
1156       // Possibility to stop the generation at this stage.
1157       if (!doHadronLevel) {
1158         processLevel.accumulate();
1159         partonLevel.accumulate();
1160
1161         // Optionally check final event for problems.
1162         if (checkEvent && !check()) {
1163           info.errorMsg("Abort from Pythia::next: "
1164             "check of event revealed problems");
1165           return false;
1166         }
1167         info.addCounter(4);
1168         if (doLHA && nPrevious < nShowLHA) lhaUpPtr->listEvent();
1169         if (nPrevious < nShowInfo) info.list();
1170         if (nPrevious < nShowProc) process.list(showSaV, showMaD);
1171         if (nPrevious < nShowEvt)  event.list(showSaV, showMaD);
1172         return true;
1173       }
1174
1175       // Hadron-level: hadronization, decays.
1176       info.addCounter(16);
1177       if ( !hadronLevel.next( event) ) {
1178         info.errorMsg("Error in Pythia::next: "
1179           "hadronLevel failed; try again");
1180         physical = false;
1181         continue;
1182       }
1183
1184       // If R-hadrons have been formed, then (optionally) let them decay.
1185       if (decayRHadrons && rHadrons.exist() && !doRHadronDecays()) {
1186         info.errorMsg("Error in Pythia::next: "
1187           "decayRHadrons failed; try again");
1188         physical = false;
1189         continue;
1190       }
1191       info.addCounter(17);
1192
1193       // Optionally check final event for problems.
1194       if (checkEvent && !check()) {
1195         info.errorMsg("Error in Pythia::next: "
1196           "check of event revealed problems");
1197         physical = false;
1198         continue;
1199       }
1200
1201       // Stop parton- and hadron-level looping if you got this far.
1202       info.addCounter(18);
1203       break;
1204     }
1205
1206     // If event vetoed then to make a new try.
1207     if (hasVetoed)  {
1208       if (abortIfVeto) return false;
1209       continue;
1210     }
1211
1212     // If event failed any other way (after ten tries) then give up.
1213     if (!physical) {
1214       info.errorMsg("Abort from Pythia::next: "
1215         "parton+hadronLevel failed; giving up");
1216       return false;
1217     }
1218
1219     // Process- and parton-level statistics. Event scale.
1220     processLevel.accumulate();
1221     partonLevel.accumulate();
1222     event.scale( process.scale() );
1223
1224     // End of outer loop over hard processes. Done with normal option.
1225     info.addCounter(13);
1226     break;
1227   }
1228
1229   // List events.
1230   if (doLHA && nPrevious < nShowLHA) lhaUpPtr->listEvent();
1231   if (nPrevious < nShowInfo) info.list();
1232   if (nPrevious < nShowProc) process.list(showSaV,showMaD);
1233   if (nPrevious < nShowEvt)  event.list(showSaV, showMaD);
1234
1235   // Done.
1236   info.addCounter(4);
1237   return true;
1238
1239 }
1240
1241 //--------------------------------------------------------------------------
1242
1243 // Generate only the hadronization/decay stage, using internal machinery.
1244 // The "event" instance should already contain a parton-level configuration.
1245
1246 bool Pythia::forceHadronLevel(bool findJunctions) {
1247
1248   // Can only generate event if initialization worked.
1249   if (!isInit) {
1250     info.errorMsg("Abort from Pythia::forceHadronLevel: "
1251       "not properly initialized so cannot generate events");
1252     return false;
1253   }
1254
1255   // Check whether any junctions in system. (Normally done in ProcessLevel.)
1256   if (findJunctions) {
1257     event.clearJunctions();
1258     processLevel.findJunctions( event);
1259   }
1260
1261   // Save spare copy of event in case of failure.
1262   Event spareEvent = event;
1263
1264   // Allow up to ten tries for hadron-level processing.
1265   bool physical = true;
1266   for (int iTry = 0; iTry < NTRY; ++ iTry) {
1267     physical = true;
1268
1269     // Check whether any resonances need to be handled at process level.
1270     if (doResDec) {
1271       process = event;
1272       processLevel.nextDecays( process);
1273       bool hasDecays = false;
1274       for (int i = 1; i < process.size(); ++i)
1275         if (process[i].status() < 0) hasDecays = true;
1276
1277       // Allow for showers if decays happened at process level.
1278       if (hasDecays) {
1279         if (doFSRinRes) {
1280           partonLevel.setupShowerSys( process, event);
1281           partonLevel.resonanceShowers( process, event, false);
1282         } else event = process;
1283       }
1284     }
1285
1286     // Hadron-level: hadronization, decays.
1287     if (hadronLevel.next( event)) break;
1288
1289     // If failure then warn, restore original configuration and try again.
1290     info.errorMsg("Error in Pythia::forceHadronLevel: "
1291       "hadronLevel failed; try again");
1292     physical = false;
1293     event    = spareEvent;
1294   }
1295
1296   // Done for simpler option.
1297   if (!physical)  {
1298     info.errorMsg("Abort from Pythia::forceHadronLevel: "
1299       "hadronLevel failed; giving up");
1300     return false;
1301   }
1302
1303   // Optionally check final event for problems.
1304   if (checkEvent && !check()) {
1305     info.errorMsg("Abort from Pythia::forceHadronLevel: "
1306       "check of event revealed problems");
1307     return false;
1308   }
1309
1310   // Done.
1311   return true;
1312
1313 }
1314
1315 //--------------------------------------------------------------------------
1316
1317 // Recalculate kinematics for each event when beam momentum has a spread.
1318
1319 void Pythia::nextKinematics() {
1320
1321   // Read out momentum shift to give current beam momenta.
1322   pAnow = pAinit + beamShapePtr->deltaPA();
1323   pAnow.e( sqrt(pAnow.pAbs2() + mA * mA) );
1324   pBnow = pBinit + beamShapePtr->deltaPB();
1325   pBnow.e( sqrt(pBnow.pAbs2() + mB * mB) );
1326
1327   // Construct CM frame kinematics.
1328   eCM   = (pAnow + pBnow).mCalc();
1329   pzAcm = 0.5 * sqrtpos( (eCM + mA + mB) * (eCM - mA - mB)
1330         * (eCM - mA + mB) * (eCM + mA - mB) ) / eCM;
1331   pzBcm = -pzAcm;
1332   eA    = sqrt(mA*mA + pzAcm*pzAcm);
1333   eB    = sqrt(mB*mB + pzBcm*pzBcm);
1334
1335   // Set relevant info for other classes to use.
1336   info.setBeamA( idA, pzAcm, eA, mA);
1337   info.setBeamB( idB, pzBcm, eB, mB);
1338   info.setECM( eCM);
1339   beamA.newPzE( pzAcm, eA);
1340   beamB.newPzE( pzBcm, eB);
1341
1342   // Set boost/rotation matrices from/to CM frame.
1343   MfromCM.reset();
1344   MfromCM.fromCMframe( pAnow, pBnow);
1345   MtoCM = MfromCM;
1346   MtoCM.invert();
1347
1348 }
1349
1350 //--------------------------------------------------------------------------
1351
1352 // Boost from CM frame to lab frame, or inverse. Set production vertex.
1353
1354 void Pythia::boostAndVertex( bool toLab, bool setVertex) {
1355
1356   // Boost process from CM frame to lab frame.
1357   if (toLab) {
1358     if      (boostType == 2) process.bst(0., 0., betaZ, gammaZ);
1359     else if (boostType == 3) process.rotbst(MfromCM);
1360
1361     // Boost nonempty event from CM frame to lab frame.
1362     if (event.size() > 0) {
1363       if      (boostType == 2) event.bst(0., 0., betaZ, gammaZ);
1364       else if (boostType == 3) event.rotbst(MfromCM);
1365     }
1366
1367   // Boost process from lab frame to CM frame.
1368   } else {
1369     if      (boostType == 2) process.bst(0., 0., -betaZ, gammaZ);
1370     else if (boostType == 3) process.rotbst(MtoCM);
1371
1372     // Boost nonempty event from lab frame to CM frame.
1373     if (event.size() > 0) {
1374       if      (boostType == 2) event.bst(0., 0., -betaZ, gammaZ);
1375       else if (boostType == 3) event.rotbst(MtoCM);
1376     }
1377   }
1378
1379   // Set production vertex; assumes particles are in lab frame and at origin.
1380   if (setVertex && doVertexSpread) {
1381     Vec4 vertex = beamShapePtr->vertex();
1382     for (int i = 0; i < process.size(); ++i) process[i].vProd( vertex);
1383     for (int i = 0; i < event.size(); ++i) event[i].vProd( vertex);
1384   }
1385
1386 }
1387
1388 //--------------------------------------------------------------------------
1389
1390 // Perform R-hadron decays, either as part of normal evolution or forced.
1391
1392 bool Pythia::doRHadronDecays( ) {
1393
1394   // Check if R-hadrons exist to be processed.
1395   if ( !rHadrons.exist() ) return true;
1396
1397   // Do the R-hadron decay itself.
1398   if ( !rHadrons.decay( event) ) return false;
1399
1400   // Perform showers in resonance decay chains.
1401   if ( !partonLevel.resonanceShowers( process, event, false) ) return false;
1402
1403   // Subsequent hadronization and decays.
1404   if ( !hadronLevel.next( event) ) return false;
1405
1406   // Done.
1407   return true;
1408
1409 }
1410
1411 //--------------------------------------------------------------------------
1412
1413 // Print statistics on event generation.
1414
1415 void Pythia::stat() {
1416
1417   // Read out settings for what to include.
1418   bool showPrL = settings.flag("Stat:showProcessLevel");
1419   bool showPaL = settings.flag("Stat:showPartonLevel");
1420   bool showErr = settings.flag("Stat:showErrors");
1421   bool reset   = settings.flag("Stat:reset");
1422
1423   // Statistics on cross section and number of events.
1424   if (doProcessLevel) {
1425     if (showPrL) processLevel.statistics(false);
1426     if (reset)   processLevel.resetStatistics();
1427   }
1428
1429   // Statistics from other classes, currently multiparton interactions.
1430   if (showPaL) partonLevel.statistics(false);
1431   if (reset)   partonLevel.resetStatistics();
1432
1433   // Warning if significant part of events considerably above the
1434   // merging scale value.
1435   if (hasMergingHooks && doMerging && mergingHooksPtr->stats() ) {
1436     string message="Warning in MergingHooks::stats: Most LH";
1437     message+=" Events significantly above Merging:TMS cut. Please check.";
1438     info.errorMsg(message);
1439   }
1440
1441   // Summary of which and how many warnings/errors encountered.
1442   if (showErr) info.errorStatistics();
1443   if (reset)   info.errorReset();
1444
1445 }
1446
1447 //--------------------------------------------------------------------------
1448
1449 // Print statistics on event generation - deprecated version.
1450
1451 void Pythia::statistics(bool all, bool reset) {
1452
1453   // Statistics on cross section and number of events.
1454   if (doProcessLevel) processLevel.statistics(reset);
1455
1456   // Statistics from other classes, e.g. multiparton interactions.
1457   if (all) partonLevel.statistics(reset);
1458
1459   // Warning if significant part of events considerably above the
1460   // merging scale value.
1461   if (hasMergingHooks && doMerging && mergingHooksPtr->stats() ) {
1462     string message="Warning in MergingHooks::stats: Most LH";
1463     message+=" Events significantly above Merging:TMS cut. Please check.";
1464     info.errorMsg(message);
1465   }
1466
1467   // Summary of which and how many warnings/errors encountered.
1468   info.errorStatistics();
1469   if (reset) info.errorReset();
1470
1471 }
1472
1473 //--------------------------------------------------------------------------
1474
1475 // Write the Pythia banner, with symbol and version information.
1476
1477 void Pythia::banner(ostream& os) {
1478
1479   // Read in version number and last date of change.
1480   double versionNumber = settings.parm("Pythia:versionNumber");
1481   int versionDate = settings.mode("Pythia:versionDate");
1482   string month[12] = {"Jan", "Feb", "Mar", "Apr", "May", "Jun",
1483     "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"};
1484
1485   // Get date and time.
1486   time_t t = time(0);
1487   char dateNow[12];
1488   strftime(dateNow,12,"%d %b %Y",localtime(&t));
1489   char timeNow[9];
1490   strftime(timeNow,9,"%H:%M:%S",localtime(&t));
1491
1492   os << "\n"
1493      << " *-------------------------------------------"
1494      << "-----------------------------------------* \n"
1495      << " |                                           "
1496      << "                                         | \n"
1497      << " |  *----------------------------------------"
1498      << "--------------------------------------*  | \n"
1499      << " |  |                                        "
1500      << "                                      |  | \n"
1501      << " |  |                                        "
1502      << "                                      |  | \n"
1503      << " |  |   PPP   Y   Y  TTTTT  H   H  III    A  "
1504      << "    Welcome to the Lund Monte Carlo!  |  | \n"
1505      << " |  |   P  P   Y Y     T    H   H   I    A A "
1506      << "    This is PYTHIA version " << fixed << setprecision(3)
1507      << setw(5) << versionNumber << "      |  | \n"
1508      << " |  |   PPP     Y      T    HHHHH   I   AAAAA"
1509      << "    Last date of change: " << setw(2) << versionDate%100
1510      << " " << month[ (versionDate/100)%100 - 1 ]
1511      << " " << setw(4) << versionDate/10000 <<  "  |  | \n"
1512      << " |  |   P       Y      T    H   H   I   A   A"
1513      << "                                      |  | \n"
1514      << " |  |   P       Y      T    H   H  III  A   A"
1515      << "    Now is " << dateNow << " at " << timeNow << "    |  | \n"
1516      << " |  |                                        "
1517      << "                                      |  | \n"
1518      << " |  |   Torbjorn Sjostrand;  Department of As"
1519      << "tronomy and Theoretical Physics,      |  | \n"
1520      << " |  |      Lund University, Solvegatan 14A, S"
1521      << "E-223 62 Lund, Sweden;                |  | \n"
1522      << " |  |      phone: + 46 - 46 - 222 48 16; e-ma"
1523      << "il: torbjorn@thep.lu.se               |  | \n"
1524      << " |  |   Stefan Ask;  Department of Physics, U"
1525      << "niversity of Cambridge,               |  | \n"
1526      << " |  |      Cavendish Laboratory, JJ Thomson A"
1527      << "ve., Cambridge CB3 0HE, UK;           |  | \n"
1528      << " |  |      phone: + 41 - 22 - 767 6707; e-mai"
1529      << "l: Stefan.Ask@cern.ch                 |  | \n"
1530      << " |  |   Richard Corke;  Niels Bohr Institute,"
1531      << " University of Copenhagen,            |  | \n"
1532      << " |  |      Blegdamsvej 17, 2100 Copenhagen, D"
1533      << "enmark;                               |  | \n"
1534      << " |  |      e-mail: richard.corke@thep.lu.se  "
1535      << "                                      |  | \n"
1536      << " |  |   Stephen Mrenna;  Computing Division, "
1537      << "Simulations Group,                    |  | \n"
1538      << " |  |      Fermi National Accelerator Laborat"
1539      << "ory, MS 234, Batavia, IL 60510, USA;  |  | \n"
1540      << " |  |      phone: + 1 - 630 - 840 - 2556; e-m"
1541      << "ail: mrenna@fnal.gov                  |  | \n"
1542      << " |  |   Stefan Prestel;  Department of Astron"
1543      << "omy and Theoretical Physics,          |  | \n"
1544      << " |  |      Lund University, Solvegatan 14A, S"
1545      << "E-223 62 Lund, Sweden;                |  | \n"
1546      << " |  |      phone: + 46 - 46 - 222 06 64; e-ma"
1547      << "il: stefan.prestel@thep.lu.se         |  | \n"
1548      << " |  |   Peter Skands;  Theoretical Physics, C"
1549      << "ERN, CH-1211 Geneva 23, Switzerland;  |  | \n"
1550      << " |  |      phone: + 41 - 22 - 767 2447; e-mai"
1551      << "l: peter.skands@cern.ch               |  | \n"
1552      << " |  |                                        "
1553      << "                                      |  | \n"
1554      << " |  |   The main program reference is the 'Br"
1555      << "ief Introduction to PYTHIA 8.1',      |  | \n"
1556      << " |  |   T. Sjostrand, S. Mrenna and P. Skands"
1557      << ", Comput. Phys. Comm. 178 (2008) 85   |  | \n"
1558      << " |  |   [arXiv:0710.3820]                    "
1559      << "                                      |  | \n"
1560      << " |  |                                        "
1561      << "                                      |  | \n"
1562      << " |  |   The main physics reference is the 'PY"
1563      << "THIA 6.4 Physics and Manual',         |  | \n"
1564      << " |  |   T. Sjostrand, S. Mrenna and P. Skands"
1565      << ", JHEP05 (2006) 026 [hep-ph/0603175]. |  | \n"
1566      << " |  |                                        "
1567      << "                                      |  | \n"
1568      << " |  |   An archive of program versions and do"
1569      << "cumentation is found on the web:      |  | \n"
1570      << " |  |   http://www.thep.lu.se/~torbjorn/Pythi"
1571      << "a.html                                |  | \n"
1572      << " |  |                                        "
1573      << "                                      |  | \n"
1574      << " |  |   This program is released under the GN"
1575      << "U General Public Licence version 2.   |  | \n"
1576      << " |  |   Please respect the MCnet Guidelines f"
1577      << "or Event Generator Authors and Users. |  | \n"
1578      << " |  |                                        "
1579      << "                                      |  | \n"
1580      << " |  |   Disclaimer: this program comes withou"
1581      << "t any guarantees.                     |  | \n"
1582      << " |  |   Beware of errors and use common sense"
1583      << " when interpreting results.           |  | \n"
1584      << " |  |                                        "
1585      << "                                      |  | \n"
1586      << " |  |   Copyright (C) 2012 Torbjorn Sjostrand"
1587      << "                                      |  | \n"
1588      << " |  |                                        "
1589      << "                                      |  | \n"
1590      << " |  |                                        "
1591      << "                                      |  | \n"
1592      << " |  *----------------------------------------"
1593      << "--------------------------------------*  | \n"
1594      << " |                                           "
1595      << "                                         | \n"
1596      << " *-------------------------------------------"
1597      << "-----------------------------------------* \n" << endl;
1598
1599 }
1600
1601 //--------------------------------------------------------------------------
1602
1603 // Check for lines in file that mark the beginning of new subrun.
1604
1605 int Pythia::readSubrun(string line, bool warn, ostream& os) {
1606
1607   // If empty line then done.
1608   int subrunLine = SUBRUNDEFAULT;
1609   if (line.find_first_not_of(" \n\t\v\b\r\f\a") == string::npos)
1610     return subrunLine;
1611
1612   // If first character is not a letter, then done.
1613   string lineNow = line;
1614   int firstChar = lineNow.find_first_not_of(" \n\t\v\b\r\f\a");
1615   if (!isalpha(lineNow[firstChar])) return subrunLine;
1616
1617   // Replace an equal sign by a blank to make parsing simpler.
1618   while (lineNow.find("=") != string::npos) {
1619     int firstEqual = lineNow.find_first_of("=");
1620     lineNow.replace(firstEqual, 1, " ");
1621   }
1622
1623   // Get first word of a line.
1624   istringstream splitLine(lineNow);
1625   string name;
1626   splitLine >> name;
1627
1628   // Replace two colons by one (:: -> :) to allow for such mistakes.
1629   while (name.find("::") != string::npos) {
1630     int firstColonColon = name.find_first_of("::");
1631     name.replace(firstColonColon, 2, ":");
1632   }
1633
1634   // Convert to lowercase.
1635   for (int i = 0; i < int(name.length()); ++i)
1636     name[i] = std::tolower(name[i]);
1637
1638   // If no match then done.
1639   if (name != "main:subrun") return subrunLine;
1640
1641   // Else find new subrun number and return it.
1642   splitLine >> subrunLine;
1643   if (!splitLine) {
1644     if (warn) os << "\n PYTHIA Warning: Main:subrun number not"
1645         << " recognized; skip:\n   " << line << endl;
1646     subrunLine = SUBRUNDEFAULT;
1647   }
1648   return subrunLine;
1649
1650 }
1651
1652 //--------------------------------------------------------------------------
1653
1654 // Check that the final event makes sense: no unknown id codes;
1655 // charge and energy-momentum conserved.
1656
1657 bool Pythia::check(ostream& os) {
1658
1659   // Reset.
1660   bool physical     = true;
1661   bool listVertices = false;
1662   bool listHistory  = false;
1663   bool listSystems  = false;
1664   bool listBeams    = false;
1665   iErrId.resize(0);
1666   iErrCol.resize(0);
1667   iErrNan.resize(0);
1668   iErrNanVtx.resize(0);
1669   Vec4 pSum;
1670   double chargeSum  = 0.;
1671
1672   // Incoming beams counted with negative momentum and charge.
1673   if (doProcessLevel) {
1674     pSum      = - (event[1].p() + event[2].p());
1675     chargeSum = - (event[1].charge() + event[2].charge());
1676
1677   // If no ProcessLevel then sum final state of process record.
1678   } else if (process.size() > 0) {
1679     pSum = - process[0].p();
1680     for (int i = 0; i < process.size(); ++i)
1681       if (process[i].isFinal()) chargeSum -= process[i].charge();
1682
1683   // If process not filled, then use outgoing primary in event.
1684   } else {
1685     pSum = - event[0].p();
1686     for (int i = 1; i < event.size(); ++i)
1687       if (event[i].statusAbs() < 10 || event[i].statusAbs() == 23)
1688         chargeSum -= event[i].charge();
1689   }
1690   double eLab = abs(pSum.e());
1691
1692   // Loop over particles in the event.
1693   for (int i = 0; i < event.size(); ++i) {
1694
1695     // Look for any unrecognized particle codes.
1696     int id = event[i].id();
1697     if (id == 0 || !particleData.isParticle(id)) {
1698       ostringstream errCode;
1699       errCode << ", i = " << i << ", id = " << id;
1700       info.errorMsg("Error in Pythia::check: "
1701         "unknown particle code", errCode.str());
1702       physical = false;
1703       iErrId.push_back(i);
1704
1705     // Check that colour assignments are the expected ones.
1706     } else {
1707       int colType = event[i].colType();
1708       int col     = event[i].col();
1709       int acol    = event[i].acol();
1710       if ( (colType ==  0 && (col  > 0 || acol  > 0))
1711         || (colType ==  1 && (col <= 0 || acol  > 0))
1712         || (colType == -1 && (col  > 0 || acol <= 0))
1713         || (colType ==  2 && (col <= 0 || acol <= 0)) ) {
1714         ostringstream errCode;
1715         errCode << ", i = " << i << ", id = " << id << " cols = " << col
1716                 << " " << acol;
1717         info.errorMsg("Error in Pythia::check: "
1718           "incorrect colours", errCode.str());
1719         physical = false;
1720         iErrCol.push_back(i);
1721       }
1722     }
1723
1724     // Look for particles with not-a-number energy/momentum/mass.
1725     if (abs(event[i].px()) >= 0. && abs(event[i].py()) >= 0.
1726       && abs(event[i].pz()) >= 0.  && abs(event[i].e()) >= 0.
1727       && abs(event[i].m()) >= 0.) ;
1728     else {
1729       info.errorMsg("Error in Pythia::check: "
1730         "not-a-number energy/momentum/mass");
1731       physical = false;
1732       iErrNan.push_back(i);
1733     }
1734
1735     // Look for particles with not-a-number vertex/lifetime.
1736     if (abs(event[i].xProd()) >= 0. && abs(event[i].yProd()) >= 0.
1737       && abs(event[i].zProd()) >= 0.  && abs(event[i].tProd()) >= 0.
1738       && abs(event[i].tau()) >= 0.) ;
1739     else {
1740       info.errorMsg("Error in Pythia::check: "
1741         "not-a-number vertex/lifetime");
1742       physical     = false;
1743       listVertices = true;
1744       iErrNanVtx.push_back(i);
1745     }
1746
1747     // Add final-state four-momentum and charge.
1748     if (event[i].isFinal()) {
1749       pSum      += event[i].p();
1750       chargeSum += event[i].charge();
1751     }
1752
1753   // End of particle loop.
1754   }
1755
1756   // Check energy-momentum/charge conservation.
1757   double epDev = abs(pSum.e()) + abs(pSum.px()) + abs(pSum.py())
1758     + abs(pSum.pz());
1759   if (epDev > epTolErr * eLab) {
1760     info.errorMsg("Error in Pythia::check: energy-momentum not conserved");
1761     physical = false;
1762   } else if (epDev > epTolWarn * eLab) {
1763     info.errorMsg("Warning in Pythia::check: "
1764       "energy-momentum not quite conserved");
1765   }
1766   if (abs(chargeSum) > 0.1) {
1767     info.errorMsg("Error in Pythia::check: charge not conserved");
1768     physical = false;
1769   }
1770
1771   // Check that beams and event records agree on incoming partons.
1772   // Only meaningful for resolved beams.
1773   if (info.isResolved())
1774   for (int iSys = 0; iSys < beamA.sizeInit(); ++iSys) {
1775     int eventANw  = partonSystems.getInA(iSys);
1776     int eventBNw  = partonSystems.getInB(iSys);
1777     int beamANw   = beamA[iSys].iPos();
1778     int beamBNw   = beamB[iSys].iPos();
1779     if (eventANw != beamANw || eventBNw != beamBNw) {
1780       info.errorMsg("Error in Pythia::check: "
1781         "event and beams records disagree");
1782       physical    = false;
1783       listSystems = true;
1784       listBeams   = true;
1785     }
1786   }
1787
1788   // Check that mother and daughter information match for each particle.
1789   vector<int> noMot;
1790   vector<int> noDau;
1791   vector< pair<int,int> > noMotDau;
1792   if (checkHistory) {
1793
1794     // Loop through the event and check that there are beam particles.
1795     bool hasBeams = false;
1796     for (int i = 0; i < event.size(); ++i) {
1797       int status = event[i].status();
1798       if (abs(status) == 12) hasBeams = true;
1799
1800       // Check that mother and daughter lists not empty where not expected to.
1801       vector<int> mList = event.motherList(i);
1802       vector<int> dList = event.daughterList(i);
1803       if (mList.size() == 0 && abs(status) != 11 && abs(status) != 12)
1804         noMot.push_back(i);
1805       if (dList.size() == 0 && status < 0 && status != -11)
1806         noDau.push_back(i);
1807
1808       // Check that the particle appears in the daughters list of each mother.
1809       for (int j = 0; j < int(mList.size()); ++j) {
1810         vector<int> dmList = event.daughterList( mList[j] );
1811         bool foundMatch = false;
1812         for (int k = 0; k < int(dmList.size()); ++k)
1813         if (dmList[k] == i) {
1814           foundMatch = true;
1815           break;
1816         }
1817         if (!hasBeams && mList.size() == 1 && mList[0] == 0) foundMatch = true;
1818         if (!foundMatch) {
1819           bool oldPair = false;
1820           for (int k = 0; k < int(noMotDau.size()); ++k)
1821           if (noMotDau[k].first == mList[j] && noMotDau[k].second == i) {
1822             oldPair = true;
1823             break;
1824           }
1825           if (!oldPair) noMotDau.push_back( make_pair( mList[j], i) );
1826         }
1827       }
1828
1829       // Check that the particle appears in the mothers list of each daughter.
1830       for (int j = 0; j < int(dList.size()); ++j) {
1831         vector<int> mdList = event.motherList( dList[j] );
1832         bool foundMatch = false;
1833         for (int k = 0; k < int(mdList.size()); ++k)
1834         if (mdList[k] == i) {
1835           foundMatch = true;
1836           break;
1837         }
1838         if (!foundMatch) {
1839           bool oldPair = false;
1840           for (int k = 0; k < int(noMotDau.size()); ++k)
1841           if (noMotDau[k].first == i && noMotDau[k].second == dList[j]) {
1842             oldPair = true;
1843             break;
1844           }
1845           if (!oldPair) noMotDau.push_back( make_pair( i, dList[j]) );
1846         }
1847       }
1848     }
1849
1850     // Warn if any errors were found.
1851     if (noMot.size() > 0 || noDau.size() > 0 || noMotDau.size() > 0) {
1852       info.errorMsg("Error in Pythia::check: "
1853         "mismatch in daughter and mother lists");
1854       physical    = false;
1855       listHistory = true;
1856     }
1857   }
1858
1859   // Done for sensible events.
1860   if (physical) return true;
1861
1862   // Print (the first few) flawed events: local info.
1863   if (nErrEvent < nErrList) {
1864     os << "\n PYTHIA erroneous event info: \n";
1865     if (iErrId.size() > 0) {
1866       os << " unknown particle codes in lines ";
1867       for (int i = 0; i < int(iErrId.size()); ++i)
1868         os << iErrId[i] << " ";
1869       os << "\n";
1870     }
1871     if (iErrCol.size() > 0) {
1872       os << " incorrect colour assignments in lines ";
1873       for (int i = 0; i < int(iErrCol.size()); ++i)
1874         os << iErrCol[i] << " ";
1875       os << "\n";
1876     }
1877     if (iErrNan.size() > 0) {
1878       os << " not-a-number energy/momentum/mass in lines ";
1879       for (int i = 0; i < int(iErrNan.size()); ++i)
1880         os << iErrNan[i] << " ";
1881       os << "\n";
1882     }
1883     if (iErrNanVtx.size() > 0) {
1884       os << " not-a-number vertex/lifetime in lines ";
1885       for (int i = 0; i < int(iErrNanVtx.size()); ++i)
1886         os << iErrNanVtx[i] << " ";
1887       os << "\n";
1888     }
1889     if (epDev > epTolErr * eLab) os << scientific << setprecision(3)
1890       << " total energy-momentum non-conservation = " << epDev << "\n";
1891     if (abs(chargeSum) > 0.1) os << fixed << setprecision(2)
1892       << " total charge non-conservation = " << chargeSum << "\n";
1893     if (noMot.size() > 0) {
1894       os << " missing mothers for particles ";
1895       for (int i = 0; i < int(noMot.size()); ++i) os << noMot[i] << " ";
1896       os << "\n";
1897     }
1898     if (noDau.size() > 0) {
1899       os << " missing daughters for particles ";
1900       for (int i = 0; i < int(noDau.size()); ++i) os << noDau[i] << " ";
1901       os << "\n";
1902     }
1903     if (noMotDau.size() > 0) {
1904       os << " inconsistent history for (mother,daughter) pairs ";
1905       for (int i = 0; i < int(noMotDau.size()); ++i)
1906         os << "(" << noMotDau[i].first << "," << noMotDau[i].second << ") ";
1907       os << "\n";
1908     }
1909
1910     // Print (the first few) flawed events: standard listings.
1911     info.list();
1912     event.list(listVertices, listHistory);
1913     if (listSystems) partonSystems.list();
1914     if (listBeams) beamA.list();
1915     if (listBeams) beamB.list();
1916   }
1917
1918   // Update error counter. Done also for flawed event.
1919   ++nErrEvent;
1920   return false;
1921
1922 }
1923
1924 //--------------------------------------------------------------------------
1925
1926 // Routine to set up a PDF pointer.
1927
1928 PDF* Pythia::getPDFPtr(int idIn, int sequence) {
1929
1930   // Temporary pointer to be returned.
1931   PDF* tempPDFPtr = 0;
1932
1933   // One option is to treat a Pomeron like a pi0.
1934   if (idIn == 990 && settings.mode("PDF:PomSet") == 2) idIn = 111;
1935
1936   // Proton beam, normal choice.
1937   if (abs(idIn) == 2212 && sequence == 1) {
1938     int  pSet      = settings.mode("PDF:pSet");
1939     bool useLHAPDF = settings.flag("PDF:useLHAPDF");
1940
1941     // Use internal sets.
1942     if (!useLHAPDF) {
1943       if      (pSet == 1) tempPDFPtr = new GRV94L(idIn);
1944       else if (pSet == 2) tempPDFPtr = new CTEQ5L(idIn);
1945       else if (pSet <= 6)
1946            tempPDFPtr = new MSTWpdf(idIn, pSet - 2, xmlPath, &info);
1947       else tempPDFPtr = new CTEQ6pdf(idIn, pSet - 6, xmlPath, &info);
1948     }
1949
1950     // Use sets from LHAPDF.
1951     else {
1952       string LHAPDFset    = settings.word("PDF:LHAPDFset");
1953       int    LHAPDFmember = settings.mode("PDF:LHAPDFmember");
1954       tempPDFPtr = new LHAPDF(idIn, LHAPDFset, LHAPDFmember, 1, &info);
1955
1956       // Optionally allow extrapolation beyond x and Q2 limits.
1957       tempPDFPtr->setExtrapolate( settings.flag("PDF:extrapolateLHAPDF") );
1958     }
1959   }
1960
1961   // Proton beam, special choice for the hard process.
1962   else if (abs(idIn) == 2212) {
1963     int  pSet      = settings.mode("PDF:pHardSet");
1964     bool useLHAPDF = settings.flag("PDF:useHardLHAPDF");
1965
1966     // Use internal sets.
1967     if (!useLHAPDF) {
1968       if      (pSet == 1) tempPDFPtr = new GRV94L(idIn);
1969       else if (pSet == 2) tempPDFPtr = new CTEQ5L(idIn);
1970       else if (pSet <= 6)
1971            tempPDFPtr = new MSTWpdf(idIn, pSet - 2, xmlPath, &info);
1972       else tempPDFPtr = new CTEQ6pdf(idIn, pSet - 6, xmlPath, &info);
1973     }
1974
1975     // Use sets from LHAPDF.
1976     else {
1977       string LHAPDFset    = settings.word("PDF:hardLHAPDFset");
1978       int    LHAPDFmember = settings.mode("PDF:hardLHAPDFmember");
1979       // May need to require LHAPDF to handle two sets simultaneously.
1980      int nSet = (settings.flag("PDF:useLHAPDF")) ? 2 : 1;
1981       tempPDFPtr = new LHAPDF(idIn, LHAPDFset, LHAPDFmember, nSet, &info);
1982
1983      // Optionally allow extrapolation beyond x and Q2 limits.
1984       tempPDFPtr->setExtrapolate( settings.flag("PDF:extrapolateLHAPDF") );
1985     }
1986   }
1987
1988   // Pion beam (or, in one option, Pomeron beam).
1989   else if (abs(idIn) == 211 || idIn == 111) {
1990     bool useLHAPDF = settings.flag("PDF:piUseLHAPDF");
1991
1992     // Use internal sets.
1993     if (!useLHAPDF) {
1994       tempPDFPtr = new GRVpiL(idIn);
1995     }
1996
1997     // Use sets from LHAPDF.
1998     else {
1999       string LHAPDFset    = settings.word("PDF:piLHAPDFset");
2000       int    LHAPDFmember = settings.mode("PDF:piLHAPDFmember");
2001       tempPDFPtr = new LHAPDF(idIn, LHAPDFset, LHAPDFmember, 1, &info);
2002
2003       // Optionally allow extrapolation beyond x and Q2 limits.
2004       tempPDFPtr->setExtrapolate( settings.flag("PDF:extrapolateLHAPDF") );
2005     }
2006   }
2007
2008   // Pomeron beam, if not treated like a pi0 beam.
2009   else if (idIn == 990) {
2010     int    pomSet  = settings.mode("PDF:PomSet");
2011     double rescale = settings.parm("PDF:PomRescale");
2012
2013     // A generic Q2-independent parametrization.
2014     if (pomSet == 1) {
2015       double gluonA      = settings.parm("PDF:PomGluonA");
2016       double gluonB      = settings.parm("PDF:PomGluonB");
2017       double quarkA      = settings.parm("PDF:PomQuarkA");
2018       double quarkB      = settings.parm("PDF:PomQuarkB");
2019       double quarkFrac   = settings.parm("PDF:PomQuarkFrac");
2020       double strangeSupp = settings.parm("PDF:PomStrangeSupp");
2021       tempPDFPtr = new PomFix( 990, gluonA, gluonB, quarkA, quarkB,
2022         quarkFrac, strangeSupp);
2023     }
2024
2025     // The H1 Q2-dependent parametrizations. Initialization requires files.
2026     else if (pomSet == 3 || pomSet == 4)
2027       tempPDFPtr = new PomH1FitAB( 990, pomSet - 2, rescale, xmlPath, &info);
2028     else if (pomSet == 5)
2029       tempPDFPtr = new PomH1Jets( 990, rescale, xmlPath, &info);
2030     else if (pomSet == 6)
2031       tempPDFPtr = new PomH1FitAB( 990, 3, rescale, xmlPath, &info);
2032   }
2033
2034   // Lepton beam: neutrino, resolved charged lepton or unresolved ditto.
2035   else if (abs(idIn) > 10 && abs(idIn) < 17) {
2036     if (abs(idIn)%2 == 0) tempPDFPtr = new NeutrinoPoint(idIn);
2037     else if (settings.flag("PDF:lepton")) tempPDFPtr = new Lepton(idIn);
2038     else tempPDFPtr = new LeptonPoint(idIn);
2039   }
2040
2041   // Failure for unrecognized particle.
2042   else tempPDFPtr = 0;
2043
2044   // Done.
2045   return tempPDFPtr;
2046 }
2047
2048 //--------------------------------------------------------------------------
2049
2050 // Initialize SUSY Les Houches Accord data.
2051
2052 bool Pythia::initSLHA() {
2053
2054   // Initial and settings values.
2055   int    ifailLHE    = 1;
2056   int    ifailSpc    = 1;
2057   int    readFrom    = settings.mode("SLHA:readFrom");
2058   string lhefFile    = settings.word("Beams:LHEF");
2059   string lhefHeader  = settings.word("Beams:LHEFheader");
2060   string slhaFile    = settings.word("SLHA:file");
2061   int    verboseSLHA = settings.mode("SLHA:verbose");
2062   bool   slhaUseDec  = settings.flag("SLHA:useDecayTable");
2063
2064   // No SUSY by default
2065   couplingsPtr->isSUSY = false;
2066
2067   // Option with no SLHA read-in at all.
2068   if (readFrom == 0) return true;
2069
2070   // First check LHEF header (if reading from LHEF)
2071   if (readFrom == 1) {
2072     if (lhefHeader != "void")
2073       ifailLHE = slha.readFile(lhefHeader, verboseSLHA, slhaUseDec );
2074     else if (lhefFile != "void")
2075       ifailLHE = slha.readFile(lhefFile, verboseSLHA, slhaUseDec );
2076   }
2077
2078   // If LHEF read successful, everything needed should already be ready
2079   if (ifailLHE == 0) {
2080     ifailSpc = 0;
2081     couplingsPtr->isSUSY = true;
2082   // If no LHEF file or no SLHA info in header, read from SLHA:file
2083   } else {
2084     lhefFile = "void";
2085     if ( settings.word("SLHA:file") == "none"
2086          || settings.word("SLHA:file") == "void"
2087          || settings.word("SLHA:file") == ""
2088          || settings.word("SLHA:file") == " ") return true;
2089     ifailSpc = slha.readFile(slhaFile,verboseSLHA, slhaUseDec);
2090   }
2091
2092   // In case of problems, print error and fail init.
2093   if (ifailSpc != 0) {
2094     info.errorMsg("Error in Pythia::initSLHA: "
2095       "problem reading SLHA file", slhaFile);
2096     return false;
2097   } else {
2098     couplingsPtr->isSUSY = true;
2099   }
2100
2101   // Check spectrum for consistency. Switch off SUSY if necessary.
2102   ifailSpc = slha.checkSpectrum();
2103
2104   // ifail >= 1 : no MODSEL found -> don't switch on SUSY
2105   if (ifailSpc == 1) {
2106     // no SUSY, but MASS ok
2107     couplingsPtr->isSUSY = false;
2108     info.errorMsg("Warning in Pythia::initSLHA: "
2109                   "No MODSEL found, keeping internal SUSY switched off");
2110   } else if (ifailSpc >= 2) {
2111     // no SUSY, but problems
2112     info.errorMsg("Warning in Pythia::initSLHA: "
2113                       "Problem with SLHA MASS or QNUMBERS.");
2114     couplingsPtr->isSUSY = false;
2115   }
2116   // ifail = 0 : MODSEL found, spectrum OK
2117   else if (ifailSpc == 0) {
2118     // Print spectrum. Done.
2119     slha.printSpectrum(0);
2120   }
2121   else if (ifailSpc < 0) {
2122     info.errorMsg("Warning in Pythia::initSLHA: "
2123                       "Problem with SLHA spectrum.",
2124                       "\n Only using masses and switching off SUSY.");
2125     settings.flag("SUSY:all", false);
2126     couplingsPtr->isSUSY = false;
2127     slha.printSpectrum(ifailSpc);
2128   }
2129
2130   // Import qnumbers
2131   if ( (ifailSpc == 1 || ifailSpc == 0) && slha.qnumbers.size() > 0) {
2132     for (int iQnum=0; iQnum < int(slha.qnumbers.size()); iQnum++) {
2133       // Always use positive id codes
2134       int id = abs(slha.qnumbers[iQnum](0));
2135       ostringstream idCode;
2136       idCode << id;
2137       if (particleData.isParticle(id)) {
2138         info.errorMsg("Warning in Pythia::initSLHA: "
2139                       "ignoring QNUMBERS", "for id = "+idCode.str()
2140                       +" (already exists)", true);
2141       } else {
2142         int qEM3    = slha.qnumbers[iQnum](1);
2143         int nSpins  = slha.qnumbers[iQnum](2);
2144         int colRep  = slha.qnumbers[iQnum](3);
2145         int hasAnti = slha.qnumbers[iQnum](4);
2146         // Translate colRep to PYTHIA colType
2147         int colType = 0;
2148         if (colRep == 3) colType = 1;
2149         else if (colRep == -3) colType = -1;
2150         else if (colRep == 8) colType = 2;
2151         else if (colRep == 6) colType = 3;
2152         else if (colRep == -6) colType = -3;
2153         // Default name: PDG code
2154         string name, antiName;
2155         ostringstream idStream;
2156         idStream<<id;
2157         name     = idStream.str();
2158         antiName = "-"+name;
2159         if (iQnum < int(slha.qnumbersName.size())) {
2160           name = slha.qnumbersName[iQnum];
2161           if (iQnum < int(slha.qnumbersAntiName.size()))
2162             antiName = slha.qnumbersAntiName[iQnum];
2163           if (antiName == "") {
2164             if (name.find("+") != string::npos) {
2165               antiName = name;
2166               antiName.replace(antiName.find("+"),1,"-");
2167             } else if (name.find("-") != string::npos) {
2168               antiName = name;
2169               antiName.replace(antiName.find("-"),1,"+");
2170             } else {
2171               antiName = name+"bar";
2172             }
2173           }
2174         }
2175         if ( hasAnti == 0) {
2176           antiName = "";
2177           particleData.addParticle(id, name, nSpins, qEM3, colType);
2178         } else {
2179           particleData.addParticle(id, name, antiName, nSpins, qEM3, colType);
2180         }
2181       }
2182     }
2183   }
2184
2185   // Import mass spectrum.
2186   bool   keepSM     = settings.flag("SLHA:keepSM");
2187   double minMassSM  = settings.parm("SLHA:minMassSM");
2188   double massMargin = settings.parm("SLHA:minDecayDeltaM");
2189   if (ifailSpc == 1 || ifailSpc == 0) {
2190
2191     // Loop through to update particle data.
2192     int    id = slha.mass.first();
2193     for (int i = 1; i <= slha.mass.size() ; i++) {
2194       double mass = abs(slha.mass(id));
2195
2196       // Ignore masses for known SM particles or particles with
2197       // default masses < minMassSM; overwrite masses for rest.
2198       if (keepSM && (id < 25 || (id > 80 && id < 1000000))) ;
2199       else if (id < 1000000 && particleData.m0(id) < minMassSM) {
2200         ostringstream idCode;
2201         idCode << id;
2202         info.errorMsg("Warning in Pythia::initSLHA: "
2203           "ignoring MASS entry", "for id = "+idCode.str()
2204           +" (m0 < SLHA:minMassSM)", true);
2205       }
2206       else particleData.m0(id,mass);
2207       id = slha.mass.next();
2208     };
2209
2210   }
2211
2212   // Update decay data.
2213   for (int iTable=0; iTable < int(slha.decays.size()); iTable++) {
2214
2215     // Pointer to this SLHA table
2216     LHdecayTable* slhaTable=&(slha.decays[iTable]);
2217
2218     // Extract ID and create pointer to corresponding particle data object
2219     int idRes     = slhaTable->getId();
2220     ParticleDataEntry* particlePtr
2221       = particleData.particleDataEntryPtr(idRes);
2222
2223     // Ignore decay channels for known SM particles or particles with
2224     // default masses < minMassSM; overwrite masses for rest.
2225     if (keepSM && (idRes < 25 || (idRes > 80 && idRes < 1000000))) continue;
2226     else if (idRes < 1000000 && particleData.m0(idRes) < minMassSM) {
2227       ostringstream idCode;
2228       idCode << idRes;
2229       info.errorMsg("Warning in Pythia::initSLHA: "
2230         "ignoring DECAY table", "for id = " + idCode.str()
2231         + " (m0 < SLHA:minMassSM)", true);
2232       continue;
2233     }
2234
2235     // Extract and store total width (absolute value, neg -> switch off)
2236     double widRes = abs(slhaTable->getWidth());
2237     particlePtr->setMWidth(widRes);
2238
2239     // Reset decay table of the particle. Allow decays, treat as resonance.
2240     if (slhaTable->size() > 0) {
2241       particlePtr->clearChannels();
2242       particleData.mayDecay(idRes,true);
2243       particleData.isResonance(idRes,true);
2244     }
2245
2246     // Reset to stable if width <= 0.0
2247     if (slhaTable->getWidth() <= 0.0) particleData.mayDecay(idRes,false);
2248
2249     // Set initial minimum mass.
2250     double brWTsum   = 0.;
2251     double massWTsum = 0.;
2252
2253     // Loop over SLHA channels, import into Pythia, treating channels
2254     // with negative branching fractions as having the equivalent positive
2255     // branching fraction, but being switched off for this run
2256     for (int iChannel=0 ; iChannel<slhaTable->size(); iChannel++) {
2257       LHdecayChannel slhaChannel = slhaTable->getChannel(iChannel);
2258       double brat      = slhaChannel.getBrat();
2259       vector<int> idDa = slhaChannel.getIdDa();
2260       if (idDa.size() >= 9) {
2261         info.errorMsg("Error in Pythia::initSLHA: "
2262                           "max number of decay products is 8.");
2263       } else if (idDa.size() <= 1) {
2264         info.errorMsg("Error in Pythia::initSLHA: "
2265                           "min number of decay products is 2.");
2266       }
2267       else {
2268         int onMode = 1;
2269         if (brat < 0.0) onMode = 0;
2270
2271         // Check phase space, including margin
2272         double massSum = massMargin;
2273         for (int jDa=0; jDa<int(idDa.size()); ++jDa)
2274           massSum += particleData.m0( idDa[jDa] );
2275         if (onMode == 1 && brat > 0.0 && massSum > particleData.m0(idRes) ) {
2276           // String containing decay name
2277           ostringstream errCode;
2278           errCode << idRes <<" ->";
2279           for (int jDa=0; jDa<int(idDa.size()); ++jDa) errCode<<" "<<idDa[jDa];
2280           info.errorMsg("Warning in Pythia::initSLHA: switching off decay",
2281             errCode.str() + " (mRes - mDa < minDecayDeltaM)"
2282             "\n       (Note: cross sections will be scaled by remaining"
2283             " open branching fractions!)" , true);
2284           onMode=0;
2285         }
2286
2287         // Branching-ratio-weighted average mass in decay.
2288         brWTsum   += abs(brat);
2289         massWTsum += abs(brat) * massSum;
2290
2291         // Add channel
2292         int id0 = idDa[0];
2293         int id1 = idDa[1];
2294         int id2 = (idDa.size() >= 3) ? idDa[2] : 0;
2295         int id3 = (idDa.size() >= 4) ? idDa[3] : 0;
2296         int id4 = (idDa.size() >= 5) ? idDa[4] : 0;
2297         int id5 = (idDa.size() >= 6) ? idDa[5] : 0;
2298         int id6 = (idDa.size() >= 7) ? idDa[6] : 0;
2299         int id7 = (idDa.size() >= 8) ? idDa[7] : 0;
2300         particlePtr->addChannel(onMode,abs(brat),101,
2301                                       id0,id1,id2,id3,id4,id5,id6,id7);
2302
2303       }
2304     }
2305
2306     // Set minimal mass, but always below nominal one.
2307     if (slhaTable->size() > 0) {
2308       double massAvg = massWTsum / brWTsum;
2309       double massMin = min( massAvg, particlePtr->m0()) - massMargin;
2310       particlePtr->setMMin(massMin);
2311     }
2312   }
2313
2314   return true;
2315
2316 }
2317
2318 //--------------------------------------------------------------------------
2319
2320 // Function to perform CKKW-L merging on the input event.
2321
2322 bool Pythia::mergeProcess() {
2323
2324   // Reset weight of the event.
2325   mergingHooksPtr->setWeightCKKWL(1.);
2326   info.setWeightCKKWL(1.);
2327   double wgt = 1.0;
2328
2329   // Store candidates for the splitting V -> qqbar'.
2330   mergingHooksPtr->storeHardProcessCandidates( process);
2331
2332   // Calculate number of clustering steps.
2333   int nSteps   = mergingHooksPtr->getNumberOfClusteringSteps( process);
2334
2335   // Check if event passes the merging scale cut.
2336   double tms   = mergingHooksPtr->tms();
2337   // Get merging scale in current event.
2338   double tnow  = 0.;
2339   // Use KT/Durham merging scale definition.
2340   if ( mergingHooksPtr->doKTMerging() || mergingHooksPtr->doMGMerging() )
2341     tnow       = mergingHooksPtr->kTms( process );
2342   // Use Lund PT merging scale definition.
2343   else if ( mergingHooksPtr->doPTLundMerging() )
2344     tnow       = mergingHooksPtr->rhoms( process, false );
2345    // Use DeltaR_{ij}, pT_i, Q_{ij} combination merging scale definition.
2346   else if ( mergingHooksPtr->doCutBasedMerging() )
2347     tnow       = mergingHooksPtr->cutbasedms( process );
2348   // Use user-defined merging scale.
2349   else
2350     tnow       = mergingHooksPtr->tmsDefinition( process );
2351
2352   bool enforceCutOnLHE  = settings.flag("Merging:enforceCutOnLHE");
2353   // Enfore merging scale cut if the event did not pass the merging scale
2354   // criterion.
2355   if ( enforceCutOnLHE && nSteps > 0 && tnow < tms ) {
2356     string message="Warning in Pythia::mergeProcess: Les Houches Event";
2357     message+=" fails merging scale cut. Cut by rejecting event.";
2358     info.errorMsg(message);
2359     return true;
2360   }
2361
2362   double TMSMISMATCHPOS = tms / 2. ;
2363   // Remember if event is significantly above the merging scale.
2364   if ( enforceCutOnLHE && nSteps > 0 && tms > 0.
2365     && tnow > tms && (tnow-tms > TMSMISMATCHPOS) ) {
2366     stringstream tmsmis;
2367     tmsmis << int(TMSMISMATCHPOS);
2368     string message="Warning in Pythia::mergeProcess: Les Houches event";
2369     message+=" more than ";
2370     message+=tmsmis.str();
2371     message+=" GeV above desired merging scale value.";
2372     info.errorMsg(message);
2373     info.addCounter(41);
2374   }
2375
2376   // For now, prefer construction of ordered histories.
2377   mergingHooksPtr->orderHistories(true);
2378
2379   // Declare pdfWeight and alpha_s weight.
2380   double RN = rndm.flat();
2381
2382   process.scale(0.0);
2383   // Generate all histories.
2384   History FullHistory( nSteps, 0.0, process, Clustering(), mergingHooksPtr,
2385             beamA, beamB, &particleData, &info, true, true, true, true,
2386             1.0, 0);
2387
2388   // Project histories onto desired branches, e.g. only ordered paths.
2389   FullHistory.projectOntoDesiredHistories();
2390
2391   // Calculate CKKWL weight:
2392   // Perform reweighting with Sudakov factors, save alpha_s ratios and
2393   // PDF ratio weights.
2394   wgt = FullHistory.weightTREE( &trialPartonLevel,
2395       mergingHooksPtr->AlphaS_FSR(), mergingHooksPtr->AlphaS_ISR(), RN);
2396
2397   // Event with production scales set for further (trial) showering
2398   // and starting conditions for the shower.
2399   FullHistory.getStartingConditions( RN, process );
2400
2401   // Allow to dampen histories in which the lowest multiplicity reclustered
2402   // state does not pass the lowest multiplicity cut of the matrix element.
2403   double dampWeight = mergingHooksPtr->dampenIfFailCuts(
2404     FullHistory.lowestMultProc(RN) );
2405
2406   // Save the weight of the event for histogramming. Only change the
2407   // event weight after trial shower on the matrix element
2408   // multiplicity event (= in doVetoStep).
2409   wgt *= dampWeight;
2410
2411   // Set QCD 2->2 starting scale different from arbitrary scale in LHEF!
2412   // --> Set to minimal mT of partons.
2413   int nFinal = 0;
2414   double muf = process[0].e();
2415   for ( int i=0; i < process.size(); ++i )
2416   if ( process[i].isFinal()
2417     && (process[i].colType() != 0 || process[i].id() == 22 ) ) {
2418     nFinal++;
2419     muf = min( muf, abs(process[i].mT()) );
2420   }
2421   // For pure QCD dijet events (only!), set the process scale to the
2422   // transverse momentum of the outgoing partons.
2423   if ( nSteps == 0 && nFinal == 2
2424     && ( mergingHooksPtr->getProcessString().compare("pp>jj") == 0
2425       || mergingHooksPtr->getProcessString().compare("pp>aj") == 0) )
2426     process.scale(muf);
2427
2428   // Save the weight of the event for histogramming.
2429   mergingHooksPtr->setWeightCKKWL(wgt);
2430   info.setWeightCKKWL(wgt);
2431
2432   // If the weight of the event is zero, do not continue evolution.
2433   if(wgt == 0.) return true;
2434
2435   // Done
2436   return false;
2437
2438 }
2439
2440 //==========================================================================
2441
2442 } // end namespace Pythia8