PYTHIA8/pythia8130/src/PartonLevel.cxx

   1 // PartonLevel.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2008 Torbjorn Sjostrand.
   3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
   4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
   5
   6 // Function definitions (not found in the header) for the 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, BeamParticle* beamAPtrIn,
  29   BeamParticle* beamBPtrIn,  SigmaTotal* sigmaTotPtr,
  30   TimeShower* timesDecPtrIn, TimeShower* timesPtrIn,
  31   SpaceShower* spacePtrIn, UserHooks* userHooksPtrIn) {
  32
  33   // Store input pointers and modes for future use.
  34   infoPtr            = infoPtrIn;
  35   beamAPtr           = beamAPtrIn;
  36   beamBPtr           = beamBPtrIn;
  37   timesDecPtr        = timesDecPtrIn;
  38   timesPtr           = timesPtrIn;
  39   spacePtr           = spacePtrIn;
  40   userHooksPtr       = userHooksPtrIn;
  41
  42   // Main flags.
  43   doMI               = Settings::flag("PartonLevel:MI");
  44   doISR              = Settings::flag("PartonLevel:ISR");
  45   bool FSR           = Settings::flag("PartonLevel:FSR");
  46   bool FSRinProcess  = Settings::flag("PartonLevel:FSRinProcess");
  47   bool interleaveFSR = Settings::flag("TimeShower:interleave");
  48   doFSRduringProcess = FSR && FSRinProcess &&  interleaveFSR;
  49   doFSRafterProcess  = FSR && FSRinProcess && !interleaveFSR;
  50   doFSRinResonances  = FSR && Settings::flag("PartonLevel:FSRinResonances");
  51   doRemnants         = true;
  52   doSecondHard       = Settings::flag("SecondHard:generate");
  53
  54   // Need MI initialization for minbias processes, even if only first MI.
  55   // But no need to initialize MI if never going to use it.
  56   doMIinit           = doMI;
  57   if (Settings::flag("SoftQCD:minBias") || Settings::flag("SoftQCD:all"))
  58     doMIinit = true;
  59   if (!Settings::flag("PartonLevel:all")) doMIinit = false;
  60
  61   // Flag if lepton beams, and if non-resolved ones. May change main flags.
  62   hasLeptonBeams  = ( beamAPtr->isLepton() || beamBPtr->isLepton() );
  63   hasPointLeptons = ( hasLeptonBeams
  64     && (beamAPtr->isUnresolved() || beamBPtr->isUnresolved() ) );
  65   if (hasLeptonBeams) {
  66     doMI       = false;
  67     doMIinit   = false;
  68   }
  69   if (hasPointLeptons) {
  70     doISR      = false;
  71     doRemnants = false;
  72   }
  73
  74   // Possibility to allow user veto during evolution.
  75   canVetoPT   = (userHooksPtr > 0) ? userHooksPtr->canVetoPT()   : false;
  76   pTvetoPT    = (canVetoPT)        ? userHooksPtr->scaleVetoPT() : -1.;
  77   canVetoStep = (userHooksPtr > 0) ? userHooksPtr->canVetoStep() : false;
  78   nVetoStep   = (canVetoStep)   ? userHooksPtr->numberVetoStep() : -1;
  79
  80   // Set info and initialize the respective program elements.
  81   timesPtr->init( beamAPtr, beamBPtr);
  82   if (doISR) spacePtr->init( beamAPtr, beamBPtr);
  83   doMI = multi.init( doMIinit, infoPtr, beamAPtr, beamBPtr, sigmaTotPtr);
  84   remnants.init( infoPtr, beamAPtr, beamBPtr);
  85
  86   // Succeeded. (Check return values from other classes??)
  87   return true;
  88 }
  89
  90 //*********
  91
  92 // Main routine to do the parton-level evolution.
  93
  94 bool PartonLevel::next( Event& process, Event& event) {
  95
  96   // Special case if unresolved = elastic/diffractive event.
  97   if (!infoPtr->isResolved()) return setupUnresolvedSys( process, event);
  98
  99   // Special case if minimum bias: do hardest interaction.
 100   if (doMIinit) multi.clear();
 101   if (infoPtr->isMinBias()) {
 102     multi.pTfirst();
 103     multi.setupFirstSys( process);
 104   }
 105
 106   // Allow up to ten tries; failure possible for beam remnants.
 107   // Main cause: inconsistent colour flow at the end of the day.
 108   bool physical = true;
 109   int  nRad     = 0;
 110   for (int iTry = 0; iTry < NTRY; ++ iTry) {
 111
 112     // Reset flag, counters and max scales.
 113     physical   = true;
 114     nMI        = (doSecondHard) ? 2 : 1;
 115     nISR       = 0;
 116     nFSRinProc = 0;
 117     nFSRinRes  = 0;
 118     nISRhard   = 0;
 119     nFSRhard   = 0;
 120     pTsaveMI   = 0.;
 121     pTsaveISR  = 0.;
 122     pTsaveFSR  = 0.;
 123
 124     // Identify hard interaction system for showers.
 125     setupHardSys( process, event);
 126
 127     // Check matching of process scale to maximum ISR/MI scales.
 128     double Q2Fac       = infoPtr->Q2Fac();
 129     double Q2Ren       = infoPtr->Q2Ren();
 130     bool limitPTmaxISR = (doISR)
 131       ? spacePtr->limitPTmax( event, Q2Fac, Q2Ren) : false;
 132     bool limitPTmaxMI  = (doMI)  ? multi.limitPTmax( event) : false;
 133
 134     // Set hard scale, maximum for showers and multiple interactions,
 135     double pTscale  = process.scale();
 136     if (doSecondHard) pTscale = max( pTscale, process.scaleSecond() );
 137     double pTmaxMI  = (limitPTmaxMI)  ? pTscale : infoPtr->eCM();
 138     double pTmaxISR = (limitPTmaxISR) ? spacePtr->enhancePTmax() * pTscale
 139                                       : infoPtr->eCM();
 140     double pTmaxFSR = timesPtr->enhancePTmax() * pTscale;
 141     double pTmax    = max( pTmaxMI, max( pTmaxISR, pTmaxFSR) );
 142     pTsaveMI        = pTmaxMI;
 143     pTsaveISR       = pTmaxISR;
 144     pTsaveFSR       = pTmaxFSR;
 145
 146     // Prepare the classes to begin the generation.
 147     if (doMI)  multi.prepare( pTmaxMI);
 148     if (doISR) spacePtr->prepare( 0, event, limitPTmaxISR);
 149     if (doFSRduringProcess) timesPtr->prepare( 0, event);
 150     if (doSecondHard && doISR) spacePtr->prepare( 1, event, limitPTmaxISR);
 151     if (doSecondHard && doFSRduringProcess) timesPtr->prepare( 1, event);
 152
 153     // Set up initial veto scale.
 154     doVeto        = false;
 155     double pTveto = pTvetoPT;
 156     typeLatest    = 0;
 157
 158     // Begin evolution down in pT from hard pT scale.
 159     do {
 160
 161       typeVetoStep = 0;
 162       nRad         =  nISR + nFSRinProc;
 163
 164       // Find next pT value for FSR, MI and ISR.
 165       // Order calls to minimize time expenditure.
 166       double pTgen = 0.;
 167       double pTtimes = (doFSRduringProcess)
 168         ? timesPtr->pTnext( event, pTmaxFSR, pTgen) : -1.;
 169       pTgen = max( pTgen, pTtimes);
 170       double pTmulti = (doMI)
 171         ? multi.pTnext( pTmaxMI, pTgen) : -1.;
 172       pTgen = max( pTgen, pTmulti);
 173       double pTspace = (doISR)
 174         ? spacePtr->pTnext( event, pTmaxISR, pTgen, nRad) : -1.;
 175
 176       // Allow a user veto. Only do it once, so remember to change pTveto.
 177       if (pTveto > 0. && pTveto > pTmulti && pTveto > pTspace
 178         && pTveto > pTtimes) {
 179         pTveto = -1.;
 180         doVeto = userHooksPtr->doVetoPT( typeLatest, event);
 181         // Abort event if vetoed.
 182         if (doVeto) return false;
 183       }
 184
 185       // Do a multiple interaction (if allowed).
 186       if (pTmulti > 0. && pTmulti > pTspace && pTmulti > pTtimes) {
 187         multi.scatter( event);
 188         typeLatest = 1;
 189         ++nMI;
 190         if (doISR)              spacePtr->prepare( nMI - 1, event);
 191         if (doFSRduringProcess) timesPtr->prepare( nMI - 1, event);
 192         pTmaxMI  = pTmulti;
 193         pTmaxISR = min( pTmulti, pTmaxISR);
 194         pTmaxFSR = min( pTmulti, pTmaxFSR);
 195         pTmax    = pTmulti;
 196       }
 197
 198       // Do an initial-state emission (if allowed).
 199       else if (pTspace > 0. && pTspace > pTtimes) {
 200         if (spacePtr->branch( event)) {
 201           typeLatest = 2;
 202           iSysNow = spacePtr->system();
 203           ++nISR;
 204           if (iSysNow == 0) ++nISRhard;
 205           if (doFSRduringProcess) timesPtr->update( iSysNow, event);
 206           if (canVetoStep && iSysNow == 0 && nISRhard <= nVetoStep)
 207             typeVetoStep = 2;
 208         }
 209         pTmaxMI  = min( pTspace, pTmaxMI);
 210         pTmaxISR = pTspace;
 211         pTmaxFSR = min( pTspace, pTmaxFSR);
 212         pTmax    = pTspace;
 213       }
 214
 215       // Do a final-state emission (if allowed).
 216       else if (pTtimes > 0.) {
 217         if (timesPtr->branch( event)) {
 218           typeLatest = 3;
 219           iSysNow = timesPtr->system();
 220           ++nFSRinProc;
 221           if (iSysNow == 0) ++nFSRhard;
 222           if (doISR) spacePtr->update( iSysNow, event);
 223           if (canVetoStep && iSysNow == 0 && nFSRhard <= nVetoStep)
 224             typeVetoStep = 3;
 225         }
 226         pTmaxMI  = min( pTtimes, pTmaxMI);
 227         pTmaxISR = min( pTtimes, pTmaxISR);
 228         pTmaxFSR = pTtimes;
 229         pTmax    = pTtimes;
 230       }
 231
 232       // If no pT scales above zero then nothing to be done.
 233       else pTmax = 0.;
 234
 235       // Optionally check for a veto after the first few emissions.
 236       if (typeVetoStep > 0) {
 237         doVeto = userHooksPtr->doVetoStep( typeVetoStep, nISRhard,
 238           nFSRhard, event);
 239         // Abort event if vetoed.
 240         if (doVeto) return false;
 241       }
 242
 243     // Keep on evolving until nothing is left to be done.
 244     } while (pTmax > 0.);
 245
 246     // Do all final-state emissions if not already considered above.
 247     if (doFSRafterProcess) {
 248
 249       // Find largest scale for final partons.
 250       pTmax = 0.;
 251       for (int i = 0; i < event.size(); ++i)
 252         if (event[i].isFinal() && event[i].scale() > pTmax)
 253           pTmax = event[i].scale();
 254       pTsaveFSR = pTmax;
 255
 256       // Prepare all subsystems for evolution.
 257       for (int iSys = 0; iSys < event.sizeSystems(); ++iSys)
 258         timesPtr->prepare( iSys, event);
 259
 260       // Set up initial veto scale.
 261       doVeto = false;
 262       pTveto = pTvetoPT;
 263
 264       // Begin evolution down in pT from hard pT scale.
 265       do {
 266         typeVetoStep = 0;
 267         double pTtimes = timesPtr->pTnext( event, pTmax, 0.);
 268
 269         // Allow a user veto. Only do it once, so remember to change pTveto.
 270         if (pTveto > 0. && pTveto > pTtimes) {
 271           pTveto = -1.;
 272           doVeto = userHooksPtr->doVetoPT( 4, event);
 273           // Abort event if vetoed.
 274           if (doVeto) return false;
 275         }
 276
 277         // Do a final-state emission (if allowed).
 278         if (pTtimes > 0.) {
 279           if (timesPtr->branch( event)) {
 280             iSysNow = timesPtr->system();
 281             ++nFSRinProc;
 282             if (iSysNow == 0) ++nFSRhard;
 283             if (canVetoStep && iSysNow == 0 && nFSRhard <= nVetoStep)
 284               typeVetoStep = 4;
 285           }
 286           pTmax = pTtimes;
 287         }
 288
 289         // If no pT scales above zero then nothing to be done.
 290         else pTmax = 0.;
 291
 292         // Optionally check for a veto after the first few emissions.
 293         if (typeVetoStep > 0) {
 294           doVeto = userHooksPtr->doVetoStep( typeVetoStep, nISRhard,
 295             nFSRhard, event);
 296           // Abort event if vetoed.
 297           if (doVeto) return false;
 298         }
 299
 300       // Keep on evolving until nothing is left to be done.
 301       } while (pTmax > 0.);
 302     }
 303
 304     // Add beam remnants, including primordial kT kick and colour tracing.
 305     if (doRemnants && !remnants.add( event)) physical = false;
 306
 307     // If no problems then done, else restore and loop.
 308     if (physical) break;
 309     event.clear();
 310     beamAPtr->clear();
 311     beamBPtr->clear();
 312
 313   // End loop over ten tries. Hopefully it worked
 314   }
 315   if (!physical) return false;
 316
 317   // Perform showers in resonance decay chains.
 318   nFSRinRes = resonanceShowers( process, event);
 319
 320   // Store event properties.
 321   infoPtr->setImpact( multi.bMI(), multi.enhanceMI());
 322   infoPtr->setEvolution( pTsaveMI, pTsaveISR, pTsaveFSR,
 323     nMI, nISR, nFSRinProc, nFSRinRes);
 324
 325   // Done.
 326   return true;
 327 }
 328
 329 //*********
 330
 331 // Set up the hard process(es), excluding subsequent resonance decays.
 332
 333 void PartonLevel::setupHardSys( Event& process, Event& event) {
 334
 335   // Incoming partons to hard process are stored in slots 3 and 4.
 336   int inP = 3;
 337   int inM = 4;
 338
 339   // If two hard interactions then find where second begins.
 340   int iBeginSecond = process.size();
 341   if (doSecondHard) {
 342     iBeginSecond = 5;
 343     while (process[iBeginSecond].status() != -21) ++iBeginSecond;
 344   }
 345
 346   // If incoming partons are massive then recalculate to put them massless.
 347   if (process[inP].m() != 0. || process[inM].m() != 0.) {
 348     double pPlus  = process[inP].pPlus() + process[inM].pPlus();
 349     double pMinus = process[inP].pMinus() + process[inM].pMinus();
 350     process[inP].pz( 0.5 * pPlus);
 351     process[inP].e(  0.5 * pPlus);
 352     process[inP].m(  0.);
 353     process[inM].pz(-0.5 * pMinus);
 354     process[inM].e(  0.5 * pMinus);
 355     process[inM].m(  0.);
 356   }
 357
 358   // Add incoming hard-scattering partons to list in beam remnants.
 359   double x1 = process[inP].pPlus() / process[0].e();
 360   beamAPtr->append( inP, process[inP].id(), x1);
 361   double x2 = process[inM].pMinus() / process[0].e();
 362   beamBPtr->append( inM, process[inM].id(), x2);
 363
 364   // Scale. Find whether incoming partons are valence or sea. Store.
 365   double scale = process.scale();
 366   beamAPtr->xfISR( 0, process[inP].id(), x1, scale*scale);
 367   int vsc1 = beamAPtr->pickValSeaComp();
 368   beamBPtr->xfISR( 0, process[inM].id(), x2, scale*scale);
 369   int vsc2 = beamBPtr->pickValSeaComp();
 370   bool isVal1 = (vsc1 == -3);
 371   bool isVal2 = (vsc2 == -3);
 372   infoPtr->setValence( isVal1, isVal2);
 373
 374   // Initialize info needed for subsequent sequential decays + showers.
 375   nHardDone = 0;
 376   iPosBefShow.resize( process.size() );
 377
 378   // Add the beam and hard subprocess partons to the event record.
 379   for (int i = 0; i < iBeginSecond; ++ i) {
 380     if (process[i].mother1() > inM) break;
 381     event.append(process[i]);
 382     iPosBefShow[i] = i;
 383
 384     // Currently outgoing ones should not count as decayed.
 385     if (event[i].status() == -22) {
 386       event[i].statusPos();
 387       event[i].daughters(0, 0);
 388     }
 389
 390     // Complete task of copying hard subsystem into event record.
 391     ++nHardDone;
 392   }
 393
 394   // Store participating partons as first set in list of all systems.
 395   event.newSystem();
 396   for (int i = inP; i < nHardDone; ++i) event.addToSystem(0, i);
 397
 398   // Identify second hard process where applicable.
 399   // Since internally generated incoming partons are guaranteed massless.
 400   if (doSecondHard) {
 401     int inP2 = iBeginSecond;
 402     int inM2 = iBeginSecond + 1;
 403
 404     // Add incoming hard-scattering partons to list in beam remnants.
 405     // Not valid if not in rest frame??
 406     x1 = process[inP2].pPlus() / process[0].e();
 407     beamAPtr->append( inP2, process[inP2].id(), x1);
 408     x2 = process[inM2].pMinus() / process[0].e();
 409     beamBPtr->append( inM2, process[inM2].id(), x2);
 410
 411     // Find whether incoming partons are valence or sea.
 412     scale = process.scaleSecond();
 413     beamAPtr->xfISR( 1, process[inP2].id(), x1, scale*scale);
 414     beamAPtr->pickValSeaComp();
 415     beamBPtr->xfISR( 1, process[inM2].id(), x2, scale*scale);
 416     beamBPtr->pickValSeaComp();
 417
 418     // Add the beam and hard subprocess partons to the event record.
 419     for (int i = inP2; i < process.size(); ++ i) {
 420       int mother = process[i].mother1();
 421       if ( (mother > 2 && mother < inP2) || mother > inM2 ) break;
 422       event.append(process[i]);
 423       iPosBefShow[i] = i;
 424
 425       // Currently outgoing ones should not count as decayed.
 426       if (event[i].status() == -22) {
 427         event[i].statusPos();
 428         event[i].daughters(0, 0);
 429       }
 430
 431       // Complete task of copying hard subsystem into event record.
 432       ++nHardDone;
 433     }
 434
 435     // Store participating partons as second set in list of all systems.
 436     event.newSystem();
 437     for (int i = inP2; i < nHardDone; ++i) event.addToSystem(1, i);
 438
 439   // End code for second hard process.
 440   }
 441
 442   // Update event colour tag to maximum in whole process.
 443   int maxColTag = 0;
 444   for (int i = 0; i < process.size(); ++ i) {
 445     if (process[i].col() > maxColTag) maxColTag = process[i].col();
 446     if (process[i].acol() > maxColTag) maxColTag = process[i].acol();
 447   }
 448   event.initColTag(maxColTag);
 449
 450   // Copy junctions from process to event.
 451   for (int i = 0; i < process.sizeJunction(); ++i)
 452     event.appendJunction( process.getJunction(i));
 453
 454   // Done.
 455 }
 456
 457 //*********
 458
 459 // Set up an unresolved process, i.e. elastic or diffractive.
 460
 461 bool PartonLevel::setupUnresolvedSys( Event& process, Event& event) {
 462
 463   // Copy particles from process to event.
 464   for (int i = 0; i < process.size(); ++ i) event.append( process[i]);
 465
 466   // Loop to find diffractively excited beams.
 467   for (int i = 0; i < 2; ++i)
 468   if ( (i == 0 && infoPtr->isDiffractiveA())
 469     || (i == 1 && infoPtr->isDiffractiveB()) ) {
 470     int iBeam = i + 3;
 471     BeamParticle* beamPtr = (i == 0) ? beamAPtr : beamBPtr;
 472
 473     // Diffractive mass. Reconstruct boost and rotation to event cm frame.
 474     double mDiff  = process[iBeam].m();
 475     double m2Diff = mDiff*mDiff;
 476     double beta   = process[iBeam].pAbs() / process[iBeam].e();
 477     double theta  = process[iBeam].theta();
 478     double phi    = process[iBeam].phi();
 479
 480     // Pick quark or gluon kicked out and flavour subdivision.
 481     bool gluonIsKicked = beamPtr->pickGluon(mDiff);
 482     int id1 = beamPtr->pickValence();
 483     int id2 = beamPtr->pickRemnant();
 484
 485     // Find flavour masses. Scale them down if too big.
 486     double m1 = ParticleDataTable::constituentMass(id1);
 487     double m2 = ParticleDataTable::constituentMass(id2);
 488     if (m1 + m2 > 0.5 * mDiff) {
 489       double reduce = 0.5 * mDiff / (m1 + m2);
 490       m1 *= reduce;
 491       m2 *= reduce;
 492     }
 493
 494     // If quark is kicked out, then trivial kinematics in rest frame.
 495     if (!gluonIsKicked) {
 496       double pAbs = sqrt( pow2(m2Diff - m1*m1 - m2*m2)
 497         - pow2(2. * m1 * m2) ) / (2. * mDiff);
 498       double e1 = (m2Diff + m1*m1 - m2*m2) / (2. * mDiff);
 499       double e2 = (m2Diff + m2*m2 - m1*m1) / (2. * mDiff);
 500       Vec4 p1(0.,0., -pAbs, e1);
 501       Vec4 p2(0.,0., pAbs, e2);
 502
 503       // Boost and rotate to event cm frame.
 504       p1.bst(0., 0., beta); p1.rot(theta, phi);
 505       p2.bst(0., 0., beta); p2.rot(theta, phi);
 506
 507       // Set colours.
 508       int col1, acol1, col2, acol2;
 509       if (ParticleDataTable::colType(id1) == 1) {
 510         col1 = event.nextColTag(); acol1 = 0;
 511         col2 = 0; acol2 = col1;
 512       } else {
 513         col1 = 0; acol1 = event.nextColTag();
 514         col2 = acol1; acol2 = 0;
 515       }
 516
 517       // Store partons of diffractive system and mark system decayed.
 518       int iDauBeg = event.append( id1, 23, iBeam, 0, 0, 0, col1, acol1,
 519         p1, m1);
 520       int iDauEnd = event.append( id2, 63, iBeam, 0, 0, 0, col2, acol2,
 521         p2, m2);
 522       event[iBeam].statusNeg();
 523       event[iBeam].daughters(iDauBeg, iDauEnd);
 524
 525
 526     // If gluon is kicked out: share momentum between two remnants.
 527     } else {
 528       double m2Sys, zSys, pxSys, pySys, mTS1, mTS2;
 529       zSys = beamPtr->zShare(mDiff, m1, m2);
 530
 531       // Provide relative pT kick in remnant. Construct (transverse) masses.
 532       pxSys = beamPtr->pxShare();
 533       pySys = beamPtr->pyShare();
 534       mTS1  = m1*m1 + pxSys*pxSys + pySys*pySys;
 535       mTS2  = m2*m2 + pxSys*pxSys + pySys*pySys;
 536       m2Sys = mTS1 / zSys + mTS2 / (1. - zSys);
 537
 538       // Momentum of kicked-out massless gluon in diffractive rest frame.
 539       double pAbs = (m2Diff - m2Sys) / (2. * mDiff);
 540       Vec4 pG(0., 0., -pAbs, pAbs);
 541       Vec4 pRem(0., 0., pAbs, mDiff - pAbs);
 542
 543       // Momenta of the two beam remnant flavours. (Lightcone p+ = m_diff!)
 544       double e1 = 0.5 * (zSys * mDiff + mTS1 / (zSys * mDiff));
 545       double pL1 = 0.5 * (zSys * mDiff - mTS1 / (zSys * mDiff));
 546       Vec4 p1(pxSys, pySys, pL1, e1);
 547       Vec4 p2 = pRem - p1;
 548
 549       // Boost and rotate to event cm frame.
 550       pG.bst(0., 0., beta); pG.rot(theta, phi);
 551       p1.bst(0., 0., beta); p1.rot(theta, phi);
 552       p2.bst(0., 0., beta); p2.rot(theta, phi);
 553
 554       // Set colours.
 555       int colG, acolG, col1, acol1, col2, acol2;
 556       if (ParticleDataTable::colType(id1) == 1) {
 557         col1 = event.nextColTag(); acol1 = 0;
 558         colG = event.nextColTag(); acolG = col1;
 559         col2 = 0; acol2 = colG;
 560       } else {
 561         col1 = 0; acol1 = event.nextColTag();
 562         colG = acol1; acolG = event.nextColTag();
 563         col2 = acolG; acol2 = 0;
 564       }
 565
 566       // Store partons of diffractive system and mark system decayed.
 567       int iDauBeg = event.append( 21, 23, iBeam, 0, 0, 0, colG, acolG,
 568         pG, m1);
 569       event.append( id1, 63, iBeam, 0, 0, 0, col1, acol1, p1, m1);
 570       int iDauEnd = event.append( id2, 63, iBeam, 0, 0, 0, col2, acol2,
 571         p2, m2);
 572       event[iBeam].statusNeg();
 573       event[iBeam].daughters(iDauBeg, iDauEnd);
 574     }
 575
 576   // End loop over beams. Done.
 577   }
 578   return true;
 579 }
 580
 581 //*********
 582
 583 // Handle showers in successive resonance decays.
 584
 585 int PartonLevel::resonanceShowers( Event& process, Event& event) {
 586
 587   // Isolate next system to be processed, if anything remains.
 588   int nFSRres = 0;
 589   while (nHardDone < process.size()) {
 590     int iBegin = nHardDone;
 591
 592     // Mother in hard process and in complete event (after shower).
 593     int iHardMother      = process[iBegin].mother1();
 594     Particle& hardMother = process[iHardMother];
 595     int iBefMother       = iPosBefShow[iHardMother];
 596     int iAftMother       = event.iBotCopyId(iBefMother);
 597     Particle& aftMother  = event[iAftMother];
 598
 599     // From now mother counts as decayed.
 600     aftMother.statusNeg();
 601
 602     // Mother can have been moved by showering (in any of previous steps),
 603     // so prepare to update colour and momentum information for system.
 604     int colBef  = hardMother.col();
 605     int acolBef = hardMother.acol();
 606     int colAft  = aftMother.col();
 607     int acolAft = aftMother.acol();
 608     RotBstMatrix M;
 609     M.bst( hardMother.p(), aftMother.p());
 610
 611     // Extract next partons from hard event into normal event record.
 612     for (int i = iBegin; i < process.size(); ++ i) {
 613       if (process[i].mother1() != iHardMother) break;
 614       int iNow = event.append( process[i] );
 615       iPosBefShow[i] = iNow;
 616       Particle& now = event.back();
 617
 618       // Currently outgoing ones should not count as decayed.
 619       if (now.status() == -22) {
 620         now.statusPos();
 621         now.daughters(0, 0);
 622       }
 623
 624       // Update daughter and mother information.
 625       if (i == iBegin) aftMother.daughter1( iNow);
 626       else aftMother.daughter2( iNow);
 627       now.mother1(iAftMother);
 628
 629       // Update colour and momentum information.
 630       if (now.col() == colBef) now.col( colAft);
 631       if (now.acol() == acolBef) now.acol( acolAft);
 632       now.rotbst( M);
 633
 634       // Complete task of copying next subsystem into event record.
 635       ++nHardDone;
 636     }
 637     int iEnd = nHardDone - 1;
 638
 639     // Do parton showers inside subsystem: maximum scale by mother mass.
 640     if (doFSRinResonances) {
 641       double pTmax = 0.5 * hardMother.m();
 642       nFSRhard     = 0;
 643
 644       // Add new system, with two empty beam slots.
 645       int iSys = event.newSystem();
 646       event.addToSystem( iSys, 0);
 647       event.addToSystem( iSys, 0);
 648
 649       // Loop over allowed range to find all final-state particles.
 650       for (int i = iPosBefShow[iBegin]; i <= iPosBefShow[iEnd]; ++i)
 651       if (event[i].isFinal()) event.addToSystem( iSys, i);
 652
 653       // Let prepare routine do the setup.
 654       timesDecPtr->prepare( iSys, event);
 655
 656       // Set up initial veto scale.
 657       doVeto        = false;
 658       double pTveto = pTvetoPT;
 659
 660       // Begin evolution down in pT from hard pT scale.
 661       do {
 662         typeVetoStep = 0;
 663         double pTtimes = timesDecPtr->pTnext( event, pTmax, 0.);
 664
 665         // Allow a user veto. Only do it once, so remember to change pTveto.
 666         if (pTveto > 0. && pTveto > pTtimes) {
 667           pTveto = -1.;
 668           doVeto = userHooksPtr->doVetoPT( 5, event);
 669           // Abort event if vetoed.
 670           if (doVeto) return false;
 671         }
 672
 673         // Do a final-state emission (if allowed).
 674         if (pTtimes > 0.) {
 675           if (timesDecPtr->branch( event)) {
 676             ++nFSRres;
 677             ++nFSRhard;
 678             if (canVetoStep && nFSRhard <= nVetoStep) typeVetoStep = 5;
 679           }
 680           pTmax = pTtimes;
 681         }
 682
 683         // If no pT scales above zero then nothing to be done.
 684         else pTmax = 0.;
 685
 686         // Optionally check for a veto after the first few emissions.
 687         if (typeVetoStep > 0) {
 688           doVeto = userHooksPtr->doVetoStep( typeVetoStep, 0, nFSRhard,
 689             event);
 690           // Abort event if vetoed.
 691           if (doVeto) return false;
 692         }
 693
 694       // Keep on evolving until nothing is left to be done.
 695       } while (pTmax > 0.);
 696
 697     }
 698
 699   // No more systems to be processed. Return total number of emissions.
 700   }
 701   return nFSRres;
 702
 703 }
 704
 705 //**************************************************************************
 706
 707 } // end namespace Pythia8