PYTHIA8/pythia8170/src/TimeShower.cxx

   1 // TimeShower.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 TimeShower class.
   7
   8 #include "TimeShower.h"
   9
  10 namespace Pythia8 {
  11
  12 //==========================================================================
  13
  14 // The TimeShower 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 // For small x approximate 1 - sqrt(1 - x) by x/2.
  22 const double TimeShower::SIMPLIFYROOT = 1e-8;
  23
  24 // Do not allow x too close to 0 or 1 in matrix element expressions.
  25 // Warning: cuts into phase space for E_CM > 2 * pTmin * sqrt(1/XMARGIN),
  26 // i.e. will become problem roughly for E_CM > 10^6 GeV.
  27 const double TimeShower::XMARGIN      = 1e-12;
  28 const double TimeShower::XMARGINCOMB  = 1e-4;
  29
  30 // Lower limit on PDF value in order to avoid division by zero.
  31 const double TimeShower::TINYPDF      = 1e-10;
  32
  33 // Big starting value in search for smallest invariant-mass pair.
  34 const double TimeShower::LARGEM2      = 1e20;
  35
  36 // In g -> q qbar or gamma -> f fbar require m2_pair > this * m2_q/f.
  37 const double TimeShower::THRESHM2      = 4.004;
  38
  39 // Never pick pT so low that alphaS is evaluated too close to Lambda_3.
  40 const double TimeShower::LAMBDA3MARGIN = 1.1;
  41
  42 // Rescatter: rescattering + ISR + FSR + primordial kT can lead to
  43 //            systems not locally conserving momentum.
  44 // Fix up momentum in intermediate systems with rescattering
  45 const bool   TimeShower::FIXRESCATTER          = true;
  46 // Veto negative energies when using FIXRESCATTER option.
  47 const bool   TimeShower::VETONEGENERGY         = false;
  48 // Do not allow too large time- or spacelike virtualities in fixing-up.
  49 const double TimeShower::MAXVIRTUALITYFRACTION = 0.5;
  50 // Do not allow too large negative spacelike energy in system rest frame.
  51 const double TimeShower::MAXNEGENERGYFRACTION  = 0.7;
  52
  53 //--------------------------------------------------------------------------
  54
  55 // Initialize alphaStrong, alphaEM and related pTmin parameters.
  56
  57 void TimeShower::init( BeamParticle* beamAPtrIn,
  58   BeamParticle* beamBPtrIn) {
  59
  60   // Store input pointers for future use.
  61   beamAPtr           = beamAPtrIn;
  62   beamBPtr           = beamBPtrIn;
  63
  64   // Main flags.
  65   doQCDshower        = settingsPtr->flag("TimeShower:QCDshower");
  66   doQEDshowerByQ     = settingsPtr->flag("TimeShower:QEDshowerByQ");
  67   doQEDshowerByL     = settingsPtr->flag("TimeShower:QEDshowerByL");
  68   doQEDshowerByGamma = settingsPtr->flag("TimeShower:QEDshowerByGamma");
  69   doMEcorrections    = settingsPtr->flag("TimeShower:MEcorrections");
  70   doMEafterFirst     = settingsPtr->flag("TimeShower:MEafterFirst");
  71   doPhiPolAsym       = settingsPtr->flag("TimeShower:phiPolAsym");
  72   doInterleave       = settingsPtr->flag("TimeShower:interleave");
  73   allowBeamRecoil    = settingsPtr->flag("TimeShower:allowBeamRecoil");
  74   dampenBeamRecoil   = settingsPtr->flag("TimeShower:dampenBeamRecoil");
  75   recoilToColoured   = settingsPtr->flag("TimeShower:recoilToColoured");
  76
  77   // Matching in pT of hard interaction or MPI to shower evolution.
  78   pTmaxMatch         = settingsPtr->mode("TimeShower:pTmaxMatch");
  79   pTdampMatch        = settingsPtr->mode("TimeShower:pTdampMatch");
  80   pTmaxFudge         = settingsPtr->parm("TimeShower:pTmaxFudge");
  81   pTmaxFudgeMPI      = settingsPtr->parm("TimeShower:pTmaxFudgeMPI");
  82   pTdampFudge        = settingsPtr->parm("TimeShower:pTdampFudge");
  83
  84   // Charm and bottom mass thresholds.
  85   mc                 = particleDataPtr->m0(4);
  86   mb                 = particleDataPtr->m0(5);
  87   m2c                = mc * mc;
  88   m2b                = mb * mb;
  89
  90   // Parameters of scale choices.
  91   renormMultFac     = settingsPtr->parm("TimeShower:renormMultFac");
  92   factorMultFac     = settingsPtr->parm("TimeShower:factorMultFac");
  93
  94   // Parameters of alphaStrong generation.
  95   alphaSvalue        = settingsPtr->parm("TimeShower:alphaSvalue");
  96   alphaSorder        = settingsPtr->mode("TimeShower:alphaSorder");
  97   alphaS2pi          = 0.5 * alphaSvalue / M_PI;
  98
  99   // Initialize alphaStrong generation.
 100   alphaS.init( alphaSvalue, alphaSorder);
 101
 102   // Lambda for 5, 4 and 3 flavours.
 103   Lambda3flav        = alphaS.Lambda3();
 104   Lambda4flav        = alphaS.Lambda4();
 105   Lambda5flav        = alphaS.Lambda5();
 106   Lambda5flav2       = pow2(Lambda5flav);
 107   Lambda4flav2       = pow2(Lambda4flav);
 108   Lambda3flav2       = pow2(Lambda3flav);
 109
 110   // Parameters of QCD evolution. Warn if pTmin must be raised.
 111   nGluonToQuark      = settingsPtr->mode("TimeShower:nGluonToQuark");
 112   pTcolCutMin        = settingsPtr->parm("TimeShower:pTmin");
 113   if (pTcolCutMin > LAMBDA3MARGIN * Lambda3flav / sqrt(renormMultFac))
 114     pTcolCut         = pTcolCutMin;
 115   else {
 116     pTcolCut         = LAMBDA3MARGIN * Lambda3flav / sqrt(renormMultFac);
 117     ostringstream newPTcolCut;
 118     newPTcolCut << fixed << setprecision(3) << pTcolCut;
 119     infoPtr->errorMsg("Warning in TimeShower::init: pTmin too low",
 120                       ", raised to " + newPTcolCut.str() );
 121     infoPtr->setTooLowPTmin(true);
 122   }
 123   pT2colCut          = pow2(pTcolCut);
 124
 125   // Parameters of alphaEM generation .
 126   alphaEMorder       = settingsPtr->mode("TimeShower:alphaEMorder");
 127
 128   // Initialize alphaEM generation.
 129   alphaEM.init( alphaEMorder, settingsPtr);
 130
 131   // Parameters of QED evolution.
 132   nGammaToQuark      = settingsPtr->mode("TimeShower:nGammaToQuark");
 133   nGammaToLepton     = settingsPtr->mode("TimeShower:nGammaToLepton");
 134   pTchgQCut          = settingsPtr->parm("TimeShower:pTminChgQ");
 135   pT2chgQCut         = pow2(pTchgQCut);
 136   pTchgLCut          = settingsPtr->parm("TimeShower:pTminChgL");
 137   pT2chgLCut         = pow2(pTchgLCut);
 138   mMaxGamma          = settingsPtr->parm("TimeShower:mMaxGamma");
 139   m2MaxGamma         = pow2(mMaxGamma);
 140
 141   // Consisteny check for gamma -> f fbar variables.
 142   if (nGammaToQuark <= 0 && nGammaToLepton <= 0) doQEDshowerByGamma = false;
 143
 144   // Possibility of a global recoil stategy, e.g. for MC@NLO.
 145   globalRecoil       = settingsPtr->flag("TimeShower:globalRecoil");
 146   nMaxGlobalRecoil   = settingsPtr->mode("TimeShower:nMaxGlobalRecoil");
 147
 148   // Fraction and colour factor of gluon emission off onium octat state.
 149   octetOniumFraction = settingsPtr->parm("TimeShower:octetOniumFraction");
 150   octetOniumColFac   = settingsPtr->parm("TimeShower:octetOniumColFac");
 151
 152   // Z0 properties needed for gamma/Z0 mixing.
 153   mZ                 = particleDataPtr->m0(23);
 154   gammaZ             = particleDataPtr->mWidth(23);
 155   thetaWRat          = 1. / (16. * coupSMPtr->sin2thetaW()
 156                        * coupSMPtr->cos2thetaW());
 157
 158   // May have to fix up recoils related to rescattering.
 159   allowRescatter     = settingsPtr->flag("PartonLevel:MPI")
 160     && settingsPtr->flag("MultipartonInteractions:allowRescatter");
 161
 162   // Hidden Valley scenario with further shower activity.
 163   doHVshower         = settingsPtr->flag("HiddenValley:FSR");
 164   nCHV               = settingsPtr->mode("HiddenValley:Ngauge");
 165   alphaHVfix         = settingsPtr->parm("HiddenValley:alphaFSR");
 166   pThvCut            = settingsPtr->parm("HiddenValley:pTminFSR");
 167   pT2hvCut           = pThvCut * pThvCut;
 168   CFHV               = (nCHV == 1) ? 1. : (nCHV * nCHV - 1.)/(2. * nCHV);
 169   idHV               = (nCHV == 1) ? 4900022 : 4900021;
 170   mHV                = particleDataPtr->m0(idHV);
 171   brokenHVsym        = (nCHV == 1 && mHV > 0.);
 172
 173   // Possibility to allow user veto of emission step.
 174   canVetoEmission    = (userHooksPtr != 0)
 175                      ? userHooksPtr->canVetoFSREmission() : false;
 176
 177 }
 178
 179 //--------------------------------------------------------------------------
 180
 181 // Find whether to limit maximum scale of emissions.
 182 // Also allow for dampening at factorization or renormalization scale.
 183
 184 bool TimeShower::limitPTmax( Event& event, double Q2Fac, double Q2Ren) {
 185
 186   // Find whether to limit pT. Begin by user-set cases.
 187   bool dopTlimit = false;
 188   if      (pTmaxMatch == 1) dopTlimit = true;
 189   else if (pTmaxMatch == 2) dopTlimit = false;
 190
 191   // Look if any quark (u, d, s, c, b), gluon or photon in final state.
 192   else {
 193     for (int i = 5; i < event.size(); ++i)
 194     if (event[i].status() != -21) {
 195       int idAbs = event[i].idAbs();
 196       if (idAbs <= 5 || idAbs == 21 || idAbs == 22) dopTlimit = true;
 197     }
 198   }
 199
 200   // Dampening at factorization or renormalization scale.
 201   dopTdamp   = false;
 202   pT2damp    = 0.;
 203   if ( !dopTlimit && (pTdampMatch == 1 || pTdampMatch == 2) ) {
 204     dopTdamp = true;
 205     pT2damp  = pow2(pTdampFudge) * ((pTdampMatch == 1) ? Q2Fac : Q2Ren);
 206   }
 207
 208   // Done.
 209   return dopTlimit;
 210
 211 }
 212
 213 //--------------------------------------------------------------------------
 214
 215 // Top-level routine to do a full time-like shower in resonance decay.
 216
 217 int TimeShower::shower( int iBeg, int iEnd, Event& event, double pTmax,
 218   int nBranchMax) {
 219
 220   // Add new system, automatically with two empty beam slots.
 221   int iSys = partonSystemsPtr->addSys();
 222
 223   // Loop over allowed range to find all final-state particles.
 224   Vec4 pSum;
 225   for (int i = iBeg; i <= iEnd; ++i) if (event[i].isFinal()) {
 226     partonSystemsPtr->addOut( iSys, i);
 227     pSum += event[i].p();
 228   }
 229   partonSystemsPtr->setSHat( iSys, pSum.m2Calc() );
 230
 231   // Let prepare routine do the setup.
 232   prepare( iSys, event, true);
 233
 234   // Begin evolution down in pT from hard pT scale.
 235   int nBranch  = 0;
 236   pTLastBranch = 0.;
 237   do {
 238     double pTtimes = pTnext( event, pTmax, 0.);
 239
 240     // Do a final-state emission (if allowed).
 241     if (pTtimes > 0.) {
 242       if (branch( event)) {
 243         ++nBranch;
 244         pTLastBranch = pTtimes;
 245       }
 246       pTmax = pTtimes;
 247     }
 248
 249     // Keep on evolving until nothing is left to be done.
 250     else pTmax = 0.;
 251   } while (pTmax > 0. && (nBranchMax <= 0 || nBranch < nBranchMax));
 252
 253   // Return number of emissions that were performed.
 254   return nBranch;
 255
 256 }
 257
 258 //--------------------------------------------------------------------------
 259
 260 // Prepare system for evolution; identify ME.
 261
 262 void TimeShower::prepare( int iSys, Event& event, bool limitPTmaxIn) {
 263
 264   // Reset dipole-ends list for first interaction and for resonance decays.
 265   int iInA = partonSystemsPtr->getInA(iSys);
 266   int iInB = partonSystemsPtr->getInB(iSys);
 267   if (iSys == 0 || iInA == 0) dipEnd.resize(0);
 268   int dipEndSizeBeg = dipEnd.size();
 269
 270   // No dipoles for 2 -> 1 processes.
 271   if (partonSystemsPtr->sizeOut(iSys) < 2) return;
 272
 273   // Loop through final state of system to find possible dipole ends.
 274   for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) {
 275     int iRad = partonSystemsPtr->getOut( iSys, i);
 276     if (event[iRad].isFinal() && event[iRad].scale() > 0.) {
 277
 278       // Identify colour octet onium state. Check whether QCD shower allowed.
 279       int idRad    = event[iRad].id();
 280       int idRadAbs = abs(idRad);
 281       bool isOctetOnium
 282         = ( idRad == 9900441 || idRad == 9900443 || idRad == 9910441
 283          || idRad == 9900551 || idRad == 9900553 || idRad == 9910551 );
 284       bool doQCD = doQCDshower;
 285       if (doQCD && isOctetOnium)
 286         doQCD = (rndmPtr->flat() < octetOniumFraction);
 287
 288       // Find dipole end formed by colour index.
 289       int colTag = event[iRad].col();
 290       if (doQCD && colTag > 0) setupQCDdip( iSys, i,  colTag,  1, event,
 291         isOctetOnium, limitPTmaxIn);
 292
 293       // Find dipole end formed by anticolour index.
 294       int acolTag = event[iRad].acol();
 295       if (doQCD && acolTag > 0) setupQCDdip( iSys, i, acolTag, -1, event,
 296         isOctetOnium, limitPTmaxIn);
 297
 298       // Find "charge-dipole" and "photon-dipole" ends.
 299       int  chgType  = event[iRad].chargeType();
 300       bool doChgDip = (chgType != 0)
 301                        && ( ( doQEDshowerByQ && event[iRad].isQuark()  )
 302                          || ( doQEDshowerByL && event[iRad].isLepton() ) );
 303       int  gamType  = (idRad == 22) ? 1 : 0;
 304       bool doGamDip = (gamType == 1) && doQEDshowerByGamma;
 305       if (doChgDip || doGamDip) setupQEDdip( iSys, i, chgType, gamType,
 306          event, limitPTmaxIn);
 307
 308       // Find Hidden Valley dipole ends.
 309       bool isHVrad =  (idRadAbs > 4900000 && idRadAbs < 4900007)
 310                    || (idRadAbs > 4900010 && idRadAbs < 4900017)
 311                    || idRadAbs == 4900101;
 312       if (doHVshower && isHVrad) setupHVdip( iSys, i, event, limitPTmaxIn);
 313
 314     // End loop over system final state. Have now found the dipole ends.
 315     }
 316   }
 317
 318   // Loop through dipole ends to find matrix element corrections.
 319   for (int iDip = dipEndSizeBeg; iDip < int(dipEnd.size()); ++iDip)
 320     findMEtype( event, dipEnd[iDip]);
 321
 322   // Update dipole list after a multiparton interactions rescattering.
 323   if (iSys > 0 && ( (iInA > 0 && event[iInA].status() == -34)
 324     || (iInB > 0 && event[iInB].status() == -34) ) )
 325     rescatterUpdate( iSys, event);
 326
 327 }
 328
 329 //--------------------------------------------------------------------------
 330
 331 // Update dipole list after a multiparton interactions rescattering.
 332 void TimeShower::rescatterUpdate( int iSys, Event& event) {
 333
 334   // Loop over two incoming partons in system; find their rescattering mother.
 335   // (iOut is outgoing from old system = incoming iIn of rescattering system.)
 336   for (int iResc = 0; iResc < 2; ++iResc) {
 337     int iIn = (iResc == 0) ? partonSystemsPtr->getInA(iSys)
 338                            : partonSystemsPtr->getInB(iSys);
 339     if (iIn == 0 || event[iIn].status() != -34) continue;
 340     int iOut = event[iIn].mother1();
 341
 342     // Loop over all dipoles.
 343     int dipEndSize = dipEnd.size();
 344     for (int iDip = 0; iDip < dipEndSize; ++iDip) {
 345       TimeDipoleEnd& dipNow = dipEnd[iDip];
 346
 347       // Kill dipoles where rescattered parton is radiator.
 348       if (dipNow.iRadiator == iOut) {
 349         dipNow.colType = 0;
 350         dipNow.chgType = 0;
 351         dipNow.gamType = 0;
 352         continue;
 353       }
 354       // No matrix element for dipoles between scatterings.
 355       if (dipNow.iMEpartner == iOut) {
 356         dipNow.MEtype     =  0;
 357         dipNow.iMEpartner = -1;
 358       }
 359
 360       // Update dipoles where outgoing rescattered parton is recoiler.
 361       if (dipNow.iRecoiler == iOut) {
 362         int iRad = dipNow.iRadiator;
 363
 364         // Colour dipole: recoil in final state, initial state or new.
 365         if (dipNow.colType > 0) {
 366           int colTag = event[iRad].col();
 367           bool done  = false;
 368           for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) {
 369             int iRecNow = partonSystemsPtr->getOut( iSys, i);
 370             if (event[iRecNow].acol() == colTag) {
 371               dipNow.iRecoiler = iRecNow;
 372               dipNow.systemRec = iSys;
 373               dipNow.MEtype    = 0;
 374               done             = true;
 375               break;
 376             }
 377           }
 378           if (!done) {
 379             int iIn2 = (iResc == 0) ? partonSystemsPtr->getInB(iSys)
 380                                     : partonSystemsPtr->getInA(iSys);
 381             if (event[iIn2].col() == colTag) {
 382               dipNow.iRecoiler = iIn2;
 383               dipNow.systemRec = iSys;
 384               dipNow.MEtype    = 0;
 385               int isrType      = event[iIn2].mother1();
 386               // This line in case mother is a rescattered parton.
 387               while (isrType > 2 + beamOffset)
 388                 isrType = event[isrType].mother1();
 389               if (isrType > 2) isrType -= beamOffset;
 390               dipNow.isrType   = isrType;
 391               done             = true;
 392             }
 393           }
 394           // If above options failed, then create new dipole.
 395           if (!done) {
 396             int iRadNow = partonSystemsPtr->getIndexOfOut(dipNow.system, iRad);
 397             if (iRadNow != -1)
 398               setupQCDdip(dipNow.system, iRadNow, event[iRad].col(), 1,
 399                           event, dipNow.isOctetOnium, true);
 400             else
 401               infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: "
 402               "failed to locate radiator in system");
 403
 404             dipNow.colType = 0;
 405             dipNow.chgType = 0;
 406             dipNow.gamType = 0;
 407
 408             infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: "
 409             "failed to locate new recoiling colour partner");
 410           }
 411
 412         // Anticolour dipole: recoil in final state, initial state or new.
 413         } else if (dipNow.colType < 0) {
 414           int  acolTag = event[iRad].acol();
 415           bool done    = false;
 416           for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) {
 417             int iRecNow = partonSystemsPtr->getOut( iSys, i);
 418             if (event[iRecNow].col() == acolTag) {
 419               dipNow.iRecoiler = iRecNow;
 420               dipNow.systemRec = iSys;
 421               dipNow.MEtype    = 0;
 422               done             = true;
 423               break;
 424             }
 425           }
 426           if (!done) {
 427             int iIn2 = (iResc == 0) ? partonSystemsPtr->getInB(iSys)
 428                                     : partonSystemsPtr->getInA(iSys);
 429             if (event[iIn2].acol() == acolTag) {
 430               dipNow.iRecoiler = iIn2;
 431               dipNow.systemRec = iSys;
 432               dipNow.MEtype    = 0;
 433               int isrType      = event[iIn2].mother1();
 434               // This line in case mother is a rescattered parton.
 435               while (isrType > 2 + beamOffset)
 436                 isrType = event[isrType].mother1();
 437               if (isrType > 2) isrType -= beamOffset;
 438               dipNow.isrType   = isrType;
 439               done             = true;
 440             }
 441           }
 442           // If above options failed, then create new dipole.
 443           if (!done) {
 444             int iRadNow = partonSystemsPtr->getIndexOfOut(dipNow.system, iRad);
 445             if (iRadNow != -1)
 446               setupQCDdip(dipNow.system, iRadNow, event[iRad].acol(), -1,
 447                           event, dipNow.isOctetOnium, true);
 448             else
 449               infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: "
 450               "failed to locate radiator in system");
 451
 452             dipNow.colType = 0;
 453             dipNow.chgType = 0;
 454             dipNow.gamType = 0;
 455
 456             infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: "
 457             "failed to locate new recoiling colour partner");
 458           }
 459
 460         // Charge or photon dipoles: same flavour in final or initial state.
 461         } else if (dipNow.chgType != 0 || dipNow.gamType != 0) {
 462           int  idTag = event[dipNow.iRecoiler].id();
 463           bool done  = false;
 464           for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) {
 465             int iRecNow = partonSystemsPtr->getOut( iSys, i);
 466             if (event[iRecNow].id() == idTag) {
 467               dipNow.iRecoiler = iRecNow;
 468               dipNow.systemRec = iSys;
 469               dipNow.MEtype    = 0;
 470               done             = true;
 471               break;
 472             }
 473           }
 474           if (!done) {
 475             int iIn2 = (iResc == 0) ? partonSystemsPtr->getInB(iSys)
 476                                     : partonSystemsPtr->getInA(iSys);
 477             if (event[iIn2].id() == -idTag) {
 478               dipNow.iRecoiler = iIn2;
 479               dipNow.systemRec = iSys;
 480               dipNow.MEtype    = 0;
 481               int isrType      = event[iIn2].mother1();
 482               // This line in case mother is a rescattered parton.
 483               while (isrType > 2 + beamOffset)
 484                 isrType = event[isrType].mother1();
 485               if (isrType > 2) isrType -= beamOffset;
 486               dipNow.isrType   = isrType;
 487               done             = true;
 488             }
 489           }
 490           // If above options failed, then create new dipole
 491           if (!done) {
 492             int iRadNow = partonSystemsPtr->getIndexOfOut(dipNow.system, iRad);
 493             if (iRadNow != -1)
 494               setupQEDdip(dipNow.system, iRadNow, dipNow.chgType,
 495                           dipNow.gamType, event, true);
 496             else
 497               infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: "
 498               "failed to locate radiator in system");
 499
 500             dipNow.colType = 0;
 501             dipNow.chgType = 0;
 502             dipNow.gamType = 0;
 503
 504             infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: "
 505             "failed to locate new recoiling charge partner");
 506           }
 507         }
 508       }
 509
 510     // End of loop over dipoles and two incoming sides.
 511     }
 512   }
 513
 514 }
 515
 516 //--------------------------------------------------------------------------
 517
 518 // Update dipole list after each ISR emission (so not used for resonances).
 519
 520 void TimeShower::update( int iSys, Event& event) {
 521
 522   // Start list of rescatterers that gave further changed systems in ISR.
 523   vector<int> iRescatterer;
 524
 525   // Find new and old positions of incoming partons in the system.
 526   vector<int> iNew, iOld;
 527   iNew.push_back( partonSystemsPtr->getInA(iSys) );
 528   iOld.push_back( event[iNew[0]].daughter2() );
 529   iNew.push_back( partonSystemsPtr->getInB(iSys) );
 530   iOld.push_back( event[iNew[1]].daughter2() );
 531
 532   // Ditto for outgoing partons, except the newly created one.
 533   int sizeOut = partonSystemsPtr->sizeOut(iSys) - 1;
 534   for (int i = 0; i < sizeOut; ++i) {
 535     int iNow = partonSystemsPtr->getOut(iSys, i);
 536     iNew.push_back( iNow );
 537     iOld.push_back( event[iNow].mother1() );
 538     // Add non-final to list of rescatterers.
 539     if (!event[iNow].isFinal()) iRescatterer.push_back( iNow );
 540   }
 541   int iNewNew = partonSystemsPtr->getOut(iSys, sizeOut);
 542
 543   // Swap beams to let 0 be side on which branching occured.
 544   if (event[iNew[0]].status() != -41) {
 545     swap( iNew[0], iNew[1]);
 546     swap( iOld[0], iOld[1]);
 547   }
 548
 549   // Loop over all dipole ends belonging to the system
 550   // or to the recoil system, if different.
 551   for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip)
 552   if (dipEnd[iDip].system == iSys || dipEnd[iDip].systemRec == iSys) {
 553     TimeDipoleEnd& dipNow = dipEnd[iDip];
 554
 555     // Replace radiator (always in final state so simple).
 556     for (int i = 2; i < 2 + sizeOut; ++i)
 557     if (dipNow.iRadiator == iOld[i]) {
 558       dipNow.iRadiator = iNew[i];
 559       break;
 560     }
 561
 562     // Replace ME partner (always in final state, if exists, so simple).
 563     for (int i = 2; i < 2 + sizeOut; ++i)
 564     if (dipNow.iMEpartner == iOld[i]) {
 565       dipNow.iMEpartner = iNew[i];
 566       break;
 567     }
 568
 569     // Recoiler: by default pick old one, only moved. Note excluded beam.
 570     int iRec = 0;
 571     if (dipNow.systemRec == iSys) {
 572       for (int i = 1; i < 2 + sizeOut; ++i)
 573       if (dipNow.iRecoiler == iOld[i]) {
 574         iRec = iNew[i];
 575         break;
 576       }
 577
 578       // QCD recoiler: check if colour hooks up with new final parton.
 579       if ( dipNow.colType > 0
 580         && event[dipNow.iRadiator].col() == event[iNewNew].acol() ) {
 581         iRec = iNewNew;
 582         dipNow.isrType = 0;
 583       }
 584       if ( dipNow.colType < 0
 585         && event[dipNow.iRadiator].acol() == event[iNewNew].col() ) {
 586         iRec = iNewNew;
 587         dipNow.isrType = 0;
 588       }
 589
 590       // QCD recoiler: check if colour hooks up with new beam parton.
 591       if ( iRec == 0 && dipNow.colType > 0
 592         && event[dipNow.iRadiator].col()  == event[iNew[0]].col() )
 593         iRec = iNew[0];
 594       if ( iRec == 0 && dipNow.colType < 0
 595         && event[dipNow.iRadiator].acol() == event[iNew[0]].acol() )
 596         iRec = iNew[0];
 597
 598       // QED/photon recoiler: either to new particle or remains to beam.
 599       if ( iRec == 0 && (dipNow.chgType != 0 || dipNow.gamType != 0) ) {
 600         if ( event[iNew[0]].chargeType() == 0 ) {
 601           iRec = iNewNew;
 602           dipNow.isrType = 0;
 603         } else {
 604           iRec = iNew[0];
 605         }
 606       }
 607
 608     // Recoiler in another system: keep it as is.
 609     } else iRec = dipNow.iRecoiler;
 610
 611     // Done. Kill dipole if failed to find new recoiler.
 612     dipNow.iRecoiler = iRec;
 613     if ( iRec == 0 && (dipNow.colType != 0 || dipNow.chgType != 0
 614       || dipNow.gamType != 0) ) {
 615       dipNow.colType = 0;
 616       dipNow.chgType = 0;
 617       dipNow.gamType = 0;
 618       infoPtr->errorMsg("Error in TimeShower::update: "
 619       "failed to locate new recoiling partner");
 620     }
 621   }
 622
 623   // Find new dipole end formed by colour index.
 624   int colTag = event[iNewNew].col();
 625   if (doQCDshower && colTag > 0)
 626     setupQCDdip( iSys, sizeOut, colTag, 1, event, false, true);
 627
 628   // Find new dipole end formed by anticolour index.
 629   int acolTag = event[iNewNew].acol();
 630   if (doQCDshower && acolTag > 0)
 631     setupQCDdip( iSys, sizeOut, acolTag, -1, event, false, true);
 632
 633   // Find new "charge-dipole" and "photon-dipole" ends.
 634   int  chgType  = event[iNewNew].chargeType();
 635   bool doChgDip = (chgType != 0)
 636                   && ( ( doQEDshowerByQ && event[iNewNew].isQuark()  )
 637                     || ( doQEDshowerByL && event[iNewNew].isLepton() ) );
 638   int  gamType  = (event[iNewNew].id() == 22) ? 1 : 0;
 639   bool doGamDip = (gamType == 1) && doQEDshowerByGamma;
 640   if (doChgDip || doGamDip)
 641     setupQEDdip( iSys, sizeOut, chgType, gamType, event, true);
 642
 643   // Start iterate over list of rescatterers - may be empty.
 644   int iRescNow = -1;
 645   while (++iRescNow < int(iRescatterer.size())) {
 646
 647     // Identify systems that rescatterers belong to.
 648     int iOutNew    = iRescatterer[iRescNow];
 649     int iInNew     = event[iOutNew].daughter1();
 650     int iSysResc   = partonSystemsPtr->getSystemOf(iInNew, true);
 651
 652     // Find new and old positions of incoming partons in the system.
 653     iNew.resize(0);
 654     iOld.resize(0);
 655     iNew.push_back( partonSystemsPtr->getInA(iSysResc) );
 656     iOld.push_back( event[iNew[0]].daughter1() );
 657     iNew.push_back( partonSystemsPtr->getInB(iSysResc) );
 658     iOld.push_back( event[iNew[1]].daughter1() );
 659
 660     // Ditto for outgoing partons.
 661     sizeOut = partonSystemsPtr->sizeOut(iSysResc);
 662     for (int i = 0; i < sizeOut; ++i) {
 663       int iNow = partonSystemsPtr->getOut(iSysResc, i);
 664       iNew.push_back( iNow );
 665       iOld.push_back( event[iNow].mother1() );
 666       // Add non-final to list of rescatterers.
 667       if (!event[iNow].isFinal()) iRescatterer.push_back( iNow );
 668     }
 669
 670     // Loop over all dipole ends belonging to the system
 671     // or to the recoil system, if different.
 672     for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip)
 673     if (dipEnd[iDip].system == iSysResc
 674       || dipEnd[iDip].systemRec == iSysResc) {
 675       TimeDipoleEnd& dipNow = dipEnd[iDip];
 676
 677       // Replace radiator (always in final state so simple).
 678       for (int i = 2; i < 2 + sizeOut; ++i)
 679       if (dipNow.iRadiator == iOld[i]) {
 680         dipNow.iRadiator = iNew[i];
 681         break;
 682       }
 683
 684       // Replace ME partner (always in final state, if exists, so simple).
 685       for (int i = 2; i < 2 + sizeOut; ++i)
 686       if (dipNow.iMEpartner == iOld[i]) {
 687         dipNow.iMEpartner = iNew[i];
 688         break;
 689       }
 690
 691       // Replace recoiler.
 692       for (int i = 0; i < 2 + sizeOut; ++i)
 693       if (dipNow.iRecoiler == iOld[i]) {
 694         dipNow.iRecoiler = iNew[i];
 695         break;
 696       }
 697     }
 698
 699   // End iterate over list of rescatterers.
 700   }
 701
 702 }
 703
 704 //--------------------------------------------------------------------------
 705
 706 // Setup a dipole end for a QCD colour charge.
 707
 708 void TimeShower::setupQCDdip( int iSys, int i, int colTag, int colSign,
 709   Event& event, bool isOctetOnium, bool limitPTmaxIn) {
 710
 711   // Initial values. Find if allowed to hook up beams.
 712   int iRad     = partonSystemsPtr->getOut(iSys, i);
 713   int iRec     = 0;
 714   int sizeAllA = partonSystemsPtr->sizeAll(iSys);
 715   int sizeOut  = partonSystemsPtr->sizeOut(iSys);
 716   int sizeAll  = ( allowBeamRecoil ) ? sizeAllA : sizeOut;
 717   int sizeIn   = sizeAll - sizeOut;
 718   int sizeInA  = sizeAllA - sizeIn - sizeOut;
 719   int iOffset  = i + sizeAllA - sizeOut;
 720   bool otherSystemRec = false;
 721   bool allowInitial   = (partonSystemsPtr->hasInAB(iSys)) ? true : false;
 722   // PS dec 2010: possibility to allow for several recoilers and each with
 723   // flexible normalization
 724   bool   isFlexible   = false;
 725   double flexFactor   = 1.0;
 726   vector<int> iRecVec(0);
 727
 728   // Colour: other end by same index in beam or opposite in final state.
 729   // Exclude rescattered incoming and not final outgoing.
 730   if (colSign > 0)
 731   for (int j = 0; j < sizeAll; ++j) if (j + sizeInA != iOffset) {
 732     int iRecNow = partonSystemsPtr->getAll(iSys, j + sizeInA);
 733     if ( ( j <  sizeIn && event[iRecNow].col()  == colTag
 734       && !event[iRecNow].isRescatteredIncoming() )
 735       || ( j >= sizeIn && event[iRecNow].acol() == colTag
 736       && event[iRecNow].isFinal() ) ) {
 737       iRec = iRecNow;
 738       break;
 739     }
 740   }
 741
 742   // Anticolour: other end by same index in beam or opposite in final state.
 743   // Exclude rescattered incoming and not final outgoing.
 744   if (colSign < 0)
 745   for (int j = 0; j < sizeAll; ++j) if (j + sizeInA != iOffset) {
 746     int iRecNow = partonSystemsPtr->getAll(iSys, j + sizeInA);
 747     if ( ( j <  sizeIn && event[iRecNow].acol()  == colTag
 748       && !event[iRecNow].isRescatteredIncoming() )
 749       || ( j >= sizeIn && event[iRecNow].col() == colTag
 750       && event[iRecNow].isFinal() ) ) {
 751       iRec = iRecNow;
 752       break;
 753     }
 754   }
 755
 756   // Resonance decays (= no instate):
 757   // other end to nearest recoiler in same system final state,
 758   // by (p_i + p_j)^2 - (m_i + m_j)^2 = 2 (p_i p_j - m_i m_j).
 759   // (junction colours more involved, so keep track if junction colour)
 760   bool hasJunction = false;
 761   if (iRec == 0 && !allowInitial) {
 762     for (int iJun = 0; iJun < event.sizeJunction(); ++ iJun) {
 763       // For types 1&2, all legs in final state
 764       // For types 3&4, two legs in final state
 765       // For types 5&6, one leg in final state
 766       int iBeg = (event.kindJunction(iJun)-1)/2;
 767       for (int iLeg = iBeg; iLeg < 3; ++iLeg)
 768         if (event.endColJunction( iJun, iLeg) == colTag) hasJunction  = true;
 769     }
 770     double ppMin = LARGEM2;
 771     for (int j = 0; j < sizeOut; ++j) if (j != i) {
 772         int iRecNow  = partonSystemsPtr->getOut(iSys, j);
 773         if (!event[iRecNow].isFinal()) continue;
 774         double ppNow = event[iRecNow].p() * event[iRad].p()
 775           - event[iRecNow].m() * event[iRad].m();
 776         if (ppNow < ppMin) {
 777           iRec  = iRecNow;
 778           ppMin = ppNow;
 779         }
 780       }
 781   }
 782
 783   // If no success then look for matching (anti)colour anywhere in final state.
 784   if ( iRec == 0 || (!doInterleave && !event[iRec].isFinal()) ) {
 785     iRec = 0;
 786     for (int j = 0; j < event.size(); ++j) if (event[j].isFinal())
 787     if ( (colSign > 0 && event[j].acol() == colTag)
 788       || (colSign < 0 && event[j].col()  == colTag) ) {
 789       iRec = j;
 790       otherSystemRec = true;
 791       break;
 792     }
 793
 794     // If no success then look for match to non-rescattered in initial state.
 795     if (iRec == 0 && allowInitial) {
 796       for (int iSysR = 0; iSysR < partonSystemsPtr->sizeSys(); ++iSysR)
 797       if (iSysR != iSys) {
 798         int j = partonSystemsPtr->getInA(iSysR);
 799         if (j > 0 && event[j].isRescatteredIncoming()) j = 0;
 800         if (j > 0 && ( (colSign > 0 && event[j].col() == colTag)
 801           || (colSign < 0 && event[j].acol()  == colTag) ) ) {
 802           iRec = j;
 803           otherSystemRec = true;
 804           break;
 805         }
 806         j = partonSystemsPtr->getInB(iSysR);
 807         if (j > 0 && event[j].isRescatteredIncoming()) j = 0;
 808         if (j > 0 && ( (colSign > 0 && event[j].col() == colTag)
 809           || (colSign < 0 && event[j].acol()  == colTag) ) ) {
 810           iRec = j;
 811           otherSystemRec = true;
 812           break;
 813         }
 814       }
 815     }
 816   }
 817
 818   // Junctions (PS&ND dec 2010)
 819   // For types 1&2: all legs in final state
 820   //                half-strength dipoles between all legs
 821   // For types 3&4, two legs in final state
 822   //                full-strength dipole between final-state legs
 823   // For types 5&6, one leg in final state
 824   //                no final-state dipole end
 825
 826   if (hasJunction) {
 827     for (int iJun = 0; iJun < event.sizeJunction(); ++ iJun) {
 828       int kindJun = event.kindJunction(iJun);
 829       int iBeg = (kindJun-1)/2;
 830       for (int iLeg = iBeg; iLeg < 3; ++iLeg) {
 831         if (event.endColJunction( iJun, iLeg) == colTag) {
 832           // For types 5&6, no other leg to recoil against. Switch off if
 833           // no other particles at all, since radiation then handled by ISR.
 834           // Example: qq -> ~t* : no radiation off ~t*
 835           // Allow radiation + recoil if unconnected partners available
 836           // Example: qq -> ~t* -> tbar ~chi0 : allow radiation off tbar,
 837           //                                    with ~chi0 as recoiler
 838           if (kindJun >= 5) {
 839             if (sizeOut == 1) return;
 840             else break;
 841           }
 842           // For junction types 3 & 4, span one full-strength dipole
 843           // (only look inside same decay system)
 844           else if (kindJun >= 3) {
 845             int iLegRec = 3-iLeg;
 846             int colTagRec = event.endColJunction( iJun, iLegRec);
 847             for (int j = 0; j < sizeOut; ++j) if (j != i) {
 848                 int iRecNow  = partonSystemsPtr->getOut(iSys, j);
 849                 if (!event[iRecNow].isFinal()) continue;
 850                 if ( (colSign > 0 && event[iRecNow].col()  == colTagRec)
 851                   || (colSign < 0 && event[iRecNow].acol() == colTagRec) ) {
 852                   // Only accept if staying inside same system
 853                   iRec = iRecNow;
 854                   break;
 855                 }
 856               }
 857           }
 858           // For junction types 1 & 2, span two half-strength dipoles
 859           // (only look inside same decay system)
 860           else {
 861             // Loop over two half-strength dipole connections
 862             for (int jLeg = 1; jLeg <= 2; jLeg++) {
 863               int iLegRec = (iLeg + jLeg) % 3;
 864               int colTagRec = event.endColJunction( iJun, iLegRec);
 865               for (int j = 0; j < sizeOut; ++j) if (j != i) {
 866                   int iRecNow  = partonSystemsPtr->getOut(iSys, j);
 867                   if (!event[iRecNow].isFinal()) continue;
 868                   if ( (colSign > 0 && event[iRecNow].col()  == colTagRec)
 869                     || (colSign < 0 && event[iRecNow].acol() == colTagRec) ) {
 870                     // Store recoilers in temporary array
 871                     iRecVec.push_back(iRecNow);
 872                     // Set iRec != 0 for checks below
 873                     iRec = iRecNow;
 874                   }
 875                 }
 876             }
 877
 878           }     // End if-then-else of junction kinds
 879
 880         }       // End if leg has right color tag
 881       }         // End of loop over junction legs
 882     }           // End loop over junctions
 883
 884   }             // End main junction if
 885
 886   // If fail, then other end to nearest recoiler in same system final state,
 887   // by (p_i + p_j)^2 - (m_i + m_j)^2 = 2 (p_i p_j - m_i m_j).
 888   if (iRec == 0) {
 889     double ppMin = LARGEM2;
 890     for (int j = 0; j < sizeOut; ++j) if (j != i) {
 891       int iRecNow  = partonSystemsPtr->getOut(iSys, j);
 892       if (!event[iRecNow].isFinal()) continue;
 893       double ppNow = event[iRecNow].p() * event[iRad].p()
 894                    - event[iRecNow].m() * event[iRad].m();
 895       if (ppNow < ppMin) {
 896         iRec  = iRecNow;
 897         ppMin = ppNow;
 898       }
 899     }
 900   }
 901
 902   // If fail, then other end to nearest recoiler in any system final state,
 903   // by (p_i + p_j)^2 - (m_i + m_j)^2 = 2 (p_i p_j - m_i m_j).
 904   if (iRec == 0) {
 905     double ppMin = LARGEM2;
 906     for (int iRecNow = 0; iRecNow < event.size(); ++iRecNow)
 907     if (iRecNow != iRad && event[iRecNow].isFinal()) {
 908       double ppNow = event[iRecNow].p() * event[iRad].p()
 909                    - event[iRecNow].m() * event[iRad].m();
 910       if (ppNow < ppMin) {
 911         iRec  = iRecNow;
 912         otherSystemRec = true;
 913         ppMin = ppNow;
 914       }
 915     }
 916   }
 917
 918   // PS dec 2010: make sure iRec is stored in iRecVec
 919   if (iRecVec.size() == 0 && iRec != 0) iRecVec.push_back(iRec);
 920
 921   // Remove any zero recoilers from normalization
 922   int nRec = iRecVec.size();
 923   for (unsigned int mRec = 0; mRec < iRecVec.size(); ++mRec)
 924     if (iRecVec[mRec] <= 0) nRec--;
 925   if (nRec >= 2) {
 926     isFlexible = true;
 927     flexFactor = 1.0/nRec;
 928   }
 929
 930   // Check for failure to locate any recoiler
 931   if ( nRec <= 0 ) {
 932     infoPtr->errorMsg("Error in TimeShower::setupQCDdip: "
 933                       "failed to locate any recoiling partner");
 934     return;
 935   }
 936
 937   // Store dipole colour end(s).
 938   for (unsigned int mRec = 0; mRec < iRecVec.size(); ++mRec) {
 939     iRec = iRecVec[mRec];
 940     if (iRec <= 0) continue;
 941     // Max scale either by parton scale or by half dipole mass.
 942     double pTmax = event[iRad].scale();
 943     if (limitPTmaxIn) {
 944       if (iSys == 0) pTmax *= pTmaxFudge;
 945       if (iSys > 0 && sizeIn > 0) pTmax *= pTmaxFudgeMPI;
 946     } else pTmax = 0.5 * m( event[iRad], event[iRec]);
 947     int colType  = (event[iRad].id() == 21) ? 2 * colSign : colSign;
 948     int isrType  = (event[iRec].isFinal()) ? 0 : event[iRec].mother1();
 949     // This line in case mother is a rescattered parton.
 950     while (isrType > 2 + beamOffset) isrType = event[isrType].mother1();
 951     if (isrType > 2) isrType -= beamOffset;
 952     dipEnd.push_back( TimeDipoleEnd( iRad, iRec, pTmax,
 953       colType, 0, 0, isrType, iSys, -1, -1, isOctetOnium) );
 954
 955     // If hooked up with other system then find which.
 956     if (otherSystemRec) {
 957       int systemRec = partonSystemsPtr->getSystemOf(iRec, true);
 958       if (systemRec >= 0) dipEnd.back().systemRec = systemRec;
 959       dipEnd.back().MEtype = 0;
 960     }
 961
 962     // PS dec 2010
 963     // If non-unity (flexible) normalization, set normalization factor
 964     if (isFlexible) {
 965       dipEnd.back().isFlexible = true;
 966       dipEnd.back().flexFactor = flexFactor;
 967     }
 968   }
 969
 970 }
 971
 972 //--------------------------------------------------------------------------
 973
 974 // Setup a dipole end for a QED colour charge or a photon.
 975 // No failsafe choice of recoiler, so gradually widen search.
 976
 977 void TimeShower::setupQEDdip( int iSys, int i, int chgType, int gamType,
 978   Event& event, bool limitPTmaxIn) {
 979
 980   // Initial values. Find if allowed to hook up beams.
 981   int iRad     = partonSystemsPtr->getOut(iSys, i);
 982   int idRad    = event[iRad].id();
 983   int iRec     = 0;
 984   int sizeAllA = partonSystemsPtr->sizeAll(iSys);
 985   int sizeOut  = partonSystemsPtr->sizeOut(iSys);
 986   int sizeAll  = ( allowBeamRecoil ) ? sizeAllA : sizeOut;
 987   int sizeIn   = sizeAll - sizeOut;
 988   int sizeInA  = sizeAllA - sizeIn - sizeOut;
 989   int iOffset  = i + sizeAllA - sizeOut;
 990   double ppMin = LARGEM2;
 991   bool hasRescattered = false;
 992   bool otherSystemRec = false;
 993
 994   // Find nearest same- (opposide-) flavour recoiler in initial (final)
 995   // state of same system, excluding rescattered (in or out) partons.
 996   // Also find if system is involved in rescattering.
 997   // Note: (p_i + p_j)2 - (m_i + m_j)2 = 2 (p_i p_j - m_i m_j).
 998   for (int j = 0; j < sizeAll; ++j) if (j + sizeInA != iOffset) {
 999     int iRecNow  = partonSystemsPtr->getAll(iSys, j + sizeInA);
1000     if ( (j <  sizeIn && !event[iRecNow].isRescatteredIncoming())
1001       || (j >= sizeIn && event[iRecNow].isFinal()) ) {
1002       if ( (j <  sizeIn && event[iRecNow].id() ==  idRad)
1003         || (j >= sizeIn && event[iRecNow].id() == -idRad) ) {
1004         double ppNow = event[iRecNow].p() * event[iRad].p()
1005                      - event[iRecNow].m() * event[iRad].m();
1006         if (ppNow < ppMin) {
1007           iRec  = iRecNow;
1008           ppMin = ppNow;
1009         }
1010       }
1011     } else hasRescattered = true;
1012   }
1013
1014   // If rescattering then find nearest opposite-flavour recoiler
1015   // anywhere in final state.
1016   if (iRec == 0 && hasRescattered) {
1017     for (int iRecNow = 0; iRecNow < event.size(); ++iRecNow)
1018     if (event[iRecNow].id() == -idRad && event[iRecNow].isFinal()) {
1019       double ppNow = event[iRecNow].p() * event[iRad].p()
1020                    - event[iRecNow].m() * event[iRad].m();
1021       if (ppNow < ppMin) {
1022         iRec  = iRecNow;
1023         ppMin = ppNow;
1024         otherSystemRec = true;
1025       }
1026     }
1027   }
1028
1029   // Find nearest recoiler in same system, charge-squared-weighted,
1030   // including initial state, but excluding rescatterer.
1031   if (iRec == 0)
1032   for (int j = 0; j < sizeAll; ++j) if (j + sizeInA != iOffset) {
1033     int iRecNow = partonSystemsPtr->getAll(iSys, j + sizeInA);
1034     int chgTypeRecNow = event[iRecNow].chargeType();
1035     if (chgTypeRecNow == 0) continue;
1036     if ( (j <  sizeIn && !event[iRecNow].isRescatteredIncoming())
1037       || (j >= sizeIn && event[iRecNow].isFinal()) ) {
1038       double ppNow = (event[iRecNow].p() * event[iRad].p()
1039                    -  event[iRecNow].m() * event[iRad].m())
1040                    / pow2(chgTypeRecNow);
1041       if (ppNow < ppMin) {
1042         iRec  = iRecNow;
1043         ppMin = ppNow;
1044       }
1045     }
1046   }
1047
1048   // If rescattering then find nearest recoiler in the final state,
1049   // charge-squared-weighted.
1050   if (iRec == 0 && hasRescattered) {
1051     for (int iRecNow = 0; iRecNow < event.size(); ++iRecNow)
1052     if (iRecNow != iRad && event[iRecNow].isFinal()) {
1053       int chgTypeRecNow = event[iRecNow].chargeType();
1054       if (chgTypeRecNow != 0 && event[iRecNow].isFinal()) {
1055         double ppNow = (event[iRecNow].p() * event[iRad].p()
1056                      -  event[iRecNow].m() * event[iRad].m())
1057                      / pow2(chgTypeRecNow);
1058         if (ppNow < ppMin) {
1059           iRec  = iRecNow;
1060           ppMin = ppNow;
1061           otherSystemRec = true;
1062         }
1063       }
1064     }
1065   }
1066
1067   // Find any nearest recoiler in final state of same system.
1068   if (iRec == 0)
1069   for (int j = 0; j < sizeOut; ++j) if (j != i) {
1070     int iRecNow  = partonSystemsPtr->getOut(iSys, j);
1071     double ppNow = event[iRecNow].p() * event[iRad].p()
1072                  - event[iRecNow].m() * event[iRad].m();
1073     if (ppNow < ppMin) {
1074       iRec  = iRecNow;
1075       ppMin = ppNow;
1076     }
1077   }
1078
1079   // Find any nearest recoiler in final state.
1080   if (iRec == 0)
1081   for (int iRecNow = 0; iRecNow < event.size(); ++iRecNow)
1082   if (iRecNow != iRad && event[iRecNow].isFinal()) {
1083     double ppNow = event[iRecNow].p() * event[iRad].p()
1084                  - event[iRecNow].m() * event[iRad].m();
1085     if (ppNow < ppMin) {
1086       iRec  = iRecNow;
1087       ppMin = ppNow;
1088       otherSystemRec = true;
1089     }
1090   }
1091
1092   // Fill charge-dipole or photon-dipole end.
1093   if (iRec > 0) {
1094     // Max scale either by parton scale or by half dipole mass.
1095     double pTmax = event[iRad].scale();
1096     if (limitPTmaxIn) {
1097       if (iSys == 0) pTmax *= pTmaxFudge;
1098       if (iSys > 0 && sizeIn > 0) pTmax *= pTmaxFudgeMPI;
1099     } else pTmax = 0.5 * m( event[iRad], event[iRec]);
1100     int isrType = (event[iRec].isFinal()) ? 0 : event[iRec].mother1();
1101     // This line in case mother is a rescattered parton.
1102     while (isrType > 2 + beamOffset) isrType = event[isrType].mother1();
1103     if (isrType > 2) isrType -= beamOffset;
1104     dipEnd.push_back( TimeDipoleEnd(iRad, iRec, pTmax,
1105       0, chgType, gamType, isrType, iSys, -1) );
1106
1107     // If hooked up with other system then find which.
1108     if (otherSystemRec) {
1109       int systemRec = partonSystemsPtr->getSystemOf(iRec);
1110       if (systemRec >= 0) dipEnd.back().systemRec = systemRec;
1111       dipEnd.back().MEtype = 0;
1112     }
1113
1114   // Failure to find other end of dipole.
1115   } else {
1116     infoPtr->errorMsg("Error in TimeShower::setupQEDdip: "
1117       "failed to locate any recoiling partner");
1118   }
1119
1120 }
1121
1122 //--------------------------------------------------------------------------
1123
1124 // Setup a dipole end for a Hidden Valley colour charge.
1125
1126 void TimeShower::setupHVdip( int iSys, int i, Event& event,
1127   bool limitPTmaxIn) {
1128
1129   // Initial values.
1130   int iRad    = partonSystemsPtr->getOut(iSys, i);
1131   int iRec    = 0;
1132   int idRad   = event[iRad].id();
1133   int sizeOut = partonSystemsPtr->sizeOut(iSys);
1134
1135   // Hidden Valley colour positive for positive id, and vice versa.
1136   // Find opposte HV colour in final state of same system.
1137   for (int j = 0; j < sizeOut; ++j) if (j != i) {
1138     int iRecNow = partonSystemsPtr->getOut(iSys, j);
1139     int idRec   = event[iRecNow].id();
1140     if ( (abs(idRec) > 4900000 && abs(idRec) < 4900017)
1141       && idRad * idRec < 0) {
1142       iRec = iRecNow;
1143       break;
1144     }
1145   }
1146
1147   // Else find heaviest other final-state in same system.
1148   // (Intended for decays; should mainly be two-body so unique.)
1149   double mMax = -sqrt(LARGEM2);
1150    if (iRec == 0)
1151   for (int j = 0; j < sizeOut; ++j) if (j != i) {
1152     int iRecNow = partonSystemsPtr->getOut(iSys, j);
1153     if (event[iRecNow].m() > mMax) {
1154       iRec = iRecNow;
1155       mMax = event[iRecNow].m();
1156     }
1157   }
1158
1159   // Set up dipole end, or report failure.
1160   if (iRec > 0) {
1161     // Max scale either by parton scale or by half dipole mass.
1162     double pTmax = event[iRad].scale();
1163     if (limitPTmaxIn) {
1164       if (iSys == 0) pTmax *= pTmaxFudge;
1165     } else pTmax = 0.5 * m( event[iRad], event[iRec]);
1166     int colvType  = (event[iRad].id() > 0) ? 1 : -1;
1167     dipEnd.push_back( TimeDipoleEnd( iRad, iRec, pTmax, 0, 0, 0, 0,
1168       iSys, -1, -1, false, true, colvType) );
1169   } else infoPtr->errorMsg("Error in TimeShower::setupHVdip: "
1170       "failed to locate any recoiling partner");
1171
1172 }
1173
1174 //--------------------------------------------------------------------------
1175
1176 // Select next pT in downwards evolution of the existing dipoles.
1177
1178 double TimeShower::pTnext( Event& event, double pTbegAll, double pTendAll) {
1179
1180   // Begin loop over all possible radiating dipole ends.
1181   dipSel  = 0;
1182   iDipSel = -1;
1183   double pT2sel = pTendAll * pTendAll;
1184   for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip) {
1185     TimeDipoleEnd& dip = dipEnd[iDip];
1186
1187     // Check if global recoil should be used.
1188     useLocalRecoilNow = !(globalRecoil && dip.system == 0
1189       && partonSystemsPtr->sizeOut(0) <= nMaxGlobalRecoil);
1190
1191     // Dipole properties; normal local recoil.
1192     dip.mRad   = event[dip.iRadiator].m();
1193     if (useLocalRecoilNow) {
1194       dip.mRec = event[dip.iRecoiler].m();
1195       dip.mDip = m( event[dip.iRadiator], event[dip.iRecoiler] );
1196
1197     // Dipole properties, alternative global recoil. Squares.
1198     } else {
1199       Vec4 pSumGlobal;
1200       for (int i = 0; i < partonSystemsPtr->sizeOut( dip.system); ++i) {
1201         int ii = partonSystemsPtr->getOut( dip.system, i);
1202         if (ii !=  dip.iRadiator) pSumGlobal += event[ii].p();
1203       }
1204       dip.mRec = pSumGlobal.mCalc();
1205       dip.mDip = m( event[dip.iRadiator].p(), pSumGlobal);
1206     }
1207     dip.m2Rad  = pow2(dip.mRad);
1208     dip.m2Rec  = pow2(dip.mRec);
1209     dip.m2Dip  = pow2(dip.mDip);
1210
1211     // Find maximum evolution scale for dipole.
1212     dip.m2DipCorr    = pow2(dip.mDip - dip.mRec) - dip.m2Rad;
1213     double pTbegDip  = min( pTbegAll, dip.pTmax );
1214     double pT2begDip = min( pow2(pTbegDip), 0.25 * dip.m2DipCorr);
1215
1216     // Do QCD, QED or HV evolution if it makes sense.
1217     dip.pT2 = 0.;
1218     if (pT2begDip > pT2sel) {
1219       if      (dip.colType != 0)
1220         pT2nextQCD(pT2begDip, pT2sel, dip, event);
1221       else if (dip.chgType != 0 || dip.gamType != 0)
1222         pT2nextQED(pT2begDip, pT2sel, dip, event);
1223       else if (dip.colvType != 0)
1224         pT2nextHV(pT2begDip, pT2sel, dip, event);
1225
1226       // Update if found larger pT than current maximum.
1227       if (dip.pT2 > pT2sel) {
1228         pT2sel  = dip.pT2;
1229         dipSel  = &dip;
1230         iDipSel = iDip;
1231       }
1232     }
1233   }
1234
1235   // Return nonvanishing value if found pT bigger than already found.
1236   return (dipSel == 0) ? 0. : sqrt(pT2sel);
1237
1238 }
1239
1240 //--------------------------------------------------------------------------
1241
1242 // Evolve a QCD dipole end.
1243
1244 void TimeShower::pT2nextQCD(double pT2begDip, double pT2sel,
1245   TimeDipoleEnd& dip, Event& event) {
1246
1247   // Lower cut for evolution. Return if no evolution range.
1248   double pT2endDip = max( pT2sel, pT2colCut );
1249   if (pT2begDip < pT2endDip) return;
1250
1251   // Upper estimate for matrix element weighting and colour factor.
1252   // Special cases for triplet recoiling against gluino and octet onia.
1253   // Note that g -> g g and g -> q qbar are split on two sides.
1254   int    colTypeAbs = abs(dip.colType);
1255   double wtPSglue   = 2.;
1256   double colFac     = (colTypeAbs == 1) ? 4./3. : 3./2.;
1257   if (dip.MEgluinoRec)  colFac  = 3.;
1258   if (dip.isOctetOnium) colFac *= 0.5 * octetOniumColFac;
1259   // PS dec 2010. Include possibility for flexible normalization,
1260   // e.g., for dipoles stretched to junctions or to switch off radiation.
1261   if (dip.isFlexible)   colFac *= dip.flexFactor;
1262   double wtPSqqbar  = (colTypeAbs == 2) ? 0.25 * nGluonToQuark : 0.;
1263
1264   // Variables used inside evolution loop. (Mainly dummy start values.)
1265   dip.pT2              = pT2begDip;
1266   int    nFlavour      = 3;
1267   double zMinAbs       = 0.5;
1268   double pT2min        = pT2endDip;
1269   double b0            = 4.5;
1270   double Lambda2       = Lambda3flav2;
1271   double emitCoefGlue  = 0.;
1272   double emitCoefQqbar = 0.;
1273   double emitCoefTot   = 0.;
1274   double wt            = 0.;
1275   bool   mustFindRange = true;
1276
1277   // Begin evolution loop towards smaller pT values.
1278   do {
1279
1280     // Initialize evolution coefficients at the beginning and
1281     // reinitialize when crossing c and b flavour thresholds.
1282     if (mustFindRange) {
1283
1284       // Determine overestimated z range; switch at c and b masses.
1285       if (dip.pT2 > m2b) {
1286         nFlavour = 5;
1287         pT2min   = max( m2b, pT2endDip);
1288         b0       = 23./6.;
1289         Lambda2  = Lambda5flav2;
1290       } else if (dip.pT2 > m2c) {
1291         nFlavour = 4;
1292         pT2min   = max( m2c, pT2endDip);
1293         b0       = 25./6.;
1294         Lambda2  = Lambda4flav2;
1295       } else {
1296         nFlavour = 3;
1297         pT2min   = pT2endDip;
1298         b0       = 27./6.;
1299         Lambda2  = Lambda3flav2;
1300       }
1301       // A change of renormalization scale expressed by a change of Lambda.
1302       Lambda2 /= renormMultFac;
1303       zMinAbs = 0.5 - sqrtpos( 0.25 - pT2min / dip.m2DipCorr );
1304       if (zMinAbs < SIMPLIFYROOT) zMinAbs = pT2min / dip.m2DipCorr;
1305
1306       // Find emission coefficients for X -> X g and g -> q qbar.
1307       emitCoefGlue = wtPSglue * colFac * log(1. / zMinAbs - 1.);
1308       emitCoefTot  = emitCoefGlue;
1309       if (colTypeAbs == 2 && event[dip.iRadiator].id() == 21) {
1310         emitCoefQqbar = wtPSqqbar * (1. - 2. * zMinAbs);
1311         emitCoefTot  += emitCoefQqbar;
1312       }
1313
1314       // Initialization done for current range.
1315       mustFindRange = false;
1316     }
1317
1318     // Pick pT2 (in overestimated z range) for fixed alpha_strong.
1319     if (alphaSorder == 0) {
1320       dip.pT2 = dip.pT2 * pow( rndmPtr->flat(),
1321         1. / (alphaS2pi * emitCoefTot) );
1322
1323     // Ditto for first-order alpha_strong.
1324     } else if (alphaSorder == 1) {
1325       dip.pT2 = Lambda2 * pow( dip.pT2 / Lambda2,
1326         pow( rndmPtr->flat(), b0 / emitCoefTot) );
1327
1328       // For second order reject by second term in alpha_strong expression.
1329     } else {
1330       do dip.pT2 = Lambda2 * pow( dip.pT2 / Lambda2,
1331         pow( rndmPtr->flat(), b0 / emitCoefTot) );
1332       while (alphaS.alphaS2OrdCorr(renormMultFac * dip.pT2) < rndmPtr->flat()
1333         && dip.pT2 > pT2min);
1334     }
1335     wt = 0.;
1336
1337     // If crossed c or b thresholds: continue evolution from threshold.
1338     if (nFlavour == 5 && dip.pT2 < m2b) {
1339       mustFindRange = true;
1340       dip.pT2       = m2b;
1341     } else if ( nFlavour == 4 && dip.pT2 < m2c) {
1342       mustFindRange = true;
1343       dip.pT2       = m2c;
1344
1345     // Abort evolution if below cutoff scale, or below another branching.
1346     } else {
1347       if ( dip.pT2 < pT2endDip) { dip.pT2 = 0.; return; }
1348
1349       // Pick kind of branching: X -> X g or g -> q qbar.
1350       dip.flavour  = 21;
1351       dip.mFlavour = 0.;
1352       if (colTypeAbs == 2 && emitCoefQqbar > rndmPtr->flat()
1353         * emitCoefTot) dip.flavour = 0;
1354
1355       // Pick z: either dz/(1-z) or flat dz.
1356       if (dip.flavour == 21) {
1357         dip.z = 1. - zMinAbs * pow( 1. / zMinAbs - 1., rndmPtr->flat() );
1358       } else {
1359         dip.z = zMinAbs + (1. - 2. * zMinAbs) * rndmPtr->flat();
1360       }
1361
1362       // Do not accept branching if outside allowed z range.
1363       double zMin = 0.5 - sqrtpos( 0.25 - dip.pT2 / dip.m2DipCorr );
1364       if (zMin < SIMPLIFYROOT) zMin = dip.pT2 / dip.m2DipCorr;
1365       dip.m2 = dip.m2Rad + dip.pT2 / (dip.z * (1. - dip.z));
1366       if (dip.z > zMin && dip.z < 1. - zMin
1367         && dip.m2 * dip.m2Dip < dip.z * (1. - dip.z)
1368         * pow2(dip.m2Dip + dip.m2 - dip.m2Rec) ) {
1369
1370         // Flavour choice for g -> q qbar.
1371         if (dip.flavour == 0) {
1372           dip.flavour  = min(5, 1 + int(nGluonToQuark * rndmPtr->flat()));
1373           dip.mFlavour = particleDataPtr->m0(dip.flavour);
1374         }
1375
1376         // No z weight, except threshold, if to do ME corrections later on.
1377         if (dip.MEtype > 0) {
1378           wt = 1.;
1379           if (dip.flavour < 10 && dip.m2 < THRESHM2 * pow2(dip.mFlavour))
1380             wt = 0.;
1381
1382         // z weight for X -> X g.
1383         } else if (dip.flavour == 21
1384           && (colTypeAbs == 1 || colTypeAbs == 3) ) {
1385           wt = (1. + pow2(dip.z)) / wtPSglue;
1386         } else if (dip.flavour == 21) {
1387           wt = (1. + pow3(dip.z)) / wtPSglue;
1388
1389         // z weight for g -> q qbar.
1390         } else {
1391           double beta  = sqrtpos( 1. - 4. * pow2(dip.mFlavour) / dip.m2 );
1392           wt = beta * ( pow2(dip.z) + pow2(1. - dip.z) );
1393         }
1394
1395         // Suppression factors for dipole to beam remnant.
1396         if (dip.isrType != 0 && useLocalRecoilNow) {
1397          BeamParticle& beam = (dip.isrType == 1) ? *beamAPtr : *beamBPtr;
1398           int iSysRec = dip.systemRec;
1399           double xOld = beam[iSysRec].x();
1400           double xNew = xOld * (1. + (dip.m2 - dip.m2Rad) /
1401             (dip.m2Dip - dip.m2Rad));
1402           double xMaxAbs = beam.xMax(iSysRec);
1403           if (xMaxAbs < 0.) {
1404             infoPtr->errorMsg("Warning in TimeShower::pT2nextQCD: "
1405             "xMaxAbs negative");
1406             return;
1407           }
1408
1409           // Firstly reduce by PDF ratio.
1410           if (xNew > xMaxAbs) wt = 0.;
1411           else {
1412             int idRec     = event[dip.iRecoiler].id();
1413             double pdfOld = max ( TINYPDF,
1414               beam.xfISR( iSysRec, idRec, xOld, factorMultFac * dip.pT2) );
1415             double pdfNew =
1416               beam.xfISR( iSysRec, idRec, xNew, factorMultFac * dip.pT2);
1417             wt *= min( 1., pdfNew / pdfOld);
1418           }
1419
1420           // Secondly optionally reduce by 4 pT2_hard / (4 pT2_hard + m2).
1421           if (dampenBeamRecoil) {
1422             double pTpT = sqrt(event[dip.iRadiator].pT2() * dip.pT2);
1423             wt *= pTpT / (pTpT + dip.m2);
1424           }
1425         }
1426       }
1427     }
1428
1429   // Iterate until acceptable pT (or have fallen below pTmin).
1430   } while (wt < rndmPtr->flat());
1431
1432 }
1433
1434 //--------------------------------------------------------------------------
1435
1436 // Evolve a QED dipole end, either charged or photon.
1437
1438 void TimeShower::pT2nextQED(double pT2begDip, double pT2sel,
1439   TimeDipoleEnd& dip, Event& event) {
1440
1441   // Lower cut for evolution. Return if no evolution range.
1442   double pT2chgCut = (dip.chgType != 0 && abs(dip.chgType) != 3)
1443     ? pT2chgQCut : pT2chgLCut;
1444   double pT2endDip = max( pT2sel, pT2chgCut );
1445   if (pT2begDip < pT2endDip) return;
1446
1447   // Emission of photon or photon branching.
1448   bool hasCharge = (dip.chgType != 0);
1449
1450   // Default values.
1451   double wtPSgam     = 0.;
1452   double chg2Sum     = 0.;
1453   double chg2SumL    = 0.;
1454   double chg2SumQ    = 0.;
1455   double zMinAbs     = 0.;
1456   double emitCoefTot = 0.;
1457
1458   // alpha_em at maximum scale provides upper estimate.
1459   double alphaEMmax  = alphaEM.alphaEM(renormMultFac * pT2begDip);
1460   double alphaEM2pi  = alphaEMmax / (2. * M_PI);
1461
1462   // Emission: upper estimate for matrix element weighting; charge factor.
1463   if (hasCharge) {
1464     wtPSgam     = 2.;
1465     double chg2 = pow2(dip.chgType / 3.);
1466
1467     // Determine overestimated z range. Find evolution coefficient.
1468     zMinAbs = 0.5 - sqrtpos( 0.25 - pT2endDip / dip.m2DipCorr );
1469     if (zMinAbs < SIMPLIFYROOT) zMinAbs = pT2endDip / dip.m2DipCorr;
1470     emitCoefTot = alphaEM2pi * chg2 * wtPSgam * log(1. / zMinAbs - 1.);
1471
1472   // Branching: sum of squared charge factors for lepton and quark daughters.
1473   } else {
1474     chg2SumL = max(0, min(3, nGammaToLepton));
1475     if      (nGammaToQuark > 4) chg2SumQ = 11. / 9.;
1476     else if (nGammaToQuark > 3) chg2SumQ = 10. / 9.;
1477     else if (nGammaToQuark > 2) chg2SumQ =  6. / 9.;
1478     else if (nGammaToQuark > 1) chg2SumQ =  5. / 9.;
1479     else if (nGammaToQuark > 0) chg2SumQ =  1. / 9.;
1480
1481     // Total sum of squared charge factors. Find evolution coefficient.
1482     chg2Sum     = chg2SumL + 3. * chg2SumQ;
1483     emitCoefTot = alphaEM2pi * chg2Sum;
1484   }
1485
1486   // Variables used inside evolution loop.
1487   dip.pT2 = pT2begDip;
1488   double wt;
1489
1490   // Begin evolution loop towards smaller pT values.
1491   do {
1492
1493     // Pick pT2 (in overestimated z range).
1494     dip.pT2 = dip.pT2 * pow(rndmPtr->flat(), 1. / emitCoefTot);
1495     wt = 0.;
1496
1497     // Abort evolution if below cutoff scale, or below another branching.
1498     if ( dip.pT2 < pT2endDip) { dip.pT2 = 0.; return; }
1499
1500     // Pick z according to dz/(1-z) or flat.
1501     if (hasCharge) dip.z = 1. - zMinAbs
1502       * pow( 1. / zMinAbs - 1., rndmPtr->flat() );
1503     else           dip.z = rndmPtr->flat();
1504
1505     // Do not accept branching if outside allowed z range.
1506     double zMin = 0.5 - sqrtpos( 0.25 - dip.pT2 / dip.m2DipCorr );
1507     if (zMin < SIMPLIFYROOT) zMin = dip.pT2 / dip.m2DipCorr;
1508     dip.m2 = dip.m2Rad + dip.pT2 / (dip.z * (1. - dip.z));
1509     if (dip.z > zMin && dip.z < 1. - zMin
1510       && dip.m2 * dip.m2Dip < dip.z * (1. - dip.z)
1511         * pow2(dip.m2Dip + dip.m2 - dip.m2Rec)
1512       // For gamma -> f fbar also impose maximum mass.
1513       && (hasCharge || dip.m2 < m2MaxGamma) ) {
1514
1515       // Photon emission: unique flavour choice.
1516       if (hasCharge) {
1517         dip.flavour = 22;
1518         dip.mFlavour = 0.;
1519
1520       // Photon branching: either lepton or quark flavour choice.
1521       } else {
1522         if (rndmPtr->flat() * chg2Sum < chg2SumL)
1523           dip.flavour  = 9 + 2 * min(3, 1 + int(chg2SumL * rndmPtr->flat()));
1524         else {
1525           double rndmQ = 9. * chg2SumQ * rndmPtr->flat();
1526           if      (rndmQ <  1.) dip.flavour = 1;
1527           else if (rndmQ <  5.) dip.flavour = 2;
1528           else if (rndmQ <  6.) dip.flavour = 3;
1529           else if (rndmQ < 10.) dip.flavour = 4;
1530           else                  dip.flavour = 5;
1531         }
1532         dip.mFlavour = particleDataPtr->m0(dip.flavour);
1533       }
1534
1535       // No z weight, except threshold, if to do ME corrections later on.
1536       if (dip.MEtype > 0) {
1537         wt = 1.;
1538         if (dip.flavour < 20 && dip.m2 < THRESHM2 * pow2(dip.mFlavour))
1539           wt = 0.;
1540
1541       // z weight for X -> X gamma.
1542       } else if (hasCharge) {
1543         wt = (1. + pow2(dip.z)) / wtPSgam;
1544
1545       // z weight for gamma -> f fbar.
1546       } else {
1547         double beta  = sqrtpos( 1. - 4. * pow2(dip.mFlavour) / dip.m2 );
1548         wt = beta * ( pow2(dip.z) + pow2(1. - dip.z) );
1549       }
1550
1551       // Correct to current value of alpha_EM.
1552       double alphaEMnow = alphaEM.alphaEM(renormMultFac * dip.pT2);
1553       wt *= (alphaEMnow / alphaEMmax);
1554
1555       // Suppression factors for dipole to beam remnant.
1556       if (dip.isrType != 0 && useLocalRecoilNow) {
1557         BeamParticle& beam = (dip.isrType == 1) ? *beamAPtr : *beamBPtr;
1558         int iSys    = dip.system;
1559         double xOld = beam[iSys].x();
1560         double xNew = xOld * (1. + (dip.m2 - dip.m2Rad) /
1561           (dip.m2Dip - dip.m2Rad));
1562         double xMaxAbs = beam.xMax(iSys);
1563         if (xMaxAbs < 0.) {
1564           infoPtr->errorMsg("Warning in TimeShower::pT2nextQED: "
1565           "xMaxAbs negative");
1566           return;
1567         }
1568
1569         // Firstly reduce by PDF ratio.
1570         if (xNew > xMaxAbs) wt = 0.;
1571         else {
1572           int idRec     = event[dip.iRecoiler].id();
1573           double pdfOld = max ( TINYPDF,
1574             beam.xfISR( iSys, idRec, xOld, factorMultFac * dip.pT2) );
1575           double pdfNew =
1576             beam.xfISR( iSys, idRec, xNew, factorMultFac * dip.pT2);
1577           wt *= min( 1., pdfNew / pdfOld);
1578         }
1579
1580         // Secondly optionally reduce by 4 pT2_hard / (4 pT2_hard + m2).
1581         if (dampenBeamRecoil) {
1582           double pT24 = 4. * event[dip.iRadiator].pT2();
1583           wt *= pT24 / (pT24 + dip.m2);
1584         }
1585       }
1586     }
1587
1588   // Iterate until acceptable pT (or have fallen below pTmin).
1589   } while (wt < rndmPtr->flat());
1590
1591 }
1592
1593 //--------------------------------------------------------------------------
1594
1595 // Evolve a Hidden Valley dipole end.
1596
1597 void TimeShower::pT2nextHV(double pT2begDip, double pT2sel,
1598   TimeDipoleEnd& dip, Event& ) {
1599
1600   // Lower cut for evolution. Return if no evolution range.
1601   double pT2endDip = max( pT2sel, pT2hvCut );
1602   if (pT2begDip < pT2endDip) return;
1603
1604   // C_F * alpha_HV/2 pi.
1605   int    colvTypeAbs = abs(dip.colvType);
1606   double colvFac     = (colvTypeAbs == 1) ? CFHV : 0.5 * nCHV;
1607   double alphaHV2pi  = colvFac * (alphaHVfix / (2. * M_PI));
1608
1609   // Determine overestimated z range. Find evolution coefficient.
1610   double zMinAbs = 0.5 - sqrtpos( 0.25 - pT2endDip / dip.m2DipCorr );
1611   if (zMinAbs < SIMPLIFYROOT) zMinAbs = pT2endDip / dip.m2DipCorr;
1612   double emitCoefTot = alphaHV2pi * 2. * log(1. / zMinAbs - 1.);
1613
1614   // Variables used inside evolution loop.
1615   dip.pT2 = pT2begDip;
1616   double wt;
1617
1618   // Begin evolution loop towards smaller pT values.
1619   do {
1620
1621     // Pick pT2 (in overestimated z range).
1622     dip.pT2 = dip.pT2 * pow(rndmPtr->flat(), 1. / emitCoefTot);
1623     wt = 0.;
1624
1625     // Abort evolution if below cutoff scale, or below another branching.
1626     if ( dip.pT2 < pT2endDip) { dip.pT2 = 0.; return; }
1627
1628     // Pick z according to dz/(1-z).
1629     dip.z = 1. - zMinAbs * pow( 1. / zMinAbs - 1., rndmPtr->flat() );
1630
1631     // Do not accept branching if outside allowed z range.
1632     double zMin = 0.5 - sqrtpos( 0.25 - dip.pT2 / dip.m2DipCorr );
1633     if (zMin < SIMPLIFYROOT) zMin = dip.pT2 / dip.m2DipCorr;
1634     dip.m2 = dip.m2Rad + dip.pT2 / (dip.z * (1. - dip.z));
1635     if (dip.z > zMin && dip.z < 1. - zMin
1636       && dip.m2 * dip.m2Dip < dip.z * (1. - dip.z)
1637         * pow2(dip.m2Dip + dip.m2 - dip.m2Rec) ) {
1638
1639       // HV gamma or gluon emission: unique flavour choice.
1640       dip.flavour  = idHV;
1641       dip.mFlavour = mHV;
1642
1643       // No z weight, except threshold, if to do ME corrections later on.
1644       if (dip.MEtype > 0) wt = 1.;
1645
1646       // z weight for X -> X g_HV.
1647       else if (colvTypeAbs == 1) wt = (1. + pow2(dip.z)) / 2.;
1648       else wt = (1. + pow3(dip.z)) / 2.;
1649     }
1650
1651   // Iterate until acceptable pT (or have fallen below pTmin).
1652   } while (wt < rndmPtr->flat());
1653
1654 }
1655
1656 //--------------------------------------------------------------------------
1657
1658 // ME corrections and kinematics that may give failure.
1659 // Notation: radBef, recBef = radiator, recoiler before emission,
1660 //           rad, rec, emt = radiator, recoiler, emitted efter emission.
1661 //           (rad, emt distinguished by colour flow for g -> q qbar.)
1662
1663 bool TimeShower::branch( Event& event, bool isInterleaved) {
1664
1665   // Check if global recoil should be used.
1666   useLocalRecoilNow = !(globalRecoil && dipSel->system == 0
1667     && partonSystemsPtr->sizeOut(0) <= nMaxGlobalRecoil);
1668
1669   // Find initial radiator and recoiler particles in dipole branching.
1670   int iRadBef      = dipSel->iRadiator;
1671   int iRecBef      = dipSel->iRecoiler;
1672   Particle& radBef = event[iRadBef];
1673   Particle& recBef = event[iRecBef];
1674
1675   // Find their momenta, with special sum for global recoil.
1676   Vec4 pRadBef     = event[iRadBef].p();
1677   Vec4 pRecBef;
1678   vector<int> iGRecBef, iGRec;
1679   if (useLocalRecoilNow) pRecBef =  event[iRecBef].p();
1680   else {
1681     for (int i = 0; i < partonSystemsPtr->sizeOut( dipSel->system); ++i) {
1682       int iG = partonSystemsPtr->getOut( dipSel->system, i);
1683       if (iG !=  dipSel->iRadiator) {
1684         iGRecBef.push_back(iG);
1685         pRecBef += event[iG].p();
1686       }
1687     }
1688   }
1689
1690   // Default flavours and colour tags for new particles in dipole branching.
1691   int idRad        = radBef.id();
1692   int idEmt        = dipSel->flavour;
1693   int colRad       = radBef.col();
1694   int acolRad      = radBef.acol();
1695   int colEmt       = 0;
1696   int acolEmt      = 0;
1697   iSysSel          = dipSel->system;
1698   int iSysSelRec   = dipSel->systemRec;
1699
1700   // Default OK for photon, photon_HV or gluon_HV emission.
1701   if (dipSel->flavour == 22 || dipSel->flavour == idHV) {
1702   // New colour tag required for gluon emission.
1703   } else if (dipSel->flavour == 21 && dipSel->colType > 0) {
1704     colEmt  = colRad;
1705     colRad  = event.nextColTag();
1706     acolEmt = colRad;
1707   } else if (dipSel->flavour == 21) {
1708     acolEmt = acolRad;
1709     acolRad = event.nextColTag();
1710     colEmt  = acolRad;
1711   // New flavours for g -> q qbar; split colours.
1712   } else if (dipSel->colType > 0) {
1713     idEmt   = dipSel->flavour ;
1714     idRad   = -idEmt;
1715     colEmt  = colRad;
1716     colRad  = 0;
1717   } else if (dipSel->colType < 0) {
1718     idEmt   = -dipSel->flavour ;
1719     idRad   = -idEmt;
1720     acolEmt = acolRad;
1721     acolRad = 0;
1722   // New flavours for gamma -> f fbar, and maybe also colours.
1723   } else if (dipSel->gamType == 1 && rndmPtr->flat() > 0.5) {
1724     idEmt   = -dipSel->flavour ;
1725     idRad   = -idEmt;
1726     if (idRad < 10) colRad = event.nextColTag();
1727     acolEmt = colRad;
1728   } else if (dipSel->gamType == 1) {
1729     idEmt   = dipSel->flavour ;
1730     idRad   = -idEmt;
1731     if (idEmt < 10) colEmt = event.nextColTag();
1732     acolRad = colEmt;
1733   }
1734
1735   // Construct kinematics in dipole rest frame:
1736   // begin simple (like g -> g g).
1737   double pTorig       = sqrt( dipSel->pT2);
1738   double eRadPlusEmt  = 0.5 * (dipSel->m2Dip + dipSel->m2 - dipSel->m2Rec)
1739     / dipSel->mDip;
1740   double e2RadPlusEmt = pow2(eRadPlusEmt);
1741   double pzRadPlusEmt = 0.5 * sqrtpos( pow2(dipSel->m2Dip - dipSel->m2
1742     - dipSel->m2Rec) - 4. * dipSel->m2 * dipSel->m2Rec ) / dipSel->mDip;
1743   double pT2corr = dipSel->m2 * (e2RadPlusEmt * dipSel->z * (1. - dipSel->z)
1744                       - 0.25 * dipSel->m2) / pow2(pzRadPlusEmt);
1745   double pTcorr       = sqrtpos( pT2corr );
1746   double pzRad        = (e2RadPlusEmt * dipSel->z - 0.5 * dipSel->m2)
1747                       / pzRadPlusEmt;
1748   double pzEmt        = (e2RadPlusEmt * (1. - dipSel->z) - 0.5 * dipSel->m2)
1749                       / pzRadPlusEmt;
1750   double mRad         = dipSel->mRad;
1751   double mEmt         = 0.;
1752
1753   // Kinematics reduction for q -> q gamma_v when m_q > 0 and m_gamma_v > 0.
1754   if ( abs(dipSel->colvType) == 1 && dipSel->mFlavour > 0.) {
1755     mEmt              = dipSel->mFlavour;
1756     if (pow2(mRad + mEmt) > dipSel->m2) return false;
1757     double m2Emt      = pow2(mEmt);
1758     double lambda     = sqrtpos( pow2(dipSel->m2 - dipSel->m2Rad - m2Emt)
1759                       - 4. * dipSel->m2Rad * m2Emt );
1760     kRad              = 0.5 * (dipSel->m2 - lambda + m2Emt - dipSel->m2Rad)
1761                       / dipSel->m2;
1762     kEmt              = 0.5 * (dipSel->m2 - lambda + dipSel->m2Rad - m2Emt)
1763                       / dipSel->m2;
1764     pTorig           *= 1. - kRad - kEmt;
1765     pTcorr           *= 1. - kRad - kEmt;
1766     double pzMove     = kRad * pzRad - kEmt * pzEmt;
1767     pzRad            -= pzMove;
1768     pzEmt            += pzMove;
1769
1770   // Kinematics reduction for q -> q g/gamma/g_HV when m_q > 0.
1771   } else if (abs(dipSel->colType) == 1 || dipSel->chgType != 0
1772     || abs(dipSel->colvType) == 1) {
1773     pTorig           *= 1. - dipSel->m2Rad / dipSel->m2;
1774     pTcorr           *= 1. - dipSel->m2Rad / dipSel->m2;
1775     pzRad            += pzEmt * dipSel->m2Rad / dipSel->m2;
1776     pzEmt            *= 1. - dipSel->m2Rad / dipSel->m2;
1777
1778   // Kinematics reduction for g -> q qbar or gamma -> f fbar when m_f > 0;
1779   } else if (abs(dipSel->flavour) < 20) {
1780     mEmt              = dipSel->mFlavour;
1781     mRad              = mEmt;
1782     double beta       = sqrtpos( 1. - 4. * pow2(mEmt) / dipSel->m2 );
1783     pTorig           *= beta;
1784     pTcorr           *= beta;
1785     pzRad             = 0.5 * ( (1. + beta) * pzRad + (1. - beta) * pzEmt );
1786     pzEmt             = pzRadPlusEmt - pzRad;
1787   }
1788
1789   // Reject g emission where mass effects have reduced pT below cutoff.
1790   if (idEmt == 21 && pTorig < pTcolCut) return false;
1791
1792   // Find rest frame and angles of original dipole.
1793   RotBstMatrix M;
1794   M.fromCMframe(pRadBef, pRecBef);
1795
1796   // Evaluate coefficient of azimuthal asymmetry from gluon polarization.
1797   findAsymPol( event, dipSel);
1798
1799   // Begin construction of new dipole kinematics: pick azimuthal angle.
1800   Vec4 pRad, pEmt, pRec;
1801   double wtPhi = 1.;
1802   do {
1803     double phi = 2. * M_PI * rndmPtr->flat();
1804
1805     // Define kinematics of branching in dipole rest frame.
1806     pRad = Vec4( pTcorr * cos(phi), pTcorr * sin(phi), pzRad,
1807       sqrt( pow2(pTcorr) + pow2(pzRad) + pow2(mRad) ) );
1808     pEmt = Vec4( -pRad.px(), -pRad.py(), pzEmt,
1809       sqrt( pow2(pTcorr) + pow2(pzEmt) + pow2(mEmt) ) );
1810     pRec = Vec4( 0., 0., -pzRadPlusEmt, sqrt( pow2(pzRadPlusEmt)
1811       + dipSel->m2Rec ) );
1812
1813     // Rotate and boost dipole products to the event frame.
1814     pRad.rotbst(M);
1815     pEmt.rotbst(M);
1816     pRec.rotbst(M);
1817
1818     // Azimuthal phi weighting: loop to new phi value if required.
1819     if (dipSel->asymPol != 0.) {
1820       Vec4 pAunt = event[dipSel->iAunt].p();
1821       double cosPhi = cosphi( pRad, pAunt, pRadBef );
1822       wtPhi = ( 1. + dipSel->asymPol * (2. * pow2(cosPhi) - 1.) )
1823         / ( 1. + abs(dipSel->asymPol) );
1824     }
1825   } while (wtPhi < rndmPtr->flat()) ;
1826
1827   // Kinematics when recoiler is initial-state parton.
1828   int isrTypeNow  = dipSel->isrType;
1829   int isrTypeSave = isrTypeNow;
1830   if (!useLocalRecoilNow) isrTypeNow = 0;
1831   if (isrTypeNow != 0) pRec = 2. * recBef.p() - pRec;
1832
1833   // PS dec 2010: check if radiator has flexible normalization
1834   bool isFlexible = dipSel->isFlexible;
1835
1836   // Define new particles from dipole branching.
1837   double pTsel = sqrt(dipSel->pT2);
1838   Particle rad = Particle(idRad, 51, iRadBef, 0, 0, 0,
1839     colRad, acolRad, pRad, mRad, pTsel);
1840   Particle emt = Particle(idEmt, 51, iRadBef, 0, 0, 0,
1841     colEmt, acolEmt, pEmt, mEmt, pTsel);
1842
1843   // Recoiler either in final or in initial state
1844   Particle rec = (isrTypeNow == 0)
1845     ? Particle(recBef.id(),  52, iRecBef, iRecBef, 0, 0,
1846       recBef.col(), recBef.acol(), pRec, dipSel->mRec, pTsel)
1847     : Particle(recBef.id(), -53, 0, 0, iRecBef, iRecBef,
1848       recBef.col(), recBef.acol(), pRec, 0., 0.);
1849
1850   // ME corrections can lead to branching being rejected.
1851   if (dipSel->MEtype > 0) {
1852     Particle& partner = (dipSel->iMEpartner == iRecBef)
1853       ? rec : event[dipSel->iMEpartner];
1854     if ( findMEcorr( dipSel, rad, partner, emt) < rndmPtr->flat() )
1855       return false;
1856   }
1857
1858   // Rescatter: if the recoiling partner is not in the same system
1859   //            as the radiator, fix up intermediate systems (can lead
1860   //            to emissions being vetoed)
1861   if (allowRescatter && FIXRESCATTER && isInterleaved
1862     && iSysSel != iSysSelRec) {
1863     Vec4 pNew = rad.p() + emt.p();
1864     if (!rescatterPropagateRecoil(event, pNew)) return false;
1865   }
1866
1867   // Save properties to be restored in case of user-hook veto of emission.
1868   int eventSizeOld = event.size();
1869   int iRadStatusV  = event[iRadBef].status();
1870   int iRadDau1V    = event[iRadBef].daughter1();
1871   int iRadDau2V    = event[iRadBef].daughter2();
1872   int iRecStatusV  = event[iRecBef].status();
1873   int iRecMot1V    = event[iRecBef].mother1();
1874   int iRecMot2V    = event[iRecBef].mother2();
1875   int iRecDau1V    = event[iRecBef].daughter1();
1876   int iRecDau2V    = event[iRecBef].daughter2();
1877   int beamOff1     = 1 + beamOffset;
1878   int beamOff2     = 2 + beamOffset;
1879   int ev1Dau1V     = event[beamOff1].daughter1();
1880   int ev2Dau1V     = event[beamOff2].daughter1();
1881
1882   // Shower may occur at a displaced vertex.
1883   if (radBef.hasVertex()) {
1884     rad.vProd( radBef.vProd() );
1885     emt.vProd( radBef.vProd() );
1886   }
1887   if (recBef.hasVertex()) rec.vProd( recBef.vProd() );
1888
1889   // Put new particles into the event record.
1890   // Mark original dipole partons as branched and set daughters/mothers.
1891   int iRad = event.append(rad);
1892   int iEmt = event.append(emt);
1893   event[iRadBef].statusNeg();
1894   event[iRadBef].daughters( iRad, iEmt);
1895   int iRec = 0;
1896   if (useLocalRecoilNow) {
1897     iRec = event.append(rec);
1898     if (isrTypeNow == 0) {
1899       event[iRecBef].statusNeg();
1900       event[iRecBef].daughters( iRec, iRec);
1901     } else {
1902       event[iRecBef].mothers( iRec, iRec);
1903       event[iRec].mothers( iRecMot1V, iRecMot2V);
1904       if (iRecMot1V == beamOff1) event[beamOff1].daughter1( iRec);
1905       if (iRecMot1V == beamOff2) event[beamOff2].daughter1( iRec);
1906     }
1907
1908   // Global recoil: need to find relevant rotation+boost for recoilers:
1909   // boost+rotate to rest frame, boost along z axis, rotate+boost back.
1910   } else {
1911     RotBstMatrix MG = M;
1912     MG.invert();
1913     double pzRecBef = -0.5 * sqrtpos( pow2(dipSel->m2Dip - dipSel->m2Rad
1914       - dipSel->m2Rec) - 4. * dipSel->m2Rad * dipSel->m2Rec ) / dipSel->mDip;
1915     double eRecBef  = sqrt( pow2(pzRecBef) + dipSel->m2Rec);
1916     double pzRecAft = -pzRadPlusEmt;
1917     double eRecAft  = sqrt( pow2(pzRecAft) + dipSel->m2Rec);
1918     MG.bst( Vec4(0., 0., pzRecBef, eRecBef), Vec4(0., 0., pzRecAft, eRecAft) );
1919     MG.rotbst( M);
1920
1921     // Global recoil: copy particles, and rotate+boost momenta (not vertices).
1922     for (int iG = 0; iG < int(iGRecBef.size()); ++iG) {
1923       iRec = event.copy( iGRecBef[iG], 52);
1924       iGRec.push_back( iRec);
1925       Vec4 pGRec = event[iRec].p();
1926       pGRec.rotbst( MG);
1927       event[iRec].p( pGRec);
1928     }
1929   }
1930
1931   // Allow veto of branching. If so restore event record to before emission.
1932   bool inResonance = (partonSystemsPtr->getInA(iSysSel) == 0) ? true : false;
1933   if ( canVetoEmission && userHooksPtr->doVetoFSREmission( eventSizeOld,
1934     event, iSysSel, inResonance) ) {
1935     event.popBack( event.size() - eventSizeOld);
1936     event[iRadBef].status( iRadStatusV);
1937     event[iRadBef].daughters( iRadDau1V, iRadDau2V);
1938     if (useLocalRecoilNow && isrTypeNow == 0) {
1939       event[iRecBef].status( iRecStatusV);
1940       event[iRecBef].daughters( iRecDau1V, iRecDau2V);
1941     } else if (useLocalRecoilNow) {
1942       event[iRecBef].mothers( iRecMot1V, iRecMot2V);
1943       if (iRecMot1V == beamOff1) event[beamOff1].daughter1( ev1Dau1V);
1944       if (iRecMot1V == beamOff2) event[beamOff2].daughter1( ev2Dau1V);
1945     } else {
1946       for (int iG = 0; iG < int(iGRecBef.size()); ++iG) {
1947         event[iGRecBef[iG]].statusPos();
1948         event[iGRecBef[iG]].daughters( 0, 0);
1949       }
1950     }
1951     return false;
1952   }
1953
1954   // For global recoil restore the one nominal recoiler, for bookkeeping.
1955   if (!useLocalRecoilNow) {
1956     iRec = iRecBef;
1957     for (int iG = 0; iG < int(iGRecBef.size()); ++iG)
1958     if (iGRecBef[iG] == iRecBef) iRec = iGRec[iG];
1959   }
1960
1961   // For initial-state recoiler also update beam and sHat info.
1962   if (isrTypeNow != 0) {
1963     BeamParticle& beamRec = (isrTypeNow == 1) ? *beamAPtr : *beamBPtr;
1964     double xOld = beamRec[iSysSelRec].x();
1965     double xRec = 2. * pRec.e() / (beamAPtr->e() + beamBPtr->e());
1966     beamRec[iSysSelRec].iPos( iRec);
1967     beamRec[iSysSelRec].x( xRec);
1968     partonSystemsPtr->setSHat( iSysSelRec,
1969     partonSystemsPtr->getSHat(iSysSelRec) * xRec / xOld);
1970   }
1971
1972   // Photon emission: update to new dipole ends; add new photon "dipole".
1973   if (dipSel->flavour == 22) {
1974     dipSel->iRadiator = iRad;
1975     dipSel->iRecoiler = iRec;
1976     // When recoiler was uncharged particle, in resonance decays,
1977     // assign recoil to emitted photons.
1978     if (recoilToColoured && inResonance && event[iRec].chargeType() == 0)
1979       dipSel->iRecoiler = iEmt;
1980     dipSel->pTmax = pTsel;
1981     if (doQEDshowerByGamma) dipEnd.push_back( TimeDipoleEnd(iEmt, iRad,
1982       pTsel, 0, 0, 1, 0, iSysSel, 0));
1983
1984   // Gluon emission: update both dipole ends and add two new ones.
1985   } else if (dipSel->flavour == 21) {
1986     dipSel->iRadiator  = iRad;
1987     dipSel->iRecoiler  = iEmt;
1988     dipSel->systemRec  = iSysSel;
1989     dipSel->isrType    = 0;
1990     dipSel->pTmax      = pTsel;
1991     // Optionally also kill ME corrections after first emission.
1992     if (!doMEafterFirst) dipSel->MEtype = 0;
1993     // PS dec 2010: check normalization of radiating dipole
1994     // Dipole corresponding to the newly created color tag has normal strength
1995     double flexFactor  = (isFlexible) ? dipSel->flexFactor : 1.0;
1996     dipSel->isFlexible = false;
1997     for (int i = 0; i < int(dipEnd.size()); ++i) {
1998       if (dipEnd[i].iRadiator == iRecBef && dipEnd[i].iRecoiler == iRadBef
1999         && dipEnd[i].colType != 0) {
2000         dipEnd[i].iRadiator = iRec;
2001         dipEnd[i].iRecoiler = iEmt;
2002         // Optionally also kill ME corrections after first emission.
2003         if (!doMEafterFirst) dipEnd[i].MEtype = 0;
2004         // Strive to match colour to anticolour inside closed system.
2005         if ( !isFlexible && dipEnd[i].colType * dipSel->colType > 0)
2006           dipEnd[i].iRecoiler = iRad;
2007         dipEnd[i].pTmax = pTsel;
2008         // PS dec 2010: if the (iRadBef,iRecBef) dipole was flexible, the
2009         // same should be true for this (opposite) end. If so, this end keeps
2010         // the modified normalization, so we shouldn't need to do anything.
2011       }
2012     }
2013     int colType = (dipSel->colType > 0) ? 2 : -2 ;
2014     // When recoiler was uncoloured particle, in resonance decays,
2015     // assign recoil to coloured particle.
2016     int iRecMod = iRec;
2017     if (recoilToColoured && inResonance && event[iRec].col() == 0
2018       && event[iRec].acol() == 0) iRecMod = iRad;
2019     dipEnd.push_back( TimeDipoleEnd(iEmt, iRecMod, pTsel,
2020        colType, 0, 0, isrTypeSave, iSysSel, 0));
2021     dipEnd.back().systemRec = iSysSelRec;
2022     // PS dec 2010: the (iEmt,iRec) dipole "inherits" flexible normalization
2023     if (isFlexible) {
2024       dipEnd.back().isFlexible = true;
2025       dipEnd.back().flexFactor = flexFactor;
2026     }
2027     dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel,
2028       -colType, 0, 0, 0, iSysSel, 0));
2029
2030   // Gluon branching to q qbar: update current dipole and other of gluon.
2031   } else if (dipSel->colType != 0) {
2032     for (int i = 0; i < int(dipEnd.size()); ++i) {
2033       // Strive to match colour to anticolour inside closed system.
2034       if ( !isFlexible && dipEnd[i].iRecoiler == iRadBef
2035         && dipEnd[i].colType * dipSel->colType < 0 )
2036         dipEnd[i].iRecoiler = iEmt;
2037       if (dipEnd[i].iRadiator == iRadBef && abs(dipEnd[i].colType) == 2) {
2038         dipEnd[i].colType /= 2;
2039         if (dipEnd[i].system != dipEnd[i].systemRec) continue;
2040
2041         // Note: gluino -> quark + squark gives a deeper radiation dip than
2042         // the more obvious alternative photon decay, so is more realistic.
2043         dipEnd[i].MEtype = 66;
2044         if (&dipEnd[i] == dipSel) dipEnd[i].iMEpartner = iRad;
2045         else                      dipEnd[i].iMEpartner = iEmt;
2046       }
2047     }
2048     dipSel->iRadiator = iEmt;
2049     dipSel->iRecoiler = iRec;
2050     dipSel->pTmax     = pTsel;
2051
2052     // Gluon branching to q qbar: also add two charge dipole ends.
2053     // Note: gluino -> quark + squark gives a deeper radiation dip than
2054     // the more obvious alternative photon decay, so is more realistic.
2055     if (doQEDshowerByQ) {
2056       int chgType = event[iRad].chargeType();
2057       dipEnd.push_back( TimeDipoleEnd(iRad, iEmt, pTsel,
2058         0,  chgType, 0, 0, iSysSel, 66, iEmt));
2059       dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel,
2060         0, -chgType, 0, 0, iSysSel, 66, iRad));
2061     }
2062
2063   // Photon branching to f fbar: inactivate photon "dipole";
2064   // optionally add new charge and colour dipole ends.
2065   } else if (dipSel->gamType != 0) {
2066     dipSel->gamType = 0;
2067     int chgType = event[iRad].chargeType();
2068     int colType = event[iRad].colType();
2069     // MEtype = 102 for charge in vector decay.
2070     if ( chgType != 0 && ( ( doQEDshowerByQ && colType != 0 )
2071       || ( doQEDshowerByL && colType == 0 ) ) ) {
2072       dipEnd.push_back( TimeDipoleEnd(iRad, iEmt, pTsel,
2073         0,  chgType, 0, 0, iSysSel, 102, iEmt));
2074       dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel,
2075         0, -chgType, 0, 0, iSysSel, 102, iRad));
2076     }
2077     // MEtype = 11 for colour in vector decay.
2078     if (colType != 0 && doQCDshower) {
2079       dipEnd.push_back( TimeDipoleEnd(iRad, iEmt, pTsel,
2080          colType, 0, 0, 0, iSysSel, 11, iEmt));
2081       dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel,
2082         -colType, 0, 0, 0, iSysSel, 11, iRad));
2083     }
2084
2085   // Photon_HV emission: update to new dipole ends.
2086   } else if (dipSel->flavour == 4900022) {
2087     dipSel->iRadiator = iRad;
2088     dipSel->iRecoiler = iRec;
2089     dipSel->pTmax = pTsel;
2090
2091   // Gluon_HV emission: update to new dipole ends.
2092   } else if (dipSel->flavour == 4900021) {
2093     dipSel->iRadiator = iRad;
2094     dipSel->iRecoiler = iEmt;
2095     dipSel->pTmax     = pTsel;
2096     for (int i = 0; i < int(dipEnd.size()); ++i)
2097     if (dipEnd[i].iRadiator == iRecBef && dipEnd[i].iRecoiler == iRadBef
2098       && dipEnd[i].isHiddenValley) {
2099       dipEnd[i].iRadiator = iRec;
2100       dipEnd[i].iRecoiler = iEmt;
2101       dipEnd[i].pTmax     = pTsel;
2102     }
2103     int colvType = (dipSel->colvType > 0) ? 2 : -2 ;
2104     dipEnd.push_back( TimeDipoleEnd(iEmt, iRec, pTsel,
2105        0, 0, 0, isrTypeSave, iSysSel, 0, -1, false, true, colvType) );
2106     dipEnd.back().systemRec = iSysSelRec;
2107     dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel,
2108       0, 0, 0, 0, iSysSel, 0, -1, false, true, -colvType) );
2109   }
2110
2111   // Copy or set lifetime for new final state.
2112   if (event[iRad].id() == event[iRadBef].id()) {
2113     event[iRad].tau( event[iRadBef].tau() );
2114   } else {
2115     event[iRad].tau( event[iRad].tau0() * rndmPtr->exp() );
2116     event[iEmt].tau( event[iEmt].tau0() * rndmPtr->exp() );
2117   }
2118   event[iRec].tau( event[iRecBef].tau() );
2119
2120   // Now update other dipoles that also involved the radiator or recoiler.
2121   for (int i = 0; i < int(dipEnd.size()); ++i) {
2122     // PS dec 2010: if radiator was flexible and now is normal, there may
2123     // be other flexible dipoles that need updating.
2124     if (isFlexible && !dipSel->isFlexible && dipEnd[i].isFlexible) {
2125       if (dipEnd[i].iRecoiler  == iRadBef) dipEnd[i].iRecoiler = iEmt;
2126       if (dipEnd[i].iRadiator  == iRadBef) {
2127         dipEnd[i].iRadiator = iEmt;
2128         if (dipEnd[i].colType == 1 && dipSel->flavour == 21)
2129           dipEnd[i].colType = 2;
2130         if (dipEnd[i].colType ==-1 && dipSel->flavour == 21)
2131           dipEnd[i].colType =-2;
2132       }
2133     }
2134     if (dipEnd[i].iRadiator  == iRadBef) dipEnd[i].iRadiator  = iRad;
2135     if (dipEnd[i].iRecoiler  == iRadBef) dipEnd[i].iRecoiler  = iRad;
2136     if (dipEnd[i].iMEpartner == iRadBef) dipEnd[i].iMEpartner = iRad;
2137     if (useLocalRecoilNow) {
2138       if (dipEnd[i].iRadiator  == iRecBef) dipEnd[i].iRadiator  = iRec;
2139       if (dipEnd[i].iRecoiler  == iRecBef) dipEnd[i].iRecoiler  = iRec;
2140       if (dipEnd[i].iMEpartner == iRecBef) dipEnd[i].iMEpartner = iRec;
2141     } else {
2142       for (int iG = 0; iG < int(iGRecBef.size()); ++iG) {
2143         if (dipEnd[i].iRadiator  == iGRecBef[iG])
2144             dipEnd[i].iRadiator  =  iGRec[iG];
2145         if (dipEnd[i].iRecoiler  == iGRecBef[iG])
2146             dipEnd[i].iRecoiler  =  iGRec[iG];
2147         if (dipEnd[i].iMEpartner == iGRecBef[iG])
2148             dipEnd[i].iMEpartner =  iGRec[iG];
2149       }
2150     }
2151   }
2152
2153   // PS Apr 2011
2154   // Update any junctions downstream of this branching (if necessary)
2155   // (This happens, e.g., via LHEF, when adding showers to intermediate
2156   //  coloured resonances whose decays involved junctions)
2157   for (int iJun = 0; iJun < event.sizeJunction(); iJun++) {
2158     // Number of incoming colour lines for this junction.
2159     int nIncoming = (event.kindJunction(iJun)-1)/2;
2160     // Check radiator colour or anticolour, depending on junction kind
2161     // (if junction, incoming = anticolours, and vice versa)
2162     int colChk = 0;
2163     colChk = ( event.kindJunction(iJun) % 2 == 0 )
2164            ? event[iRadBef].col() : event[iRadBef].acol();
2165     // Loop over incoming junction ends
2166     for (int iCol = 0; iCol < nIncoming; iCol++) {
2167       int colJun = event.colJunction( iJun, iCol);
2168       // If match, update junction end with new upstream (anti)colour
2169       if (colJun == colChk) {
2170         int colNew = 0;
2171         if ( event.kindJunction(iJun) % 2 == 0 ) colNew = colRad;
2172         else colNew = acolRad;
2173         event.colJunction( iJun, iCol, colNew );
2174       }
2175     }
2176   }
2177
2178   // Finally update the list of all partons in all systems.
2179   partonSystemsPtr->replace(iSysSel, iRadBef, iRad);
2180   partonSystemsPtr->addOut(iSysSel, iEmt);
2181   if (useLocalRecoilNow)
2182     partonSystemsPtr->replace(iSysSelRec, iRecBef, iRec);
2183   else {
2184     for (int iG = 0; iG < int(iGRecBef.size()); ++iG)
2185     partonSystemsPtr->replace(iSysSel, iGRecBef[iG], iGRec[iG]);
2186   }
2187
2188   // Done.
2189   return true;
2190
2191 }
2192
2193 //--------------------------------------------------------------------------
2194
2195 // Rescatter: If a dipole stretches between two different systems, those
2196 //            systems will no longer locally conserve momentum. These
2197 //            imbalances become problematic when ISR or primordial kT
2198 //            is switched on as these steps involve Lorentz boosts.
2199 //
2200 //            'rescatterPropagateRecoil' tries to fix momentum in all
2201 //            systems by propogating recoil momentum through all
2202 //            intermediate systems. As the momentum transfer is already
2203 //            defined, this can lead to internal lines gaining a
2204 //            virtuality.
2205
2206 // Useful definitions for a pair of integers and a vector of pairs
2207 typedef pair < int, int >  pairInt;
2208 typedef vector < pairInt > vectorPairInt;
2209
2210 //--------------------------------------------------------------------------
2211
2212 // findParentSystems
2213 //  Utility routine to find all parent systems of a given system
2214 //  Returns a vector of pairs of integers with:
2215 //   a) The system index, including the starting system (negative
2216 //      if (b) points to a parent system, positive if (b) points
2217 //      to a daughter system
2218 //   b) The event record index that is the path out of the system
2219 //      (if forwards == false, this is an incoming parton to the
2220 //      system, and is +ve if side A or -ve if side B,
2221 //      if forwards == true, this is an outgoing parton from the
2222 //      system).
2223 //  Returns as empty vector on failure
2224 //  Note: this assumes single rescattering only and therefore only
2225 //        one possible parent system
2226
2227 inline vectorPairInt findParentSystems(const int sys,
2228   Event& event, PartonSystems* partonSystemsPtr, bool forwards) {
2229
2230   vectorPairInt parentSystems;
2231   parentSystems.reserve(10);
2232
2233   int iSysCur = sys;
2234   while (true) {
2235     // Get two incoming partons
2236     int iInA = partonSystemsPtr->getInA(iSysCur);
2237     int iInB = partonSystemsPtr->getInB(iSysCur);
2238
2239     // Check if either of these links to another system
2240     int iIn = 0;
2241     if (event[iInA].isRescatteredIncoming()) iIn =  iInA;
2242     if (event[iInB].isRescatteredIncoming()) iIn = -iInB;
2243
2244     // Save the current system to the vector
2245     parentSystems.push_back( pairInt(-iSysCur, iIn) );
2246     if (iIn == 0) break;
2247
2248     int iInAbs  = abs(iIn);
2249     int iMother = event[iInAbs].mother1();
2250     iSysCur     = partonSystemsPtr->getSystemOf(iMother);
2251     if (iSysCur == -1) {
2252       parentSystems.clear();
2253       break;
2254     }
2255   } // while (true)
2256
2257   // If forwards is set, change all event record indices to go to daughter
2258   // systems rather than parent systems
2259   if (forwards) {
2260     vectorPairInt::reverse_iterator rit;
2261     for (rit = parentSystems.rbegin(); rit < (parentSystems.rend() - 1);
2262          ++rit) {
2263       pairInt &cur  = *rit;
2264       pairInt &next = *(rit + 1);
2265       cur.first     = -cur.first;
2266       cur.second    = (next.second < 0) ? -event[abs(next.second)].mother1() :
2267                                            event[abs(next.second)].mother1();
2268     }
2269   }
2270
2271   return parentSystems;
2272 }
2273
2274 //--------------------------------------------------------------------------
2275
2276 // rescatterPropagateRecoil
2277 //  Fix up momentum in all intermediate systems when radiator and recoiler
2278 //  systems are different. The strategy is to look at all parent systems
2279 //  from the radiator system and the recoiler system and find where they
2280 //  intersect.
2281
2282 bool TimeShower::rescatterPropagateRecoil( Event& event, Vec4& pNew) {
2283
2284   // Some useful variables for later
2285   int  iRadBef    = dipSel->iRadiator;
2286        iSysSel    = dipSel->system;
2287   int  iSysSelRec = dipSel->systemRec;
2288   Vec4 pImbal     = pNew - event[iRadBef].p();
2289
2290   // Store changes locally at first in case we veto the branching
2291   // eventMod stores index into the event record and new 4-vector
2292   vector < pair < int, Vec4 > > eventMod;
2293   eventMod.reserve(10);
2294   // systemMod stores system index (iSys) and system-parton index (iMem)
2295   //   iMem >=  0 - index into outgoing partons (iOut)
2296   //   iMem == -1 - incoming A
2297   //   iMem == -2 - incoming B
2298   vectorPairInt systemMod;
2299   systemMod.reserve(10);
2300
2301   // Find all parent systems from radiating and recoiling systems
2302   vectorPairInt radParent = findParentSystems(iSysSel, event,
2303                                               partonSystemsPtr, false);
2304   vectorPairInt recParent = findParentSystems(iSysSelRec, event,
2305                                               partonSystemsPtr, true);
2306   if (radParent.size() == 0 || recParent.size() == 0) {
2307     // This should never happen
2308     infoPtr->errorMsg("Error in TimeShower::rescatterPropagateRecoil: "
2309       "couldn't find parent system; branching vetoed");
2310     return false;
2311   }
2312   // Find the system that connects radiating and recoiling system
2313   bool foundPath = false;
2314   unsigned int iRadP = 0;
2315   unsigned int iRecP = 0;
2316   for (iRadP = 0; iRadP < radParent.size(); iRadP++) {
2317     for (iRecP = 0; iRecP < recParent.size(); iRecP++)
2318       if (abs(radParent[iRadP].first) == abs(recParent[iRecP].first)) {
2319         foundPath = true;
2320         break;
2321       }
2322     if (foundPath) break;
2323   }
2324   if (!foundPath) {
2325     // Can fail e.g. for QED dipoles where there is no connection
2326     // between radiator and recoiler systems
2327     infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: "
2328       "couldn't find recoil path; branching vetoed");
2329     return false;
2330   }
2331
2332   // Join together to form complete path from radiating system
2333   // to recoiling system
2334   vectorPairInt path;
2335   if (radParent.size() > 1)
2336     path.assign(radParent.begin(), radParent.begin() + iRadP);
2337   if (recParent.size() > 1)
2338     path.insert(path.end(), recParent.rend() - iRecP - 1,
2339                 recParent.rend() - 1);
2340
2341   // Follow the path fixing up momenta as we go
2342   for (unsigned int i = 0; i < path.size(); i++) {
2343     // Line out of the current system
2344     bool isIncoming  = (path[i].first < 0) ? true : false;
2345     int  iSysCur     = abs(path[i].first);
2346     bool isIncomingA = (path[i].second > 0) ? true : false;
2347     int  iLink       = abs(path[i].second);
2348
2349     int iMemCur;
2350     if (isIncoming) iMemCur = (isIncomingA) ? -1 : -2;
2351     else {
2352       iMemCur = -1;
2353       for (int j = 0; j < partonSystemsPtr->sizeOut(iSysCur); j++)
2354         if (partonSystemsPtr->getOut(iSysCur, j) == iLink) {
2355           iMemCur = j;
2356           break;
2357         }
2358       if (iMemCur == -1) {
2359         // This should never happen
2360         infoPtr->errorMsg("Error in TimeShower::rescatterPropagateRecoil: "
2361           "couldn't find parton system; branching vetoed");
2362         return false;
2363       }
2364     }
2365
2366     Vec4 pMod = (isIncoming) ? event[iLink].p() + pImbal :
2367                                event[iLink].p() - pImbal;
2368     eventMod.push_back(pair <int, Vec4> (iLink, pMod));
2369     systemMod.push_back(pairInt(iSysCur, iMemCur));
2370
2371     // Calculate sHat of iSysCur
2372     int  iInCurA = partonSystemsPtr->getInA(iSysCur);
2373     int  iInCurB = partonSystemsPtr->getInB(iSysCur);
2374     Vec4 pTotCur = event[iInCurA].p() + event[iInCurB].p();
2375
2376     // If iMemCur is -1 or -2, then we must have changed the sHat of iSysCur
2377     if (iMemCur < 0) pTotCur += (isIncoming) ? pImbal : -pImbal;
2378     double sHatCur = pTotCur.m2Calc();
2379
2380     // The fixed-up incoming and outgoing partons should not have
2381     // too large a virtuality in relation to the system mass-square.
2382     if (abs(pMod.m2Calc()) > MAXVIRTUALITYFRACTION * sHatCur) {
2383       infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: "
2384         "virtuality much larger than sHat; branching vetoed");
2385       return false;
2386     }
2387
2388     // Outgoing ones should also not have too large negative energy
2389     // in the rest frame of the system.
2390     if (!isIncoming && pMod * pTotCur < -MAXNEGENERGYFRACTION * sHatCur) {
2391       infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: "
2392         "rest frame energy too negative; branching vetoed");
2393       return false;
2394     }
2395
2396     // Veto negative sHat.
2397     if (sHatCur < 0.0) {
2398       infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: "
2399         "sHat became negative; branching vetoed");
2400       return false;
2401     }
2402
2403     // Line into the new current system
2404     iLink   = (isIncoming) ? event[iLink].mother1()  :
2405                              event[iLink].daughter1();
2406     iSysCur = partonSystemsPtr->getSystemOf(iLink, true);
2407
2408     if (!isIncoming) iMemCur = (isIncomingA) ? -1 : -2;
2409     else {
2410       iMemCur = -1;
2411       for (int j = 0; j < partonSystemsPtr->sizeOut(iSysCur); j++)
2412         if (partonSystemsPtr->getOut(iSysCur, j) == iLink) {
2413           iMemCur = j;
2414           break;
2415         }
2416       if (iMemCur == -1) {
2417         // This should never happen
2418         infoPtr->errorMsg("Error in TimeShower::rescatterPropagateRecoil: "
2419           "couldn't find parton system; branching vetoed");
2420         return false;
2421       }
2422     }
2423
2424     pMod = (isIncoming) ? event[iLink].p() + pImbal :
2425                           event[iLink].p() - pImbal;
2426     eventMod.push_back(pair <int, Vec4> (iLink, pMod));
2427     systemMod.push_back(pairInt(iSysCur, iMemCur));
2428
2429     // Calculate sHat of iSysCur
2430     iInCurA = partonSystemsPtr->getInA(iSysCur);
2431     iInCurB = partonSystemsPtr->getInB(iSysCur);
2432     pTotCur = event[iInCurA].p() + event[iInCurB].p();
2433
2434     // If iMemCur is -1 or -2, then we must have changed the sHat of iSysCur
2435     if (iMemCur < 0) pTotCur += (isIncoming) ? pImbal : -pImbal;
2436     sHatCur = pTotCur.m2Calc();
2437
2438     // The fixed-up incoming and outgoing partons should not have
2439     // too large a virtuality in relation to the system mass-square.
2440     if (abs(pMod.m2Calc()) > MAXVIRTUALITYFRACTION * sHatCur) {
2441       infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: "
2442         "virtuality much larger than sHat; branching vetoed");
2443       return false;
2444     }
2445
2446     // Outgoing ones should also not have too large negative energy
2447     // in the rest frame of the system.
2448     if (!isIncoming && pMod * pTotCur < -MAXNEGENERGYFRACTION * sHatCur) {
2449       infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: "
2450         "rest frame energy too negative; branching vetoed");
2451       return false;
2452     }
2453
2454     // Veto negative sHat
2455     if (sHatCur < 0.0) {
2456       infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: "
2457         "sHat became negative; branching vetoed");
2458       return false;
2459     }
2460
2461     // Do negative energy veto
2462     if (VETONEGENERGY && pMod.e() < 0.0) {
2463       infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: "
2464         "energy became negative; branching vetoed");
2465       return false;
2466     }
2467
2468   } // for (unsigned int i = 0; i < path.size(); i++)
2469
2470   // If no vetos by this point, apply the changes to the event record
2471   // An incoming particle with changed momentum is given status code -54,
2472   // an outgoing particle with changed momentum is given status code -55
2473   for (unsigned int i = 0; i < eventMod.size(); i++) {
2474     int idx    = eventMod[i].first;
2475     Vec4 &pMod = eventMod[i].second;
2476     int iSys   = systemMod[i].first;
2477     int iMem   = systemMod[i].second;
2478
2479     // If incoming to a process then set the copy to be the mother
2480     if (event[idx].isRescatteredIncoming()) {
2481       int mother1 = event[idx].mother1();
2482       idx = event.copy(idx, -54);
2483       event[mother1].daughters(idx, idx);
2484
2485       // Update beam information if necessary
2486       double eCM = sqrt(m2( beamAPtr->p(), beamBPtr->p()));
2487       if        (iMem == -1) {
2488         partonSystemsPtr->setInA(iSys, idx);
2489         (*beamAPtr)[iSys].x((pMod.e() + pMod.pz()) / eCM);
2490         (*beamAPtr)[iSys].m(pMod.mCalc());
2491         (*beamAPtr)[iSys].p(pMod);
2492         (*beamAPtr)[iSys].iPos(idx);
2493       } else if (iMem == -2) {
2494         partonSystemsPtr->setInB(iSys, idx);
2495         (*beamBPtr)[iSys].x((pMod.e() - pMod.pz()) / eCM);
2496         (*beamBPtr)[iSys].m(pMod.mCalc());
2497         (*beamBPtr)[iSys].p(pMod);
2498         (*beamBPtr)[iSys].iPos(idx);
2499       } else {
2500         // This should never happen
2501         infoPtr->errorMsg("Error in TimeShower::rescatterPropagateRecoil: "
2502         "internal bookeeping error");
2503       }
2504
2505     // Otherwise set the new event record entry to be the daughter
2506     } else {
2507       int daughter1 = event[idx].daughter1();
2508       idx = event.copy(idx, 55);
2509       event[idx].statusNeg();
2510       event[daughter1].mothers(idx, idx);
2511
2512       partonSystemsPtr->setOut(iSys, iMem, idx);
2513     }
2514
2515     event[idx].p( eventMod[i].second );
2516     event[idx].m( event[idx].mCalc() );
2517   }
2518
2519   return true;
2520 }
2521
2522
2523 //--------------------------------------------------------------------------
2524
2525 // Find class of QCD ME correction.
2526 // MEtype classification follow codes in Norrbin article,
2527 // additionally -1 = try to find type, 0 = no ME corrections.
2528 // Warning: not yet tried out to do a correct assignment in
2529 // arbitrary multiparton configurations! ??
2530
2531 void TimeShower::findMEtype( Event& event, TimeDipoleEnd& dip) {
2532
2533   // Initial value. Mark if no ME corrections to be applied.
2534   bool setME = true;
2535   if (!doMEcorrections) setME = false;
2536   int iMother  = event[dip.iRadiator].mother1();
2537   int iMother2 = event[dip.iRadiator].mother2();
2538
2539   // Allow ME corrections for Hidden Valley pair in 2 -> 2.
2540   if (dip.isHiddenValley && event[dip.iRecoiler].id()
2541     == -event[dip.iRadiator].id());
2542
2543   // Else no ME corrections in 2 -> n processes.
2544   else {
2545     if (iMother2 != iMother && iMother2 != 0) setME = false;
2546     if (event[dip.iRecoiler].mother1() != iMother)  setME = false;
2547     if (event[dip.iRecoiler].mother2() != iMother2) setME = false;
2548   }
2549
2550   // No ME corrections for recoiler in initial state.
2551   if (event[dip.iRecoiler].status() < 0) setME = false;
2552
2553   // No ME corrections for recoiler not in same system
2554   if (dip.system != dip.systemRec) setME = false;
2555
2556   // Done if no ME to be set.
2557   if (!setME) {
2558     dip.MEtype = 0;
2559     return;
2560   }
2561
2562   // If no ME partner set, assume it is the recoiler.
2563   if (dip.iMEpartner < 0) dip.iMEpartner = dip.iRecoiler;
2564
2565   // Now begin processing of colour dipole, including Hidden Valley.
2566   if (dip.colType != 0 || dip.colvType != 0) {
2567     bool isHiddenColour = (dip.colvType != 0);
2568
2569     // Find daughter types (may or may not be used later on).
2570     int idDau1      = event[dip.iRadiator].id();
2571     int idDau2      = event[dip.iMEpartner].id();
2572     int dau1Type    = findMEparticle(idDau1, isHiddenColour);
2573     int dau2Type    = findMEparticle(idDau2, isHiddenColour);
2574     int minDauType  = min(dau1Type, dau2Type);
2575     int maxDauType  = max(dau1Type, dau2Type);
2576
2577     // Reorder dipole ends in kinematics. Split ME expression in two sides.
2578     dip.MEorder     = (dau2Type >= dau1Type);
2579     dip.MEsplit     = (maxDauType <= 6);
2580     dip.MEgluinoRec = false;
2581
2582     // If type already set (or set not to have) then done.
2583     if (minDauType == 0 && dip.MEtype < 0) dip.MEtype = 0;
2584     if (dip.MEtype >= 0) return;
2585     dip.MEtype = 0;
2586
2587     // For H -> gg -> ggg we found that DGLAP kernels do better than eikonal.
2588     if (dau1Type == 4 && dau2Type == 4) return;
2589
2590     // Find mother type.
2591     int idMother = 0;
2592     if ( event[dip.iRecoiler].mother1() == iMother && iMother >= 0)
2593       idMother = event[iMother].id();
2594     int motherType = (idMother != 0)
2595       ? findMEparticle(idMother, isHiddenColour) : 0;
2596
2597     // When a mother if not known then use colour and spin content to guess.
2598     if (motherType == 0) {
2599       int col1  = event[dip.iRadiator].col();
2600       int acol1 = event[dip.iRadiator].acol();
2601       int col2  = event[dip.iMEpartner].col();
2602       int acol2 = event[dip.iMEpartner].acol();
2603       // spinT = 0/1 = integer or half-integer.
2604       int spinT = ( event[dip.iRadiator].spinType()
2605                 + event[dip.iMEpartner].spinType() )%2;
2606       // Colour singlet mother.
2607       if ( col1 == acol2 && acol1 == col2 )
2608         motherType = (spinT == 0) ? 7 : 9;
2609       // Colour octet mother.
2610       else if ( (col1 == acol2 && acol1 != 0 && col2 != 0)
2611         || (acol1 == col2 && col1 != 0 && acol2 != 0) )
2612         motherType = (spinT == 0) ? 4 : 5;
2613       // Colour triplet mother.
2614       else if ( (col1 == acol2 && acol1 != col2)
2615         || (acol1 == col2 && col1 != acol2) )
2616         motherType = (spinT == 0) ? 2 : 1;
2617       // If no colours are matched then cannot have common mother, so done.
2618       else return;
2619     }
2620
2621     // Now start from default, which is eikonal ME corrections,
2622     // and try to find matching ME cases below.
2623     int MEkind = 0;
2624     int MEcombi = 4;
2625     dip.MEmix = 0.5;
2626
2627     // Hidden Valley with massive gamma_v covered by two special cases.
2628     if (isHiddenColour && brokenHVsym) {
2629       MEkind = (dau2Type == 0 || dau2Type > 6) ? 30 : 31;
2630       dip.MEtype = 5 * MEkind + 1;
2631       return;
2632     }
2633
2634     // Triplet recoiling against gluino needs enhanced radiation
2635     // to match to matrix elements.
2636     dip.MEgluinoRec = (dau1Type >= 1 && dau1Type <= 3 && dau2Type == 5);
2637
2638     // Vector/axial vector -> q + qbar.
2639     if (minDauType == 1 && maxDauType == 1 &&
2640       (motherType == 4 || motherType == 7) ) {
2641       MEkind = 2;
2642       if (idMother == 21 || idMother == 22) MEcombi = 1;
2643       else if (idMother == 23 || idDau1 + idDau2 == 0) {
2644         MEcombi = 3;
2645         dip.MEmix = gammaZmix( event, iMother, dip.iRadiator, dip.iRecoiler );
2646       }
2647       else if (idMother == 24) MEcombi = 4;
2648     }
2649     // For chi -> chi q qbar, use V/A -> q qbar as first approximation.
2650     else if (minDauType == 1 && maxDauType == 1 && motherType == 9)
2651       MEkind = 2;
2652
2653     // q -> q + V.
2654     else if (minDauType == 1 && maxDauType == 7 && motherType == 1)
2655       MEkind = 3;
2656       if (idDau1 == 22 || idDau2 == 22) MEcombi = 1;
2657
2658     // Scalar/pseudoscalar -> q + qbar; q -> q + S.
2659     else if (minDauType == 1 && maxDauType == 1 && motherType == 8) {
2660       MEkind = 4;
2661       if (idMother == 25 || idMother == 35 || idMother == 37) MEcombi = 1;
2662       else if (idMother == 36) MEcombi = 2;
2663     }
2664     else if (minDauType == 1 && maxDauType == 8 && motherType == 1)
2665       MEkind = 5;
2666
2667     // V -> ~q + ~qbar; ~q -> ~q + V; S -> ~q + ~qbar; ~q -> ~q + S.
2668     else if (minDauType == 2 && maxDauType == 2 && (motherType == 4
2669       || motherType == 7) ) MEkind = 6;
2670     else if (minDauType == 2 && (maxDauType == 4 || maxDauType == 7)
2671       && motherType == 2) MEkind = 7;
2672     else if (minDauType == 2 && maxDauType == 2 && motherType == 8)
2673       MEkind = 8;
2674     else if (minDauType == 2 && maxDauType == 8 && motherType == 2)
2675       MEkind = 9;
2676
2677     // chi -> q + ~qbar; ~q -> q + chi; q -> ~q + chi.
2678     else if (minDauType == 1 && maxDauType == 2 && motherType == 9)
2679       MEkind = 10;
2680     else if (minDauType == 1 && maxDauType == 9 && motherType == 2)
2681       MEkind = 11;
2682     else if (minDauType == 2 && maxDauType == 9 && motherType == 1)
2683       MEkind = 12;
2684
2685     // ~g -> q + ~qbar; ~q -> q + ~g; q -> ~q + ~g.
2686     else if (minDauType == 1 && maxDauType == 2 && motherType == 5)
2687       MEkind = 13;
2688     else if (minDauType == 1 && maxDauType == 5 && motherType == 2)
2689       MEkind = 14;
2690     else if (minDauType == 2 && maxDauType == 5 && motherType == 1)
2691       MEkind = 15;
2692
2693     // In cases where coloured spin 1 particle involved use spin 0.
2694     // V_coloured -> q + l.
2695     else if (minDauType == 1 && maxDauType == 9 && motherType == 3)
2696       MEkind = 11;
2697     // q -> V_coloured + l;
2698     else if (minDauType == 3 && maxDauType == 9 && motherType == 1)
2699       MEkind = 12;
2700
2701     // g (+V, S) -> ~g + ~g (eikonal approximation).
2702     else if (minDauType == 5 && maxDauType == 5) MEkind = 16;
2703
2704     // Save ME type and gamma_5 admixture.
2705     dip.MEtype = 5 * MEkind + MEcombi;
2706
2707   // Now begin processing of charge dipole - still primitive.
2708   } else if (dip.chgType != 0) {
2709
2710     // Set defaults for QED case; then possibly done.
2711     dip.MEorder = true;
2712     dip.MEsplit = true;
2713     if (dip.MEtype >= 0) return;
2714
2715     // So far only ME corrections for q qbar or l lbar.
2716     int idDau1 = event[dip.iRadiator].id();
2717     int idDau2 = event[dip.iMEpartner].id();
2718     if (abs(idDau1) < 9 && abs(idDau2) < 9 && idDau1 * idDau2 < 0) ;
2719     else if (abs(idDau1) > 10 && abs(idDau1) < 19 && abs(idDau2) > 10
2720       && abs(idDau2) < 19 && idDau1 * idDau2 < 0) ;
2721     else { dip.MEtype = 0; return; }
2722
2723     // Distinguish charge sum != 0 or = 0; in latter assume vector source.
2724     dip.MEtype = 101;
2725     if (idDau1 + idDau2 == 0) dip.MEtype = 102;
2726     dip.MEmix = 1.;
2727   }
2728
2729 }
2730
2731 //--------------------------------------------------------------------------
2732
2733 // Find type of particle for ME type: 0 = unknown, 1 = quark, 2 = squark,
2734 // 3 = spare triplet, 4 = gluon, 5 = gluino, 6 = spare octet,
2735 // 7 = vector boson, 8 = colourless scalar, 9 = colourless spin 1/2.
2736
2737 int TimeShower::findMEparticle( int id, bool isHiddenColour) {
2738
2739   // find colour and spin of particle.
2740   int type = 0;
2741   int colType = abs(particleDataPtr->colType(id));
2742   int spinType = particleDataPtr->spinType(id);
2743
2744   // For hidden valley particle treat HV colour as normal one.
2745   // Note: no need to assign gv/gammav since not in ME.
2746   if (isHiddenColour) {
2747     colType = 0;
2748     int idAbs = abs(id);
2749     if (  (idAbs > 4900000 && idAbs < 4900007)
2750        || (idAbs > 4900010 && idAbs < 4900017)
2751        || idAbs == 4900101) colType = 1;
2752   }
2753
2754   // Find particle type from colour and spin.
2755   if      (colType == 1 && spinType == 2) type = 1;
2756   else if (colType == 1 && spinType == 1) type = 2;
2757   else if (colType == 1)                  type = 3;
2758   else if (colType == 2 && spinType == 3) type = 4;
2759   else if (colType == 2 && spinType == 2) type = 5;
2760   else if (colType == 2)                  type = 6;
2761   else if (colType == 0 && spinType == 3) type = 7;
2762   else if (colType == 0 && spinType == 1) type = 8;
2763   else if (colType == 0 && spinType == 2) type = 9;
2764
2765   // Done.
2766   return type;
2767
2768 }
2769
2770 //--------------------------------------------------------------------------
2771
2772 // Find mixture of V and A in gamma/Z: energy- and flavour-dependent.
2773
2774 double TimeShower::gammaZmix( Event& event, int iRes, int iDau1, int iDau2) {
2775
2776   // Try to identify initial flavours; use e+e- as default.
2777   int idIn1 = -11;
2778   int idIn2 = 11;
2779   int iIn1  = (iRes >= 0) ? event[iRes].mother1() : -1;
2780   int iIn2  = (iRes >= 0) ? event[iRes].mother2() : -1;
2781   if (iIn1 >=0) idIn1 = event[iIn1].id();
2782   if (iIn2 >=0) idIn2 = event[iIn1].id();
2783
2784   // In processes f + g/gamma -> f + Z only need find one fermion.
2785   if (idIn1 == 21 || idIn1 == 22) idIn1 = -idIn2;
2786   if (idIn2 == 21 || idIn2 == 22) idIn2 = -idIn1;
2787
2788   // Initial flavours and couplings; return if don't make sense.
2789   if (idIn1 + idIn2 != 0 ) return 0.5;
2790   int idInAbs = abs(idIn1);
2791   if (idInAbs == 0 || idInAbs > 18 ) return 0.5;
2792   double ei = coupSMPtr->ef(idInAbs);
2793   double vi = coupSMPtr->vf(idInAbs);
2794   double ai = coupSMPtr->af(idInAbs);
2795
2796   // Final flavours and couplings; return if don't make sense.
2797   if (event[iDau1].id() + event[iDau2].id() != 0) return 0.5;
2798   int idOutAbs = abs(event[iDau1].id());
2799   if (idOutAbs == 0 || idOutAbs >18 ) return 0.5;
2800   double ef = coupSMPtr->ef(idOutAbs);
2801   double vf = coupSMPtr->vf(idOutAbs);
2802   double af = coupSMPtr->af(idOutAbs);
2803
2804   // Calculate prefactors for interference and resonance part.
2805   Vec4 psum = event[iDau1].p() + event[iDau2].p();
2806   double sH = psum.m2Calc();
2807   double intNorm = 2. * thetaWRat * sH * (sH - mZ*mZ)
2808     / ( pow2(sH - mZ*mZ) + pow2(sH * gammaZ / mZ) );
2809   double resNorm = pow2(thetaWRat * sH)
2810     / ( pow2(sH - mZ*mZ) + pow2(sH * gammaZ / mZ) );
2811
2812   // Calculate vector and axial expressions and find mix.
2813   double vect = ei*ei * ef*ef + ei*vi * intNorm * ef*vf
2814     + (vi*vi + ai*ai) * resNorm * vf*vf;
2815   double axiv = (vi*vi + ai*ai) * resNorm * af*af;
2816   return vect / (vect + axiv);
2817 }
2818
2819 //--------------------------------------------------------------------------
2820
2821 // Set up to calculate QCD ME correction with calcMEcorr.
2822 // Normally for primary particles, but also from g/gamma -> f fbar.
2823
2824 double TimeShower::findMEcorr(TimeDipoleEnd* dip, Particle& rad,
2825   Particle& partner, Particle& emt) {
2826
2827   // Initial values and matrix element kind.
2828   double wtME    = 1.;
2829   double wtPS    = 1.;
2830   int    MEkind  = dip->MEtype / 5;
2831   int    MEcombi = dip->MEtype % 5;
2832
2833   // Construct ME variables.
2834   Vec4   sum     = rad.p() + partner.p() + emt.p();
2835   double eCMME   = sum.mCalc();
2836   double x1      = 2. * (sum * rad.p()) / pow2(eCMME);
2837   double x2      = 2. * (sum * partner.p()) / pow2(eCMME);
2838   double r1      = rad.m() / eCMME;
2839   double r2      = partner.m() / eCMME;
2840   double r3      = 0.;
2841
2842   // Evaluate kinematics for Hidden Valley with massive gamma_v.
2843   double gammavCorr = 1.;
2844   if (dip->colvType != 0 && brokenHVsym) {
2845     r3              = emt.m() / eCMME;
2846     double x3Tmp    = 2. - x1 - x2;
2847     gammavCorr      = x3Tmp / (x3Tmp - kRad * (x1 + x3Tmp));
2848     // For Q_v Qbar_v pair correct kinematics to common average mass.
2849     if (MEkind == 31) {
2850       double m2Pair = (rad.p() + partner.p()).m2Calc();
2851       double m2Avg  = 0.5 * (rad.m2() + partner.m2())
2852                     - 0.25 * pow2(rad.m2() - partner.m2()) / m2Pair;
2853       r1            = sqrt(m2Avg) / eCMME;
2854       r2            = r1;
2855       double xShift = 0.5 * (x1 + x2) * (partner.m2() - rad.m2()) / m2Pair;
2856       x1           += xShift;
2857       x2           -= xShift;
2858     }
2859   }
2860
2861   // Derived ME variables, suitably protected.
2862   double x1minus = max(XMARGIN, 1. + r1*r1 - r2*r2 - x1);
2863   double x2minus = max(XMARGIN, 1. + r2*r2 - r1*r1 - x2) ;
2864   double x3      = max(XMARGIN, 2. - x1 - x2);
2865
2866   // Begin processing of QCD dipoles.
2867   if (dip->colType !=0 || dip->colvType != 0) {
2868
2869     // Evaluate normal ME, for proper order of particles.
2870     if (dip->MEorder)
2871          wtME = calcMEcorr(MEkind, MEcombi, dip->MEmix, x1, x2, r1, r2, r3);
2872     else wtME = calcMEcorr(MEkind, MEcombi, dip->MEmix, x2, x1, r2, r1, r3);
2873
2874     // Split up total ME when two radiating particles.
2875     if (dip->MEsplit) wtME = wtME * x1minus / x3;
2876
2877     // Evaluate shower rate to be compared with.
2878     wtPS = 2. / ( x3 * x2minus );
2879     if (dip->MEgluinoRec) wtPS *= 9./4.;
2880     if (dip->colvType != 0 && brokenHVsym) wtPS *= gammavCorr;
2881
2882   // For generic charge combination currently only massless expression.
2883   // (Masses included only to respect phase space boundaries.)
2884   } else if (dip->chgType !=0 && dip->MEtype == 101) {
2885     double chg1 = particleDataPtr->charge(rad.id());
2886     double chg2 = particleDataPtr->charge(partner.id());
2887     wtME = (x1*x1 + x2*x2) * pow2( chg1 * x1minus / x3
2888       - chg2 * x2minus / x3 );
2889     wtPS = 2. * ( chg1*chg1 * x1minus / x3 + chg2*chg2 * x2minus / x3 );
2890
2891   // For flavour neutral system assume vector source and include masses.
2892   } else if (dip->chgType !=0 && dip->MEtype == 102) {
2893     wtME = calcMEcorr(2, 1, dip->MEmix, x1, x2, r1, r2) * x1minus / x3;
2894     wtPS = 2. / ( x3 * x2minus );
2895   }
2896   if (wtME > wtPS) infoPtr->errorMsg("Warning in TimeShower::findMEcorr: "
2897     "ME weight above PS one");
2898
2899   // Return ratio of actual ME to assumed PS rate of emission.
2900   return wtME / wtPS;
2901 }
2902
2903 //--------------------------------------------------------------------------
2904
2905 // Matrix elements for gluon (or photon) emission from
2906 // a two-body state; to be used by the parton shower routine.
2907 // Here x_i = 2 E_i/E_cm, r_i = m_i/E_cm and
2908 // 1/sigma_0 d(sigma)/d(x_1)d(x_2) = (alpha-strong/2 pi) * C_F * (this),
2909 // i.e. normalization is such that one recovers the familiar
2910 // (x_1^2 + x_2^2)/((1-x_1)*(1-x_2)) for the massless case.
2911 // Coupling structure:
2912 // kind =  1 : eikonal soft-gluon expression (spin-independent)
2913 //      =  2 : V -> q qbar (V = vector/axial vector colour singlet)
2914 //      =  3 : q -> q V
2915 //      =  4 : S -> q qbar (S = scalar/pseudoscalar colour singlet)
2916 //      =  5 : q -> q S
2917 //      =  6 : V -> ~q ~qbar (~q = squark)
2918 //      =  7 : ~q -> ~q V
2919 //      =  8 : S -> ~q ~qbar
2920 //      =  9 : ~q -> ~q S
2921 //      = 10 : chi -> q ~qbar (chi = neutralino/chargino)
2922 //      = 11 : ~q -> q chi
2923 //      = 12 : q -> ~q chi
2924 //      = 13 : ~g -> q ~qbar
2925 //      = 14 : ~q -> q ~g
2926 //      = 15 : q -> ~q ~g
2927 //      = 16 : (9/4)*(eikonal) for gg -> ~g ~g
2928 //      = 30 : Dv -> d qv     (Dv= hidden valley fermion, qv= valley scalar)
2929 //      = 31 : S  -> Dv Dvbar (S=scalar color singlet)
2930 // Note that the order of the decay products is important.
2931 // combi = 1 : pure non-gamma5, i.e. vector/scalar/...
2932 //       = 2 : pure gamma5, i.e. axial vector/pseudoscalar/....
2933 //       = 3 : mixture mix*(combi=1) + (1-mix)*(combi=2)
2934 //       = 4 : mixture (combi=1) +- (combi=2)
2935
2936 double TimeShower::calcMEcorr( int kind, int combiIn, double mixIn,
2937   double x1, double x2, double r1, double r2, double r3) {
2938
2939   // Frequent variable combinations.
2940   double x3     = 2. - x1 - x2;
2941   double x1s    = x1 * x1;
2942   double x2s    = x2 * x2;
2943   double x3s    = x3 * x3;
2944   double x1c    = x1 * x1s;
2945   double x2c    = x2 * x2s;
2946   double x3c    = x3 * x3s;
2947   double r1s    = r1 * r1;
2948   double r2s    = r2 * r2;
2949   double r1c    = r1 * r1s;
2950   double r2c    = r2 * r2s;
2951   double r1q    = r1s * r1s;
2952   double r2q    = r2s * r2s;
2953   double prop1  = 1. + r1s - r2s - x1;
2954   double prop2  = 1. + r2s - r1s - x2;
2955   double prop1s = prop1 * prop1;
2956   double prop2s = prop2 * prop2;
2957   double prop12 = prop1 * prop2;
2958   double prop13 = prop1 * x3;
2959   double prop23 = prop2 * x3;
2960
2961   // Special case: Hidden-Valley massive photon.
2962   double r3s    = r3 * r3;
2963   double prop3  = r3s - x3;
2964   double prop3s = prop3 * prop3;
2965   if (kind == 30) prop13 = prop1 * prop3;
2966
2967   // Check input values. Return zero outside allowed phase space.
2968   if (x1 - 2.*r1 < XMARGIN || prop1 < XMARGIN) return 0.;
2969   if (x2 - 2.*r2 < XMARGIN || prop2 < XMARGIN) return 0.;
2970   // Limits not worked out for r3 > 0.
2971   if (kind != 30 && kind != 31) {
2972     if (x1 + x2 - 1. - pow2(r1+r2) < XMARGIN) return 0.;
2973     // Note: equivalent rewritten form 4. * ( (1. - x1) * (1. - x2)
2974     // * (1. - r1s - r2s - x3) + r1s * (1. - x2s - x3) + r2s
2975     // * (1. - x1s - x3) - pow2(r1s - r2s) ) gives about same result.
2976     if ( (x1s - 4.*r1s) * (x2s - 4.*r2s)
2977       - pow2( 2. * (1. - x1 - x2 + r1s + r2s) + x1*x2 )
2978       < XMARGIN * (XMARGINCOMB + r1 + r2) ) return 0.;
2979   }
2980
2981   // Initial values; phase space.
2982   int combi   = max(1, min(4, combiIn) );
2983   double mix  = max(0., min(1., mixIn) );
2984   bool isSet1 = false;
2985   bool isSet2 = false;
2986   bool isSet4 = false;
2987   double ps = sqrtpos( pow2(1. - r1*r1 - r2*r2) - pow2(2. * r1 * r2) );
2988   double rLO = 0., rFO = 0., rLO1 = 0., rFO1 = 0., rLO2 = 0.,
2989     rFO2 = 0., rLO4 = 0., rFO4 = 0.;
2990   double offset = 0;
2991
2992   // Select which kind of ME to use.
2993   switch (kind) {
2994
2995     // case 1 is equal to default, i.e. eikonal expression.
2996
2997     // V -> q qbar (V = gamma*/Z0/W+-/...).
2998     case 2:
2999       if (combi == 1 || combi == 3) {
3000         rLO1 = ps*(2.-r1s-r1q+6.*r1*r2-r2s+2.*r1s*r2s-r2q)/2.;
3001         rFO1 = -(3.+6.*r1s+r1q-6.*r1*r2+6.*r1c*r2-2.*r2s-6.*r1s*r2s
3002         +6.*r1*r2c+r2q-3.*x1+6.*r1*r2*x1+2.*r2s*x1+x1s-2.*r1s*x1s
3003         +3.*r1s*x3+6.*r1*r2*x3-r2s*x3-2.*x1*x3-5.*r1s*x1*x3
3004         +r2s*x1*x3+x1s*x3-3.*x3s-3.*r1s*x3s+r2s*x3s
3005         +2.*x1*x3s+x3c-x2)
3006         /prop2s
3007         -2.*(-3.+r1s-6.*r1*r2+6.*r1c*r2+3.*r2s-4.*r1s*r2s
3008         +6.*r1*r2c+2.*x1+3.*r1s*x1+r2s*x1-x1s-r1s*x1s
3009         -r2s*x1s+4.*x3+2.*r1s*x3+3.*r1*r2*x3-r2s*x3-3.*x1*x3
3010         -2.*r1s*x1*x3+x1s*x3-x3s-r1s*x3s+r1*r2*x3s+x1*x3s)
3011         /prop12
3012         -(-1.+2.*r1s+r1q+6.*r1*r2+6.*r1c*r2-2.*r2s-6.*r1s*r2s
3013         +6.*r1*r2c+r2q-x1-2.*r1s*x1-6.*r1*r2*x1+8.*r2s*x1+x1s
3014         -2.*r2s*x1s-r1s*x3+r2s*x3-r1s*x1*x3+r2s*x1*x3+x1s*x3+x2)
3015         /prop1s;
3016         rFO1 = rFO1/2.;
3017         isSet1 = true;
3018       }
3019       if (combi == 2 || combi == 3) {
3020         rLO2 = ps*(2.-r1s-r1q-6.*r1*r2-r2s+2.*r1s*r2s-r2q)/2.;
3021         rFO2 = -(3.+6.*r1s+r1q+6.*r1*r2-6.*r1c*r2-2.*r2s-6.*r1s*r2s
3022         -6.*r1*r2c+r2q-3.*x1-6.*r1*r2*x1+2.*r2s*x1+x1s-2.*r1s*x1s
3023         +3.*r1s*x3-6.*r1*r2*x3-r2s*x3-2.*x1*x3-5.*r1s*x1*x3
3024         +r2s*x1*x3+x1s*x3-3.*x3s-3.*r1s*x3s+r2s*x3s+2.*x1*x3s+x3c-x2)
3025         /prop2s
3026         -2.*(-3+r1s+6.*r1*r2-6.*r1c*r2+3.*r2s-4.*r1s*r2s-6.*r1*r2c
3027         +2.*x1+3.*r1s*x1+r2s*x1-x1s-r1s*x1s-r2s*x1s+4.*x3+2.*r1s*x3
3028         -3.*r1*r2*x3-r2s*x3-3.*x1*x3-2.*r1s*x1*x3+x1s*x3-x3s-r1s*x3s
3029         -r1*r2*x3s+x1*x3s)
3030         /prop12
3031         -(-1.+2.*r1s+r1q-6.*r1*r2-6.*r1c*r2-2.*r2s-6.*r1s*r2s
3032         -6.*r1*r2c+r2q-x1-2.*r1s*x1+6.*r1*r2*x1+8.*r2s*x1+x1s
3033         -2.*r2s*x1s-r1s*x3+r2s*x3-r1s*x1*x3+r2s*x1*x3+x1s*x3+x2)
3034         /prop1s;
3035         rFO2 = rFO2/2.;
3036         isSet2 = true;
3037       }
3038       if (combi == 4) {
3039         rLO4 = ps*(2.-r1s-r1q-r2s+2.*r1s*r2s-r2q)/2.;
3040         rFO4 = (1.-r1q+6.*r1s*r2s-r2q+x1+3.*r1s*x1-9.*r2s*x1-3.*x1s
3041         -r1s*x1s+3.*r2s*x1s+x1c-x2-r1s*x2+r2s*x2-r1s*x1*x2+r2s*x1*x2
3042         +x1s*x2)
3043         /prop1s
3044         -2.*(1.+r1s+r2s-4.*r1s*r2s+r1s*x1+2.*r2s*x1-x1s-r2s*x1s
3045         +2.*r1s*x2+r2s*x2-3.*x1*x2+x1s*x2-x2s-r1s*x2s+x1*x2s)
3046         /prop12
3047         +(1.-r1q+6.*r1s*r2s-r2q-x1+r1s*x1-r2s*x1+x2-9.*r1s*x2
3048         +3.*r2s*x2+r1s*x1*x2-r2s*x1*x2-3.*x2s+3.*r1s*x2s-r2s*x2s
3049         +x1*x2s+x2c)
3050         /prop2s;
3051         rFO4 = rFO4/2.;
3052         isSet4 = true;
3053       }
3054       break;
3055
3056     // q -> q V.
3057     case 3:
3058       if (combi == 1 || combi == 3) {
3059         rLO1 = ps*(1.-2.*r1s+r1q+r2s-6.*r1*r2s+r1s*r2s-2.*r2q);
3060         rFO1 = -2.*(-1.+r1-2.*r1s+2.*r1c-r1q+pow5(r1)-r2s+r1*r2s
3061         -5.*r1s*r2s+r1c*r2s-2.*r1*r2q+2.*x1-2.*r1*x1+2.*r1s*x1
3062         -2.*r1c*x1+2.*r2s*x1+5.*r1*r2s*x1+r1s*r2s*x1+2.*r2q*x1
3063         -x1s+r1*x1s-r2s*x1s+3.*x2+4.*r1s*x2+r1q*x2+2.*r2s*x2
3064         +2.*r1s*r2s*x2-4.*x1*x2-2.*r1s*x1*x2-r2s*x1*x2+x1s*x2
3065         -2.*x2s-2.*r1s*x2s+x1*x2s)
3066         /prop23
3067         +(2.*r2s+6.*r1*r2s-6.*r1s*r2s+6.*r1c*r2s+2.*r2q+6.*r1*r2q
3068         -r2s*x1+r1s*r2s*x1-r2q*x1+x2-r1q*x2-3.*r2s*x2-6.*r1*r2s*x2
3069         +9.*r1s*r2s*x2-2.*r2q*x2-x1*x2+r1s*x1*x2-x2s-3.*r1s*x2s
3070         +2.*r2s*x2s+x1*x2s)
3071         /prop2s
3072         +(-4.-8.*r1s-4.*r1q+4.*r2s-4.*r1s*r2s+8.*r2q+9.*x1+10.*r1s*x1
3073         +r1q*x1-3.*r2s*x1+6.*r1*r2s*x1+r1s*r2s*x1-2.*r2q*x1-6.*x1s-
3074         2.*r1s*x1s+x1c+7.*x2+8.*r1s*x2+r1q*x2-7.*r2s*x2+6.*r1*r2s*x2
3075         +r1s*r2s*x2-2.*r2q*x2-9.*x1*x2-3.*r1s*x1*x2+2.*r2s*x1*x2
3076         +2.*x1s*x2-3.*x2s-r1s*x2s+2.*r2s*x2s+x1*x2s)
3077         /x3s;
3078         isSet1 = true;
3079       }
3080       if (combi == 2 || combi == 3) {
3081         rLO2 = ps*(1.-2.*r1s+r1q+r2s+6.*r1*r2s+r1s*r2s-2.*r2q);
3082         rFO2 = 2*(1.+r1+2.*r1s+2.*r1c+r1q+pow5(r1)+r2s+r1*r2s
3083         +5.*r1s*r2s+r1c*r2s-2.*r1*r2q-2.*x1-2.*r1*x1-2.*r1s*x1
3084         -2.*r1c*x1-2.*r2s*x1+5.*r1*r2s*x1-r1s*r2s*x1-2.*r2q*x1+x1s
3085         +r1*x1s+r2s*x1s-3.*x2-4.*r1s*x2-r1q*x2-2.*r2s*x2
3086         -2.*r1s*r2s*x2+4.*x1*x2+2.*r1s*x1*x2+r2s*x1*x2-x1s*x2
3087         +2.*x2s+2.*r1s*x2s-x1*x2s)
3088         /prop23
3089         +(2.*r2s-6.*r1*r2s-6.*r1s*r2s-6.*r1c*r2s+2.*r2q-6.*r1*r2q
3090         -r2s*x1+r1s*r2s*x1-r2q*x1+x2-r1q*x2-3.*r2s*x2+6.*r1*r2s*x2
3091         +9.*r1s*r2s*x2-2.*r2q*x2-x1*x2+r1s*x1*x2-x2s-3.*r1s*x2s
3092         +2.*r2s*x2s+x1*x2s)
3093         /prop2s
3094         +(-4.-8.*r1s-4.*r1q+4.*r2s-4.*r1s*r2s+8.*r2q+9.*x1+10.*r1s*x1
3095         +r1q*x1-3.*r2s*x1-6.*r1*r2s*x1+r1s*r2s*x1-2.*r2q*x1-6.*x1s
3096         -2.*r1s*x1s+x1c+7.*x2+8.*r1s*x2+r1q*x2-7.*r2s*x2-6.*r1*r2s*x2
3097         +r1s*r2s*x2-2.*r2q*x2-9.*x1*x2-3.*r1s*x1*x2+2.*r2s*x1*x2
3098         +2.*x1s*x2-3.*x2s-r1s*x2s+2.*r2s*x2s+x1*x2s)
3099         /x3s;
3100         isSet2 = true;
3101       }
3102       if (combi == 4) {
3103         rLO4 = ps*(1.-2.*r1s+r1q+r2s+r1s*r2s-2.*r2q);
3104         rFO4 = 2*(1.+2.*r1s+r1q+r2s+5.*r1s*r2s-2.*x1-2.*r1s*x1
3105         -2.*r2s*x1-r1s*r2s*x1-2.*r2q*x1+x1s+r2s*x1s-3.*x2-4.*r1s*x2
3106         -r1q*x2-2.*r2s*x2-2.*r1s*r2s*x2+4.*x1*x2+2.*r1s*x1*x2+r2s*x1*x2
3107         -x1s*x2+2.*x2s+2.*r1s*x2s-x1*x2s)
3108         /prop23
3109         +(2.*r2s-6.*r1s*r2s+2.*r2q-r2s*x1+r1s*r2s*x1-r2q*x1+x2-r1q*x2
3110         -3.*r2s*x2+9.*r1s*r2s*x2-2.*r2q*x2-x1*x2+r1s*x1*x2-x2s-3.*r1s*x2s
3111         +2.*r2s*x2s+x1*x2s)
3112         /prop2s
3113         +(-4.-8.*r1s-4.*r1q+4.*r2s-4.*r1s*r2s+8.*r2q+9.*x1+10.*r1s*x1
3114         +r1q*x1-3.*r2s*x1+r1s*r2s*x1-2.*r2q*x1-6.*x1s-2.*r1s*x1s+x1c
3115         +7.*x2+8.*r1s*x2+r1q*x2-7.*r2s*x2+r1s*r2s*x2-2.*r2q*x2-9.*x1*x2
3116         -3.*r1s*x1*x2+2.*r2s*x1*x2+2.*x1s*x2-3.*x2s-r1s*x2s+2.*r2s*x2s
3117         +x1*x2s)
3118         /x3s;
3119         isSet4 = true;
3120       }
3121       break;
3122
3123     // S -> q qbar    (S = h0/H0/A0/H+-/...).
3124     case 4:
3125       if (combi == 1 || combi == 3) {
3126         rLO1 = ps*(1.-r1s-r2s-2.*r1*r2);
3127         rFO1 = -(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q+x1
3128         -r1s*x1+2.*r1*r2*x1+3.*r2s*x1+x2+r1s*x2-r2s*x2-x1*x2)
3129         /prop1s
3130         -2.*(r1s+r1q-2.*r1c*r2+r2s-6.*r1s*r2s-2.*r1*r2c+r2q-r1s*x1
3131         +r1*r2*x1+2.*r2s*x1+2.*r1s*x2+r1*r2*x2-r2s*x2-x1*x2)
3132         /prop12
3133         -(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q+x1-r1s*x1
3134         +r2s*x1+x2+3.*r1s*x2+2.*r1*r2*x2-r2s*x2-x1*x2)
3135         /prop2s;
3136         isSet1 = true;
3137       }
3138       if (combi == 2 || combi == 3) {
3139         rLO2 = ps*(1.-r1s-r2s+2.*r1*r2);
3140         rFO2 = -(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1
3141         -r1s*x1-2.*r1*r2*x1+3.*r2s*x1+x2+r1s*x2-r2s*x2-x1*x2)
3142         /prop1s
3143         -(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1
3144         -r1s*x1+r2s*x1+x2+3.*r1s*x2-2.*r1*r2*x2-r2s*x2-x1*x2)
3145         /prop2s
3146         +2.*(-r1s-r1q-2.*r1c*r2-r2s+6.*r1s*r2s-2.*r1*r2c-r2q+r1s*x1
3147         +r1*r2*x1-2.*r2s*x1-2.*r1s*x2+r1*r2*x2+r2s*x2+x1*x2)
3148         /prop12;
3149         isSet2 = true;
3150       }
3151       if (combi == 4) {
3152         rLO4 = ps*(1.-r1s-r2s);
3153         rFO4 = -(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+3.*r2s*x1+x2
3154         +r1s*x2-r2s*x2-x1*x2)
3155         /prop1s
3156         -2.*(r1s+r1q+r2s-6.*r1s*r2s+r2q-r1s*x1
3157         +2.*r2s*x1+2.*r1s*x2-r2s*x2-x1*x2)
3158         /prop12
3159         -(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+r2s*x1
3160         +x2+3.*r1s*x2-r2s*x2-x1*x2)
3161         /prop2s;
3162         isSet4 = true;
3163       }
3164       break;
3165
3166     // q -> q S.
3167     case 5:
3168       if (combi == 1 || combi == 3) {
3169         rLO1 = ps*(1.+r1s-r2s+2.*r1);
3170         rFO1 = (4.-4.*r1s+4.*r2s-3.*x1-2.*r1*x1+r1s*x1-r2s*x1-5.*x2
3171         -2.*r1*x2+r1s*x2-r2s*x2+x1*x2+x2s)
3172         /x3s
3173         -2.*(3.-r1-5.*r1s-r1c+3.*r2s+r1*r2s-2.*x1-r1*x1
3174         +r1s*x1-4.*x2+2.*r1s*x2-r2s*x2+x1*x2+x2s)
3175         /prop23
3176         +(2.-2.*r1-6.*r1s-2.*r1c+2.*r2s-2.*r1*r2s-x1+r1s*x1
3177         -r2s*x1-3.*x2+2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s)
3178         /prop2s;
3179         isSet1 = true;
3180       }
3181       if (combi == 2 || combi == 3) {
3182         rLO2 = ps*(1.+r1s-r2s-2.*r1);
3183         rFO2 = (4.-4.*r1s+4.*r2s-3.*x1+2.*r1*x1+r1s*x1-r2s*x1-5.*x2
3184         +2.*r1*x2+r1s*x2-r2s*x2+x1*x2+x2s)
3185         /x3s
3186         -2.*(3.+r1-5.*r1s+r1c+3.*r2s-r1*r2s-2.*x1+r1*x1
3187         +r1s*x1-4.*x2+2.*r1s*x2-r2s*x2+x1*x2+x2s)
3188         /prop23
3189         +(2.+2.*r1-6.*r1s+2.*r1c+2.*r2s+2.*r1*r2s-x1+r1s*x1
3190         -r2s*x1-3.*x2-2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s)
3191         /prop2s;
3192         isSet2 = true;
3193       }
3194       if (combi == 4) {
3195         rLO4 = ps*(1.+r1s-r2s);
3196         rFO4 = (4.-4.*r1s+4.*r2s-3.*x1+r1s*x1-r2s*x1-5.*x2+r1s*x2
3197         -r2s*x2+x1*x2+x2s)
3198         /x3s
3199         -2.*(3.-5.*r1s+3.*r2s-2.*x1+r1s*x1-4.*x2+2.*r1s*x2
3200         -r2s*x2+x1*x2+x2s)
3201         /prop23
3202         +(2.-6.*r1s+2.*r2s-x1+r1s*x1-r2s*x1-3.*x2+3.*r1s*x2
3203         -r2s*x2+x1*x2+x2s)
3204         /prop2s;
3205         isSet4 = true;
3206       }
3207       break;
3208
3209     // V -> ~q ~qbar  (~q = squark).
3210     case 6:
3211       rLO1 = ps*(1.-2.*r1s+r1q-2.*r2s-2.*r1s*r2s+r2q);
3212       rFO1 = 2.*3.+(1.+r1s+r2s-x1)*(4.*r1s-x1s)
3213       /prop1s
3214       +2.*(-1.-3.*r1s-r2s+x1+x1s*0.5+x2-x1*x2*0.5)
3215       /prop1
3216       +(1.+r1s+r2s-x2)*(4.*r2s-x2s)
3217       /prop2s
3218       +2.*(-1.-r1s-3.*r2s+x1+x2-x1*x2*0.5+x2s*0.5)
3219       /prop2
3220       -(-4.*r1s-4.*r1q-4.*r2s-8.*r1s*r2s-4.*r2q+2.*x1+6.*r1s*x1
3221       +6.*r2s*x1-2.*x1s+2.*x2+6.*r1s*x2+6.*r2s*x2-4.*x1*x2
3222       -2.*r1s*x1*x2-2.*r2s*x1*x2+x1s*x2-2.*x2s+x1*x2s)
3223       /prop12;
3224       isSet1 = true;
3225       break;
3226
3227     // ~q -> ~q V.
3228     case 7:
3229       rLO1 = ps*(1.-2.*r1s+r1q-2.*r2s-2.*r1s*r2s+r2q);
3230       rFO1 = 16.*r2s-8.*(4.*r2s+2.*r2s*x1+x2+r1s*x2+r2s*x2-x1*x2
3231       -2.*x2s)
3232       /(3.*prop2)
3233       +8.*(1.+r1s+r2s-x2)*(4.*r2s-x2s)
3234       /(3.*prop2s)
3235       +8.*(x1+x2)*(-1.-2.*r1s-r1q-2.*r2s+2.*r1s*r2s-r2q+2.*x1
3236       +2.*r1s*x1+2.*r2s*x1-x1s+2.*x2+2.*r1s*x2+2.*r2s*x2-2.*x1*x2-x2s)
3237       /(3.*x3s)
3238       +8.*(-1.-r1s+r2s-x1)*(2.*r2s*x1+x2+r1s*x2+r2s*x2-x1*x2-x2s)
3239       /(3.*prop2*x3)
3240       -8.*(1.+2.*r1s+r1q+2.*r2s-2.*r1s*r2s+r2q-2.*x1-2.*r1s*x1
3241       -4.*r2s*x1+x1s-3.*x2-3.*r1s*x2-3.*r2s*x2+3.*x1*x2+2.*x2s)
3242       /(3.*x3);
3243       rFO1 = 3.*rFO1/8.;
3244       isSet1 = true;
3245       break;
3246
3247     // S -> ~q ~qbar.
3248     case 8:
3249       rLO1 = ps;
3250       rFO1 = (-1.-2.*r1s-r1q-2.*r2s+2.*r1s*r2s-r2q+2.*x1+2.*r1s*x1
3251       +2.*r2s*x1-x1s-r2s*x1s+2.*x2+2.*r1s*x2+2.*r2s*x2-3.*x1*x2
3252       -r1s*x1*x2-r2s*x1*x2+x1s*x2-x2s-r1s*x2s+x1*x2s)
3253       /(prop1s*prop2s);
3254       rFO1 = 2.*rFO1;
3255       isSet1 = true;
3256       break;
3257
3258     // ~q -> ~q S.
3259     case 9:
3260       rLO1 = ps;
3261       rFO1 = (-1.-r1s-r2s+x2)
3262       /prop2s
3263       +(1.+r1s-r2s+x1)
3264       /prop23
3265       -(x1+x2)
3266       /x3s;
3267       isSet1 = true;
3268       break;
3269
3270     // chi -> q ~qbar   (chi = neutralino/chargino).
3271     case 10:
3272       if (combi == 1 || combi == 3) {
3273         rLO1 = ps*(1.+r1s-r2s+2.*r1);
3274         rFO1 = (2.*r1+x1)*(-1.-r1s-r2s+x1)
3275         /prop1s
3276         +2.*(-1.-r1s-2.*r1c-r2s-2.*r1*r2s+3.*x1*0.5+r1*x1
3277         -r1s*x1*0.5-r2s*x1*0.5+x2+r1*x2+r1s*x2-x1*x2*0.5)
3278         /prop12
3279         +(2.-2.*r1-6.*r1s-2.*r1c+2.*r2s-2.*r1*r2s-x1+r1s*x1
3280         -r2s*x1-3.*x2+2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s)
3281         /prop2s;
3282         isSet1 = true;
3283       }
3284       if (combi == 2 || combi == 3) {
3285         rLO2 = ps*(1.-2.*r1+r1s-r2s);
3286         rFO2 = (2.*r1-x1)*(1.+r1s+r2s-x1)
3287         /prop1s
3288         +2.*(-1.-r1s+2.*r1c-r2s+2.*r1*r2s+3.*x1*0.5-r1*x1
3289         -r1s*x1*0.5-r2s*x1*0.5+x2-r1*x2+r1s*x2-x1*x2*0.5)
3290         /prop12
3291         +(2.+2.*r1-6.*r1s+2.*r1c+2.*r2s+2.*r1*r2s-x1+r1s*x1
3292         -r2s*x1-3.*x2-2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s)/
3293         prop2s;
3294         isSet2 = true;
3295       }
3296       if (combi == 4) {
3297         rLO4 = ps*(1.+r1s-r2s);
3298         rFO4 = x1*(-1.-r1s-r2s+x1)
3299         /prop1s
3300         +2.*(-1.-r1s-r2s+3.*x1*0.5-r1s*x1*0.5-r2s*x1*0.5
3301         +x2+r1s*x2-x1*x2*0.5)
3302         /prop12
3303         +(2.-6.*r1s+2.*r2s-x1+r1s*x1-r2s*x1-3.*x2+3.*r1s*x2
3304         -r2s*x2+x1*x2+x2s)
3305         /prop2s;
3306         isSet4 = true;
3307       }
3308       break;
3309
3310     // ~q -> q chi.
3311     case 11:
3312       if (combi == 1 || combi == 3) {
3313         rLO1 = ps*(1.-pow2(r1+r2));
3314         rFO1 = (1.+r1s+2.*r1*r2+r2s-x1-x2)*(x1+x2)
3315         /x3s
3316         -(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q+x1
3317         -r1s*x1+r2s*x1+x2+3.*r1s*x2+2.*r1*r2*x2-r2s*x2-x1*x2)
3318         /prop2s
3319         +(-1.-2.*r1s-r1q-2.*r1*r2-2.*r1c*r2+2.*r1*r2c+r2q+x1+r1s*x1
3320         -2.*r1*r2*x1-3.*r2s*x1+2.*r1s*x2-2.*r2s*x2+x1*x2+x2s)
3321         /prop23;
3322         isSet1 = true;
3323       }
3324       if (combi == 2 || combi == 3) {
3325         rLO2 = ps*(1.-pow2(r1-r2));
3326         rFO2 = (1.+r1s-2.*r1*r2+r2s-x1-x2)*(x1+x2)
3327         /x3s
3328         -(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1
3329         -r1s*x1+r2s*x1+x2+3.*r1s*x2-2.*r1*r2*x2-r2s*x2-x1*x2)
3330         /prop2s
3331         +(-1.-2.*r1s-r1q+2.*r1*r2+2.*r1c*r2-2.*r1*r2c+r2q+x1+r1s*x1
3332         +2.*r1*r2*x1-3.*r2s*x1+2.*r1s*x2-2.*r2s*x2+x1*x2+x2s)
3333         /prop23;
3334         isSet2 = true;
3335       }
3336       if (combi == 4) {
3337         rLO4 = ps*(1.-r1s-r2s);
3338         rFO4 = (1.+r1s+r2s-x1-x2)*(x1+x2)
3339         /x3s
3340         -(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+r2s*x1+x2
3341         +3.*r1s*x2-r2s*x2-x1*x2)
3342         /prop2s
3343         +(-1.-2.*r1s-r1q+r2q+x1+r1s*x1-3.*r2s*x1
3344         +2.*r1s*x2-2.*r2s*x2+x1*x2+x2s)
3345         /prop23;
3346         isSet4 = true;
3347       }
3348       break;
3349
3350     // q -> ~q chi.
3351     case 12:
3352       if (combi == 1 || combi == 3) {
3353         rLO1 = ps*(1.-r1s+r2s+2.*r2);
3354         rFO1 = (2.*r2+x2)*(-1.-r1s-r2s+x2)
3355         /prop2s
3356         +(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1-2.*r2*x1+r2s*x1+x1s
3357         -3.*x2-r1s*x2-2.*r2*x2+r2s*x2+x1*x2)
3358         /x3s
3359         +2.*(-1.-r1s+r2+r1s*r2-r2s-r2c+x1+r2*x1+r2s*x1+2.*x2
3360         +r1s*x2-x1*x2*0.5-x2s*0.5)
3361         /prop23;
3362         isSet1 = true;
3363       }
3364       if (combi == 2 || combi == 3) {
3365         rLO2 = ps*(1.-r1s+r2s-2.*r2);
3366         rFO2 = (2.*r2-x2)*(1.+r1s+r2s-x2)
3367         /prop2s
3368         +(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1+2.*r2*x1+r2s*x1+x1s
3369         -3.*x2-r1s*x2+2.*r2*x2+r2s*x2+x1*x2)
3370         /x3s
3371         +2.*(-1.-r1s-r2-r1s*r2-r2s+r2c+x1-r2*x1+r2s*x1+2.*x2
3372         +r1s*x2-x1*x2*0.5-x2s*0.5)
3373         /prop23;
3374         isSet2 = true;
3375       }
3376       if (combi == 4) {
3377         rLO4 = ps*(1.-r1s+r2s);
3378         rFO4 = x2*(-1.-r1s-r2s+x2)
3379         /prop2s
3380         +(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1+r2s*x1+x1s
3381         -3.*x2-r1s*x2+r2s*x2+x1*x2)
3382         /x3s
3383         +2.*(-1.-r1s-r2s+x1+r2s*x1+2.*x2
3384         +r1s*x2-x1*x2*0.5-x2s*0.5)
3385         /prop23;
3386         isSet4 = true;
3387       }
3388       break;
3389
3390     // ~g -> q ~qbar.
3391     case 13:
3392       if (combi == 1 || combi == 3) {
3393         rLO1 = ps*(1.+r1s-r2s+2.*r1);
3394         rFO1 = 4.*(2.*r1+x1)*(-1.-r1s-r2s+x1)
3395         /(3.*prop1s)
3396         -(-1.-r1s-2.*r1c-r2s-2.*r1*r2s+3.*x1*0.5+r1*x1-r1s*x1*0.5
3397         -r2s*x1*0.5+x2+r1*x2+r1s*x2-x1*x2*0.5)
3398         /(3.*prop12)
3399         +3.*(-1.+r1-r1s-r1c-r2s+r1*r2s+2.*x1+r2s*x1-x1s*0.5+x2+r1*x2
3400         +r1s*x2-x1*x2*0.5)
3401         /prop13
3402         +3.*(4.-4.*r1s+4.*r2s-3.*x1-2.*r1*x1+r1s*x1-r2s*x1-5.*x2
3403         -2.*r1*x2+r1s*x2-r2s*x2+x1*x2+x2s)
3404         /x3s
3405         -3.*(3.-r1-5.*r1s-r1c+3.*r2s+r1*r2s-2.*x1-r1*x1+r1s*x1
3406         -4.*x2+2.*r1s*x2-r2s*x2+x1*x2+x2s)
3407         /prop23
3408         +4.*(2.-2.*r1-6.*r1s-2.*r1c+2.*r2s-2.*r1*r2s-x1+r1s*x1-r2s*x1
3409         -3.*x2+2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s)
3410         /(3.*prop2s);
3411         rFO1 = 3.*rFO1/4.;
3412         isSet1 = true;
3413       }
3414       if (combi == 2 || combi == 3) {
3415         rLO2 = ps*(1.+r1s-r2s-2.*r1);
3416         rFO2 = 4.*(2.*r1-x1)*(1.+r1s+r2s-x1)
3417         /(3.*prop1s)
3418         +3.*(-1.-r1-r1s+r1c-r2s-r1*r2s+2.*x1+r2s*x1-x1s*0.5
3419         +x2-r1*x2+r1s*x2-x1*x2*0.5)
3420         /prop13
3421         +(2.+2.*r1s-4.*r1c+2.*r2s-4.*r1*r2s-3.*x1+2.*r1*x1
3422         +r1s*x1+r2s*x1-2.*x2+2.*r1*x2-2.*r1s*x2+x1*x2)
3423         /(6.*prop12)
3424         +3.*(4.-4.*r1s+4.*r2s-3.*x1+2.*r1*x1+r1s*x1-r2s*x1-5.*x2
3425         +2.*r1*x2+r1s*x2-r2s*x2+x1*x2+x2s)
3426         /x3s
3427         -3.*(3.+r1-5.*r1s+r1c+3.*r2s-r1*r2s-2.*x1+r1*x1+r1s*x1-4.*x2
3428         +2.*r1s*x2-r2s*x2+x1*x2+x2s)
3429         /prop23
3430         +4.*(2.+2.*r1-6.*r1s+2.*r1c+2.*r2s+2.*r1*r2s-x1+r1s*x1-r2s*x1
3431         -3.*x2-2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s)
3432         /(3.*prop2s);
3433         rFO2 = 3.*rFO2/4.;
3434         isSet2 = true;
3435       }
3436       if (combi == 4) {
3437         rLO4 = ps*(1.+r1s-r2s);
3438         rFO4 = 8.*x1*(-1.-r1s-r2s+x1)
3439         /(3.*prop1s)
3440         +6.*(-1-r1s-r2s+2.*x1+r2s*x1-x1s*0.5+x2+r1s*x2-x1*x2*0.5)
3441         /prop13
3442         +(2.+2.*r1s+2.*r2s-3.*x1+r1s*x1+r2s*x1-2.*x2-2.*r1s*x2+x1*x2)
3443         /(3.*prop12)
3444         +6.*(4.-4.*r1s+4.*r2s-3.*x1+r1s*x1-r2s*x1-5.*x2+r1s*x2-r2s*x2
3445         +x1*x2+x2s)
3446         /x3s
3447         -6.*(3.-5.*r1s+3.*r2s-2.*x1+r1s*x1-4.*x2+2.*r1s*x2-r2s*x2+x1*x2+x2s)
3448         /prop23
3449         +8.*(2.-6.*r1s+2.*r2s-x1+r1s*x1-r2s*x1-3.*x2+3.*r1s*x2-r2s*x2
3450         +x1*x2+x2s)
3451         /(3.*prop2s);
3452         rFO4 = 3.*rFO4/8.;
3453         isSet4 = true;
3454       }
3455       break;
3456
3457     // ~q -> q ~g.
3458     case 14:
3459       if (combi == 1 || combi == 3) {
3460         rLO1 = ps*(1.-r1s-r2s-2.*r1*r2);
3461         rFO1 = 64.*(1.+r1s+2.*r1*r2+r2s-x1-x2)*(x1+x2)
3462         /(9.*x3s)
3463         -16.*(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q
3464         +x1-r1s*x1+2.*r1*r2*x1+3.*r2s*x1+x2+r1s*x2-r2s*x2-x1*x2)
3465         /prop1s
3466         -16.*(r1s+r1q-2.*r1c*r2+r2s-6.*r1s*r2s-2.*r1*r2c+r2q-r1s*x1
3467         +r1*r2*x1+2.*r2s*x1+2.*r1s*x2+r1*r2*x2-r2s*x2-x1*x2)
3468         /prop12
3469         -64.*(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q+x1
3470         -r1s*x1+r2s*x1+x2+3.*r1s*x2+2.*r1*r2*x2-r2s*x2-x1*x2)
3471         /(9.*prop2s)
3472         +8.*(-1.+r1q-2.*r1*r2+2.*r1c*r2-2.*r2s-2.*r1*r2c-r2q-2.*r1s*x1
3473         +2.*r2s*x1+x1s+x2-3.*r1s*x2-2.*r1*r2*x2+r2s*x2+x1*x2)
3474         /prop13
3475         -8.*(-1.-2.*r1s-r1q-2.*r1*r2-2.*r1c*r2+2.*r1*r2c+r2q+x1+r1s*x1
3476         -2.*r1*r2*x1-3.*r2s*x1+2.*r1s*x2-2.*r2s*x2+x1*x2+x2s)
3477         /(9.*prop23);
3478         rFO1 = 9.*rFO1/64.;
3479         isSet1 = true;
3480       }
3481       if (combi == 2 || combi == 3) {
3482         rLO2 = ps*(1.-r1s-r2s+2.*r1*r2);
3483         rFO2 = 64.*(1.+r1s-2.*r1*r2+r2s-x1-x2)*(x1+x2)
3484         /(9.*x3s)
3485         -16.*(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1
3486         -r1s*x1-2.*r1*r2*x1+3.*r2s*x1+x2+r1s*x2-r2s*x2-x1*x2)
3487         /prop1s
3488         -64.*(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1
3489         -r1s*x1+r2s*x1+x2+3.*r1s*x2-2.*r1*r2*x2-r2s*x2-x1*x2)
3490         /(9.*prop2s)
3491         +16.*(-r1s-r1q-2.*r1c*r2-r2s+6.*r1s*r2s-2.*r1*r2c-r2q+r1s*x1
3492         +r1*r2*x1-2.*r2s*x1-2.*r1s*x2+r1*r2*x2+r2s*x2+x1*x2)
3493         /prop12
3494         +8.*(-1.+r1q+2.*r1*r2-2.*r1c*r2-2.*r2s+2.*r1*r2c-r2q-2.*r1s*x1
3495         +2.*r2s*x1+x1s+x2-3.*r1s*x2+2.*r1*r2*x2+r2s*x2+x1*x2)
3496         /prop13
3497         -8.*(-1.-2.*r1s-r1q+2.*r1*r2+2.*r1c*r2-2.*r1*r2c+r2q+x1+r1s*x1+
3498         2.*r1*r2*x1-3.*r2s*x1+2.*r1s*x2-2.*r2s*x2+x1*x2+x2s)
3499         /(9.*prop23);
3500         rFO2 = 9.*rFO2/64.;
3501         isSet2 = true;
3502       }
3503       if (combi == 4) {
3504         rLO4 = ps*(1.-r1s-r2s);
3505         rFO4 = 128.*(1.+r1s+r2s-x1-x2)*(x1+x2)
3506         /(9.*x3s)
3507         -32*(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+3.*r2s*x1+x2
3508         +r1s*x2-r2s*x2-x1*x2)
3509         /prop1s
3510         -32.*(r1s+r1q+r2s-6.*r1s*r2s+r2q-r1s*x1+2.*r2s*x1+2.*r1s*x2
3511         -r2s*x2-x1*x2)
3512         /prop12
3513         -128.*(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+r2s*x1+x2+3.*r1s*x2
3514         -r2s*x2-x1*x2)
3515         /(9.*prop2s)
3516         +16.*(-1.+r1q-2.*r2s-r2q-2.*r1s*x1+2.*r2s*x1+x1s
3517         +x2-3.*r1s*x2+r2s*x2+x1*x2)
3518         /prop13
3519         -16.*(-1.-2.*r1s-r1q+r2q+x1+r1s*x1-3.*r2s*x1
3520         +2.*r1s*x2-2.*r2s*x2+x1*x2+x2s)
3521         /(9.*prop23);
3522         rFO4 = 9.*rFO4/128.;
3523         isSet4 = true;
3524       }
3525       break;
3526
3527     // q -> ~q ~g.
3528     case 15:
3529       if (combi == 1 || combi == 3) {
3530         rLO1 = ps*(1.-r1s+r2s+2.*r2);
3531         rFO1 = 32*(2.*r2+x2)*(-1.-r1s-r2s+x2)
3532         /(9.*prop2s)
3533         +8.*(-1.-r1s-2.*r1s*r2-r2s-2.*r2c+x1+r2*x1+r2s*x1
3534         +3.*x2*0.5-r1s*x2*0.5+r2*x2-r2s*x2*0.5-x1*x2*0.5)
3535         /prop12
3536         +8.*(2.+2.*r1s-2.*r2-2.*r1s*r2-6.*r2s-2.*r2c-3.*x1-r1s*x1
3537         +2.*r2*x1+3.*r2s*x1+x1s-x2-r1s*x2+r2s*x2+x1*x2)
3538         /prop1s
3539         +32.*(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1-2.*r2*x1+r2s*x1+x1s
3540         -3.*x2-r1s*x2-2.*r2*x2+r2s*x2+x1*x2)
3541         /(9.*x3s)
3542         -8.*(3.+3.*r1s-r2+r1s*r2-5.*r2s-r2c-4.*x1-r1s*x1
3543         +2.*r2s*x1+x1s-2.*x2-r2*x2+r2s*x2+x1*x2)
3544         /prop13
3545         -8.*(-1.-r1s+r2+r1s*r2-r2s-r2c+x1+r2*x1+r2s*x1+2.*x2+r1s*x2
3546         -x1*x2*0.5-x2s*0.5)
3547         /(9.*prop23);
3548         rFO1 = 9.*rFO1/32.;
3549         isSet1 = true;
3550       }
3551       if (combi == 2 || combi == 3) {
3552         rLO2 = ps*(1.-r1s+r2s-2.*r2);
3553         rFO2 = 32*(2.*r2-x2)*(1.+r1s+r2s-x2)
3554         /(9.*prop2s)
3555         +8.*(-1.-r1s+2.*r1s*r2-r2s+2.*r2c+x1-r2*x1+r2s*x1
3556         +3.*x2*0.5-r1s*x2*0.5-r2*x2-r2s*x2*0.5-x1*x2*0.5)
3557         /prop12
3558         +8.*(2.+2.*r1s+2.*r2+2.*r1s*r2-6.*r2s+2.*r2c-3.*x1-r1s*x1
3559         -2.*r2*x1+3.*r2s*x1+x1s-x2-r1s*x2+r2s*x2+x1*x2)
3560         /prop1s
3561         -8.*(3.+3.*r1s+r2-r1s*r2-5.*r2s+r2c-4.*x1-r1s*x1+2.*r2s*x1+x1s
3562         -2.*x2+r2*x2+r2s*x2+x1*x2)
3563         /prop13
3564         +32*(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1+2.*r2*x1+r2s*x1
3565         +x1s-3.*x2-r1s*x2+2.*r2*x2+r2s*x2+x1*x2)
3566         /(9.*x3s)
3567         -8.*(-1.-r1s-r2-r1s*r2-r2s+r2c+x1-r2*x1+r2s*x1+2.*x2+r1s*x2
3568         -x1*x2*0.5-x2s*0.5)
3569         /(9.*prop23);
3570         rFO2 = 9.*rFO2/32.;
3571         isSet2 = true;
3572       }
3573       if (combi == 4) {
3574         rLO4 = ps*(1.-r1s+r2s);
3575         rFO4 = 64.*x2*(-1.-r1s-r2s+x2)
3576         /(9.*prop2s)
3577         +16.*(-1.-r1s-r2s+x1+r2s*x1+3.*x2*0.5-r1s*x2*0.5
3578         -r2s*x2*0.5-x1*x2*0.5)
3579         /prop12
3580         -16.*(3.+3.*r1s-5.*r2s-4.*x1-r1s*x1+2.*r2s*x1+x1s-2.*x2+r2s*x2
3581         +x1*x2)
3582         /prop13
3583         +64.*(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1+r2s*x1+x1s-3.*x2
3584         -r1s*x2+r2s*x2+x1*x2)
3585         /(9.*x3s)
3586         +16.*(2.+2.*r1s-6.*r2s-3.*x1-r1s*x1+3.*r2s*x1+x1s
3587         -x2-r1s*x2+r2s*x2+x1*x2)
3588         /prop1s
3589         -16.*(-1.-r1s-r2s+x1+r2s*x1+2.*x2+r1s*x2-x1*x2*0.5-x2s*0.5)
3590         /(9.*prop23);
3591         rFO4 = 9.*rFO4/64.;
3592         isSet4 = true;
3593       }
3594       break;
3595
3596     // g -> ~g ~g. Use (9/4)*eikonal. May be changed in the future.
3597     case 16:
3598       rLO = ps;
3599       if      (combi == 2) offset = x3s;
3600       else if (combi == 3) offset = mix * x3s;
3601       else if (combi == 4) offset = 0.5 * x3s;
3602       rFO = ps * 4.5 * ( (x1+x2-1.+offset-r1s-r2s)/prop12
3603       - r1s/prop2s - r2s/prop1s );
3604       break;
3605
3606     // Dv -> qv d.
3607     case 30:
3608       rLO = ps*(1.-r1s+r2s+2.*r2);
3609       rFO = ( 0.5*r3s + 2.*r1q + 0.5*r2s*r3s + r2*r3s - 2.*r1s
3610              - 0.5*r1s*r3s - 2.*r1s*r2s - 4.*r1s*r2 ) / prop2s
3611           + ( -2. + 2.*r2q + 2.*r1q + 2.*r2s*r3s - 4.*r2 + 2.*r2*r3s
3612              + 4.*r2*r2s - 4.*r1s*r2s - 4.*r1s*r2 ) /prop23
3613           + ( -2. - 0.5*r3s - 2.*r2s - 4.*r2 + 2.*r1s ) / prop2
3614           + ( -2. - r3s - 2.*r2s - r2s*r3s - 4.*r2 - 2.*r2*r3s
3615              + 2.*r1s + r1s*r3s ) / prop3s
3616           + ( -1. - r3s - r2s - 4.*r2 + r1s - x2 ) / prop3
3617           + 1.;
3618       break;
3619
3620     // S -> Dv Dvbar
3621     case 31:
3622       rLO = ps*(1.-4.*r1s);
3623       rFO = (r3s + 2.*r1s) * (-1. + 4.*r1s) * (1./prop1s + 1./prop2s)
3624           + (-1. + 8.*r1s - x2) / prop1
3625           + (-1. + 8.*r1s - x1) / prop2
3626           + 2. * (1. - 6.*r1s + 8.*r1q + 4.*r3s*r1s) / prop12
3627           + 2.;
3628       break;
3629
3630     // Eikonal expression for kind == 1; also acts as default.
3631     default:
3632       rLO = ps;
3633       if      (combi == 2) offset = x3s;
3634       else if (combi == 3) offset = mix * x3s;
3635       else if (combi == 4) offset = 0.5 * x3s;
3636       rFO = ps * 2. * ( (x1+x2-1.+offset-r1s-r2s)/prop12
3637       - r1s/prop2s - r2s/prop1s );
3638       break;
3639
3640   // End of ME cases.
3641   }
3642
3643   // Find relevant leading and first order expressions.
3644   if      (combi == 1 && isSet1) {
3645     rLO = rLO1;
3646     rFO = rFO1; }
3647   else if (combi == 2 && isSet2) {
3648     rLO = rLO2;
3649     rFO = rFO2; }
3650   else if (combi == 3 && isSet1 && isSet2) {
3651     rLO = mix * rLO1 + (1.-mix) * rLO2;
3652     rFO = mix * rFO1 + (1.-mix) * rFO2; }
3653   else if (isSet4) {
3654     rLO = rLO4;
3655     rFO = rFO4; }
3656   else if (combi == 4 && isSet1 && isSet2) {
3657     rLO = 0.5 * (rLO1 + rLO2);
3658     rFO = 0.5 * (rFO1 + rFO2); }
3659   else if (isSet1) {
3660     rLO = rLO1;
3661     rFO = rFO1; }
3662
3663   // Return ratio of first to leading order cross section.
3664   return rFO / rLO;
3665 }
3666
3667 //--------------------------------------------------------------------------
3668
3669 // Find coefficient of azimuthal asymmetry from gluon polarization.
3670
3671 void TimeShower::findAsymPol( Event& event, TimeDipoleEnd* dip) {
3672
3673   // Default is no asymmetry. Only gluons are studied.
3674   dip->asymPol = 0.;
3675   dip->iAunt = 0;
3676   int iRad = dip->iRadiator;
3677   if (!doPhiPolAsym || event[iRad].id() != 21) return;
3678
3679   // Trace grandmother via possibly intermediate recoil copies.
3680   int iMother = event.iTopCopy(iRad);
3681   int iGrandM = event[iMother].mother1();
3682
3683   // If grandmother in initial state of hard scattering,
3684   // then only keep gg and qq initial states.
3685   int statusGrandM = event[iGrandM].status();
3686   bool isHardProc  = (statusGrandM == -21 || statusGrandM == -31);
3687   if (isHardProc) {
3688     if (event[iGrandM + 1].status() != statusGrandM) return;
3689     if (event[iGrandM].isGluon() && event[iGrandM + 1].isGluon());
3690     else if (event[iGrandM].isQuark() && event[iGrandM + 1].isQuark());
3691     else return;
3692   }
3693
3694   // Set aunt by history or, for hard scattering, by colour flow.
3695   if (isHardProc) dip->iAunt = dip->iRecoiler;
3696   else dip->iAunt = (event[iGrandM].daughter1() == iMother)
3697     ? event[iGrandM].daughter2() : event[iGrandM].daughter1();
3698
3699   // Coefficient from gluon production (approximate z by energy).
3700   // For hard process arbitrarily put z = 1/2.
3701   double zProd = (isHardProc) ? 0.5 : event[iRad].e()
3702     / (event[iRad].e() + event[dip->iAunt].e());
3703   if (event[iGrandM].isGluon()) dip->asymPol = pow2( (1. - zProd)
3704     / (1. - zProd * (1. - zProd) ) );
3705   else dip->asymPol = 2. * (1. - zProd) / (1. + pow2(1. - zProd) );
3706
3707   // Coefficients from gluon decay.
3708   if (dip->flavour == 21) dip->asymPol *= pow2( (1. - dip->z)
3709     / (1. - dip->z * (1. - dip->z) ) );
3710   else  dip->asymPol *= -2. * dip->z *( 1. - dip->z )
3711     / (1. - 2. * dip->z * (1. - dip->z) );
3712
3713 }
3714
3715 //--------------------------------------------------------------------------
3716
3717 // Print the list of dipoles.
3718
3719 void TimeShower::list(ostream& os) const {
3720
3721   // Header.
3722   os << "\n --------  PYTHIA TimeShower Dipole Listing  ----------------"
3723      << "--------------------------------------------- \n \n    i    rad"
3724      << "    rec       pTmax  col  chg  gam  oni   hv  isr  sys sysR typ"
3725      << "e  MErec     mix  ord  spl  ~gR \n" << fixed << setprecision(3);
3726
3727   // Loop over dipole list and print it.
3728   for (int i = 0; i < int(dipEnd.size()); ++i)
3729   os << setw(5) << i                     << setw(7) << dipEnd[i].iRadiator
3730      << setw(7) << dipEnd[i].iRecoiler   << setw(12) << dipEnd[i].pTmax
3731      << setw(5) << dipEnd[i].colType     << setw(5) << dipEnd[i].chgType
3732      << setw(5) << dipEnd[i].gamType     << setw(5) << dipEnd[i].isOctetOnium
3733      << setw(5) << dipEnd[i].isHiddenValley << setw(5) << dipEnd[i].isrType
3734      << setw(5) << dipEnd[i].system      << setw(5) << dipEnd[i].systemRec
3735      << setw(5) << dipEnd[i].MEtype      << setw(7) << dipEnd[i].iMEpartner
3736      << setw(8) << dipEnd[i].MEmix       << setw(5) << dipEnd[i].MEorder
3737      << setw(5) << dipEnd[i].MEsplit     << setw(5) << dipEnd[i].MEgluinoRec
3738      << "\n";
3739
3740   // Done.
3741   os << "\n --------  End PYTHIA TimeShower Dipole Listing  ------------"
3742      << "---------------------------------------------" << endl;
3743
3744 }
3745
3746 //==========================================================================
3747
3748 } // end namespace Pythia8
3749