PYTHIA8/pythia8145/src/SpaceShower.cxx

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