PYTHIA8/pythia8145/src/PartonLevel.cxx

   1 // PartonLevel.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 Torbjorn Sjostrand.
   3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
   4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
   5
   6 // Function definitions (not found in the header) for the PartonLevel class.
   7
   8 #include "PartonLevel.h"
   9
  10 namespace Pythia8 {
  11
  12 //==========================================================================
  13
  14 // The PartonLevel class.
  15
  16 //--------------------------------------------------------------------------
  17
  18 // Constants: could be changed here if desired, but normally should not.
  19 // These are of technical nature, as described for each.
  20
  21 // Maximum number of tries to produce parton level from given input.
  22 const int PartonLevel::NTRY = 10;
  23
  24 //--------------------------------------------------------------------------
  25
  26 // Main routine to initialize the parton-level generation process.
  27
  28 bool PartonLevel::init( Info* infoPtrIn, Settings& settings,
  29   ParticleData* particleDataPtrIn, Rndm* rndmPtrIn,
  30   BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
  31   BeamParticle* beamPomAPtrIn, BeamParticle* beamPomBPtrIn,
  32   Couplings* couplingsPtrIn, PartonSystems* partonSystemsPtrIn,
  33   SigmaTotal* sigmaTotPtr, TimeShower* timesDecPtrIn, TimeShower* timesPtrIn,
  34   SpaceShower* spacePtrIn, UserHooks* userHooksPtrIn) {
  35
  36   // Store input pointers and modes for future use.
  37   infoPtr            = infoPtrIn;
  38   particleDataPtr    = particleDataPtrIn;
  39   rndmPtr            = rndmPtrIn;
  40   beamAPtr           = beamAPtrIn;
  41   beamBPtr           = beamBPtrIn;
  42   beamHadAPtr        = beamAPtr;
  43   beamHadBPtr        = beamBPtr;
  44   beamPomAPtr        = beamPomAPtrIn;
  45   beamPomBPtr        = beamPomBPtrIn;
  46   couplingsPtr       = couplingsPtrIn;
  47   partonSystemsPtr   = partonSystemsPtrIn;
  48   timesDecPtr        = timesDecPtrIn;
  49   timesPtr           = timesPtrIn;
  50   spacePtr           = spacePtrIn;
  51   userHooksPtr       = userHooksPtrIn;
  52
  53   // Min bias and single diffraction processes need special treatment.
  54   doMinBias          =  settings.flag("SoftQCD:all")
  55                      || settings.flag("SoftQCD:minBias");
  56   doDiffraction      =  settings.flag("SoftQCD:all")
  57                      || settings.flag("SoftQCD:singleDiffractive")
  58                      || settings.flag("SoftQCD:doubleDiffractive");
  59
  60   // Separate low-mass (unresolved) and high-mass (perturbative) diffraction.
  61   mMinDiff           = settings.parm("Diffraction:mMinPert");
  62   mWidthDiff         = settings.parm("Diffraction:mWidthPert");
  63   if (mMinDiff + mWidthDiff > infoPtr->eCM()) doDiffraction = false;
  64
  65   // Need MI initialization for soft QCD processes, even if only first MI.
  66   // But no need to initialize MI if never going to use it.
  67   doMI               = settings.flag("PartonLevel:MI");
  68   doMIMB             = doMI;
  69   doMISDA            = doMI;
  70   doMISDB            = doMI;
  71   doMIinit           = doMI;
  72   if (doMinBias || doDiffraction) doMIinit = true;
  73   if (!settings.flag("PartonLevel:all")) doMIinit = false;
  74
  75   // Flags for showers: ISR and FSR.
  76   doISR              = settings.flag("PartonLevel:ISR");
  77   bool FSR           = settings.flag("PartonLevel:FSR");
  78   bool FSRinProcess  = settings.flag("PartonLevel:FSRinProcess");
  79   bool interleaveFSR = settings.flag("TimeShower:interleave");
  80   doFSRduringProcess = FSR && FSRinProcess &&  interleaveFSR;
  81   doFSRafterProcess  = FSR && FSRinProcess && !interleaveFSR;
  82   doFSRinResonances  = FSR && settings.flag("PartonLevel:FSRinResonances");
  83
  84   // Some other flags.
  85   doRemnants         = settings.flag("PartonLevel:Remnants");
  86   doSecondHard       = settings.flag("SecondHard:generate");
  87
  88   // Flag if lepton beams, and if non-resolved ones. May change main flags.
  89   hasLeptonBeams  = ( beamAPtr->isLepton() || beamBPtr->isLepton() );
  90   hasPointLeptons = ( hasLeptonBeams
  91     && (beamAPtr->isUnresolved() || beamBPtr->isUnresolved() ) );
  92   if (hasLeptonBeams) {
  93     doMIMB     = false;
  94     doMISDA    = false;
  95     doMISDB    = false;
  96     doMIinit   = false;
  97   }
  98   if (hasPointLeptons) {
  99     doISR      = false;
 100     doRemnants = false;
 101   }
 102
 103   // Possibility to allow user veto during evolution.
 104   canVetoPT   = (userHooksPtr > 0) ? userHooksPtr->canVetoPT()   : false;
 105   pTvetoPT    = (canVetoPT)        ? userHooksPtr->scaleVetoPT() : -1.;
 106   canVetoStep = (userHooksPtr > 0) ? userHooksPtr->canVetoStep() : false;
 107   nVetoStep   = (canVetoStep)   ? userHooksPtr->numberVetoStep() : -1;
 108   canVetoMIStep = (userHooksPtr > 0) ? userHooksPtr->canVetoMIStep() : false;
 109   nVetoMIStep = (canVetoStep)   ? userHooksPtr->numberVetoMIStep() : -1;
 110
 111   // Possibility to set maximal shower scale in resonance decays.
 112   canSetScale = (userHooksPtr > 0) ? userHooksPtr->canSetResonanceScale()
 113                                    : false;
 114
 115   // Set info and initialize the respective program elements.
 116   timesPtr->init( beamAPtr, beamBPtr);
 117   if (doISR) spacePtr->init( beamAPtr, beamBPtr);
 118   doMIMB  =  multiMB.init( doMIinit, 0, infoPtr, settings, particleDataPtr,
 119     rndmPtr, beamAPtr, beamBPtr, couplingsPtr, partonSystemsPtr, sigmaTotPtr);
 120   if (doDiffraction) doMISDA = multiSDA.init( doMIinit, 1, infoPtr,
 121     settings, particleDataPtr, rndmPtr, beamPomBPtr, beamAPtr, couplingsPtr,
 122     partonSystemsPtr, sigmaTotPtr);
 123   if (doDiffraction) doMISDB = multiSDB.init( doMIinit, 2, infoPtr,
 124     settings, particleDataPtr, rndmPtr, beamPomAPtr, beamBPtr, couplingsPtr,
 125     partonSystemsPtr, sigmaTotPtr);
 126   remnants.init( infoPtr, settings, rndmPtr, beamAPtr, beamBPtr,
 127     partonSystemsPtr);
 128
 129   // Succeeded, or not.
 130   multiPtr       = &multiMB;
 131   if (doMIinit && !doMIMB) return false;
 132   if (doMIinit && doDiffraction && (!doMISDA || !doMISDB)) return false;
 133   if (!doMIMB || !doMISDA || !doMISDB) doMI = false;
 134   return true;
 135 }
 136
 137 //--------------------------------------------------------------------------
 138
 139 // Main routine to do the parton-level evolution.
 140
 141 bool PartonLevel::next( Event& process, Event& event) {
 142
 143   // Current event classification.
 144   isResolved     = infoPtr->isResolved();
 145   isResolvedA    = isResolved;
 146   isResolvedB    = isResolved;
 147   isDiffA        = infoPtr->isDiffractiveA();
 148   isDiffB        = infoPtr->isDiffractiveB();
 149   isDiff         = isDiffA || isDiffB;
 150   isDoubleDiff   = isDiffA && isDiffB;
 151   isSingleDiff   = isDiff && !isDoubleDiff;
 152   isMinBias      = infoPtr->isMinBias();
 153
 154   // nHardLoop counts how many hard-scattering subsystems are to be processed.
 155   // Almost always 1, but elastic and low-mass diffraction gives 0, while
 156   // double diffraction can give up to 2. Not to be confused with SecondHard.
 157   int nHardLoop  = 1;
 158   if (!isResolved) nHardLoop = (isDiff) ? decideResolvedDiff( process) : 0;
 159
 160   // Handle unresolved subsystems. Done if no resolved ones.
 161   sizeProcess    = 0;
 162   sizeEvent      = 0;
 163   if (!isResolvedA || !isResolvedB) {
 164     bool physical = setupUnresolvedSys( process, event);
 165     if (!physical || nHardLoop == 0) return physical;
 166     sizeProcess  = process.size();
 167     sizeEvent    = event.size();
 168   }
 169
 170   // Big outer loop to handle up to two systems (in double diffraction),
 171   // but normally one. (Not indented in following, but end clearly marked.)
 172   for (int iHardLoop = 1; iHardLoop <= nHardLoop; ++iHardLoop) {
 173     infoPtr->setCounter(20, iHardLoop);
 174     infoPtr->setCounter(21);
 175
 176   // Process and event records can be out of step for diffraction.
 177   if (iHardLoop == 2) {
 178     sizeProcess = process.size();
 179     sizeEvent   = event.size();
 180     partonSystemsPtr->clear();
 181   }
 182
 183   // If you need to restore then do not throw existing diffractive system.
 184   if (isDiff) {
 185     event.saveSize();
 186     event.saveJunctionSize();
 187   }
 188
 189   // Allow special treatment of diffractive systems.
 190   if (isDiff) setupResolvedDiff(iHardLoop, process);
 191
 192   // Prepare to do multiple interactions; at new mass for diffraction.
 193   if (doMIinit) multiPtr->reset();
 194
 195   // Special case if minimum bias: do hardest interaction.
 196   if (isMinBias || isDiff) {
 197     multiPtr->pTfirst();
 198     multiPtr->setupFirstSys( process);
 199   }
 200
 201   // Allow up to ten tries; failure possible for beam remnants.
 202   // Main cause: inconsistent colour flow at the end of the day.
 203   bool physical = true;
 204   int  nRad     = 0;
 205   for (int iTry = 0; iTry < NTRY; ++ iTry) {
 206     infoPtr->addCounter(21);
 207     for (int i = 22; i < 32; ++i) infoPtr->setCounter(i);
 208
 209     // Reset flag, counters and max scales.
 210     physical   = true;
 211     nMI        = (doSecondHard) ? 2 : 1;
 212     nISR       = 0;
 213     nFSRinProc = 0;
 214     nFSRinRes  = 0;
 215     nISRhard   = 0;
 216     nFSRhard   = 0;
 217     pTsaveMI   = 0.;
 218     pTsaveISR  = 0.;
 219     pTsaveFSR  = 0.;
 220
 221     // Identify hard interaction system for showers.
 222     setupHardSys( iHardLoop, process, event);
 223
 224     // Optionally check for a veto after the hardest interaction.
 225     if (canVetoMIStep) {
 226       doVeto = userHooksPtr->doVetoMIStep( 1, event);
 227       // Abort event if vetoed.
 228       if (doVeto) {
 229         if (isDiff) leaveResolvedDiff( iHardLoop, event);
 230         return false;
 231       }
 232     }
 233
 234     // Check matching of process scale to maximum ISR/FSR/MI scales.
 235     double Q2Fac       = infoPtr->Q2Fac();
 236     double Q2Ren       = infoPtr->Q2Ren();
 237     bool limitPTmaxISR = (doISR)
 238       ? spacePtr->limitPTmax( event, Q2Fac, Q2Ren) : false;
 239     bool limitPTmaxFSR = (doFSRduringProcess)
 240       ? timesPtr->limitPTmax( event, Q2Fac, Q2Ren) : false;
 241     bool limitPTmaxMI  = (doMI)  ? multiPtr->limitPTmax( event) : false;
 242
 243     // Set hard scale, maximum for showers and multiple interactions,
 244     double pTscale  = process.scale();
 245     if (doSecondHard) pTscale = max( pTscale, process.scaleSecond() );
 246     double pTmaxMI  = (limitPTmaxMI)  ? pTscale : infoPtr->eCM();
 247     double pTmaxISR = (limitPTmaxISR) ? spacePtr->enhancePTmax() * pTscale
 248                                       : infoPtr->eCM();
 249     double pTmaxFSR = (limitPTmaxFSR) ? timesPtr->enhancePTmax() * pTscale
 250                                       : infoPtr->eCM();
 251     double pTmax    = max( pTmaxMI, max( pTmaxISR, pTmaxFSR) );
 252     pTsaveMI        = pTmaxMI;
 253     pTsaveISR       = pTmaxISR;
 254     pTsaveFSR       = pTmaxFSR;
 255
 256     // Prepare the classes to begin the generation.
 257     if (doMI)  multiPtr->prepare( pTmaxMI);
 258     if (doISR) spacePtr->prepare( 0, event, limitPTmaxISR);
 259     if (doFSRduringProcess) timesPtr->prepare( 0, event, limitPTmaxFSR);
 260     if (doSecondHard && doISR) spacePtr->prepare( 1, event, limitPTmaxISR);
 261     if (doSecondHard && doFSRduringProcess)
 262       timesPtr->prepare( 1, event, limitPTmaxFSR);
 263
 264     // Set up initial veto scale.
 265     doVeto        = false;
 266     double pTveto = pTvetoPT;
 267     typeLatest    = 0;
 268
 269     // Begin evolution down in pT from hard pT scale.
 270     do {
 271       infoPtr->addCounter(22);
 272       typeVetoStep = 0;
 273       nRad         =  nISR + nFSRinProc;
 274
 275       // Find next pT value for FSR, MI and ISR.
 276       // Order calls to minimize time expenditure.
 277       double pTgen = 0.;
 278       double pTtimes = (doFSRduringProcess)
 279         ? timesPtr->pTnext( event, pTmaxFSR, pTgen) : -1.;
 280       pTgen = max( pTgen, pTtimes);
 281       double pTmulti = (doMI)
 282         ? multiPtr->pTnext( pTmaxMI, pTgen, event) : -1.;
 283       pTgen = max( pTgen, pTmulti);
 284       double pTspace = (doISR)
 285         ? spacePtr->pTnext( event, pTmaxISR, pTgen, nRad) : -1.;
 286       double pTnow = max( pTtimes, max( pTmulti, pTspace));
 287       infoPtr->setPTnow( pTnow);
 288
 289       // Allow a user veto. Only do it once, so remember to change pTveto.
 290       if (pTveto > 0. && pTveto > pTnow) {
 291         pTveto = -1.;
 292         doVeto = userHooksPtr->doVetoPT( typeLatest, event);
 293         // Abort event if vetoed.
 294         if (doVeto) {
 295           if (isDiff) leaveResolvedDiff( iHardLoop, event);
 296           return false;
 297         }
 298       }
 299
 300       // Do a multiple interaction (if allowed).
 301       if (pTmulti > 0. && pTmulti > pTspace && pTmulti > pTtimes) {
 302         infoPtr->addCounter(23);
 303         multiPtr->scatter( event);
 304         typeLatest = 1;
 305         ++nMI;
 306         if (canVetoMIStep && nMI <= nVetoMIStep) typeVetoStep = 1;
 307
 308         // Update ISR and FSR dipoles.
 309         if (doISR)              spacePtr->prepare( nMI - 1, event);
 310         if (doFSRduringProcess) timesPtr->prepare( nMI - 1, event);
 311
 312         // Set maximal scales for next pT to pick.
 313         pTmaxMI  = pTmulti;
 314         pTmaxISR = min( pTmulti, pTmaxISR);
 315         pTmaxFSR = min( pTmulti, pTmaxFSR);
 316         pTmax    = pTmulti;
 317       }
 318
 319       // Do an initial-state emission (if allowed).
 320       else if (pTspace > 0. && pTspace > pTtimes) {
 321         infoPtr->addCounter(24);
 322         if (spacePtr->branch( event)) {
 323           typeLatest = 2;
 324           iSysNow = spacePtr->system();
 325           ++nISR;
 326           if (iSysNow == 0) ++nISRhard;
 327           if (canVetoStep && iSysNow == 0 && nISRhard <= nVetoStep)
 328             typeVetoStep = 2;
 329
 330           // Update FSR dipoles.
 331           if (doFSRduringProcess) timesPtr->update( iSysNow, event);
 332
 333         // Rescatter: it is possible for kinematics to fail, in which
 334         //            case we need to restart the parton level processing.
 335         } else if (spacePtr->doRestart()) {
 336           physical = false;
 337           break;
 338         }
 339
 340         // Set maximal scales for next pT to pick.
 341         pTmaxMI  = min( pTspace, pTmaxMI);
 342         pTmaxISR = pTspace;
 343         pTmaxFSR = min( pTspace, pTmaxFSR);
 344         pTmax    = pTspace;
 345       }
 346
 347       // Do a final-state emission (if allowed).
 348       else if (pTtimes > 0.) {
 349         infoPtr->addCounter(25);
 350         if (timesPtr->branch( event, true)) {
 351           typeLatest = 3;
 352           iSysNow = timesPtr->system();
 353           ++nFSRinProc;
 354           if (iSysNow == 0) ++nFSRhard;
 355           if (canVetoStep && iSysNow == 0 && nFSRhard <= nVetoStep)
 356             typeVetoStep = 3;
 357
 358           // Update ISR dipoles.
 359           if (doISR) spacePtr->update( iSysNow, event);
 360         }
 361
 362         // Set maximal scales for next pT to pick.
 363         pTmaxMI  = min( pTtimes, pTmaxMI);
 364         pTmaxISR = min( pTtimes, pTmaxISR);
 365         pTmaxFSR = pTtimes;
 366         pTmax    = pTtimes;
 367       }
 368
 369       // If no pT scales above zero then nothing to be done.
 370       else pTmax = 0.;
 371
 372       // Optionally check for a veto after the first few interactions,
 373       // or after the first few emissions, ISR or FSR, in the hardest system.
 374       if (typeVetoStep == 1) {
 375         doVeto = userHooksPtr->doVetoMIStep( nMI, event);
 376       } else if (typeVetoStep > 1) {
 377         doVeto = userHooksPtr->doVetoStep( typeVetoStep, nISRhard,
 378           nFSRhard, event);
 379       }
 380
 381       // Abort event if vetoed.
 382       if (doVeto) {
 383         if (isDiff) leaveResolvedDiff( iHardLoop, event);
 384         return false;
 385       }
 386
 387       // Keep on evolving until nothing is left to be done.
 388       if (typeLatest > 0 && typeLatest < 4)
 389         infoPtr->addCounter(25 + typeLatest);
 390       infoPtr->setPartEvolved( nMI, nISR);
 391     } while (pTmax > 0.);
 392
 393     // Do all final-state emissions if not already considered above.
 394     if (doFSRafterProcess) {
 395
 396       // Find largest scale for final partons.
 397       pTmax = 0.;
 398       for (int i = 0; i < event.size(); ++i)
 399         if (event[i].isFinal() && event[i].scale() > pTmax)
 400           pTmax = event[i].scale();
 401       pTsaveFSR = pTmax;
 402
 403       // Prepare all subsystems for evolution.
 404       for (int iSys = 0; iSys < partonSystemsPtr->sizeSys(); ++iSys)
 405         timesPtr->prepare( iSys, event);
 406
 407       // Set up initial veto scale.
 408       doVeto        = false;
 409       pTveto = pTvetoPT;
 410
 411       // Begin evolution down in pT from hard pT scale.
 412       do {
 413         infoPtr->addCounter(29);
 414         typeVetoStep = 0;
 415         double pTtimes = timesPtr->pTnext( event, pTmax, 0.);
 416         infoPtr->setPTnow( pTtimes);
 417
 418         // Allow a user veto. Only do it once, so remember to change pTveto.
 419         if (pTveto > 0. && pTveto > pTtimes) {
 420           pTveto = -1.;
 421           doVeto = userHooksPtr->doVetoPT( 4, event);
 422           // Abort event if vetoed.
 423           if (doVeto) {
 424             if (isDiff) leaveResolvedDiff( iHardLoop, event);
 425             return false;
 426           }
 427         }
 428
 429         // Do a final-state emission (if allowed).
 430         if (pTtimes > 0.) {
 431           infoPtr->addCounter(30);
 432           if (timesPtr->branch( event, true)) {
 433             iSysNow = timesPtr->system();
 434             ++nFSRinProc;
 435             if (iSysNow == 0) ++nFSRhard;
 436             if (canVetoStep && iSysNow == 0 && nFSRhard <= nVetoStep)
 437             typeVetoStep = 4;
 438           }
 439           pTmax = pTtimes;
 440         }
 441
 442         // If no pT scales above zero then nothing to be done.
 443         else pTmax = 0.;
 444
 445         // Optionally check for a veto after the first few emissions.
 446         if (typeVetoStep > 0) {
 447           doVeto = userHooksPtr->doVetoStep( typeVetoStep, nISRhard,
 448             nFSRhard, event);
 449           // Abort event if vetoed.
 450           if (doVeto) {
 451             if (isDiff) leaveResolvedDiff( iHardLoop, event);
 452             return false;
 453           }
 454         }
 455
 456       // Keep on evolving until nothing is left to be done.
 457         infoPtr->addCounter(31);
 458       } while (pTmax > 0.);
 459     }
 460
 461     // Add beam remnants, including primordial kT kick and colour tracing.
 462     if (physical && doRemnants && !remnants.add( event)) physical = false;
 463
 464     // If no problems then done.
 465     if (physical) break;
 466
 467     // Else restore and loop, but do not throw existing diffractive system.
 468     if (!isDiff) event.clear();
 469     else {
 470       event.restoreSize();
 471       event.restoreJunctionSize();
 472     }
 473     beamAPtr->clear();
 474     beamBPtr->clear();
 475     partonSystemsPtr->clear();
 476
 477   // End loop over ten tries. Restore from diffraction. Hopefully it worked.
 478   }
 479   if (isDiff) leaveResolvedDiff( iHardLoop, event);
 480   if (!physical) return false;
 481
 482   // End big outer loop to handle two systems in double diffraction.
 483   }
 484
 485   // Perform showers in resonance decay chains.
 486   doVeto = !resonanceShowers( process, event);
 487   // Abort event if vetoed.
 488   if (doVeto) return false;
 489
 490   // Store event properties.
 491   infoPtr->setImpact( multiPtr->bMI(), multiPtr->enhanceMI());
 492   infoPtr->setEvolution( pTsaveMI, pTsaveISR, pTsaveFSR,
 493     nMI, nISR, nFSRinProc, nFSRinRes);
 494
 495   // Done.
 496   return true;
 497 }
 498
 499 //--------------------------------------------------------------------------
 500
 501 // Decide which diffractive subsystems are resolved (= perturbative).
 502
 503 int PartonLevel::decideResolvedDiff( Event& process) {
 504
 505   // Loop over two systems.
 506   int nHighMass = 0;
 507   for (int iDiffSys = 1; iDiffSys <= 2; ++iDiffSys) {
 508     int iDiffMot = iDiffSys + 2;
 509
 510     // Only high-mass diffractive systems should be resolved.
 511     double mDiff = process[iDiffMot].m();
 512     bool isHighMass = ( mDiff > mMinDiff
 513       && rndmPtr->flat() < 1. - exp( -(mDiff - mMinDiff) / mWidthDiff ) );
 514
 515     // Set outcome and done.
 516     if (isHighMass) ++nHighMass;
 517     if (iDiffSys == 1) isResolvedA = isHighMass;
 518     if (iDiffSys == 2) isResolvedB = isHighMass;
 519   }
 520   return nHighMass;
 521
 522 }
 523
 524 //--------------------------------------------------------------------------
 525
 526 // Set up an unresolved process, i.e. elastic or diffractive.
 527
 528 bool PartonLevel::setupUnresolvedSys( Event& process, Event& event) {
 529
 530   // No hard scale in event.
 531   process.scale( 0.);
 532
 533   // Copy particles from process to event.
 534   for (int i = 0; i < process.size(); ++ i) event.append( process[i]);
 535
 536   // Loop to find diffractively excited beams.
 537   for (int i = 0; i < 2; ++i)
 538   if ( (i == 0 && isDiffA && !isResolvedA)
 539     || (i == 1 && isDiffB && !isResolvedB) ) {
 540     int iBeam = i + 3;
 541     BeamParticle* beamPtr = (i == 0) ? beamAPtr : beamBPtr;
 542
 543     // Diffractive mass. Reconstruct boost and rotation to event cm frame.
 544     double mDiff  = process[iBeam].m();
 545     double m2Diff = mDiff*mDiff;
 546     double beta   = process[iBeam].pAbs() / process[iBeam].e();
 547     double theta  = process[iBeam].theta();
 548     double phi    = process[iBeam].phi();
 549
 550     // Pick quark or gluon kicked out and flavour subdivision.
 551     bool gluonIsKicked = beamPtr->pickGluon(mDiff);
 552     int id1 = beamPtr->pickValence();
 553     int id2 = beamPtr->pickRemnant();
 554
 555     // Find flavour masses. Scale them down if too big.
 556     double m1 = particleDataPtr->constituentMass(id1);
 557     double m2 = particleDataPtr->constituentMass(id2);
 558     if (m1 + m2 > 0.5 * mDiff) {
 559       double reduce = 0.5 * mDiff / (m1 + m2);
 560       m1 *= reduce;
 561       m2 *= reduce;
 562     }
 563
 564     // If quark is kicked out, then trivial kinematics in rest frame.
 565     if (!gluonIsKicked) {
 566       double pAbs = sqrt( pow2(m2Diff - m1*m1 - m2*m2)
 567         - pow2(2. * m1 * m2) ) / (2. * mDiff);
 568       double e1 = (m2Diff + m1*m1 - m2*m2) / (2. * mDiff);
 569       double e2 = (m2Diff + m2*m2 - m1*m1) / (2. * mDiff);
 570       Vec4 p1(0.,0., -pAbs, e1);
 571       Vec4 p2(0.,0., pAbs, e2);
 572
 573       // Boost and rotate to event cm frame.
 574       p1.bst(0., 0., beta); p1.rot(theta, phi);
 575       p2.bst(0., 0., beta); p2.rot(theta, phi);
 576
 577       // Set colours.
 578       int col1, acol1, col2, acol2;
 579       if (particleDataPtr->colType(id1) == 1) {
 580         col1 = event.nextColTag(); acol1 = 0;
 581         col2 = 0; acol2 = col1;
 582       } else {
 583         col1 = 0; acol1 = event.nextColTag();
 584         col2 = acol1; acol2 = 0;
 585       }
 586      // Update process colours to stay in step.
 587       process.nextColTag();
 588
 589       // Store partons of diffractive system and mark system decayed.
 590       int iDauBeg = event.append( id1, 23, iBeam, 0, 0, 0, col1, acol1,
 591         p1, m1);
 592       int iDauEnd = event.append( id2, 63, iBeam, 0, 0, 0, col2, acol2,
 593         p2, m2);
 594       event[iBeam].statusNeg();
 595       event[iBeam].daughters(iDauBeg, iDauEnd);
 596
 597
 598     // If gluon is kicked out: share momentum between two remnants.
 599     } else {
 600       double m2Sys, zSys, pxSys, pySys, mTS1, mTS2;
 601       zSys = beamPtr->zShare(mDiff, m1, m2);
 602
 603       // Provide relative pT kick in remnant. Construct (transverse) masses.
 604       pxSys = beamPtr->pxShare();
 605       pySys = beamPtr->pyShare();
 606       mTS1  = m1*m1 + pxSys*pxSys + pySys*pySys;
 607       mTS2  = m2*m2 + pxSys*pxSys + pySys*pySys;
 608       m2Sys = mTS1 / zSys + mTS2 / (1. - zSys);
 609
 610       // Momentum of kicked-out massless gluon in diffractive rest frame.
 611       double pAbs = (m2Diff - m2Sys) / (2. * mDiff);
 612       Vec4 pG(0., 0., -pAbs, pAbs);
 613       Vec4 pRem(0., 0., pAbs, mDiff - pAbs);
 614
 615       // Momenta of the two beam remnant flavours. (Lightcone p+ = m_diff!)
 616       double e1 = 0.5 * (zSys * mDiff + mTS1 / (zSys * mDiff));
 617       double pL1 = 0.5 * (zSys * mDiff - mTS1 / (zSys * mDiff));
 618       Vec4 p1(pxSys, pySys, pL1, e1);
 619       Vec4 p2 = pRem - p1;
 620
 621       // Boost and rotate to event cm frame.
 622       pG.bst(0., 0., beta); pG.rot(theta, phi);
 623       p1.bst(0., 0., beta); p1.rot(theta, phi);
 624       p2.bst(0., 0., beta); p2.rot(theta, phi);
 625
 626       // Set colours.
 627       int colG, acolG, col1, acol1, col2, acol2;
 628       if (particleDataPtr->colType(id1) == 1) {
 629         col1 = event.nextColTag(); acol1 = 0;
 630         colG = event.nextColTag(); acolG = col1;
 631         col2 = 0; acol2 = colG;
 632       } else {
 633         col1 = 0; acol1 = event.nextColTag();
 634         colG = acol1; acolG = event.nextColTag();
 635         col2 = acolG; acol2 = 0;
 636       }
 637       // Update process colours to stay in step.
 638       process.nextColTag();
 639       process.nextColTag();
 640
 641       // Store partons of diffractive system and mark system decayed.
 642       int iDauBeg = event.append( 21, 23, iBeam, 0, 0, 0, colG, acolG,
 643         pG, m1);
 644       event.append( id1, 63, iBeam, 0, 0, 0, col1, acol1, p1, m1);
 645       int iDauEnd = event.append( id2, 63, iBeam, 0, 0, 0, col2, acol2,
 646         p2, m2);
 647       event[iBeam].statusNeg();
 648       event[iBeam].daughters(iDauBeg, iDauEnd);
 649     }
 650
 651   // End loop over beams. Done.
 652   }
 653   return true;
 654
 655 }
 656
 657 //--------------------------------------------------------------------------
 658
 659 // Set up the hard process(es), excluding subsequent resonance decays.
 660
 661   void PartonLevel::setupHardSys( int iHardLoop, Event& process,
 662     Event& event) {
 663
 664   // Incoming partons to hard process are stored in slots 3 and 4.
 665   int inS = 0;
 666   int inP = 3;
 667   int inM = 4;
 668
 669   // Identify any diffractive system, mother, last entry. Offset.
 670   int iDiffSys =  (iHardLoop == 2 || !isResolvedA) ? 2 : 1;
 671   int iDiffMot = iDiffSys + 2;
 672   int iDiffDau = process.size() - 1;
 673   int nOffset  = sizeEvent - sizeProcess;
 674
 675   // Resolved diffraction means more entries.
 676    if (isDiff) {
 677     inS   = iDiffMot;
 678     inP   = iDiffDau - 3;
 679     inM   = iDiffDau - 2;
 680
 681     // Diffractively excited particle described as Pomeron-hadron beams.
 682     event[inS].statusNeg();
 683     event[inS].daughters( inP - 2 + nOffset, inM - 2 + nOffset);
 684   }
 685
 686   // If two hard interactions then find where second begins.
 687   int iBeginSecond = process.size();
 688   if (doSecondHard) {
 689     iBeginSecond = 5;
 690     while (process[iBeginSecond].status() != -21) ++iBeginSecond;
 691   }
 692
 693   // If incoming partons are massive then recalculate to put them massless.
 694   if (process[inP].m() != 0. || process[inM].m() != 0.) {
 695     double pPos = process[inP].pPos() + process[inM].pPos();
 696     double pNeg = process[inP].pNeg() + process[inM].pNeg();
 697     process[inP].pz( 0.5 * pPos);
 698     process[inP].e(  0.5 * pPos);
 699     process[inP].m(  0.);
 700     process[inM].pz(-0.5 * pNeg);
 701     process[inM].e(  0.5 * pNeg);
 702     process[inM].m(  0.);
 703   }
 704
 705   // Add incoming hard-scattering partons to list in beam remnants.
 706   double x1 = process[inP].pPos() / process[inS].m();
 707   beamAPtr->append( inP + nOffset, process[inP].id(), x1);
 708   double x2 = process[inM].pNeg() / process[inS].m();
 709   beamBPtr->append( inM + nOffset, process[inM].id(), x2);
 710
 711   // Scale. Find whether incoming partons are valence or sea. Store.
 712   double scale = process.scale();
 713   beamAPtr->xfISR( 0, process[inP].id(), x1, scale*scale);
 714   int vsc1 = beamAPtr->pickValSeaComp();
 715   beamBPtr->xfISR( 0, process[inM].id(), x2, scale*scale);
 716   int vsc2 = beamBPtr->pickValSeaComp();
 717   bool isVal1 = (vsc1 == -3);
 718   bool isVal2 = (vsc2 == -3);
 719   infoPtr->setValence( isVal1, isVal2);
 720
 721   // Initialize info needed for subsequent sequential decays + showers.
 722   nHardDone = sizeProcess;
 723   iPosBefShow.resize( process.size() );
 724
 725   // Add the beam and hard subprocess partons to the event record.
 726   for (int i = sizeProcess; i < iBeginSecond; ++ i) {
 727     if (process[i].mother1() > inM) break;
 728     int j = event.append(process[i]);
 729     iPosBefShow[i] = i;
 730
 731     // Offset history if process and event not in step.
 732     if (nOffset != 0) {
 733       int iOrd = i - iBeginSecond + 7;
 734       if (iOrd == 1 || iOrd == 2)
 735         event[j].offsetHistory( 0, 0, 0, nOffset);
 736       else if (iOrd == 3 || iOrd == 4)
 737         event[j].offsetHistory( 0, nOffset, 0, nOffset);
 738       else if (iOrd == 5 || iOrd == 6)
 739         event[j].offsetHistory( 0, nOffset, 0, 0);
 740     }
 741
 742     // Currently outgoing ones should not count as decayed.
 743     if (event[j].status() == -22) {
 744       event[j].statusPos();
 745       event[j].daughters(0, 0);
 746     }
 747
 748     // Complete task of copying hard subsystem into event record.
 749     ++nHardDone;
 750   }
 751
 752   // Store participating partons as first set in list of all systems.
 753   partonSystemsPtr->addSys();
 754   partonSystemsPtr->setInA(0, inP + nOffset);
 755   partonSystemsPtr->setInB(0, inM + nOffset);
 756   for (int i = inM + 1; i < nHardDone; ++i)
 757     partonSystemsPtr->addOut(0, i + nOffset);
 758   partonSystemsPtr->setSHat( 0,
 759     (event[inP + nOffset].p() + event[inM + nOffset].p()).m2Calc() );
 760
 761   // Identify second hard process where applicable.
 762   // Since internally generated incoming partons are guaranteed massless.
 763   if (doSecondHard) {
 764     int inP2 = iBeginSecond;
 765     int inM2 = iBeginSecond + 1;
 766
 767     // Add incoming hard-scattering partons to list in beam remnants.
 768     // Not valid if not in rest frame??
 769     x1 = process[inP2].pPos() / process[0].e();
 770     beamAPtr->append( inP2, process[inP2].id(), x1);
 771     x2 = process[inM2].pNeg() / process[0].e();
 772     beamBPtr->append( inM2, process[inM2].id(), x2);
 773
 774     // Find whether incoming partons are valence or sea.
 775     scale = process.scaleSecond();
 776     beamAPtr->xfISR( 1, process[inP2].id(), x1, scale*scale);
 777     beamAPtr->pickValSeaComp();
 778     beamBPtr->xfISR( 1, process[inM2].id(), x2, scale*scale);
 779     beamBPtr->pickValSeaComp();
 780
 781     // Add the beam and hard subprocess partons to the event record.
 782     for (int i = inP2; i < process.size(); ++ i) {
 783       int mother = process[i].mother1();
 784       if ( (mother > 2 && mother < inP2) || mother > inM2 ) break;
 785       event.append(process[i]);
 786       iPosBefShow[i] = i;
 787
 788       // Currently outgoing ones should not count as decayed.
 789       if (event[i].status() == -22) {
 790         event[i].statusPos();
 791         event[i].daughters(0, 0);
 792       }
 793
 794       // Complete task of copying hard subsystem into event record.
 795       ++nHardDone;
 796     }
 797
 798     // Store participating partons as second set in list of all systems.
 799     partonSystemsPtr->addSys();
 800     partonSystemsPtr->setInA(1, inP2);
 801     partonSystemsPtr->setInB(1, inM2);
 802     for (int i = inM2 + 1; i < nHardDone; ++i)
 803       partonSystemsPtr->addOut(1, i);
 804     partonSystemsPtr->setSHat( 1,
 805       (event[inP2].p() + event[inM2].p()).m2Calc() );
 806
 807   // End code for second hard process.
 808   }
 809
 810   // Update event colour tag to maximum in whole process.
 811   int maxColTag = 0;
 812   for (int i = 0; i < process.size(); ++ i) {
 813     if (process[i].col() > maxColTag) maxColTag = process[i].col();
 814     if (process[i].acol() > maxColTag) maxColTag = process[i].acol();
 815   }
 816   event.initColTag(maxColTag);
 817
 818   // Copy junctions from process to event.
 819   for (int i = 0; i < process.sizeJunction(); ++i)
 820     event.appendJunction( process.getJunction(i));
 821
 822   // Done.
 823 }
 824
 825 //--------------------------------------------------------------------------
 826
 827 // Resolved diffraction: replace full event with diffractive subsystem.
 828
 829 void PartonLevel::setupResolvedDiff(int iHardLoop, Event& process) {
 830
 831   // Identify diffractive system, mother, last entry.
 832   int iDiffSys =  (iHardLoop == 2 || !isResolvedA) ? 2 : 1;
 833   int iDiffMot = iDiffSys + 2;
 834   int iDiffDau = process.size() - 1;
 835
 836   // Diffractive system mass.
 837   double mDiff = process[iDiffMot].m();
 838   double m2Diff = mDiff * mDiff;
 839
 840   // Diffractively excited particle described as Pomeron-hadron beams.
 841   process[iDiffMot].statusNeg();
 842   process[iDiffMot].daughters( iDiffDau + 1, iDiffDau + 2);
 843
 844   // Set up Pomeron-proton system as if it were the complete collision.
 845   int    idHad = process[iDiffSys].id();
 846   double mHad  = process[iDiffSys].m();
 847   double m2Had = mHad * mHad;
 848   double m2Pom = (process[2].p() - process[4].p()).m2Calc();
 849   double mPom  = (m2Pom >= 0.) ? sqrt(m2Pom) : -sqrt(-m2Pom);
 850   double eHad  = 0.5 * (m2Diff + m2Had - m2Pom) / mDiff;
 851   double ePom  = 0.5 * (m2Diff + m2Pom - m2Had) / mDiff;
 852   double pzHP  = 0.5 * sqrtpos( pow2(m2Diff - m2Had - m2Pom)
 853                - 4. * m2Had * m2Pom ) / mDiff;
 854   process.append(   990, 13, iDiffMot, 0, 0, 0, 0, 0,
 855     0., 0.,  pzHP, ePom, mPom);
 856   process.append( idHad, 13, iDiffMot, 0, 0, 0, 0, 0,
 857     0., 0., -pzHP, eHad, mHad);
 858
 859   // Reassign multiple interactions pointer to right object.
 860   multiPtr = (iDiffSys == 1) ? &multiSDA : &multiSDB;
 861
 862   // Reassign one beam pointer to refer to incoming Pomeron.
 863   if (iDiffSys == 1) {
 864     beamAPtr = beamPomBPtr;
 865     beamBPtr = beamHadAPtr;
 866   } else {
 867     beamAPtr = beamPomAPtr;
 868     beamBPtr = beamHadBPtr;
 869   }
 870
 871   // Beams not found in normal slots 1 and 2.
 872   int beamOffset = (sizeEvent > 0) ? sizeEvent - 1 : 4;
 873
 874   // Reassign beam pointers in other classes.
 875   timesPtr->reassignBeamPtrs( beamAPtr, beamBPtr, beamOffset);
 876   spacePtr->reassignBeamPtrs( beamAPtr, beamBPtr);
 877   remnants.reassignBeamPtrs( beamAPtr, beamBPtr);
 878
 879   // Pretend that the diffractive system is the whole collision.
 880   eCMsave = infoPtr->eCM();
 881   infoPtr->setECM( mDiff);
 882   beamAPtr->newPzE(  pzHP, ePom);
 883   beamBPtr->newPzE( -pzHP, eHad);
 884
 885 }
 886
 887 //--------------------------------------------------------------------------
 888
 889 // Resolved diffraction: restore to original behaviour.
 890
 891 void PartonLevel::leaveResolvedDiff( int iHardLoop, Event& event) {
 892
 893   // Identify diffractive system.
 894   int iDiffSys =  (iHardLoop == 2 || !isResolvedA) ? 2 : 1;
 895
 896   // Reconstruct boost and rotation to event cm frame.
 897   Vec4 pPom = event[3 - iDiffSys].p() - event[5 - iDiffSys].p();
 898   Vec4 pHad =  event[iDiffSys].p();
 899   RotBstMatrix MtoCM;
 900   MtoCM.fromCMframe( pPom, pHad);
 901
 902   // Perform rotation and boost on diffractive system.
 903   int iFirst = (iHardLoop == 1) ? 5 + sizeEvent - sizeProcess : sizeEvent;
 904   for (int i = iFirst; i < event.size(); ++i)
 905     event[i].rotbst( MtoCM);
 906
 907   // Restore multiple interactions pointer to default object.
 908   multiPtr = &multiMB;
 909
 910   // Restore beam pointers to incoming hadrons.
 911   beamAPtr = beamHadAPtr;
 912   beamBPtr = beamHadBPtr;
 913
 914   // Reassign beam pointers in other classes.
 915   timesPtr->reassignBeamPtrs( beamAPtr, beamBPtr, 0);
 916   spacePtr->reassignBeamPtrs( beamAPtr, beamBPtr);
 917   remnants.reassignBeamPtrs( beamAPtr, beamBPtr);
 918
 919   // Restore cm energy.
 920   infoPtr->setECM( eCMsave);
 921   beamAPtr->newPzE( event[1].pz(), event[1].e());
 922   beamBPtr->newPzE( event[2].pz(), event[2].e());
 923
 924 }
 925
 926 //--------------------------------------------------------------------------
 927
 928 // Handle showers in successive resonance decays.
 929
 930 bool PartonLevel::resonanceShowers( Event& process, Event& event) {
 931
 932   // Isolate next system to be processed, if anything remains.
 933   int nRes    = 0;
 934   int nFSRres = 0;
 935   while (nHardDone < process.size()) {
 936     ++nRes;
 937     int iBegin = nHardDone;
 938
 939     // Mother in hard process and in complete event (after shower).
 940     int iHardMother      = process[iBegin].mother1();
 941     Particle& hardMother = process[iHardMother];
 942     int iBefMother       = iPosBefShow[iHardMother];
 943     int iAftMother       = event.iBotCopyId(iBefMother);
 944     Particle& aftMother  = event[iAftMother];
 945
 946     // From now on mother counts as decayed.
 947     aftMother.statusNeg();
 948
 949     // Mother can have been moved by showering (in any of previous steps),
 950     // so prepare to update colour and momentum information for system.
 951     int colBef  = hardMother.col();
 952     int acolBef = hardMother.acol();
 953     int colAft  = aftMother.col();
 954     int acolAft = aftMother.acol();
 955     RotBstMatrix M;
 956     M.bst( hardMother.p(), aftMother.p());
 957
 958     // Extract next partons from hard event into normal event record.
 959     for (int i = iBegin; i < process.size(); ++ i) {
 960       if (process[i].mother1() != iHardMother) break;
 961       int iNow = event.append( process[i] );
 962       iPosBefShow[i] = iNow;
 963       Particle& now = event.back();
 964
 965       // Currently outgoing ones should not count as decayed.
 966       if (now.status() == -22) {
 967         now.statusPos();
 968         now.daughters(0, 0);
 969       }
 970
 971       // Update daughter and mother information.
 972       if (i == iBegin) event[iAftMother].daughter1( iNow);
 973       else             event[iAftMother].daughter2( iNow);
 974       now.mother1(iAftMother);
 975
 976       // Update colour and momentum information.
 977       if (now.col() == colBef) now.col( colAft);
 978       if (now.acol() == acolBef) now.acol( acolAft);
 979       now.rotbst( M);
 980
 981       // Complete task of copying next subsystem into event record.
 982       ++nHardDone;
 983     }
 984     int iEnd = nHardDone - 1;
 985
 986     // Do parton showers inside subsystem: maximum scale by mother mass.
 987     if (doFSRinResonances) {
 988       double pTmax = 0.5 * hardMother.m();
 989       if (canSetScale) pTmax = userHooksPtr->scaleResonance( iAftMother, event);
 990       nFSRhard     = 0;
 991
 992       // Add new system, automatically with two empty beam slots.
 993       int iSys = partonSystemsPtr->addSys();
 994       partonSystemsPtr->setSHat(iSys, pow2(hardMother.m()) );
 995
 996       // Loop over allowed range to find all final-state particles.
 997       for (int i = iPosBefShow[iBegin]; i <= iPosBefShow[iEnd]; ++i)
 998       if (event[i].isFinal()) partonSystemsPtr->addOut( iSys, i);
 999
1000       // Let prepare routine do the setup.
1001       timesDecPtr->prepare( iSys, event);
1002
1003       // Set up initial veto scale.
1004       doVeto        = false;
1005       double pTveto = pTvetoPT;
1006
1007       // Begin evolution down in pT from hard pT scale.
1008       do {
1009         typeVetoStep = 0;
1010         double pTtimes = timesDecPtr->pTnext( event, pTmax, 0.);
1011
1012         // Allow a user veto. Only do it once, so remember to change pTveto.
1013         if (pTveto > 0. && pTveto > pTtimes) {
1014           pTveto = -1.;
1015           doVeto = userHooksPtr->doVetoPT( 5, event);
1016           // Abort event if vetoed.
1017           if (doVeto) return false;
1018         }
1019
1020         // Do a final-state emission (if allowed).
1021         if (pTtimes > 0.) {
1022           if (timesDecPtr->branch( event)) {
1023             ++nFSRres;
1024             ++nFSRhard;
1025             if (canVetoStep && nFSRhard <= nVetoStep) typeVetoStep = 5;
1026           }
1027           pTmax = pTtimes;
1028         }
1029
1030         // If no pT scales above zero then nothing to be done.
1031         else pTmax = 0.;
1032
1033         // Optionally check for a veto after the first few emissions.
1034         if (typeVetoStep > 0) {
1035           doVeto = userHooksPtr->doVetoStep( typeVetoStep, 0, nFSRhard,
1036             event);
1037           // Abort event if vetoed.
1038           if (doVeto) return false;
1039         }
1040
1041       // Keep on evolving until nothing is left to be done.
1042       } while (pTmax > 0.);
1043
1044     }
1045
1046   // No more systems to be processed. Set total number of emissions.
1047   }
1048   nFSRinRes = nFSRres;
1049   return true;
1050
1051 }
1052
1053 //==========================================================================
1054
1055 } // end namespace Pythia8