PYTHIA8/pythia8170/src/SpaceShower.cxx

   1 // SpaceShower.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
   7 // SpaceShower class.
   8
   9 #include "SpaceShower.h"
  10
  11 namespace Pythia8 {
  12
  13 //==========================================================================
  14
  15 // The SpaceShower class.
  16
  17 //--------------------------------------------------------------------------
  18
  19 // Constants: could be changed here if desired, but normally should not.
  20 // These are of technical nature, as described for each.
  21
  22 // Leftover companion can give PDF > 0 at small Q2 where other PDF's = 0,
  23 // and then one can end in infinite loop of impossible kinematics.
  24 const int    SpaceShower::MAXLOOPTINYPDF = 10;
  25
  26 // Switch to alternative (but equivalent) backwards evolution for
  27 // g -> Q Qbar (Q = c or b) when below QTHRESHOLD * mQ2.
  28 const double SpaceShower::CTHRESHOLD     = 2.0;
  29 const double SpaceShower::BTHRESHOLD     = 2.0;
  30
  31 // Renew evaluation of PDF's when the pT2 step is bigger than this
  32 // (in addition to initial scale and c and b thresholds.)
  33 const double SpaceShower::EVALPDFSTEP    = 0.1;
  34
  35 // Lower limit on PDF value in order to avoid division by zero.
  36 const double SpaceShower::TINYPDF        = 1e-10;
  37
  38 // Lower limit on estimated evolution rate, below which stop.
  39 const double SpaceShower::TINYKERNELPDF  = 1e-6;
  40
  41 // Lower limit on pT2, below which branching is rejected.
  42 const double SpaceShower::TINYPT2        = 0.25e-6;
  43
  44 // No attempt to do backwards evolution of a heavy (c or b) quark
  45 // if evolution starts at a scale pT2 < HEAVYPT2EVOL * mQ2.
  46 const double SpaceShower::HEAVYPT2EVOL   = 1.1;
  47
  48 // No attempt to do backwards evolution of a heavy (c or b) quark
  49 // if evolution starts at a  x > HEAVYXEVOL * x_max, where
  50 // x_max is the largest possible x value for a g -> Q Qbar branching.
  51 const double SpaceShower::HEAVYXEVOL     = 0.9;
  52
  53 // When backwards evolution Q -> g + Q creates a heavy quark Q,
  54 // an earlier branching g -> Q + Qbar will restrict kinematics
  55 // to  M_{Q Qbar}^2 > EXTRASPACEQ * 4 m_Q^2. (Smarter to be found??)
  56 const double SpaceShower::EXTRASPACEQ    = 2.0;
  57
  58 // Never pick pT so low that alphaS is evaluated too close to Lambda_3.
  59 const double SpaceShower::LAMBDA3MARGIN  = 1.1;
  60
  61 // Do not warn for large PDF ratios at small pT2 scales.
  62 const double SpaceShower::PT2MINWARN = 1.;
  63
  64 // Cutoff for f_e^e at x < 1 - 10^{-10} to be used in z selection.
  65 // Note: the x_min quantity come from 1 - x_max.
  66 const double SpaceShower::LEPTONXMIN     = 1e-10;
  67 const double SpaceShower::LEPTONXMAX     = 1. - 1e-10;
  68
  69 // Stop l -> l gamma evolution slightly above m2l.
  70 const double SpaceShower::LEPTONPT2MIN   = 1.2;
  71
  72 // Enhancement of l -> l gamma trial rate to compensate imperfect modelling.
  73 const double SpaceShower::LEPTONFUDGE    = 10.;
  74
  75 //--------------------------------------------------------------------------
  76
  77 // Initialize alphaStrong, alphaEM and related pTmin parameters.
  78
  79 void SpaceShower::init( BeamParticle* beamAPtrIn,
  80   BeamParticle* beamBPtrIn) {
  81
  82   // Store input pointers for future use.
  83   beamAPtr        = beamAPtrIn;
  84   beamBPtr        = beamBPtrIn;
  85
  86   // Main flags to switch on and off branchings.
  87   doQCDshower     = settingsPtr->flag("SpaceShower:QCDshower");
  88   doQEDshowerByQ  = settingsPtr->flag("SpaceShower:QEDshowerByQ");
  89   doQEDshowerByL  = settingsPtr->flag("SpaceShower:QEDshowerByL");
  90
  91   // Matching in pT of hard interaction to shower evolution.
  92   pTmaxMatch      = settingsPtr->mode("SpaceShower:pTmaxMatch");
  93   pTdampMatch     = settingsPtr->mode("SpaceShower:pTdampMatch");
  94   pTmaxFudge      = settingsPtr->parm("SpaceShower:pTmaxFudge");
  95   pTmaxFudgeMPI   = settingsPtr->parm("SpaceShower:pTmaxFudgeMPI");
  96   pTdampFudge     = settingsPtr->parm("SpaceShower:pTdampFudge");
  97
  98   // Optionally force emissions to be ordered in rapidity/angle.
  99   doRapidityOrder = settingsPtr->flag("SpaceShower:rapidityOrder");
 100
 101   // Charm, bottom and lepton mass thresholds.
 102   mc              = particleDataPtr->m0(4);
 103   mb              = particleDataPtr->m0(5);
 104   m2c             = pow2(mc);
 105   m2b             = pow2(mb);
 106
 107   // Parameters of scale choices.
 108   renormMultFac     = settingsPtr->parm("SpaceShower:renormMultFac");
 109   factorMultFac     = settingsPtr->parm("SpaceShower:factorMultFac");
 110
 111   // Parameters of alphaStrong generation.
 112   alphaSvalue     = settingsPtr->parm("SpaceShower:alphaSvalue");
 113   alphaSorder     = settingsPtr->mode("SpaceShower:alphaSorder");
 114   alphaS2pi       = 0.5 * alphaSvalue / M_PI;
 115
 116   // Initialize alpha_strong generation.
 117   alphaS.init( alphaSvalue, alphaSorder);
 118
 119   // Lambda for 5, 4 and 3 flavours.
 120   Lambda5flav     = alphaS.Lambda5();
 121   Lambda4flav     = alphaS.Lambda4();
 122   Lambda3flav     = alphaS.Lambda3();
 123   Lambda5flav2    = pow2(Lambda5flav);
 124   Lambda4flav2    = pow2(Lambda4flav);
 125   Lambda3flav2    = pow2(Lambda3flav);
 126
 127   // Regularization of QCD evolution for pT -> 0. Can be taken
 128   // same as for multiparton interactions, or be set separately.
 129   useSamePTasMPI  = settingsPtr->flag("SpaceShower:samePTasMPI");
 130   if (useSamePTasMPI) {
 131     pT0Ref        = settingsPtr->parm("MultipartonInteractions:pT0Ref");
 132     ecmRef        = settingsPtr->parm("MultipartonInteractions:ecmRef");
 133     ecmPow        = settingsPtr->parm("MultipartonInteractions:ecmPow");
 134     pTmin         = settingsPtr->parm("MultipartonInteractions:pTmin");
 135   } else {
 136     pT0Ref        = settingsPtr->parm("SpaceShower:pT0Ref");
 137     ecmRef        = settingsPtr->parm("SpaceShower:ecmRef");
 138     ecmPow        = settingsPtr->parm("SpaceShower:ecmPow");
 139     pTmin         = settingsPtr->parm("SpaceShower:pTmin");
 140   }
 141
 142   // Calculate nominal invariant mass of events. Set current pT0 scale.
 143   sCM             = m2( beamAPtr->p(), beamBPtr->p());
 144   eCM             = sqrt(sCM);
 145   pT0             = pT0Ref * pow(eCM / ecmRef, ecmPow);
 146
 147   // Restrict pTmin to ensure that alpha_s(pTmin^2 + pT_0^2) does not blow up.
 148   double pTminAbs = sqrtpos(pow2(LAMBDA3MARGIN) * Lambda3flav2 / renormMultFac
 149                   - pT0*pT0);
 150   if (pTmin < pTminAbs) {
 151     pTmin         = pTminAbs;
 152     ostringstream newPTmin;
 153     newPTmin << fixed << setprecision(3) << pTmin;
 154     infoPtr->errorMsg("Warning in SpaceShower::init: pTmin too low",
 155                       ", raised to " + newPTmin.str() );
 156     infoPtr->setTooLowPTmin(true);
 157   }
 158
 159   // Parameters of alphaEM generation.
 160   alphaEMorder    = settingsPtr->mode("SpaceShower:alphaEMorder");
 161
 162   // Initialize alphaEM generation.
 163   alphaEM.init( alphaEMorder, settingsPtr);
 164
 165   // Parameters of QED evolution.
 166   pTminChgQ       = settingsPtr->parm("SpaceShower:pTminchgQ");
 167   pTminChgL       = settingsPtr->parm("SpaceShower:pTminchgL");
 168
 169   // Derived parameters of QCD evolution.
 170   pT20            = pow2(pT0);
 171   pT2min          = pow2(pTmin);
 172   pT2minChgQ      = pow2(pTminChgQ);
 173   pT2minChgL      = pow2(pTminChgL);
 174
 175   // Various other parameters.
 176   doMEcorrections = settingsPtr->flag("SpaceShower:MEcorrections");
 177   doMEafterFirst  = settingsPtr->flag("SpaceShower:MEafterFirst");
 178   doPhiPolAsym    = settingsPtr->flag("SpaceShower:phiPolAsym");
 179   doPhiIntAsym    = settingsPtr->flag("SpaceShower:phiIntAsym");
 180   strengthIntAsym = settingsPtr->parm("SpaceShower:strengthIntAsym");
 181   nQuarkIn        = settingsPtr->mode("SpaceShower:nQuarkIn");
 182
 183   // Optional dampening at small pT's when large multiplicities.
 184   enhanceScreening
 185     = settingsPtr->mode("MultipartonInteractions:enhanceScreening");
 186   if (!useSamePTasMPI) enhanceScreening = 0;
 187
 188   // Possibility to allow user veto of emission step.
 189   canVetoEmission = (userHooksPtr != 0)
 190                   ? userHooksPtr->canVetoISREmission() : false;
 191
 192 }
 193
 194 //--------------------------------------------------------------------------
 195
 196 // Find whether to limit maximum scale of emissions.
 197 // Also allow for dampening at factorization or renormalization scale.
 198
 199 bool SpaceShower::limitPTmax( Event& event, double Q2Fac, double Q2Ren) {
 200
 201   // Find whether to limit pT. Begin by user-set cases.
 202   bool dopTlimit = false;
 203   if      (pTmaxMatch == 1) dopTlimit = true;
 204   else if (pTmaxMatch == 2) dopTlimit = false;
 205
 206   // Look if any quark (u, d, s, c, b), gluon or photon in final state.
 207   else {
 208     for (int i = 5; i < event.size(); ++i)
 209     if (event[i].status() != -21) {
 210       int idAbs = event[i].idAbs();
 211       if (idAbs <= 5 || idAbs == 21 || idAbs == 22) dopTlimit = true;
 212     }
 213   }
 214
 215   // Dampening at factorization or renormalization scale.
 216   dopTdamp   = false;
 217   pT2damp    = 0.;
 218   if ( !dopTlimit && (pTdampMatch == 1 || pTdampMatch == 2) ) {
 219     dopTdamp = true;
 220     pT2damp  = pow2(pTdampFudge) * ((pTdampMatch == 1) ? Q2Fac : Q2Ren);
 221   }
 222
 223   // Done.
 224   return dopTlimit;
 225
 226 }
 227
 228 //--------------------------------------------------------------------------
 229
 230 // Prepare system for evolution; identify ME.
 231 // Routine may be called after multiparton interactions, for a new subystem.
 232
 233 void SpaceShower::prepare( int iSys, Event& event, bool limitPTmaxIn) {
 234
 235   // Find positions of incoming colliding partons.
 236   int in1 = partonSystemsPtr->getInA(iSys);
 237   int in2 = partonSystemsPtr->getInB(iSys);
 238
 239   // Rescattered partons cannot radiate.
 240   bool canRadiate1 = !(event[in1].isRescatteredIncoming());
 241   bool canRadiate2 = !(event[in2].isRescatteredIncoming());
 242
 243   // Reset dipole-ends list for first interaction. Also resonances.
 244   if (iSys == 0) dipEnd.resize(0);
 245   if (iSys == 0) idResFirst  = 0;
 246   if (iSys == 1) idResSecond = 0;
 247
 248   // Find matrix element corrections for system.
 249   int MEtype = findMEtype( iSys, event);
 250
 251   // Maximum pT scale for dipole ends.
 252   double pTmax1 = (limitPTmaxIn) ? event[in1].scale() : eCM;
 253   double pTmax2 = (limitPTmaxIn) ? event[in2].scale() : eCM;
 254   if (iSys == 0 && limitPTmaxIn) {
 255     pTmax1 *= pTmaxFudge;
 256     pTmax2 *= pTmaxFudge;
 257   } else if (iSys > 0 && limitPTmaxIn) {
 258     pTmax1 *= pTmaxFudgeMPI;
 259     pTmax2 *= pTmaxFudgeMPI;
 260   }
 261
 262   // Find dipole ends for QCD radiation.
 263   // Note: colour type can change during evolution, so book also if zero.
 264   if (doQCDshower) {
 265     int colType1 = event[in1].colType();
 266     if (canRadiate1) dipEnd.push_back( SpaceDipoleEnd( iSys,  1,
 267       in1, in2, pTmax1, colType1, 0, MEtype, canRadiate2) );
 268     int colType2 = event[in2].colType();
 269     if (canRadiate2) dipEnd.push_back( SpaceDipoleEnd( iSys,  2,
 270       in2, in1, pTmax2, colType2, 0, MEtype, canRadiate1) );
 271   }
 272
 273   // Find dipole ends for QED radiation.
 274   // Note: charge type can change during evolution, so book also if zero.
 275   if (doQEDshowerByQ || doQEDshowerByL) {
 276     int chgType1 = ( (event[in1].isQuark() && doQEDshowerByQ)
 277       || (event[in1].isLepton() && doQEDshowerByL) )
 278       ? event[in1].chargeType() : 0;
 279     // Special: photons have charge zero, but can evolve (only off Q for now)
 280     if (event[in1].id() == 22 && doQEDshowerByQ) chgType1 = 22 ;
 281     if (canRadiate1) dipEnd.push_back( SpaceDipoleEnd( iSys, -1,
 282       in1, in2, pTmax1, 0, chgType1, MEtype, canRadiate2) );
 283     int chgType2 = ( (event[in2].isQuark() && doQEDshowerByQ)
 284       || (event[in2].isLepton() && doQEDshowerByL) )
 285       ? event[in2].chargeType() : 0;
 286     // Special: photons have charge zero, but can evolve (only off Q for now)
 287     if (event[in2].id() == 22 && doQEDshowerByQ) chgType2 = 22 ;
 288     if (canRadiate2) dipEnd.push_back( SpaceDipoleEnd( iSys, -2,
 289       in2, in1, pTmax2, 0, chgType2, MEtype, canRadiate1) );
 290   }
 291
 292   // Store the z and pT2 values of the last previous splitting
 293   // when an event history has already been constructed.
 294   if (iSys == 0 && infoPtr->hasHistory()) {
 295     double zNow   = infoPtr->zNowISR();
 296     double pT2Now = infoPtr->pT2NowISR();
 297     for (int iDipEnd = 0; iDipEnd < int(dipEnd.size()); ++iDipEnd) {
 298       dipEnd[iDipEnd].zOld = zNow;
 299       dipEnd[iDipEnd].pT2Old = pT2Now;
 300       ++dipEnd[iDipEnd].nBranch;
 301     }
 302   }
 303
 304 }
 305
 306 //--------------------------------------------------------------------------
 307
 308 // Select next pT in downwards evolution of the existing dipoles.
 309
 310 double SpaceShower::pTnext( Event& event, double pTbegAll, double pTendAll,
 311   int nRadIn) {
 312
 313   // Current cm energy, in case it varies between events.
 314   sCM           = m2( beamAPtr->p(), beamBPtr->p());
 315   eCM           = sqrt(sCM);
 316   pTbegRef      = pTbegAll;
 317
 318   // Starting values: no radiating dipole found.
 319   nRad          = nRadIn;
 320   double pT2sel = pow2(pTendAll);
 321   iDipSel       = 0;
 322   iSysSel       = 0;
 323   dipEndSel     = 0;
 324
 325   // Loop over all possible dipole ends.
 326   for (int iDipEnd = 0; iDipEnd < int(dipEnd.size()); ++iDipEnd) {
 327     iDipNow        = iDipEnd;
 328     dipEndNow      = &dipEnd[iDipEnd];
 329     iSysNow        = dipEndNow->system;
 330     dipEndNow->pT2 = 0.;
 331
 332     // Check whether dipole end should be allowed to shower.
 333     double pT2begDip = pow2( min( pTbegAll, dipEndNow->pTmax ));
 334     if (pT2begDip > pT2sel
 335       && ( dipEndNow->colType != 0 || dipEndNow->chgType != 0 ) ) {
 336       double pT2endDip = 0.;
 337
 338       // Determine lower cut for evolution, for QCD or QED (q or l).
 339       if (dipEndNow->colType != 0) pT2endDip = max( pT2sel, pT2min );
 340       else if (abs(dipEndNow->chgType) != 3) pT2endDip
 341         = max( pT2sel, pT2minChgQ );
 342       else pT2endDip = max( pT2sel, pT2minChgL );
 343
 344       // Find properties of dipole and radiating dipole end.
 345       sideA        = ( abs(dipEndNow->side) == 1 );
 346       BeamParticle& beamNow = (sideA) ? *beamAPtr : *beamBPtr;
 347       BeamParticle& beamRec = (sideA) ? *beamBPtr : *beamAPtr;
 348       iNow         = beamNow[iSysNow].iPos();
 349       iRec         = beamRec[iSysNow].iPos();
 350       idDaughter   = beamNow[iSysNow].id();
 351       xDaughter    = beamNow[iSysNow].x();
 352       x1Now        = (sideA) ? xDaughter : beamRec[iSysNow].x();
 353       x2Now        = (sideA) ? beamRec[iSysNow].x() : xDaughter;
 354
 355       // Note dipole mass correction when recoiler is a rescatter.
 356       m2Rec        = (dipEndNow->normalRecoil) ? 0. : event[iRec].m2();
 357       m2Dip        = x1Now * x2Now * sCM + m2Rec;
 358
 359       // Now do evolution in pT2, for QCD or QED
 360       if (pT2begDip > pT2endDip) {
 361         if (dipEndNow->colType != 0) pT2nextQCD( pT2begDip, pT2endDip);
 362         else                         pT2nextQED( pT2begDip, pT2endDip);
 363       }
 364
 365       // Update if found larger pT than current maximum.
 366       if (dipEndNow->pT2 > pT2sel) {
 367         pT2sel    = dipEndNow->pT2;
 368         iDipSel   = iDipNow;
 369         iSysSel   = iSysNow;
 370         dipEndSel = dipEndNow;
 371       }
 372
 373     // End loop over dipole ends.
 374     }
 375   }
 376
 377   // Return nonvanishing value if found pT is bigger than already found.
 378   return (dipEndSel == 0) ? 0. : sqrt(pT2sel);
 379 }
 380
 381 //--------------------------------------------------------------------------
 382
 383 // Evolve a QCD dipole end.
 384
 385 void SpaceShower::pT2nextQCD( double pT2begDip, double pT2endDip) {
 386
 387   // Some properties and kinematical starting values.
 388   BeamParticle& beam = (sideA) ? *beamAPtr : *beamBPtr;
 389   bool   isGluon     = (idDaughter == 21);
 390   bool   isValence   = beam[iSysNow].isValence();
 391   int    MEtype      = dipEndNow->MEtype;
 392   double pT2         = pT2begDip;
 393   double xMaxAbs     = beam.xMax(iSysNow);
 394   double zMinAbs     = xDaughter / xMaxAbs;
 395   if (xMaxAbs < 0.) {
 396     infoPtr->errorMsg("Warning in SpaceShower::pT2nextQCD: "
 397     "xMaxAbs negative");
 398     return;
 399   }
 400
 401   // Starting values for handling of massive quarks (c/b), if any.
 402   double idMassive   = 0;
 403   if ( abs(idDaughter) == 4 ) idMassive = 4;
 404   if ( abs(idDaughter) == 5 ) idMassive = 5;
 405   bool   isMassive   = (idMassive > 0);
 406   double m2Massive   = 0.;
 407   double mRatio      = 0.;
 408   double zMaxMassive = 1.;
 409   double m2Threshold = pT2;
 410
 411   // Evolution below scale of massive quark or at large x is impossible.
 412   if (isMassive) {
 413     m2Massive = (idMassive == 4) ? m2c : m2b;
 414     if (pT2 < HEAVYPT2EVOL * m2Massive) return;
 415     mRatio = sqrt( m2Massive / m2Dip );
 416     zMaxMassive = (1. -  mRatio) / ( 1. +  mRatio * (1. -  mRatio) );
 417     if (xDaughter > HEAVYXEVOL * zMaxMassive * xMaxAbs) return;
 418
 419     // Find threshold scale below which only g -> Q + Qbar will be allowed.
 420     m2Threshold = (idMassive == 4) ? min( pT2, CTHRESHOLD * m2c)
 421       : min( pT2, BTHRESHOLD * m2b);
 422   }
 423
 424   // Variables used inside evolution loop. (Mainly dummy starting values.)
 425   int    nFlavour       = 3;
 426   double b0             = 4.5;
 427   double Lambda2        = Lambda3flav2;
 428   double pT2minNow      = pT2endDip;
 429   int    idMother       = 0;
 430   int    idSister       = 0;
 431   double z              = 0.;
 432   double zMaxAbs        = 0.;
 433   double zRootMax       = 0.;
 434   double zRootMin       = 0.;
 435   double g2gInt         = 0.;
 436   double q2gInt         = 0.;
 437   double q2qInt         = 0.;
 438   double g2qInt         = 0.;
 439   double g2Qenhance     = 0.;
 440   double xPDFdaughter   = 0.;
 441   double xPDFmother[21] = {0.};
 442   double xPDFgMother    = 0.;
 443   double xPDFmotherSum  = 0.;
 444   double kernelPDF      = 0.;
 445   double xMother        = 0.;
 446   double wt             = 0.;
 447   double Q2             = 0.;
 448   double mSister        = 0.;
 449   double m2Sister       = 0.;
 450   double pT2corr        = 0.;
 451   double pT2PDF         = pT2;
 452   bool   needNewPDF     = true;
 453
 454   // Begin evolution loop towards smaller pT values.
 455   int    loopTinyPDFdau = 0;
 456   bool   hasTinyPDFdau  = false;
 457   do {
 458     wt = 0.;
 459
 460     // Bad sign if repeated looping with small daughter PDF, so fail.
 461     // (Example: if all PDF's = 0 below Q_0, except for c/b companion.)
 462     if (hasTinyPDFdau) ++loopTinyPDFdau;
 463     if (loopTinyPDFdau > MAXLOOPTINYPDF) {
 464       infoPtr->errorMsg("Warning in SpaceShower::pT2nextQCD: "
 465       "small daughter PDF");
 466       return;
 467     }
 468
 469     // Initialize integrals of splitting kernels and evaluate parton
 470     // densities at the beginning. Reinitialize after long evolution
 471     // in pT2 or when crossing c and b flavour thresholds.
 472     if (needNewPDF || pT2 < EVALPDFSTEP * pT2PDF) {
 473       pT2PDF        = pT2;
 474       hasTinyPDFdau = false;
 475
 476       // Determine overestimated z range; switch at c and b masses.
 477       if (pT2 > m2b) {
 478         nFlavour  = 5;
 479         pT2minNow = max( m2b, pT2endDip);
 480         b0        = 23./6.;
 481         Lambda2   = Lambda5flav2;
 482       } else if (pT2 > m2c) {
 483         nFlavour  = 4;
 484         pT2minNow = max( m2c, pT2endDip);
 485         b0        = 25./6.;
 486         Lambda2   = Lambda4flav2;
 487       } else {
 488         nFlavour  = 3;
 489         pT2minNow = pT2endDip;
 490         b0        = 27./6.;
 491         Lambda2   = Lambda3flav2;
 492       }
 493       // A change of renormalization scale expressed by a change of Lambda.
 494       Lambda2    /= renormMultFac;
 495       zMaxAbs     = 1. - 0.5 * (pT2minNow / m2Dip) *
 496         ( sqrt( 1. + 4. * m2Dip / pT2minNow ) - 1. );
 497       if (isMassive) zMaxAbs = min( zMaxAbs, zMaxMassive);
 498
 499       // Go to another z range with lower mass scale if current is closed.
 500       if (zMinAbs > zMaxAbs) {
 501         if (nFlavour == 3 || (idMassive == 4 && nFlavour == 4)
 502           || idMassive == 5) return;
 503         pT2 = (nFlavour == 4) ? m2c : m2b;
 504         continue;
 505       }
 506
 507       // Parton density of daughter at current scale.
 508       xPDFdaughter = beam.xfISR(iSysNow, idDaughter, xDaughter,
 509         factorMultFac * pT2);
 510       if (xPDFdaughter < TINYPDF) {
 511         xPDFdaughter  = TINYPDF;
 512         hasTinyPDFdau = true;
 513       }
 514
 515       // Integrals of splitting kernels for gluons: g -> g, q -> g.
 516       if (isGluon) {
 517         g2gInt = 6. * log(zMaxAbs * (1.-zMinAbs)
 518           / (zMinAbs * (1.-zMaxAbs)));
 519         if (doMEcorrections) g2gInt *= calcMEmax(MEtype, 21, 21);
 520         q2gInt = (16./3.) * (1./sqrt(zMinAbs) - 1./sqrt(zMaxAbs));
 521         if (doMEcorrections) q2gInt *= calcMEmax(MEtype, 1, 21);
 522
 523         // Parton density of potential quark mothers to a g.
 524         xPDFmotherSum = 0.;
 525         for (int i = -nQuarkIn; i <= nQuarkIn; ++i) {
 526           if (i == 0) {
 527             xPDFmother[10] = 0.;
 528           } else {
 529             xPDFmother[i+10] = beam.xfISR(iSysNow, i, xDaughter,
 530               factorMultFac * pT2);
 531             xPDFmotherSum += xPDFmother[i+10];
 532           }
 533         }
 534
 535         // Total QCD evolution coefficient for a gluon.
 536         kernelPDF = g2gInt + q2gInt * xPDFmotherSum / xPDFdaughter;
 537
 538       // For valence quark only need consider q -> q g branchings.
 539       // Introduce an extra factor sqrt(z) to smooth bumps.
 540       } else if (isValence) {
 541         zRootMin = (1. + sqrt(zMinAbs)) / (1. - sqrt(zMinAbs));
 542         zRootMax = (1. + sqrt(zMaxAbs)) / (1. - sqrt(zMaxAbs));
 543         q2qInt = (8./3.) * log( zRootMax / zRootMin );
 544         if (doMEcorrections) q2qInt *= calcMEmax(MEtype, 1, 1);
 545         kernelPDF = q2qInt;
 546
 547       // Integrals of splitting kernels for quarks: q -> q, g -> q.
 548       } else {
 549         q2qInt = (8./3.) * log( (1. - zMinAbs) / (1. - zMaxAbs) );
 550         if (doMEcorrections) q2qInt *= calcMEmax(MEtype, 1, 1);
 551         g2qInt = 0.5 * (zMaxAbs - zMinAbs);
 552         if (doMEcorrections) g2qInt *= calcMEmax(MEtype, 21, 1);
 553
 554         // Increase estimated upper weight for g -> Q + Qbar.
 555         if (isMassive) {
 556           if (alphaSorder == 0) g2Qenhance = log(pT2/m2Massive)
 557             / log(m2Threshold/m2Massive);
 558           else {
 559             double m2log = log( m2Massive / Lambda2);
 560             g2Qenhance = log( log(pT2/Lambda2) / m2log )
 561               / log( log(m2Threshold/Lambda2) / m2log );
 562           }
 563           g2qInt *= g2Qenhance;
 564         }
 565
 566         // Parton density of a potential gluon mother to a q.
 567         xPDFgMother = beam.xfISR(iSysNow, 21, xDaughter, factorMultFac * pT2);
 568
 569         // Total QCD evolution coefficient for a quark.
 570         kernelPDF = q2qInt + g2qInt * xPDFgMother / xPDFdaughter;
 571       }
 572
 573       // End evaluation of splitting kernels and parton densities.
 574       needNewPDF = false;
 575     }
 576     if (kernelPDF < TINYKERNELPDF) return;
 577
 578     // Pick pT2 (in overestimated z range), for one of three different cases.
 579     // Assume form alphas(pT0^2 + pT^2) * dpT^2/(pT0^2 + pT^2).
 580     double Q2alphaS;
 581
 582     // Fixed alpha_strong.
 583     if (alphaSorder == 0) {
 584       pT2 = (pT2 + pT20) * pow( rndmPtr->flat(),
 585         1. / (alphaS2pi * kernelPDF)) - pT20;
 586
 587     // First-order alpha_strong.
 588     } else if (alphaSorder == 1) {
 589       pT2 = Lambda2 * pow( (pT2 + pT20) / Lambda2,
 590         pow(rndmPtr->flat(), b0 / kernelPDF) ) - pT20;
 591
 592     // For second order reject by second term in alpha_strong expression.
 593     } else {
 594       do {
 595         pT2 = Lambda2 * pow( (pT2 + pT20) / Lambda2,
 596           pow(rndmPtr->flat(), b0 / kernelPDF) ) - pT20;
 597         Q2alphaS = renormMultFac * max( pT2 + pT20,
 598           pow2(LAMBDA3MARGIN) * Lambda3flav2);
 599       } while (alphaS.alphaS2OrdCorr(Q2alphaS) < rndmPtr->flat()
 600         && pT2 > pT2minNow);
 601     }
 602
 603     // Check for pT2 values that prompt special action.
 604
 605     // If fallen into b threshold region, force g -> b + bbar.
 606     if (idMassive == 5 && pT2 < m2Threshold) {
 607       pT2nearQCDthreshold( beam, m2Massive, m2Threshold, xMaxAbs,
 608         zMinAbs, zMaxMassive );
 609       return;
 610
 611     // If crossed b threshold, continue evolution from this threshold.
 612     } else if (nFlavour == 5 && pT2 < m2b) {
 613       needNewPDF = true;
 614       pT2 = m2b;
 615       continue;
 616
 617     // If fallen into c threshold region, force g -> c + cbar.
 618     } else if (idMassive == 4 && pT2 < m2Threshold) {
 619       pT2nearQCDthreshold( beam, m2Massive, m2Threshold, xMaxAbs,
 620         zMinAbs, zMaxMassive );
 621       return;
 622
 623     // If crossed c threshold, continue evolution from this threshold.
 624     } else if (nFlavour == 4 && pT2 < m2c) {
 625       needNewPDF = true;
 626       pT2 = m2c;
 627       continue;
 628
 629     // Abort evolution if below cutoff scale, or below another branching.
 630     } else if (pT2 < pT2endDip) return;
 631
 632     // Select z value of branching to g, and corrective weight.
 633     if (isGluon) {
 634       // g -> g (+ g).
 635       if (rndmPtr->flat() * kernelPDF < g2gInt) {
 636         idMother = 21;
 637         idSister = 21;
 638         z = 1. / ( 1. + ((1. - zMinAbs) / zMinAbs) * pow( (zMinAbs *
 639           (1. - zMaxAbs)) / (zMaxAbs * (1. - zMinAbs)), rndmPtr->flat() ) );
 640         wt = pow2( 1. - z * (1. - z));
 641       } else {
 642       // q -> g (+ q): also select flavour.
 643         double temp = xPDFmotherSum * rndmPtr->flat();
 644         idMother = -nQuarkIn - 1;
 645         do { temp -= xPDFmother[(++idMother) + 10]; }
 646         while (temp > 0. && idMother < nQuarkIn);
 647         idSister = idMother;
 648         z = (zMinAbs * zMaxAbs) / pow2( sqrt(zMinAbs) + rndmPtr->flat()
 649           * ( sqrt(zMaxAbs)- sqrt(zMinAbs) ));
 650         wt = 0.5 * (1. + pow2(1. - z)) * sqrt(z)
 651           * xPDFdaughter / xPDFmother[idMother + 10];
 652       }
 653
 654     // Select z value of branching to q, and corrective weight.
 655     // Include massive kernel corrections for c and b quarks.
 656     } else {
 657       // q -> q (+ g).
 658       if (isValence || rndmPtr->flat() * kernelPDF < q2qInt) {
 659         idMother = idDaughter;
 660         idSister = 21;
 661         // Valence more peaked at large z.
 662         if (isValence) {
 663           double zTmp = zRootMin * pow(zRootMax / zRootMin, rndmPtr->flat() );
 664           z = pow2( (1. - zTmp) / (1. + zTmp) );
 665         } else {
 666           z = 1. - (1. - zMinAbs) * pow( (1. - zMaxAbs) / (1. - zMinAbs),
 667             rndmPtr->flat() );
 668         }
 669         if (!isMassive) {
 670           wt = 0.5 * (1. + pow2(z));
 671         } else {
 672         //?? Bug?? should be 2 more for massive part??
 673         //  wt = 0.5 * (1. + pow2(z) - z * pow2(1.-z) * m2Massive / pT2);
 674           wt = 0.5 * (1. + pow2(z)) - z * pow2(1.-z) * m2Massive / pT2;
 675         }
 676         if (isValence) wt *= sqrt(z);
 677       // g -> q (+ qbar).
 678       } else {
 679         idMother = 21;
 680         idSister = - idDaughter;
 681         z = zMinAbs + rndmPtr->flat() * (zMaxAbs - zMinAbs);
 682         if (!isMassive) {
 683           wt = (pow2(z) + pow2(1.-z)) * xPDFdaughter / xPDFgMother ;
 684         } else {
 685           wt = (pow2(z) + pow2(1.-z) + 2. * z * (1.-z) * m2Massive / pT2)
 686             * xPDFdaughter / (xPDFgMother * g2Qenhance) ;
 687         }
 688       }
 689     }
 690
 691     // Derive Q2 and x of mother from pT2 and z.
 692     Q2      = pT2 / (1. - z);
 693     xMother = xDaughter / z;
 694     // Correction to x for massive recoiler from rescattering.
 695     if (!dipEndNow->normalRecoil) {
 696       if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
 697       else       xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
 698     }
 699     if(xMother > xMaxAbs) { wt = 0.; continue; }
 700
 701     // Forbidden emission if outside allowed z range for given pT2.
 702     mSister = particleDataPtr->m0(idSister);
 703     m2Sister = pow2(mSister);
 704     pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
 705     if(pT2corr < TINYPT2) { wt = 0.; continue; }
 706
 707     // Optionally veto emissions not ordered in rapidity (= angle).
 708     if ( doRapidityOrder && dipEndNow->nBranch > 0
 709       && pT2 > pow2( (1. - z) / (z * (1. - dipEndNow->zOld)) )
 710       * dipEndNow->pT2Old ) { wt = 0.; continue; }
 711
 712     // If creating heavy quark by Q -> g + Q then next need g -> Q + Qbar.
 713     // So minimum total mass2 is 4 * m2Sister, but use more to be safe.
 714     if ( isGluon && ( abs(idMother) == 4 || abs(idMother) == 5 )) {
 715       double m2QQsister =  EXTRASPACEQ * 4. * m2Sister;
 716       double pT2QQcorr = Q2 - z * (m2Dip + Q2) * (Q2 + m2QQsister) / m2Dip;
 717       if(pT2QQcorr < TINYPT2) { wt = 0.; continue; }
 718     }
 719
 720     // Evaluation of ME correction.
 721     if (doMEcorrections) wt *= calcMEcorr(MEtype, idMother, idDaughter,
 722       m2Dip, z, Q2) / calcMEmax(MEtype, idMother, idDaughter);
 723
 724     // Optional dampening of large pT values in first radiation.
 725     if (dopTdamp && iSysNow == 0 && MEtype == 0 && nRad == 0)
 726       wt *= pT2damp / (pT2 + pT2damp);
 727
 728     // Idea suggested by Gosta Gustafson: increased screening in events
 729     // with large activity can be simulated by pT0_eff = sqrt(n) * pT0.
 730     if (enhanceScreening == 2) {
 731       int nSysNow     = infoPtr->nMPI() + infoPtr->nISR() + 1;
 732       double WTscreen = pow2( (pT2 + pT20) / (pT2 + nSysNow * pT20) );
 733       wt             *= WTscreen;
 734     }
 735
 736     // Evaluation of new daughter and mother PDF's.
 737     double xPDFdaughterNew = max ( TINYPDF,
 738       beam.xfISR(iSysNow, idDaughter, xDaughter, factorMultFac * pT2) );
 739     double xPDFmotherNew =
 740       beam.xfISR(iSysNow, idMother, xMother, factorMultFac * pT2);
 741     wt *= xPDFmotherNew / xPDFdaughterNew;
 742
 743     // Check that valence step does not cause problem.
 744     if (wt > 1. && pT2 > PT2MINWARN) infoPtr->errorMsg("Warning in "
 745       "SpaceShower::pT2nextQCD: weight above unity");
 746
 747   // Iterate until acceptable pT (or have fallen below pTmin).
 748   } while (wt < rndmPtr->flat()) ;
 749
 750   // Save values for (so far) acceptable branching.
 751   dipEndNow->store( idDaughter,idMother, idSister, x1Now, x2Now, m2Dip,
 752     pT2, z, xMother, Q2, mSister, m2Sister, pT2corr);
 753
 754 }
 755
 756 //--------------------------------------------------------------------------
 757
 758 // Evolve a QCD dipole end near threshold, with g -> Q + Qbar enforced.
 759 // Note: No explicit Sudakov factor formalism here. Instead use that
 760 // df_Q(x, pT2) = (alpha_s/2pi) * (dT2/pT2) * ((gluon) * (splitting)).
 761 // This implies that effects of Q -> Q + g are neglected in this range.
 762
 763 void SpaceShower::pT2nearQCDthreshold( BeamParticle& beam,
 764   double m2Massive, double m2Threshold, double xMaxAbs,
 765   double zMinAbs, double zMaxMassive) {
 766
 767   // Initial values, to be used in kinematics and weighting.
 768   double Lambda2       = (abs(idDaughter) == 4) ? Lambda4flav2 : Lambda5flav2;
 769   Lambda2             /= renormMultFac;
 770   double logM2Lambda2  = log( m2Massive / Lambda2 );
 771   double xPDFmotherOld = beam.xfISR(iSysNow, 21, xDaughter,
 772     factorMultFac * m2Threshold);
 773
 774   // Variables used inside evolution loop. (Mainly dummy start values.)
 775   int    loop    = 0;
 776   double wt      = 0.;
 777   double pT2     = 0.;
 778   double z       = 0.;
 779   double Q2      = 0.;
 780   double pT2corr = 0.;
 781   double xMother = 0.;
 782
 783   // Begin loop over tries to find acceptable g -> Q + Qbar branching.
 784   do {
 785     wt = 0.;
 786
 787     // Check that not caught in infinite loop with impossible kinematics.
 788     if (++loop > 100) {
 789       infoPtr->errorMsg("Error in SpaceShower::pT2nearQCDthreshold: "
 790         "stuck in loop");
 791       return;
 792     }
 793
 794     // Pick dpT2/pT2 in range [m2Massive,thresholdRatio * m2Massive].
 795     pT2 = m2Massive * pow( m2Threshold / m2Massive, rndmPtr->flat() );
 796
 797     // Pick z flat in allowed range.
 798     z = zMinAbs + rndmPtr->flat() * (zMaxMassive - zMinAbs);
 799
 800     // Check that kinematically possible choice.
 801     Q2 = pT2 / (1.-z) - m2Massive;
 802     pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Massive) / m2Dip;
 803     if(pT2corr < TINYPT2) continue;
 804
 805     // Correction factor for running alpha_s.  ??
 806     wt = logM2Lambda2 / log( pT2 / Lambda2 );
 807
 808     // Correction factor for splitting kernel.
 809     wt *= pow2(z) + pow2(1.-z) + 2. * z * (1.-z) * m2Massive / pT2;
 810
 811     // x, including correction for massive recoiler from rescattering.
 812     xMother = xDaughter / z;
 813     if (!dipEndNow->normalRecoil) {
 814       if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
 815       else       xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
 816     }
 817     if (xMother > xMaxAbs) { wt = 0.; continue; }
 818
 819     // Correction factor for gluon density.
 820     double xPDFmotherNew = beam.xfISR(iSysNow, 21, xMother,
 821       factorMultFac * pT2);
 822     wt *= xPDFmotherNew / xPDFmotherOld;
 823
 824   // Iterate until acceptable pT and z.
 825   } while (wt < rndmPtr->flat()) ;
 826
 827   // Save values for (so far) acceptable branching.
 828   double mSister = (abs(idDaughter) == 4) ? mc : mb;
 829   dipEndNow->store( idDaughter, 21, -idDaughter, x1Now, x2Now, m2Dip,
 830     pT2, z, xMother, Q2, mSister, pow2(mSister), pT2corr);
 831
 832 }
 833
 834 //--------------------------------------------------------------------------
 835
 836 // Evolve a QED dipole end.
 837
 838 void SpaceShower::pT2nextQED( double pT2begDip, double pT2endDip) {
 839
 840   // Type of dipole and starting values.
 841   BeamParticle& beam  = (sideA) ? *beamAPtr : *beamBPtr;
 842   bool   isLeptonBeam = beam.isLepton();
 843   int    MEtype       = dipEndNow->MEtype;
 844   bool   isPhoton     = (idDaughter == 22);
 845   double pT2          = pT2begDip;
 846   double m2Lepton     = (isLeptonBeam) ? pow2(beam.m()) : 0.;
 847   if (isLeptonBeam && pT2begDip < m2Lepton) return;
 848
 849   // Currently no f -> gamma branching implemented for lepton beams
 850   if (isPhoton && isLeptonBeam) return;
 851
 852   // alpha_em at maximum scale provides upper estimate.
 853   double alphaEMmax  = alphaEM.alphaEM(renormMultFac * pT2begDip);
 854   double alphaEM2pi  = alphaEMmax / (2. * M_PI);
 855
 856   // Maximum x of mother implies minimum z = xDaughter / xMother.
 857   double xMaxAbs  = (isLeptonBeam) ? LEPTONXMAX : beam.xMax(iSysNow);
 858   double zMinAbs  = xDaughter / xMaxAbs;
 859   if (xMaxAbs < 0.) {
 860     infoPtr->errorMsg("Warning in SpaceShower::pT2nextQED: "
 861     "xMaxAbs negative");
 862     return;
 863   }
 864
 865   // Maximum z from minimum pT and, for lepton, from minimum x_gamma.
 866   double zMaxAbs = 1. - 0.5 * (pT2endDip / m2Dip) *
 867     ( sqrt( 1. + 4. * m2Dip / pT2endDip ) - 1. );
 868   if (isLeptonBeam) {
 869     double zMaxLepton = xDaughter / (xDaughter + LEPTONXMIN);
 870     if (zMaxLepton < zMaxAbs) zMaxAbs = zMaxLepton;
 871   }
 872   if (zMaxAbs < zMinAbs) return;
 873
 874   // Variables used inside evolution loop. (Mainly dummy start values.)
 875   int    idMother = 0;
 876   int    idSister =22;
 877   double z        = 0.;
 878   double xMother  = 0.;
 879   double wt       = 0.;
 880   double Q2       = 0.;
 881   double mSister  = 0.;
 882   double m2Sister = 0.;
 883   double pT2corr  = 0.;
 884
 885   // QED evolution of fermions
 886   if (!isPhoton) {
 887
 888     // Integrals of splitting kernels for fermions: f -> f. Use 1 + z^2 < 2.
 889     // Ansatz f(z) = 2 / (1 - z), with + 2 / (z - xDaughter) for lepton.
 890     double f2fInt  = 0.;
 891     double f2fIntA = 2. * log( (1. - zMinAbs) / (1. - zMaxAbs) );
 892     double f2fIntB = 0.;
 893     if (isLeptonBeam) {
 894       f2fIntB      = 2. * log( (zMaxAbs - xDaughter) / (zMinAbs - xDaughter) );
 895       f2fInt       = f2fIntA + f2fIntB;
 896     } else f2fInt  = pow2(dipEndNow->chgType / 3.) * f2fIntA;
 897
 898     // Upper estimate for evolution equation, including fudge factor.
 899     if (doMEcorrections) f2fInt *= calcMEmax(MEtype, 1, 1);
 900     double kernelPDF = alphaEM2pi * f2fInt;
 901     double fudge = (isLeptonBeam) ? LEPTONFUDGE * log(m2Dip/m2Lepton) : 1.;
 902     kernelPDF *= fudge;
 903     if (kernelPDF < TINYKERNELPDF) return;
 904
 905     // Begin evolution loop towards smaller pT values.
 906     do {
 907
 908       // Pick pT2 (in overestimated z range).
 909       // For l -> l gamma include extrafactor 1 / ln(pT2 / m2l) in evolution.
 910       double shift = pow(rndmPtr->flat(), 1. / kernelPDF);
 911       if (isLeptonBeam) pT2 = m2Lepton * pow( pT2 / m2Lepton, shift);
 912       else              pT2 = pT2 * shift;
 913
 914       // Abort evolution if below cutoff scale, or below another branching.
 915       if (pT2 < pT2endDip) return;
 916       if (isLeptonBeam && pT2 < LEPTONPT2MIN * m2Lepton) return;
 917
 918       // Select z value of branching f -> f + gamma, and corrective weight.
 919       idMother = idDaughter;
 920       wt = 0.5 * (1. + pow2(z));
 921       if (isLeptonBeam) {
 922         if (f2fIntA > rndmPtr->flat() * (f2fIntA + f2fIntB)) {
 923           z = 1. - (1. - zMinAbs)
 924             * pow( (1. - zMaxAbs) / (1. - zMinAbs), rndmPtr->flat() );
 925         } else {
 926           z = xDaughter + (zMinAbs - xDaughter)
 927             * pow( (zMaxAbs - xDaughter) / (zMinAbs - xDaughter),
 928                   rndmPtr->flat() );
 929         }
 930         wt *= (z - xDaughter) / (1. - xDaughter);
 931       } else {
 932         z = 1. - (1. - zMinAbs)
 933           * pow( (1. - zMaxAbs) / (1. - zMinAbs), rndmPtr->flat() );
 934       }
 935
 936       // Derive Q2 and x of mother from pT2 and z.
 937       Q2      = pT2 / (1. - z);
 938       xMother = xDaughter / z;
 939       // Correction to x for massive recoiler from rescattering.
 940       if (!dipEndNow->normalRecoil) {
 941         if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
 942         else       xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
 943       }
 944       if(xMother > xMaxAbs) { wt = 0.; continue; }
 945
 946       // Forbidden emission if outside allowed z range for given pT2.
 947       mSister  = 0.;
 948       m2Sister = 0.;
 949       pT2corr  = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
 950       if(pT2corr < TINYPT2) { wt = 0.; continue; }
 951
 952       // Correct by ln(pT2 / m2l) and fudge factor.
 953       if (isLeptonBeam) wt *= log(pT2 / m2Lepton) / fudge;
 954
 955       // Evaluation of ME correction.
 956       if (doMEcorrections) wt *= calcMEcorr(MEtype, idMother, idDaughter,
 957          m2Dip, z, Q2) / calcMEmax(MEtype, idMother, idDaughter);
 958
 959       // Extra QED correction for f fbar -> W+- gamma. Debug??
 960       if (doMEcorrections && MEtype == 1 && idDaughter == idMother
 961         && ( (iSysNow == 0 && idResFirst  == 24)
 962           || (iSysNow == 1 && idResSecond == 24) ) ) {
 963         double tHe  = -Q2;
 964         double uHe  = Q2 - m2Dip * (1. - z) / z;
 965         double chg1 = abs(dipEndNow->chgType / 3.);
 966         double chg2 = 1. - chg1;
 967         wt *= pow2(chg1 * uHe - chg2 * tHe)
 968           / ( (tHe + uHe) * (pow2(chg1) * uHe + pow2(chg2) * tHe) );
 969       }
 970
 971       // Optional dampening of large pT values in first radiation.
 972       if (dopTdamp && iSysNow == 0 && MEtype == 0 && nRad == 0)
 973         wt *= pT2damp / (pT2 + pT2damp);
 974
 975       // Correct to current value of alpha_EM.
 976       double alphaEMnow = alphaEM.alphaEM(renormMultFac * pT2);
 977       wt *= (alphaEMnow / alphaEMmax);
 978
 979       // Evaluation of new daughter and mother PDF's.
 980       double xPDFdaughterNew = max ( TINYPDF,
 981          beam.xfISR(iSysNow, idDaughter, xDaughter, factorMultFac * pT2) );
 982       double xPDFmotherNew   =
 983          beam.xfISR(iSysNow, idMother, xMother, factorMultFac * pT2);
 984       wt *= xPDFmotherNew / xPDFdaughterNew;
 985
 986     // Iterate until acceptable pT (or have fallen below pTmin).
 987     } while (wt < rndmPtr->flat()) ;
 988   }
 989
 990   // QED evolution of photons (so far only for hadron beams).
 991   else {
 992
 993     // Initial values
 994     int    nFlavour       = 3;
 995     double kernelPDF      = 0.0;
 996     double xPDFdaughter   = 0.;
 997     double xPDFmother[21] = {0.};
 998     double xPDFmotherSum  = 0.0;
 999     double pT2PDF         = pT2;
1000     double pT2minNow      = pT2endDip;
1001     bool   needNewPDF     = true;
1002
1003     // Begin evolution loop towards smaller pT values.
1004     int    loopTinyPDFdau = 0;
1005     bool   hasTinyPDFdau  = false;
1006     do {
1007       wt = 0.;
1008
1009       // Bad sign if repeated looping with small daughter PDF, so fail.
1010       if (hasTinyPDFdau) ++loopTinyPDFdau;
1011       if (loopTinyPDFdau > MAXLOOPTINYPDF) {
1012         infoPtr->errorMsg("Warning in SpaceShower::pT2nextQED: "
1013                           "small daughter PDF");
1014         return;
1015       }
1016
1017       // Initialize integrals of splitting kernels and evaluate parton
1018       // densities at the beginning. Reinitialize after long evolution
1019       // in pT2 or when crossing c and b flavour thresholds.
1020       if (needNewPDF || pT2 < EVALPDFSTEP * pT2PDF) {
1021
1022         pT2PDF        = pT2;
1023         hasTinyPDFdau = false;
1024
1025         // Determine overestimated z range; switch at c and b masses.
1026         if (pT2 > m2b && nQuarkIn >= 5) {
1027           nFlavour  = 5;
1028           pT2minNow = max( m2b, pT2endDip);
1029         } else if (pT2 > m2c && nQuarkIn >= 4) {
1030           nFlavour  = 4;
1031           pT2minNow = max( m2c, pT2endDip);
1032         } else {
1033           nFlavour  = 3;
1034           pT2minNow = pT2endDip;
1035         }
1036
1037         // Compute upper z limit
1038         zMaxAbs = 1. - 0.5 * (pT2minNow / m2Dip) *
1039           ( sqrt( 1. + 4. * m2Dip / pT2minNow ) - 1. );
1040
1041         // Parton density of daughter at current scale.
1042         xPDFdaughter = beam.xfISR(iSysNow, idDaughter, xDaughter,
1043           factorMultFac * pT2);
1044         if (xPDFdaughter < TINYPDF) {
1045           xPDFdaughter  = TINYPDF;
1046           hasTinyPDFdau = true;
1047         }
1048
1049         // Integral over f -> gamma f splitting kernel.
1050         // Normalized so: 4/3 aS/2pi P(z) -> eq^2 * aEM/2pi P(z).
1051         // (Charge-weighting happens below.)
1052         double q2gInt = 4. * (1./sqrt(zMinAbs) - 1./sqrt(zMaxAbs));
1053
1054         // Charge-weighted Parton density of potential quark mothers.
1055         xPDFmotherSum = 0.;
1056         for (int i = -nFlavour; i <= nFlavour; ++i) {
1057           if (i == 0) {
1058             xPDFmother[10] = 0.;
1059           } else {
1060             xPDFmother[i+10] = pow2((abs(i+1) % 2 + 1)/3.0)
1061               * beam.xfISR(iSysNow, i, xDaughter, factorMultFac * pT2);
1062             xPDFmotherSum += xPDFmother[i+10];
1063           }
1064         }
1065
1066         // Total QED evolution coefficient for a photon.
1067         kernelPDF = q2gInt * xPDFmotherSum / xPDFdaughter;
1068
1069         // End evaluation of splitting kernels and parton densities.
1070         needNewPDF = false;
1071       }
1072       if (kernelPDF < TINYKERNELPDF) return;
1073
1074       // Select pT2 for next trial branching
1075       pT2 *= pow( rndmPtr->flat(), 1. / (alphaEM2pi * kernelPDF));
1076
1077       // If crossed b threshold, continue evolution from this threshold.
1078       if (nFlavour == 5 && pT2 < m2b) {
1079         needNewPDF = true;
1080         pT2 = m2b;
1081         continue;
1082       }
1083
1084       // If crossed c threshold, continue evolution from this threshold.
1085       else if (nFlavour == 4 && pT2 < m2c) {
1086         needNewPDF = true;
1087         pT2 = m2c;
1088         continue;
1089       }
1090
1091       // Abort evolution if below cutoff scale, or below another branching.
1092       else if (pT2 < pT2endDip) return;
1093
1094       // Select flavour for trial branching
1095       double temp = xPDFmotherSum * rndmPtr->flat();
1096       idMother = -nQuarkIn - 1;
1097       do {
1098         temp -= xPDFmother[(++idMother) + 10];
1099       } while (temp > 0. && idMother < nQuarkIn);
1100
1101       // Sister is same as mother, but can have m2 > 0
1102       idSister = idMother;
1103       mSister = particleDataPtr->m0(idSister);
1104       m2Sister = pow2(mSister);
1105
1106       // Select z for trial branching
1107       z = (zMinAbs * zMaxAbs) / pow2( sqrt(zMinAbs) + rndmPtr->flat()
1108                                       * ( sqrt(zMaxAbs)- sqrt(zMinAbs) ));
1109
1110       // Trial weight: splitting kernel
1111       wt = 0.5 * (1. + pow2(1. - z)) * sqrt(z);
1112
1113       // Trial weight: running alpha_EM
1114       double alphaEMnow = alphaEM.alphaEM(renormMultFac * pT2);
1115       wt *= (alphaEMnow / alphaEMmax);
1116
1117       // Derive Q2 and x of mother from pT2 and z
1118       Q2      = pT2 / (1. - z);
1119       xMother = xDaughter / z;
1120       // Correction to x for massive recoiler from rescattering.
1121       if (!dipEndNow->normalRecoil) {
1122         if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
1123         else       xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
1124       }
1125
1126       // Compute pT2corr
1127       pT2corr  = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
1128       if(pT2corr < TINYPT2) { wt = 0.; continue; }
1129
1130       // If creating heavy quark by Q -> gamma + Q then next g -> Q + Qbar.
1131       // So minimum total mass2 is 4 * m2Sister, but use more to be safe.
1132       if ( abs(idMother) == 4 || abs(idMother) == 5 ) {
1133         double m2QQsister =  EXTRASPACEQ * 4. * m2Sister;
1134         double pT2QQcorr = Q2 - z * (m2Dip + Q2) * (Q2 + m2QQsister) / m2Dip;
1135         if(pT2QQcorr < TINYPT2) { wt = 0.; continue; }
1136       }
1137
1138       // Optional dampening of large pT values in first radiation.
1139       if (dopTdamp && iSysNow == 0 && MEtype == 0 && nRad == 0)
1140         wt *= pT2damp / (pT2 + pT2damp);
1141
1142       // Evaluation of new daughter PDF
1143       double xPDFdaughterNew = beam.xfISR(iSysNow, idDaughter, xDaughter,
1144         factorMultFac * pT2);
1145       if (xPDFdaughterNew < TINYPDF) {
1146         xPDFdaughterNew = TINYPDF;
1147       }
1148
1149       // Evaluation of new charge-weighted mother PDF
1150       double xPDFmotherNew = pow2( (abs(idMother+1) % 2 + 1)/3.0 )
1151         * beam.xfISR(iSysNow, idMother, xMother, factorMultFac * pT2);
1152
1153       // Trial weight: divide out old pdf ratio
1154       wt *= xPDFdaughter / xPDFmother[idMother + 10];
1155
1156       // Trial weight: new pdf ratio
1157       wt *= xPDFmotherNew / xPDFdaughterNew;
1158
1159     // Iterate until acceptable pT (or have fallen below pTmin).
1160     } while (wt < rndmPtr->flat()) ;
1161   }
1162
1163   // Save values for (so far) acceptable branching.
1164   dipEndNow->store( idDaughter, idMother, idSister, x1Now, x2Now, m2Dip,
1165     pT2, z, xMother, Q2, mSister, m2Sister, pT2corr);
1166
1167 }
1168
1169 //--------------------------------------------------------------------------
1170
1171 // Kinematics of branching.
1172 // Construct mother -> daughter + sister, with recoiler on other side.
1173
1174 bool SpaceShower::branch( Event& event) {
1175
1176   // Side on which branching occured.
1177   int side          = abs(dipEndSel->side);
1178   double sideSign   = (side == 1) ? 1. : -1.;
1179
1180   // Read in flavour and colour variables.
1181   int iDaughter     = partonSystemsPtr->getInA(iSysSel);
1182   int iRecoiler     = partonSystemsPtr->getInB(iSysSel);
1183   if (side == 2) swap(iDaughter, iRecoiler);
1184   int idDaughterNow = dipEndSel->idDaughter;
1185   int idMother      = dipEndSel->idMother;
1186   int idSister      = dipEndSel->idSister;
1187   int colDaughter   = event[iDaughter].col();
1188   int acolDaughter  = event[iDaughter].acol();
1189
1190   // Recoil parton may be rescatterer, requiring special processing.
1191   bool normalRecoil = dipEndSel->normalRecoil;
1192   int iRecoilMother = event[iRecoiler].mother1();
1193
1194   // Read in kinematical variables.
1195   double x1         = dipEndSel->x1;
1196   double x2         = dipEndSel->x2;
1197   double xMo        = dipEndSel->xMo;
1198   double m2         = dipEndSel->m2Dip;
1199   double m          = sqrt(m2);
1200   double pT2        = dipEndSel->pT2;
1201   double z          = dipEndSel->z;
1202   double Q2         = dipEndSel->Q2;
1203   double mSister    = dipEndSel->mSister;
1204   double m2Sister   = dipEndSel->m2Sister;
1205   double pT2corr    = dipEndSel->pT2corr;
1206   double x1New      = (side == 1) ? xMo : x1;
1207   double x2New      = (side == 2) ? xMo : x2;
1208
1209   // Rescatter: kinematics may fail; use the rescatterFail flag to tell
1210   //            parton level to try again.
1211   rescatterFail     = false;
1212
1213   // Construct kinematics of mother, sister and recoiler in old rest frame.
1214   // Normally both mother and recoiler are taken massless.
1215   double eNewRec, pzNewRec, pTbranch, pzMother;
1216   if (normalRecoil) {
1217     eNewRec         = 0.5 * (m2 + Q2) / m;
1218     pzNewRec        = -sideSign * eNewRec;
1219     pTbranch        = sqrt(pT2corr) * m2 / ( z * (m2 + Q2) );
1220     pzMother        = sideSign * 0.5 * m * ( (m2 - Q2) / ( z * (m2 + Q2) )
1221                     + (Q2 + m2Sister) / m2 );
1222   // More complicated kinematics when recoiler not massless. May fail.
1223   } else {
1224     m2Rec           = event[iRecoiler].m2();
1225     double s1Tmp    = m2 + Q2 - m2Rec;
1226     double s3Tmp    = m2 / z - m2Rec;
1227     double r1Tmp    = sqrt(s1Tmp * s1Tmp + 4. * Q2 * m2Rec);
1228     eNewRec         = 0.5 * (m2 + m2Rec + Q2) / m;
1229     pzNewRec        = -sideSign * 0.5 * r1Tmp / m;
1230     double pT2br    = Q2 * s3Tmp * (m2 / z - m2 - Q2)
1231       - m2Sister * s1Tmp * s3Tmp - m2Rec * pow2(Q2 + m2Sister);
1232     if (pT2br <= 0.) return false;
1233     pTbranch        = sqrt(pT2br) / r1Tmp;
1234     pzMother        = sideSign * (m * s3Tmp
1235       - eNewRec * (m2 / z - Q2 - m2Rec - m2Sister)) / r1Tmp;
1236   }
1237   // Common final kinematics steps for both normal and rescattering.
1238   double eMother    = sqrt( pow2(pTbranch) + pow2(pzMother) );
1239   double pzSister   = pzMother + pzNewRec;
1240   double eSister    = sqrt( pow2(pTbranch) + pow2(pzSister) + m2Sister );
1241   Vec4 pMother( pTbranch, 0., pzMother, eMother );
1242   Vec4 pSister( pTbranch, 0., pzSister, eSister );
1243   Vec4 pNewRec(       0., 0., pzNewRec, eNewRec );
1244
1245   // Current event and subsystem size.
1246   int eventSizeOld  = event.size();
1247   int systemSizeOld = partonSystemsPtr->sizeAll(iSysSel);
1248
1249   // Save properties to be restored in case of user-hook veto of emission.
1250   int beamOff1 = 1 + beamOffset;
1251   int beamOff2 = 2 + beamOffset;
1252   int ev1Dau1V = event[beamOff1].daughter1();
1253   int ev2Dau1V = event[beamOff2].daughter1();
1254   vector<int> statusV, mother1V, mother2V, daughter1V, daughter2V;
1255   if (canVetoEmission) {
1256     for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
1257       int iOldCopy    = partonSystemsPtr->getAll(iSysSel, iCopy);
1258       statusV.push_back( event[iOldCopy].status());
1259       mother1V.push_back( event[iOldCopy].mother1());
1260       mother2V.push_back( event[iOldCopy].mother2());
1261       daughter1V.push_back( event[iOldCopy].daughter1());
1262       daughter2V.push_back( event[iOldCopy].daughter2());
1263     }
1264   }
1265
1266   // Take copy of existing system, to be given modified kinematics.
1267   // Incoming negative status. Rescattered also negative, but after copy.
1268   for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
1269     int iOldCopy    = partonSystemsPtr->getAll(iSysSel, iCopy);
1270     int statusOld   = event[iOldCopy].status();
1271     int statusNew   = (iOldCopy == iDaughter
1272       || iOldCopy == iRecoiler) ? statusOld : 44;
1273     int iNewCopy    = event.copy(iOldCopy, statusNew);
1274     if (statusOld < 0) event[iNewCopy].statusNeg();
1275   }
1276
1277   // Define colour flow in branching.
1278   // Default corresponds to f -> f + gamma.
1279   int colMother     = colDaughter;
1280   int acolMother    = acolDaughter;
1281   int colSister     = 0;
1282   int acolSister    = 0;
1283   if (idSister == 22) ;
1284   // q -> q + g and 50% of g -> g + g; need new colour.
1285   else if (idSister == 21 && ( (idMother > 0 && idMother < 9)
1286   || (idMother == 21 && rndmPtr->flat() < 0.5) ) ) {
1287     colMother       = event.nextColTag();
1288     colSister       = colMother;
1289     acolSister      = colDaughter;
1290   // qbar -> qbar + g and other 50% of g -> g + g; need new colour.
1291   } else if (idSister == 21) {
1292     acolMother      = event.nextColTag();
1293     acolSister      = acolMother;
1294     colSister       = acolDaughter;
1295   // q -> g + q.
1296   } else if (idDaughterNow == 21 && idMother > 0) {
1297     colMother       = colDaughter;
1298     acolMother      = 0;
1299     colSister       = acolDaughter;
1300   // qbar -> g + qbar
1301   } else if (idDaughterNow == 21) {
1302     acolMother      = acolDaughter;
1303     colMother       = 0;
1304     acolSister      = colDaughter;
1305   // g -> q + qbar.
1306   } else if (idDaughterNow > 0 && idDaughterNow < 9) {
1307     acolMother      = event.nextColTag();
1308     acolSister      = acolMother;
1309   // g -> qbar + q.
1310   } else if (idDaughterNow < 0 && idDaughterNow > -9) {
1311     colMother       = event.nextColTag();
1312     colSister       = colMother;
1313   // q -> gamma + q.
1314   } else if (idDaughterNow == 22 && idMother > 0) {
1315     colMother       = event.nextColTag();
1316     colSister       = colMother;
1317    // qbar -> gamma + qbar.
1318   } else if (idDaughterNow == 22) {
1319     acolMother      = event.nextColTag();
1320     acolSister      = acolMother;
1321   }
1322
1323   // Indices of partons involved. Add new sister.
1324   int iMother       = eventSizeOld + side - 1;
1325   int iNewRecoiler  = eventSizeOld + 2 - side;
1326   int iSister       = event.append( idSister, 43, iMother, 0, 0, 0,
1327      colSister, acolSister, pSister, mSister, sqrt(pT2) );
1328
1329   // References to the partons involved.
1330   Particle& daughter    = event[iDaughter];
1331   Particle& mother      = event[iMother];
1332   Particle& newRecoiler = event[iNewRecoiler];
1333   Particle& sister      = event.back();
1334
1335   // Replace old by new mother; update new recoiler.
1336   mother.id( idMother );
1337   mother.status( -41);
1338   mother.cols( colMother, acolMother);
1339   mother.p( pMother);
1340   newRecoiler.status( (normalRecoil) ? -42 : -46 );
1341   newRecoiler.p( pNewRec);
1342   if (!normalRecoil) newRecoiler.m( event[iRecoiler].m() );
1343
1344   // Update mother and daughter pointers; also for beams.
1345   daughter.mothers( iMother, 0);
1346   mother.daughters( iSister, iDaughter);
1347   if (iSysSel == 0) {
1348     event[beamOff1].daughter1( (side == 1) ? iMother : iNewRecoiler );
1349     event[beamOff2].daughter1( (side == 2) ? iMother : iNewRecoiler );
1350   }
1351
1352   // Find boost to old rest frame.
1353   RotBstMatrix Mtot;
1354   if (normalRecoil) Mtot.bst(0., 0., (x2 - x1) / (x1 + x2) );
1355   else if (side == 1)
1356        Mtot.toCMframe( event[iDaughter].p(), event[iRecoiler].p() );
1357   else Mtot.toCMframe( event[iRecoiler].p(), event[iDaughter].p() );
1358
1359   // Initially select phi angle of branching at random.
1360   double phi = 2. * M_PI * rndmPtr->flat();
1361
1362   // Evaluate coefficient of azimuthal asymmetry from gluon polarization.
1363   findAsymPol( event, dipEndSel);
1364   int    iFinPol = dipEndSel->iFinPol;
1365   double cPol    = dipEndSel->asymPol;
1366   double phiPol  = (iFinPol == 0) ? 0. : event[iFinPol].phi();
1367
1368   // If interference: try to match sister (anti)colour to final state.
1369   int    iFinInt = 0;
1370   double cInt    = 0.;
1371   double phiInt  = 0.;
1372   if (doPhiIntAsym) {
1373     for (int i = 0; i < partonSystemsPtr->sizeOut(iSysSel); ++ i) {
1374       int iOut = partonSystemsPtr->getOut(iSysSel, i);
1375       if ( (acolSister != 0 && event[iOut].col() == acolSister)
1376         || (colSister != 0 && event[iOut].acol() == colSister) )
1377         iFinInt = iOut;
1378     }
1379     if (iFinInt != 0) {
1380       // Boost final-state parton to current frame of new kinematics.
1381       Vec4 pFinTmp = event[iFinInt].p();
1382       pFinTmp.rotbst(Mtot);
1383       double theFin = pFinTmp.theta();
1384       if (side == 2) theFin = M_PI - theFin;
1385       double theSis = pSister.theta();
1386       if (side == 2) theSis = M_PI - theSis;
1387       cInt = strengthIntAsym * 2. * theSis * theFin
1388            / (pow2(theSis) + pow2(theFin));
1389       phiInt = event[iFinInt].phi();
1390     }
1391   }
1392
1393   // Bias phi distribution for polarization and interference.
1394   if (iFinPol != 0 || iFinInt != 0) {
1395     double cPhiPol, cPhiInt, weight;
1396     do {
1397       phi     = 2. * M_PI * rndmPtr->flat();
1398       weight  = 1.;
1399       if (iFinPol !=0 ) {
1400         cPhiPol = cos(phi - phiPol);
1401         weight *= ( 1. + cPol * (2. * pow2(cPhiPol) - 1.) )
1402           / ( 1. + abs(cPol) );
1403       }
1404       if (iFinInt !=0 ) {
1405         cPhiInt = cos(phi - phiInt);
1406         weight *= (1. - cInt) * (1. - cInt * cPhiInt)
1407           / (1. + pow2(cInt) - 2. * cInt * cPhiInt);
1408       }
1409     } while (weight < rndmPtr->flat());
1410   }
1411
1412   // Include rotation -phi on boost to old rest frame.
1413   Mtot.rot(0., -phi);
1414
1415   // Find boost from old rest frame to event cm frame.
1416   RotBstMatrix MfromRest;
1417   // The boost to the new rest frame.
1418   Vec4 sumNew       = pMother + pNewRec;
1419   double betaX      = sumNew.px() / sumNew.e();
1420   double betaZ      = sumNew.pz() / sumNew.e();
1421   MfromRest.bst( -betaX, 0., -betaZ);
1422   // Alignment of  radiator + recoiler to +- z axis, and rotation +phi.
1423   // Note: with spacelike (E < 0) recoiler p'_x_mother < 0 can happen!
1424   pMother.rotbst(MfromRest);
1425   double theta = pMother.theta();
1426   if (pMother.px() < 0.) theta = -theta;
1427   if (side == 2) theta += M_PI;
1428   MfromRest.rot(-theta, phi);
1429   // Boost to radiator + recoiler in event cm frame.
1430   if (normalRecoil) {
1431     MfromRest.bst( 0., 0., (x1New - x2New) / (x1New + x2New) );
1432   } else if (side == 1) {
1433     Vec4 pMotherWanted( 0., 0.,  0.5 * eCM * x1New, 0.5 * eCM * x1New);
1434     MfromRest.fromCMframe( pMotherWanted, event[iRecoiler].p() );
1435
1436   } else {
1437     Vec4 pMotherWanted( 0., 0., -0.5 * eCM * x2New, 0.5 * eCM * x2New);
1438     MfromRest.fromCMframe( event[iRecoiler].p(), pMotherWanted );
1439   }
1440   Mtot.rotbst(MfromRest);
1441
1442   // Perform cumulative rotation/boost operation.
1443   // Mother, recoiler and sister from old rest frame to event cm frame.
1444   mother.rotbst(MfromRest);
1445   newRecoiler.rotbst(MfromRest);
1446   sister.rotbst(MfromRest);
1447   // The rest from (and to) event cm frame.
1448   for ( int i = eventSizeOld + 2; i < eventSizeOld + systemSizeOld; ++i)
1449     event[i].rotbst(Mtot);
1450
1451   // Allow veto of branching. If so restore event record to before emission.
1452   if ( canVetoEmission
1453     && userHooksPtr->doVetoISREmission(eventSizeOld, event, iSysSel) ) {
1454     event.popBack( event.size() - eventSizeOld);
1455     event[beamOff1].daughter1( ev1Dau1V);
1456     event[beamOff2].daughter1( ev2Dau1V);
1457     for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
1458       int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
1459       event[iOldCopy].status( statusV[iCopy]);
1460       event[iOldCopy].mothers( mother1V[iCopy], mother2V[iCopy]);
1461       event[iOldCopy].daughters( daughter1V[iCopy], daughter2V[iCopy]);
1462     }
1463     return false;
1464   }
1465
1466   // Update list of partons in system; adding newly produced one.
1467   partonSystemsPtr->setInA(iSysSel, eventSizeOld);
1468   partonSystemsPtr->setInB(iSysSel, eventSizeOld + 1);
1469   for (int iCopy = 2; iCopy < systemSizeOld; ++iCopy)
1470     partonSystemsPtr->setOut(iSysSel, iCopy - 2, eventSizeOld + iCopy);
1471   partonSystemsPtr->addOut(iSysSel, eventSizeOld + systemSizeOld);
1472   partonSystemsPtr->setSHat(iSysSel, m2 / z);
1473
1474   // Update info on radiating dipole ends (QCD or QED).
1475   for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip)
1476   if ( dipEnd[iDip].system == iSysSel) {
1477     if (abs(dipEnd[iDip].side) == side) {
1478       dipEnd[iDip].iRadiator = iMother;
1479       dipEnd[iDip].iRecoiler = iNewRecoiler;
1480       if (dipEnd[iDip].side > 0)
1481         dipEnd[iDip].colType = mother.colType();
1482       else {
1483         dipEnd[iDip].chgType = 0;
1484         if ( (mother.isQuark() && doQEDshowerByQ)
1485           || (mother.isLepton() && doQEDshowerByL) )
1486           dipEnd[iDip].chgType = mother.chargeType();
1487       }
1488       // Kill ME corrections after first emission.
1489       dipEnd[iDip].MEtype = 0;
1490
1491     // Update info on recoiling dipole ends (QCD or QED).
1492     } else {
1493       dipEnd[iDip].iRadiator = iNewRecoiler;
1494       dipEnd[iDip].iRecoiler = iMother;
1495       // Optionally also kill recoiler ME corrections after first emission.
1496       if (!doMEafterFirst) dipEnd[iDip].MEtype = 0;
1497     }
1498   }
1499
1500   // Update info on beam remnants.
1501   BeamParticle& beamNow = (side == 1) ? *beamAPtr : *beamBPtr;
1502   double xNew = (side == 1) ? x1New : x2New;
1503   beamNow[iSysSel].update( iMother, idMother, xNew);
1504   // Redo choice of companion kind whenever new flavour.
1505   if (idMother != idDaughterNow) {
1506     beamNow.xfISR( iSysSel, idMother, xNew, factorMultFac * pT2);
1507     beamNow.pickValSeaComp();
1508   }
1509   BeamParticle& beamRec = (side == 1) ? *beamBPtr : *beamAPtr;
1510   beamRec[iSysSel].iPos( iNewRecoiler);
1511
1512   // Store branching values of current dipole. (For rapidity ordering.)
1513   ++dipEndSel->nBranch;
1514   dipEndSel->pT2Old = pT2;
1515   dipEndSel->zOld   = z;
1516
1517   // Update history if recoiler rescatters.
1518   if (!normalRecoil)
1519     event[iRecoilMother].daughters( iNewRecoiler, iNewRecoiler);
1520
1521   // Start list of rescatterers that force changed kinematics.
1522   vector<int> iRescatterer;
1523   for ( int i = 0; i < systemSizeOld - 2; ++i) {
1524     int iOutNew = partonSystemsPtr->getOut( iSysSel, i);
1525     if (!event[iOutNew].isFinal()) iRescatterer.push_back(iOutNew);
1526   }
1527
1528   // Start iterate over list of such rescatterers.
1529   int iRescNow = -1;
1530   while (++iRescNow < int(iRescatterer.size())) {
1531
1532     // Identify partons that induce or are affected by rescatter shift.
1533     // In following Old is before change of kinematics, New after,
1534     // Out scatterer in outstate and In in instate of another system.
1535     // Daughter sequence is (iOutOld ->) iOutNew -> iInNew -> iInOld.
1536     int iOutNew    = iRescatterer[iRescNow];
1537     int iInOld     = event[iOutNew].daughter1();
1538     int iSysResc   = partonSystemsPtr->getSystemOf(iInOld, true);
1539
1540     // Copy incoming partons of rescattered system and hook them up.
1541     int iOldA      = partonSystemsPtr->getInA(iSysResc);
1542     int iOldB      = partonSystemsPtr->getInB(iSysResc);
1543     bool rescSideA = event[iOldA].isRescatteredIncoming();
1544     int statusNewA = (rescSideA) ? -45 : -42;
1545     int statusNewB = (rescSideA) ? -42 : -45;
1546     int iNewA      = event.copy(iOldA, statusNewA);
1547     int iNewB      = event.copy(iOldB, statusNewB);
1548
1549     // Copy outgoing partons of rescattered system and hook them up.
1550     int eventSize  = event.size();
1551     int sizeOutAB  = partonSystemsPtr->sizeOut(iSysResc);
1552     int iOldAB, statusOldAB, iNewAB;
1553     for (int iOutAB = 0; iOutAB < sizeOutAB; ++iOutAB) {
1554       iOldAB       = partonSystemsPtr->getOut(iSysResc, iOutAB);
1555       statusOldAB  = event[iOldAB].status();
1556       iNewAB       = event.copy(iOldAB, 44);
1557       // Status could be negative for parton that rescatters in its turn.
1558       if (statusOldAB < 0) {
1559         event[iNewAB].statusNeg();
1560         iRescatterer.push_back(iNewAB);
1561       }
1562     }
1563
1564     // Hook up new outgoing with new incoming parton.
1565     int iInNew     = (rescSideA) ? iNewA : iNewB;
1566     event[iOutNew].daughters( iInNew, iInNew);
1567     event[iInNew].mothers( iOutNew, iOutNew);
1568
1569     // Rescale recoiling incoming parton for correct invariant mass.
1570     event[iInNew].p( event[iOutNew].p() );
1571     double momFac  = (rescSideA)
1572                    ? event[iInOld].pPos() / event[iInNew].pPos()
1573                    : event[iInOld].pNeg() / event[iInNew].pNeg();
1574     int iInRec     = (rescSideA) ? iNewB : iNewA;
1575
1576     // Rescatter: A previous boost may cause the light cone momentum of a
1577     //            rescattered parton to change sign. If this happens, tell
1578     //            parton level to try again.
1579     if (momFac < 0.0) {
1580       infoPtr->errorMsg("Warning in SpaceShower::branch: "
1581       "change in lightcone momentum sign; retrying parton level");
1582       rescatterFail = true;
1583       return false;
1584     }
1585
1586     event[iInRec].rescale4( momFac);
1587
1588     // Boost outgoing partons to new frame of incoming.
1589     RotBstMatrix MmodResc;
1590     MmodResc.toCMframe(  event[iOldA].p(), event[iOldB].p());
1591     MmodResc.fromCMframe(event[iNewA].p(), event[iNewB].p());
1592     for (int iOutAB = 0; iOutAB < sizeOutAB; ++iOutAB)
1593       event[eventSize + iOutAB].rotbst(MmodResc);
1594
1595     // Update list of partons in system.
1596     partonSystemsPtr->setInA(iSysResc, iNewA);
1597     partonSystemsPtr->setInB(iSysResc, iNewB);
1598     for (int iCopy = 0; iCopy < sizeOutAB; ++iCopy)
1599       partonSystemsPtr->setOut(iSysResc, iCopy, eventSize + iCopy);
1600
1601     // Update info on radiating dipole ends (QCD or QED).
1602     for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip)
1603     if ( dipEnd[iDip].system == iSysResc) {
1604       bool sideAnow = (abs(dipEnd[iDip].side) == 1);
1605       dipEnd[iDip].iRadiator = (sideAnow) ? iNewA : iNewB;
1606       dipEnd[iDip].iRecoiler = (sideAnow) ? iNewB : iNewA;
1607     }
1608
1609     // Update info on beam remnants.
1610     BeamParticle& beamResc = (rescSideA) ? *beamAPtr : *beamBPtr;
1611     beamResc[iSysResc].iPos( iInNew);
1612     beamResc[iSysResc].p( event[iInNew].p() );
1613     beamResc[iSysResc].scaleX( 1. / momFac  );
1614     BeamParticle& beamReco = (rescSideA) ? *beamBPtr : *beamAPtr;
1615     beamReco[iSysResc].iPos( iInRec);
1616     beamReco[iSysResc].scaleX( momFac);
1617
1618   // End iterate over list of rescatterers.
1619   }
1620
1621   // Check that beam momentum not used up by rescattered-system boosts.
1622     if (beamAPtr->xMax(-1) < 0.0 || beamBPtr->xMax(-1) < 0.0) {
1623       infoPtr->errorMsg("Warning in SpaceShower::branch: "
1624       "used up beam momentum; retrying parton level");
1625       rescatterFail = true;
1626       return false;
1627     }
1628
1629   // Done without any errors.
1630   return true;
1631
1632 }
1633
1634 //--------------------------------------------------------------------------
1635
1636 // Find class of ME correction.
1637
1638 int SpaceShower::findMEtype( int iSys, Event& event) {
1639
1640   // Default values and no action.
1641   int MEtype = 0;
1642   if (!doMEcorrections) ;
1643
1644   // Identify systems producing a single resonance.
1645   else if (partonSystemsPtr->sizeOut( iSys) == 1) {
1646     int idIn1 = event[partonSystemsPtr->getInA(iSys)].id();
1647     int idIn2 = event[partonSystemsPtr->getInA(iSys)].id();
1648     int idRes = event[partonSystemsPtr->getOut(iSys, 0)].id();
1649     if (iSys == 0) idResFirst  = abs(idRes);
1650     if (iSys == 1) idResSecond = abs(idRes);
1651
1652     // f + fbar -> vector boson.
1653     if ( (idRes == 23 || abs(idRes) == 24 || idRes == 32
1654       || idRes == 33 || abs(idRes) == 34 || abs(idRes) == 41)
1655       && abs(idIn1) < 20 && abs(idIn2) < 20 ) MEtype = 1;
1656
1657     // g + g, gamma + gamma  -> Higgs boson.
1658     if ( (idRes == 25 || idRes == 35 || idRes == 36)
1659        && ( ( idIn1 == 21 && idIn2 == 21 )
1660        || ( idIn1 == 22 && idIn2 == 22 ) ) ) MEtype = 2;
1661
1662     // f + fbar  -> Higgs boson.
1663     if ( (idRes == 25 || idRes == 35 || idRes == 36)
1664        && abs(idIn1) < 20 && abs(idIn2) < 20 ) MEtype = 3;
1665   }
1666
1667   // Done.
1668   return MEtype;
1669
1670 }
1671
1672 //--------------------------------------------------------------------------
1673
1674 // Provide maximum of expected ME weight; for preweighting of evolution.
1675
1676 double SpaceShower::calcMEmax( int MEtype, int idMother, int idDaughterIn) {
1677
1678   // Currently only one non-unity case: g(gamma) f -> V f'.
1679   if (MEtype == 1 && idMother > 20 && idDaughterIn < 20) return 3.;
1680   return 1.;
1681
1682 }
1683
1684 //--------------------------------------------------------------------------
1685
1686 // Provide actual ME weight for current branching.
1687 // Note: currently ME corrections are only allowed for first branching
1688 // on each side, so idDaughter is essentially known and checks overkill.
1689
1690 double SpaceShower::calcMEcorr(int MEtype, int idMother, int idDaughterIn,
1691   double M2, double z, double Q2) {
1692
1693   // Convert to Mandelstam variables. Sometimes may need to swap later.
1694   double sH = M2 / z;
1695   double tH = -Q2;
1696   double uH = Q2 - M2 * (1. - z) / z;
1697   int idMabs = abs(idMother);
1698   int idDabs = abs(idDaughterIn);
1699
1700   // Corrections for f + fbar -> s-channel vector boson.
1701   if (MEtype == 1) {
1702     if (idMabs < 20 && idDabs < 20) {
1703       return (tH*tH + uH*uH + 2. * M2 * sH) / (sH*sH + M2*M2);
1704     } else if (idDabs < 20) {
1705       // g(gamma) f -> V f': -Q2 = (p_g - p_f')^2 in PS while
1706       // tHat = (p_f - p_f')^2 in ME so need to swap tHat <-> uHat.
1707       swap( tH, uH);
1708       return (sH*sH + uH*uH + 2. * M2 * tH) / (pow2(sH - M2) + M2*M2);
1709     }
1710
1711   // Corrections for g + g -> Higgs boson.
1712   } else if (MEtype == 2) {
1713     if (idMabs < 20 && idDabs > 20) {
1714       return (sH*sH + uH*uH) / (sH*sH + pow2(sH - M2));
1715     } else if (idDabs > 20) {
1716       return 0.5 * (pow4(sH) + pow4(tH) + pow4(uH) + pow4(M2))
1717         / pow2(sH*sH - M2 * (sH - M2));
1718     }
1719
1720   // Corrections for f + fbar -> Higgs boson (f = b mainly).
1721   } else if (MEtype == 3) {
1722     if (idMabs < 20 && idDabs < 20) {
1723       // The PS and ME answers agree for f fbar -> H g/gamma.
1724       return 1.;
1725     } else if (idDabs < 20) {
1726       // Need to swap tHat <-> uHat, cf. vector-boson production above.
1727       swap( tH, uH);
1728       return (sH*sH + uH*uH + 2. * (M2 - uH) * (M2 - sH))
1729              / (pow2(sH - M2) + M2*M2);
1730     }
1731   }
1732
1733   return 1.;
1734
1735 }
1736
1737 //--------------------------------------------------------------------------
1738
1739 // Find coefficient of azimuthal asymmetry from gluon polarization.
1740
1741 void SpaceShower::findAsymPol( Event& event, SpaceDipoleEnd* dip) {
1742
1743   // Default is no asymmetry. Only gluons are studied.
1744   dip->iFinPol   = 0;
1745   dip->asymPol   = 0.;
1746   int iRad       = dip->iRadiator;
1747   if (!doPhiPolAsym || dip->idDaughter != 21) return;
1748
1749   // At least two particles in final state, whereof at least one coloured.
1750   int systemSizeOut = partonSystemsPtr->sizeOut( iSysSel);
1751   if (systemSizeOut < 2) return;
1752   bool foundColOut  = false;
1753   for (int ii = 0; ii < systemSizeOut; ++ii) {
1754     int i = partonSystemsPtr->getOut( iSysSel, ii);
1755     if (event[i].col() != 0 || event[i].acol() != 0) foundColOut = true;
1756   }
1757   if (!foundColOut) return;
1758
1759   // Check if granddaughter in final state of hard scattering.
1760   // (May need to trace across carbon copies to find granddaughters.)
1761   // If so, only accept 2 -> 2 scatterings with gg or qq in final state.
1762   int iGrandD1 = event[iRad].daughter1();
1763   int iGrandD2 = event[iRad].daughter2();
1764   bool traceCopy = false;
1765   do {
1766     traceCopy = false;
1767     if (iGrandD1 > 0 && iGrandD2 == iGrandD1) {
1768       iGrandD1 = event[iGrandD2].daughter1();
1769       iGrandD2 = event[iGrandD2].daughter2();
1770       traceCopy = true;
1771     }
1772   } while (traceCopy);
1773   int statusGrandD1 = event[ iGrandD1 ].statusAbs();
1774   bool isHardProc  = (statusGrandD1 == 23 || statusGrandD1 == 33);
1775   if (isHardProc) {
1776     if (iGrandD2 != iGrandD1 + 1) return;
1777     if (event[iGrandD1].isGluon() && event[iGrandD2].isGluon());
1778     else if (event[iGrandD1].isQuark() && event[iGrandD2].isQuark());
1779     else return;
1780   }
1781   dip->iFinPol = iGrandD1;
1782
1783   // Coefficient from gluon production.
1784   if (dip->idMother == 21) dip->asymPol = pow2( (1. - dip->z)
1785     / (1. - dip->z * (1. - dip->z) ) );
1786   else dip->asymPol = 2. * (1. - dip->z) / (1. + pow2(1. - dip->z) );
1787
1788   // Coefficients from gluon decay. Put z = 1/2 for hard process.
1789   double zDau  = (isHardProc) ? 0.5 : dip->zOld;
1790   if (event[iGrandD1].isGluon()) dip->asymPol *= pow2( (1. - zDau)
1791     / (1. - zDau * (1. - zDau) ) );
1792   else  dip->asymPol *= -2. * zDau *( 1. - zDau )
1793     / (1. - 2. * zDau * (1. - zDau) );
1794
1795 }
1796
1797 //--------------------------------------------------------------------------
1798
1799 // Print the list of dipoles.
1800
1801 void SpaceShower::list(ostream& os) const {
1802
1803   // Header.
1804   os << "\n --------  PYTHIA SpaceShower Dipole Listing  -------------- \n"
1805      << "\n    i  syst  side   rad   rec       pTmax  col  chg   ME rec \n"
1806      << fixed << setprecision(3);
1807
1808   // Loop over dipole list and print it.
1809   for (int i = 0; i < int(dipEnd.size()); ++i)
1810   os << setw(5) << i << setw(6) << dipEnd[i].system
1811      << setw(6) << dipEnd[i].side << setw(6) << dipEnd[i].iRadiator
1812      << setw(6) << dipEnd[i].iRecoiler << setw(12) << dipEnd[i].pTmax
1813      << setw(5) << dipEnd[i].colType << setw(5) << dipEnd[i].chgType
1814      << setw(5) << dipEnd[i].MEtype << setw(4)
1815      << dipEnd[i].normalRecoil << "\n";
1816
1817   // Done.
1818   os << "\n --------  End PYTHIA SpaceShower Dipole Listing  ----------"
1819      << endl;
1820
1821
1822 }
1823
1824 //==========================================================================
1825
1826 } // end namespace Pythia8
1827