--- /dev/null
+// SpaceShower.cc is a part of the PYTHIA event generator.
+// Copyright (C) 2012 Torbjorn Sjostrand.
+// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
+// Please respect the MCnet Guidelines, see GUIDELINES for details.
+
+// Function definitions (not found in the header) for the
+// SpaceShower class.
+
+#include "SpaceShower.h"
+
+namespace Pythia8 {
+
+//==========================================================================
+
+// The SpaceShower class.
+
+//--------------------------------------------------------------------------
+
+// Constants: could be changed here if desired, but normally should not.
+// These are of technical nature, as described for each.
+
+// Leftover companion can give PDF > 0 at small Q2 where other PDF's = 0,
+// and then one can end in infinite loop of impossible kinematics.
+const int SpaceShower::MAXLOOPTINYPDF = 10;
+
+// Switch to alternative (but equivalent) backwards evolution for
+// g -> Q Qbar (Q = c or b) when below QTHRESHOLD * mQ2.
+const double SpaceShower::CTHRESHOLD = 2.0;
+const double SpaceShower::BTHRESHOLD = 2.0;
+
+// Renew evaluation of PDF's when the pT2 step is bigger than this
+// (in addition to initial scale and c and b thresholds.)
+const double SpaceShower::EVALPDFSTEP = 0.1;
+
+// Lower limit on PDF value in order to avoid division by zero.
+const double SpaceShower::TINYPDF = 1e-10;
+
+// Lower limit on estimated evolution rate, below which stop.
+const double SpaceShower::TINYKERNELPDF = 1e-6;
+
+// Lower limit on pT2, below which branching is rejected.
+const double SpaceShower::TINYPT2 = 0.25e-6;
+
+// No attempt to do backwards evolution of a heavy (c or b) quark
+// if evolution starts at a scale pT2 < HEAVYPT2EVOL * mQ2.
+const double SpaceShower::HEAVYPT2EVOL = 1.1;
+
+// No attempt to do backwards evolution of a heavy (c or b) quark
+// if evolution starts at a x > HEAVYXEVOL * x_max, where
+// x_max is the largest possible x value for a g -> Q Qbar branching.
+const double SpaceShower::HEAVYXEVOL = 0.9;
+
+// When backwards evolution Q -> g + Q creates a heavy quark Q,
+// an earlier branching g -> Q + Qbar will restrict kinematics
+// to M_{Q Qbar}^2 > EXTRASPACEQ * 4 m_Q^2. (Smarter to be found??)
+const double SpaceShower::EXTRASPACEQ = 2.0;
+
+// Never pick pT so low that alphaS is evaluated too close to Lambda_3.
+const double SpaceShower::LAMBDA3MARGIN = 1.1;
+
+// Do not warn for large PDF ratios at small pT2 scales.
+const double SpaceShower::PT2MINWARN = 1.;
+
+// Cutoff for f_e^e at x < 1 - 10^{-10} to be used in z selection.
+// Note: the x_min quantity come from 1 - x_max.
+const double SpaceShower::LEPTONXMIN = 1e-10;
+const double SpaceShower::LEPTONXMAX = 1. - 1e-10;
+
+// Stop l -> l gamma evolution slightly above m2l.
+const double SpaceShower::LEPTONPT2MIN = 1.2;
+
+// Enhancement of l -> l gamma trial rate to compensate imperfect modelling.
+const double SpaceShower::LEPTONFUDGE = 10.;
+
+//--------------------------------------------------------------------------
+
+// Initialize alphaStrong, alphaEM and related pTmin parameters.
+
+void SpaceShower::init( BeamParticle* beamAPtrIn,
+ BeamParticle* beamBPtrIn) {
+
+ // Store input pointers for future use.
+ beamAPtr = beamAPtrIn;
+ beamBPtr = beamBPtrIn;
+
+ // Main flags to switch on and off branchings.
+ doQCDshower = settingsPtr->flag("SpaceShower:QCDshower");
+ doQEDshowerByQ = settingsPtr->flag("SpaceShower:QEDshowerByQ");
+ doQEDshowerByL = settingsPtr->flag("SpaceShower:QEDshowerByL");
+
+ // Matching in pT of hard interaction to shower evolution.
+ pTmaxMatch = settingsPtr->mode("SpaceShower:pTmaxMatch");
+ pTdampMatch = settingsPtr->mode("SpaceShower:pTdampMatch");
+ pTmaxFudge = settingsPtr->parm("SpaceShower:pTmaxFudge");
+ pTmaxFudgeMPI = settingsPtr->parm("SpaceShower:pTmaxFudgeMPI");
+ pTdampFudge = settingsPtr->parm("SpaceShower:pTdampFudge");
+
+ // Optionally force emissions to be ordered in rapidity/angle.
+ doRapidityOrder = settingsPtr->flag("SpaceShower:rapidityOrder");
+
+ // Charm, bottom and lepton mass thresholds.
+ mc = particleDataPtr->m0(4);
+ mb = particleDataPtr->m0(5);
+ m2c = pow2(mc);
+ m2b = pow2(mb);
+
+ // Parameters of scale choices.
+ renormMultFac = settingsPtr->parm("SpaceShower:renormMultFac");
+ factorMultFac = settingsPtr->parm("SpaceShower:factorMultFac");
+
+ // Parameters of alphaStrong generation.
+ alphaSvalue = settingsPtr->parm("SpaceShower:alphaSvalue");
+ alphaSorder = settingsPtr->mode("SpaceShower:alphaSorder");
+ alphaS2pi = 0.5 * alphaSvalue / M_PI;
+
+ // Initialize alpha_strong generation.
+ alphaS.init( alphaSvalue, alphaSorder);
+
+ // Lambda for 5, 4 and 3 flavours.
+ Lambda5flav = alphaS.Lambda5();
+ Lambda4flav = alphaS.Lambda4();
+ Lambda3flav = alphaS.Lambda3();
+ Lambda5flav2 = pow2(Lambda5flav);
+ Lambda4flav2 = pow2(Lambda4flav);
+ Lambda3flav2 = pow2(Lambda3flav);
+
+ // Regularization of QCD evolution for pT -> 0. Can be taken
+ // same as for multiparton interactions, or be set separately.
+ useSamePTasMPI = settingsPtr->flag("SpaceShower:samePTasMPI");
+ if (useSamePTasMPI) {
+ pT0Ref = settingsPtr->parm("MultipartonInteractions:pT0Ref");
+ ecmRef = settingsPtr->parm("MultipartonInteractions:ecmRef");
+ ecmPow = settingsPtr->parm("MultipartonInteractions:ecmPow");
+ pTmin = settingsPtr->parm("MultipartonInteractions:pTmin");
+ } else {
+ pT0Ref = settingsPtr->parm("SpaceShower:pT0Ref");
+ ecmRef = settingsPtr->parm("SpaceShower:ecmRef");
+ ecmPow = settingsPtr->parm("SpaceShower:ecmPow");
+ pTmin = settingsPtr->parm("SpaceShower:pTmin");
+ }
+
+ // Calculate nominal invariant mass of events. Set current pT0 scale.
+ sCM = m2( beamAPtr->p(), beamBPtr->p());
+ eCM = sqrt(sCM);
+ pT0 = pT0Ref * pow(eCM / ecmRef, ecmPow);
+
+ // Restrict pTmin to ensure that alpha_s(pTmin^2 + pT_0^2) does not blow up.
+ double pTminAbs = sqrtpos(pow2(LAMBDA3MARGIN) * Lambda3flav2 / renormMultFac
+ - pT0*pT0);
+ if (pTmin < pTminAbs) {
+ pTmin = pTminAbs;
+ ostringstream newPTmin;
+ newPTmin << fixed << setprecision(3) << pTmin;
+ infoPtr->errorMsg("Warning in SpaceShower::init: pTmin too low",
+ ", raised to " + newPTmin.str() );
+ infoPtr->setTooLowPTmin(true);
+ }
+
+ // Parameters of alphaEM generation.
+ alphaEMorder = settingsPtr->mode("SpaceShower:alphaEMorder");
+
+ // Initialize alphaEM generation.
+ alphaEM.init( alphaEMorder, settingsPtr);
+
+ // Parameters of QED evolution.
+ pTminChgQ = settingsPtr->parm("SpaceShower:pTminchgQ");
+ pTminChgL = settingsPtr->parm("SpaceShower:pTminchgL");
+
+ // Derived parameters of QCD evolution.
+ pT20 = pow2(pT0);
+ pT2min = pow2(pTmin);
+ pT2minChgQ = pow2(pTminChgQ);
+ pT2minChgL = pow2(pTminChgL);
+
+ // Various other parameters.
+ doMEcorrections = settingsPtr->flag("SpaceShower:MEcorrections");
+ doMEafterFirst = settingsPtr->flag("SpaceShower:MEafterFirst");
+ doPhiPolAsym = settingsPtr->flag("SpaceShower:phiPolAsym");
+ doPhiIntAsym = settingsPtr->flag("SpaceShower:phiIntAsym");
+ strengthIntAsym = settingsPtr->parm("SpaceShower:strengthIntAsym");
+ nQuarkIn = settingsPtr->mode("SpaceShower:nQuarkIn");
+
+ // Optional dampening at small pT's when large multiplicities.
+ enhanceScreening
+ = settingsPtr->mode("MultipartonInteractions:enhanceScreening");
+ if (!useSamePTasMPI) enhanceScreening = 0;
+
+ // Possibility to allow user veto of emission step.
+ canVetoEmission = (userHooksPtr != 0)
+ ? userHooksPtr->canVetoISREmission() : false;
+
+}
+
+//--------------------------------------------------------------------------
+
+// Find whether to limit maximum scale of emissions.
+// Also allow for dampening at factorization or renormalization scale.
+
+bool SpaceShower::limitPTmax( Event& event, double Q2Fac, double Q2Ren) {
+
+ // Find whether to limit pT. Begin by user-set cases.
+ bool dopTlimit = false;
+ if (pTmaxMatch == 1) dopTlimit = true;
+ else if (pTmaxMatch == 2) dopTlimit = false;
+
+ // Look if any quark (u, d, s, c, b), gluon or photon in final state.
+ else {
+ for (int i = 5; i < event.size(); ++i)
+ if (event[i].status() != -21) {
+ int idAbs = event[i].idAbs();
+ if (idAbs <= 5 || idAbs == 21 || idAbs == 22) dopTlimit = true;
+ }
+ }
+
+ // Dampening at factorization or renormalization scale.
+ dopTdamp = false;
+ pT2damp = 0.;
+ if ( !dopTlimit && (pTdampMatch == 1 || pTdampMatch == 2) ) {
+ dopTdamp = true;
+ pT2damp = pow2(pTdampFudge) * ((pTdampMatch == 1) ? Q2Fac : Q2Ren);
+ }
+
+ // Done.
+ return dopTlimit;
+
+}
+
+//--------------------------------------------------------------------------
+
+// Prepare system for evolution; identify ME.
+// Routine may be called after multiparton interactions, for a new subystem.
+
+void SpaceShower::prepare( int iSys, Event& event, bool limitPTmaxIn) {
+
+ // Find positions of incoming colliding partons.
+ int in1 = partonSystemsPtr->getInA(iSys);
+ int in2 = partonSystemsPtr->getInB(iSys);
+
+ // Rescattered partons cannot radiate.
+ bool canRadiate1 = !(event[in1].isRescatteredIncoming());
+ bool canRadiate2 = !(event[in2].isRescatteredIncoming());
+
+ // Reset dipole-ends list for first interaction. Also resonances.
+ if (iSys == 0) dipEnd.resize(0);
+ if (iSys == 0) idResFirst = 0;
+ if (iSys == 1) idResSecond = 0;
+
+ // Find matrix element corrections for system.
+ int MEtype = findMEtype( iSys, event);
+
+ // Maximum pT scale for dipole ends.
+ double pTmax1 = (limitPTmaxIn) ? event[in1].scale() : eCM;
+ double pTmax2 = (limitPTmaxIn) ? event[in2].scale() : eCM;
+ if (iSys == 0 && limitPTmaxIn) {
+ pTmax1 *= pTmaxFudge;
+ pTmax2 *= pTmaxFudge;
+ } else if (iSys > 0 && limitPTmaxIn) {
+ pTmax1 *= pTmaxFudgeMPI;
+ pTmax2 *= pTmaxFudgeMPI;
+ }
+
+ // Find dipole ends for QCD radiation.
+ // Note: colour type can change during evolution, so book also if zero.
+ if (doQCDshower) {
+ int colType1 = event[in1].colType();
+ if (canRadiate1) dipEnd.push_back( SpaceDipoleEnd( iSys, 1,
+ in1, in2, pTmax1, colType1, 0, MEtype, canRadiate2) );
+ int colType2 = event[in2].colType();
+ if (canRadiate2) dipEnd.push_back( SpaceDipoleEnd( iSys, 2,
+ in2, in1, pTmax2, colType2, 0, MEtype, canRadiate1) );
+ }
+
+ // Find dipole ends for QED radiation.
+ // Note: charge type can change during evolution, so book also if zero.
+ if (doQEDshowerByQ || doQEDshowerByL) {
+ int chgType1 = ( (event[in1].isQuark() && doQEDshowerByQ)
+ || (event[in1].isLepton() && doQEDshowerByL) )
+ ? event[in1].chargeType() : 0;
+ // Special: photons have charge zero, but can evolve (only off Q for now)
+ if (event[in1].id() == 22 && doQEDshowerByQ) chgType1 = 22 ;
+ if (canRadiate1) dipEnd.push_back( SpaceDipoleEnd( iSys, -1,
+ in1, in2, pTmax1, 0, chgType1, MEtype, canRadiate2) );
+ int chgType2 = ( (event[in2].isQuark() && doQEDshowerByQ)
+ || (event[in2].isLepton() && doQEDshowerByL) )
+ ? event[in2].chargeType() : 0;
+ // Special: photons have charge zero, but can evolve (only off Q for now)
+ if (event[in2].id() == 22 && doQEDshowerByQ) chgType2 = 22 ;
+ if (canRadiate2) dipEnd.push_back( SpaceDipoleEnd( iSys, -2,
+ in2, in1, pTmax2, 0, chgType2, MEtype, canRadiate1) );
+ }
+
+ // Store the z and pT2 values of the last previous splitting
+ // when an event history has already been constructed.
+ if (iSys == 0 && infoPtr->hasHistory()) {
+ double zNow = infoPtr->zNowISR();
+ double pT2Now = infoPtr->pT2NowISR();
+ for (int iDipEnd = 0; iDipEnd < int(dipEnd.size()); ++iDipEnd) {
+ dipEnd[iDipEnd].zOld = zNow;
+ dipEnd[iDipEnd].pT2Old = pT2Now;
+ ++dipEnd[iDipEnd].nBranch;
+ }
+ }
+
+}
+
+//--------------------------------------------------------------------------
+
+// Select next pT in downwards evolution of the existing dipoles.
+
+double SpaceShower::pTnext( Event& event, double pTbegAll, double pTendAll,
+ int nRadIn) {
+
+ // Current cm energy, in case it varies between events.
+ sCM = m2( beamAPtr->p(), beamBPtr->p());
+ eCM = sqrt(sCM);
+ pTbegRef = pTbegAll;
+
+ // Starting values: no radiating dipole found.
+ nRad = nRadIn;
+ double pT2sel = pow2(pTendAll);
+ iDipSel = 0;
+ iSysSel = 0;
+ dipEndSel = 0;
+
+ // Loop over all possible dipole ends.
+ for (int iDipEnd = 0; iDipEnd < int(dipEnd.size()); ++iDipEnd) {
+ iDipNow = iDipEnd;
+ dipEndNow = &dipEnd[iDipEnd];
+ iSysNow = dipEndNow->system;
+ dipEndNow->pT2 = 0.;
+
+ // Check whether dipole end should be allowed to shower.
+ double pT2begDip = pow2( min( pTbegAll, dipEndNow->pTmax ));
+ if (pT2begDip > pT2sel
+ && ( dipEndNow->colType != 0 || dipEndNow->chgType != 0 ) ) {
+ double pT2endDip = 0.;
+
+ // Determine lower cut for evolution, for QCD or QED (q or l).
+ if (dipEndNow->colType != 0) pT2endDip = max( pT2sel, pT2min );
+ else if (abs(dipEndNow->chgType) != 3) pT2endDip
+ = max( pT2sel, pT2minChgQ );
+ else pT2endDip = max( pT2sel, pT2minChgL );
+
+ // Find properties of dipole and radiating dipole end.
+ sideA = ( abs(dipEndNow->side) == 1 );
+ BeamParticle& beamNow = (sideA) ? *beamAPtr : *beamBPtr;
+ BeamParticle& beamRec = (sideA) ? *beamBPtr : *beamAPtr;
+ iNow = beamNow[iSysNow].iPos();
+ iRec = beamRec[iSysNow].iPos();
+ idDaughter = beamNow[iSysNow].id();
+ xDaughter = beamNow[iSysNow].x();
+ x1Now = (sideA) ? xDaughter : beamRec[iSysNow].x();
+ x2Now = (sideA) ? beamRec[iSysNow].x() : xDaughter;
+
+ // Note dipole mass correction when recoiler is a rescatter.
+ m2Rec = (dipEndNow->normalRecoil) ? 0. : event[iRec].m2();
+ m2Dip = x1Now * x2Now * sCM + m2Rec;
+
+ // Now do evolution in pT2, for QCD or QED
+ if (pT2begDip > pT2endDip) {
+ if (dipEndNow->colType != 0) pT2nextQCD( pT2begDip, pT2endDip);
+ else pT2nextQED( pT2begDip, pT2endDip);
+ }
+
+ // Update if found larger pT than current maximum.
+ if (dipEndNow->pT2 > pT2sel) {
+ pT2sel = dipEndNow->pT2;
+ iDipSel = iDipNow;
+ iSysSel = iSysNow;
+ dipEndSel = dipEndNow;
+ }
+
+ // End loop over dipole ends.
+ }
+ }
+
+ // Return nonvanishing value if found pT is bigger than already found.
+ return (dipEndSel == 0) ? 0. : sqrt(pT2sel);
+}
+
+//--------------------------------------------------------------------------
+
+// Evolve a QCD dipole end.
+
+void SpaceShower::pT2nextQCD( double pT2begDip, double pT2endDip) {
+
+ // Some properties and kinematical starting values.
+ BeamParticle& beam = (sideA) ? *beamAPtr : *beamBPtr;
+ bool isGluon = (idDaughter == 21);
+ bool isValence = beam[iSysNow].isValence();
+ int MEtype = dipEndNow->MEtype;
+ double pT2 = pT2begDip;
+ double xMaxAbs = beam.xMax(iSysNow);
+ double zMinAbs = xDaughter / xMaxAbs;
+ if (xMaxAbs < 0.) {
+ infoPtr->errorMsg("Warning in SpaceShower::pT2nextQCD: "
+ "xMaxAbs negative");
+ return;
+ }
+
+ // Starting values for handling of massive quarks (c/b), if any.
+ double idMassive = 0;
+ if ( abs(idDaughter) == 4 ) idMassive = 4;
+ if ( abs(idDaughter) == 5 ) idMassive = 5;
+ bool isMassive = (idMassive > 0);
+ double m2Massive = 0.;
+ double mRatio = 0.;
+ double zMaxMassive = 1.;
+ double m2Threshold = pT2;
+
+ // Evolution below scale of massive quark or at large x is impossible.
+ if (isMassive) {
+ m2Massive = (idMassive == 4) ? m2c : m2b;
+ if (pT2 < HEAVYPT2EVOL * m2Massive) return;
+ mRatio = sqrt( m2Massive / m2Dip );
+ zMaxMassive = (1. - mRatio) / ( 1. + mRatio * (1. - mRatio) );
+ if (xDaughter > HEAVYXEVOL * zMaxMassive * xMaxAbs) return;
+
+ // Find threshold scale below which only g -> Q + Qbar will be allowed.
+ m2Threshold = (idMassive == 4) ? min( pT2, CTHRESHOLD * m2c)
+ : min( pT2, BTHRESHOLD * m2b);
+ }
+
+ // Variables used inside evolution loop. (Mainly dummy starting values.)
+ int nFlavour = 3;
+ double b0 = 4.5;
+ double Lambda2 = Lambda3flav2;
+ double pT2minNow = pT2endDip;
+ int idMother = 0;
+ int idSister = 0;
+ double z = 0.;
+ double zMaxAbs = 0.;
+ double zRootMax = 0.;
+ double zRootMin = 0.;
+ double g2gInt = 0.;
+ double q2gInt = 0.;
+ double q2qInt = 0.;
+ double g2qInt = 0.;
+ double g2Qenhance = 0.;
+ double xPDFdaughter = 0.;
+ double xPDFmother[21] = {0.};
+ double xPDFgMother = 0.;
+ double xPDFmotherSum = 0.;
+ double kernelPDF = 0.;
+ double xMother = 0.;
+ double wt = 0.;
+ double Q2 = 0.;
+ double mSister = 0.;
+ double m2Sister = 0.;
+ double pT2corr = 0.;
+ double pT2PDF = pT2;
+ bool needNewPDF = true;
+
+ // Begin evolution loop towards smaller pT values.
+ int loopTinyPDFdau = 0;
+ bool hasTinyPDFdau = false;
+ do {
+ wt = 0.;
+
+ // Bad sign if repeated looping with small daughter PDF, so fail.
+ // (Example: if all PDF's = 0 below Q_0, except for c/b companion.)
+ if (hasTinyPDFdau) ++loopTinyPDFdau;
+ if (loopTinyPDFdau > MAXLOOPTINYPDF) {
+ infoPtr->errorMsg("Warning in SpaceShower::pT2nextQCD: "
+ "small daughter PDF");
+ return;
+ }
+
+ // Initialize integrals of splitting kernels and evaluate parton
+ // densities at the beginning. Reinitialize after long evolution
+ // in pT2 or when crossing c and b flavour thresholds.
+ if (needNewPDF || pT2 < EVALPDFSTEP * pT2PDF) {
+ pT2PDF = pT2;
+ hasTinyPDFdau = false;
+
+ // Determine overestimated z range; switch at c and b masses.
+ if (pT2 > m2b) {
+ nFlavour = 5;
+ pT2minNow = max( m2b, pT2endDip);
+ b0 = 23./6.;
+ Lambda2 = Lambda5flav2;
+ } else if (pT2 > m2c) {
+ nFlavour = 4;
+ pT2minNow = max( m2c, pT2endDip);
+ b0 = 25./6.;
+ Lambda2 = Lambda4flav2;
+ } else {
+ nFlavour = 3;
+ pT2minNow = pT2endDip;
+ b0 = 27./6.;
+ Lambda2 = Lambda3flav2;
+ }
+ // A change of renormalization scale expressed by a change of Lambda.
+ Lambda2 /= renormMultFac;
+ zMaxAbs = 1. - 0.5 * (pT2minNow / m2Dip) *
+ ( sqrt( 1. + 4. * m2Dip / pT2minNow ) - 1. );
+ if (isMassive) zMaxAbs = min( zMaxAbs, zMaxMassive);
+
+ // Go to another z range with lower mass scale if current is closed.
+ if (zMinAbs > zMaxAbs) {
+ if (nFlavour == 3 || (idMassive == 4 && nFlavour == 4)
+ || idMassive == 5) return;
+ pT2 = (nFlavour == 4) ? m2c : m2b;
+ continue;
+ }
+
+ // Parton density of daughter at current scale.
+ xPDFdaughter = beam.xfISR(iSysNow, idDaughter, xDaughter,
+ factorMultFac * pT2);
+ if (xPDFdaughter < TINYPDF) {
+ xPDFdaughter = TINYPDF;
+ hasTinyPDFdau = true;
+ }
+
+ // Integrals of splitting kernels for gluons: g -> g, q -> g.
+ if (isGluon) {
+ g2gInt = 6. * log(zMaxAbs * (1.-zMinAbs)
+ / (zMinAbs * (1.-zMaxAbs)));
+ if (doMEcorrections) g2gInt *= calcMEmax(MEtype, 21, 21);
+ q2gInt = (16./3.) * (1./sqrt(zMinAbs) - 1./sqrt(zMaxAbs));
+ if (doMEcorrections) q2gInt *= calcMEmax(MEtype, 1, 21);
+
+ // Parton density of potential quark mothers to a g.
+ xPDFmotherSum = 0.;
+ for (int i = -nQuarkIn; i <= nQuarkIn; ++i) {
+ if (i == 0) {
+ xPDFmother[10] = 0.;
+ } else {
+ xPDFmother[i+10] = beam.xfISR(iSysNow, i, xDaughter,
+ factorMultFac * pT2);
+ xPDFmotherSum += xPDFmother[i+10];
+ }
+ }
+
+ // Total QCD evolution coefficient for a gluon.
+ kernelPDF = g2gInt + q2gInt * xPDFmotherSum / xPDFdaughter;
+
+ // For valence quark only need consider q -> q g branchings.
+ // Introduce an extra factor sqrt(z) to smooth bumps.
+ } else if (isValence) {
+ zRootMin = (1. + sqrt(zMinAbs)) / (1. - sqrt(zMinAbs));
+ zRootMax = (1. + sqrt(zMaxAbs)) / (1. - sqrt(zMaxAbs));
+ q2qInt = (8./3.) * log( zRootMax / zRootMin );
+ if (doMEcorrections) q2qInt *= calcMEmax(MEtype, 1, 1);
+ kernelPDF = q2qInt;
+
+ // Integrals of splitting kernels for quarks: q -> q, g -> q.
+ } else {
+ q2qInt = (8./3.) * log( (1. - zMinAbs) / (1. - zMaxAbs) );
+ if (doMEcorrections) q2qInt *= calcMEmax(MEtype, 1, 1);
+ g2qInt = 0.5 * (zMaxAbs - zMinAbs);
+ if (doMEcorrections) g2qInt *= calcMEmax(MEtype, 21, 1);
+
+ // Increase estimated upper weight for g -> Q + Qbar.
+ if (isMassive) {
+ if (alphaSorder == 0) g2Qenhance = log(pT2/m2Massive)
+ / log(m2Threshold/m2Massive);
+ else {
+ double m2log = log( m2Massive / Lambda2);
+ g2Qenhance = log( log(pT2/Lambda2) / m2log )
+ / log( log(m2Threshold/Lambda2) / m2log );
+ }
+ g2qInt *= g2Qenhance;
+ }
+
+ // Parton density of a potential gluon mother to a q.
+ xPDFgMother = beam.xfISR(iSysNow, 21, xDaughter, factorMultFac * pT2);
+
+ // Total QCD evolution coefficient for a quark.
+ kernelPDF = q2qInt + g2qInt * xPDFgMother / xPDFdaughter;
+ }
+
+ // End evaluation of splitting kernels and parton densities.
+ needNewPDF = false;
+ }
+ if (kernelPDF < TINYKERNELPDF) return;
+
+ // Pick pT2 (in overestimated z range), for one of three different cases.
+ // Assume form alphas(pT0^2 + pT^2) * dpT^2/(pT0^2 + pT^2).
+ double Q2alphaS;
+
+ // Fixed alpha_strong.
+ if (alphaSorder == 0) {
+ pT2 = (pT2 + pT20) * pow( rndmPtr->flat(),
+ 1. / (alphaS2pi * kernelPDF)) - pT20;
+
+ // First-order alpha_strong.
+ } else if (alphaSorder == 1) {
+ pT2 = Lambda2 * pow( (pT2 + pT20) / Lambda2,
+ pow(rndmPtr->flat(), b0 / kernelPDF) ) - pT20;
+
+ // For second order reject by second term in alpha_strong expression.
+ } else {
+ do {
+ pT2 = Lambda2 * pow( (pT2 + pT20) / Lambda2,
+ pow(rndmPtr->flat(), b0 / kernelPDF) ) - pT20;
+ Q2alphaS = renormMultFac * max( pT2 + pT20,
+ pow2(LAMBDA3MARGIN) * Lambda3flav2);
+ } while (alphaS.alphaS2OrdCorr(Q2alphaS) < rndmPtr->flat()
+ && pT2 > pT2minNow);
+ }
+
+ // Check for pT2 values that prompt special action.
+
+ // If fallen into b threshold region, force g -> b + bbar.
+ if (idMassive == 5 && pT2 < m2Threshold) {
+ pT2nearQCDthreshold( beam, m2Massive, m2Threshold, xMaxAbs,
+ zMinAbs, zMaxMassive );
+ return;
+
+ // If crossed b threshold, continue evolution from this threshold.
+ } else if (nFlavour == 5 && pT2 < m2b) {
+ needNewPDF = true;
+ pT2 = m2b;
+ continue;
+
+ // If fallen into c threshold region, force g -> c + cbar.
+ } else if (idMassive == 4 && pT2 < m2Threshold) {
+ pT2nearQCDthreshold( beam, m2Massive, m2Threshold, xMaxAbs,
+ zMinAbs, zMaxMassive );
+ return;
+
+ // If crossed c threshold, continue evolution from this threshold.
+ } else if (nFlavour == 4 && pT2 < m2c) {
+ needNewPDF = true;
+ pT2 = m2c;
+ continue;
+
+ // Abort evolution if below cutoff scale, or below another branching.
+ } else if (pT2 < pT2endDip) return;
+
+ // Select z value of branching to g, and corrective weight.
+ if (isGluon) {
+ // g -> g (+ g).
+ if (rndmPtr->flat() * kernelPDF < g2gInt) {
+ idMother = 21;
+ idSister = 21;
+ z = 1. / ( 1. + ((1. - zMinAbs) / zMinAbs) * pow( (zMinAbs *
+ (1. - zMaxAbs)) / (zMaxAbs * (1. - zMinAbs)), rndmPtr->flat() ) );
+ wt = pow2( 1. - z * (1. - z));
+ } else {
+ // q -> g (+ q): also select flavour.
+ double temp = xPDFmotherSum * rndmPtr->flat();
+ idMother = -nQuarkIn - 1;
+ do { temp -= xPDFmother[(++idMother) + 10]; }
+ while (temp > 0. && idMother < nQuarkIn);
+ idSister = idMother;
+ z = (zMinAbs * zMaxAbs) / pow2( sqrt(zMinAbs) + rndmPtr->flat()
+ * ( sqrt(zMaxAbs)- sqrt(zMinAbs) ));
+ wt = 0.5 * (1. + pow2(1. - z)) * sqrt(z)
+ * xPDFdaughter / xPDFmother[idMother + 10];
+ }
+
+ // Select z value of branching to q, and corrective weight.
+ // Include massive kernel corrections for c and b quarks.
+ } else {
+ // q -> q (+ g).
+ if (isValence || rndmPtr->flat() * kernelPDF < q2qInt) {
+ idMother = idDaughter;
+ idSister = 21;
+ // Valence more peaked at large z.
+ if (isValence) {
+ double zTmp = zRootMin * pow(zRootMax / zRootMin, rndmPtr->flat() );
+ z = pow2( (1. - zTmp) / (1. + zTmp) );
+ } else {
+ z = 1. - (1. - zMinAbs) * pow( (1. - zMaxAbs) / (1. - zMinAbs),
+ rndmPtr->flat() );
+ }
+ if (!isMassive) {
+ wt = 0.5 * (1. + pow2(z));
+ } else {
+ //?? Bug?? should be 2 more for massive part??
+ // wt = 0.5 * (1. + pow2(z) - z * pow2(1.-z) * m2Massive / pT2);
+ wt = 0.5 * (1. + pow2(z)) - z * pow2(1.-z) * m2Massive / pT2;
+ }
+ if (isValence) wt *= sqrt(z);
+ // g -> q (+ qbar).
+ } else {
+ idMother = 21;
+ idSister = - idDaughter;
+ z = zMinAbs + rndmPtr->flat() * (zMaxAbs - zMinAbs);
+ if (!isMassive) {
+ wt = (pow2(z) + pow2(1.-z)) * xPDFdaughter / xPDFgMother ;
+ } else {
+ wt = (pow2(z) + pow2(1.-z) + 2. * z * (1.-z) * m2Massive / pT2)
+ * xPDFdaughter / (xPDFgMother * g2Qenhance) ;
+ }
+ }
+ }
+
+ // Derive Q2 and x of mother from pT2 and z.
+ Q2 = pT2 / (1. - z);
+ xMother = xDaughter / z;
+ // Correction to x for massive recoiler from rescattering.
+ if (!dipEndNow->normalRecoil) {
+ if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
+ else xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
+ }
+ if(xMother > xMaxAbs) { wt = 0.; continue; }
+
+ // Forbidden emission if outside allowed z range for given pT2.
+ mSister = particleDataPtr->m0(idSister);
+ m2Sister = pow2(mSister);
+ pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
+ if(pT2corr < TINYPT2) { wt = 0.; continue; }
+
+ // Optionally veto emissions not ordered in rapidity (= angle).
+ if ( doRapidityOrder && dipEndNow->nBranch > 0
+ && pT2 > pow2( (1. - z) / (z * (1. - dipEndNow->zOld)) )
+ * dipEndNow->pT2Old ) { wt = 0.; continue; }
+
+ // If creating heavy quark by Q -> g + Q then next need g -> Q + Qbar.
+ // So minimum total mass2 is 4 * m2Sister, but use more to be safe.
+ if ( isGluon && ( abs(idMother) == 4 || abs(idMother) == 5 )) {
+ double m2QQsister = EXTRASPACEQ * 4. * m2Sister;
+ double pT2QQcorr = Q2 - z * (m2Dip + Q2) * (Q2 + m2QQsister) / m2Dip;
+ if(pT2QQcorr < TINYPT2) { wt = 0.; continue; }
+ }
+
+ // Evaluation of ME correction.
+ if (doMEcorrections) wt *= calcMEcorr(MEtype, idMother, idDaughter,
+ m2Dip, z, Q2) / calcMEmax(MEtype, idMother, idDaughter);
+
+ // Optional dampening of large pT values in first radiation.
+ if (dopTdamp && iSysNow == 0 && MEtype == 0 && nRad == 0)
+ wt *= pT2damp / (pT2 + pT2damp);
+
+ // Idea suggested by Gosta Gustafson: increased screening in events
+ // with large activity can be simulated by pT0_eff = sqrt(n) * pT0.
+ if (enhanceScreening == 2) {
+ int nSysNow = infoPtr->nMPI() + infoPtr->nISR() + 1;
+ double WTscreen = pow2( (pT2 + pT20) / (pT2 + nSysNow * pT20) );
+ wt *= WTscreen;
+ }
+
+ // Evaluation of new daughter and mother PDF's.
+ double xPDFdaughterNew = max ( TINYPDF,
+ beam.xfISR(iSysNow, idDaughter, xDaughter, factorMultFac * pT2) );
+ double xPDFmotherNew =
+ beam.xfISR(iSysNow, idMother, xMother, factorMultFac * pT2);
+ wt *= xPDFmotherNew / xPDFdaughterNew;
+
+ // Check that valence step does not cause problem.
+ if (wt > 1. && pT2 > PT2MINWARN) infoPtr->errorMsg("Warning in "
+ "SpaceShower::pT2nextQCD: weight above unity");
+
+ // Iterate until acceptable pT (or have fallen below pTmin).
+ } while (wt < rndmPtr->flat()) ;
+
+ // Save values for (so far) acceptable branching.
+ dipEndNow->store( idDaughter,idMother, idSister, x1Now, x2Now, m2Dip,
+ pT2, z, xMother, Q2, mSister, m2Sister, pT2corr);
+
+}
+
+//--------------------------------------------------------------------------
+
+// Evolve a QCD dipole end near threshold, with g -> Q + Qbar enforced.
+// Note: No explicit Sudakov factor formalism here. Instead use that
+// df_Q(x, pT2) = (alpha_s/2pi) * (dT2/pT2) * ((gluon) * (splitting)).
+// This implies that effects of Q -> Q + g are neglected in this range.
+
+void SpaceShower::pT2nearQCDthreshold( BeamParticle& beam,
+ double m2Massive, double m2Threshold, double xMaxAbs,
+ double zMinAbs, double zMaxMassive) {
+
+ // Initial values, to be used in kinematics and weighting.
+ double Lambda2 = (abs(idDaughter) == 4) ? Lambda4flav2 : Lambda5flav2;
+ Lambda2 /= renormMultFac;
+ double logM2Lambda2 = log( m2Massive / Lambda2 );
+ double xPDFmotherOld = beam.xfISR(iSysNow, 21, xDaughter,
+ factorMultFac * m2Threshold);
+
+ // Variables used inside evolution loop. (Mainly dummy start values.)
+ int loop = 0;
+ double wt = 0.;
+ double pT2 = 0.;
+ double z = 0.;
+ double Q2 = 0.;
+ double pT2corr = 0.;
+ double xMother = 0.;
+
+ // Begin loop over tries to find acceptable g -> Q + Qbar branching.
+ do {
+ wt = 0.;
+
+ // Check that not caught in infinite loop with impossible kinematics.
+ if (++loop > 100) {
+ infoPtr->errorMsg("Error in SpaceShower::pT2nearQCDthreshold: "
+ "stuck in loop");
+ return;
+ }
+
+ // Pick dpT2/pT2 in range [m2Massive,thresholdRatio * m2Massive].
+ pT2 = m2Massive * pow( m2Threshold / m2Massive, rndmPtr->flat() );
+
+ // Pick z flat in allowed range.
+ z = zMinAbs + rndmPtr->flat() * (zMaxMassive - zMinAbs);
+
+ // Check that kinematically possible choice.
+ Q2 = pT2 / (1.-z) - m2Massive;
+ pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Massive) / m2Dip;
+ if(pT2corr < TINYPT2) continue;
+
+ // Correction factor for running alpha_s. ??
+ wt = logM2Lambda2 / log( pT2 / Lambda2 );
+
+ // Correction factor for splitting kernel.
+ wt *= pow2(z) + pow2(1.-z) + 2. * z * (1.-z) * m2Massive / pT2;
+
+ // x, including correction for massive recoiler from rescattering.
+ xMother = xDaughter / z;
+ if (!dipEndNow->normalRecoil) {
+ if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
+ else xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
+ }
+ if (xMother > xMaxAbs) { wt = 0.; continue; }
+
+ // Correction factor for gluon density.
+ double xPDFmotherNew = beam.xfISR(iSysNow, 21, xMother,
+ factorMultFac * pT2);
+ wt *= xPDFmotherNew / xPDFmotherOld;
+
+ // Iterate until acceptable pT and z.
+ } while (wt < rndmPtr->flat()) ;
+
+ // Save values for (so far) acceptable branching.
+ double mSister = (abs(idDaughter) == 4) ? mc : mb;
+ dipEndNow->store( idDaughter, 21, -idDaughter, x1Now, x2Now, m2Dip,
+ pT2, z, xMother, Q2, mSister, pow2(mSister), pT2corr);
+
+}
+
+//--------------------------------------------------------------------------
+
+// Evolve a QED dipole end.
+
+void SpaceShower::pT2nextQED( double pT2begDip, double pT2endDip) {
+
+ // Type of dipole and starting values.
+ BeamParticle& beam = (sideA) ? *beamAPtr : *beamBPtr;
+ bool isLeptonBeam = beam.isLepton();
+ int MEtype = dipEndNow->MEtype;
+ bool isPhoton = (idDaughter == 22);
+ double pT2 = pT2begDip;
+ double m2Lepton = (isLeptonBeam) ? pow2(beam.m()) : 0.;
+ if (isLeptonBeam && pT2begDip < m2Lepton) return;
+
+ // Currently no f -> gamma branching implemented for lepton beams
+ if (isPhoton && isLeptonBeam) return;
+
+ // alpha_em at maximum scale provides upper estimate.
+ double alphaEMmax = alphaEM.alphaEM(renormMultFac * pT2begDip);
+ double alphaEM2pi = alphaEMmax / (2. * M_PI);
+
+ // Maximum x of mother implies minimum z = xDaughter / xMother.
+ double xMaxAbs = (isLeptonBeam) ? LEPTONXMAX : beam.xMax(iSysNow);
+ double zMinAbs = xDaughter / xMaxAbs;
+ if (xMaxAbs < 0.) {
+ infoPtr->errorMsg("Warning in SpaceShower::pT2nextQED: "
+ "xMaxAbs negative");
+ return;
+ }
+
+ // Maximum z from minimum pT and, for lepton, from minimum x_gamma.
+ double zMaxAbs = 1. - 0.5 * (pT2endDip / m2Dip) *
+ ( sqrt( 1. + 4. * m2Dip / pT2endDip ) - 1. );
+ if (isLeptonBeam) {
+ double zMaxLepton = xDaughter / (xDaughter + LEPTONXMIN);
+ if (zMaxLepton < zMaxAbs) zMaxAbs = zMaxLepton;
+ }
+ if (zMaxAbs < zMinAbs) return;
+
+ // Variables used inside evolution loop. (Mainly dummy start values.)
+ int idMother = 0;
+ int idSister =22;
+ double z = 0.;
+ double xMother = 0.;
+ double wt = 0.;
+ double Q2 = 0.;
+ double mSister = 0.;
+ double m2Sister = 0.;
+ double pT2corr = 0.;
+
+ // QED evolution of fermions
+ if (!isPhoton) {
+
+ // Integrals of splitting kernels for fermions: f -> f. Use 1 + z^2 < 2.
+ // Ansatz f(z) = 2 / (1 - z), with + 2 / (z - xDaughter) for lepton.
+ double f2fInt = 0.;
+ double f2fIntA = 2. * log( (1. - zMinAbs) / (1. - zMaxAbs) );
+ double f2fIntB = 0.;
+ if (isLeptonBeam) {
+ f2fIntB = 2. * log( (zMaxAbs - xDaughter) / (zMinAbs - xDaughter) );
+ f2fInt = f2fIntA + f2fIntB;
+ } else f2fInt = pow2(dipEndNow->chgType / 3.) * f2fIntA;
+
+ // Upper estimate for evolution equation, including fudge factor.
+ if (doMEcorrections) f2fInt *= calcMEmax(MEtype, 1, 1);
+ double kernelPDF = alphaEM2pi * f2fInt;
+ double fudge = (isLeptonBeam) ? LEPTONFUDGE * log(m2Dip/m2Lepton) : 1.;
+ kernelPDF *= fudge;
+ if (kernelPDF < TINYKERNELPDF) return;
+
+ // Begin evolution loop towards smaller pT values.
+ do {
+
+ // Pick pT2 (in overestimated z range).
+ // For l -> l gamma include extrafactor 1 / ln(pT2 / m2l) in evolution.
+ double shift = pow(rndmPtr->flat(), 1. / kernelPDF);
+ if (isLeptonBeam) pT2 = m2Lepton * pow( pT2 / m2Lepton, shift);
+ else pT2 = pT2 * shift;
+
+ // Abort evolution if below cutoff scale, or below another branching.
+ if (pT2 < pT2endDip) return;
+ if (isLeptonBeam && pT2 < LEPTONPT2MIN * m2Lepton) return;
+
+ // Select z value of branching f -> f + gamma, and corrective weight.
+ idMother = idDaughter;
+ wt = 0.5 * (1. + pow2(z));
+ if (isLeptonBeam) {
+ if (f2fIntA > rndmPtr->flat() * (f2fIntA + f2fIntB)) {
+ z = 1. - (1. - zMinAbs)
+ * pow( (1. - zMaxAbs) / (1. - zMinAbs), rndmPtr->flat() );
+ } else {
+ z = xDaughter + (zMinAbs - xDaughter)
+ * pow( (zMaxAbs - xDaughter) / (zMinAbs - xDaughter),
+ rndmPtr->flat() );
+ }
+ wt *= (z - xDaughter) / (1. - xDaughter);
+ } else {
+ z = 1. - (1. - zMinAbs)
+ * pow( (1. - zMaxAbs) / (1. - zMinAbs), rndmPtr->flat() );
+ }
+
+ // Derive Q2 and x of mother from pT2 and z.
+ Q2 = pT2 / (1. - z);
+ xMother = xDaughter / z;
+ // Correction to x for massive recoiler from rescattering.
+ if (!dipEndNow->normalRecoil) {
+ if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
+ else xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
+ }
+ if(xMother > xMaxAbs) { wt = 0.; continue; }
+
+ // Forbidden emission if outside allowed z range for given pT2.
+ mSister = 0.;
+ m2Sister = 0.;
+ pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
+ if(pT2corr < TINYPT2) { wt = 0.; continue; }
+
+ // Correct by ln(pT2 / m2l) and fudge factor.
+ if (isLeptonBeam) wt *= log(pT2 / m2Lepton) / fudge;
+
+ // Evaluation of ME correction.
+ if (doMEcorrections) wt *= calcMEcorr(MEtype, idMother, idDaughter,
+ m2Dip, z, Q2) / calcMEmax(MEtype, idMother, idDaughter);
+
+ // Extra QED correction for f fbar -> W+- gamma. Debug??
+ if (doMEcorrections && MEtype == 1 && idDaughter == idMother
+ && ( (iSysNow == 0 && idResFirst == 24)
+ || (iSysNow == 1 && idResSecond == 24) ) ) {
+ double tHe = -Q2;
+ double uHe = Q2 - m2Dip * (1. - z) / z;
+ double chg1 = abs(dipEndNow->chgType / 3.);
+ double chg2 = 1. - chg1;
+ wt *= pow2(chg1 * uHe - chg2 * tHe)
+ / ( (tHe + uHe) * (pow2(chg1) * uHe + pow2(chg2) * tHe) );
+ }
+
+ // Optional dampening of large pT values in first radiation.
+ if (dopTdamp && iSysNow == 0 && MEtype == 0 && nRad == 0)
+ wt *= pT2damp / (pT2 + pT2damp);
+
+ // Correct to current value of alpha_EM.
+ double alphaEMnow = alphaEM.alphaEM(renormMultFac * pT2);
+ wt *= (alphaEMnow / alphaEMmax);
+
+ // Evaluation of new daughter and mother PDF's.
+ double xPDFdaughterNew = max ( TINYPDF,
+ beam.xfISR(iSysNow, idDaughter, xDaughter, factorMultFac * pT2) );
+ double xPDFmotherNew =
+ beam.xfISR(iSysNow, idMother, xMother, factorMultFac * pT2);
+ wt *= xPDFmotherNew / xPDFdaughterNew;
+
+ // Iterate until acceptable pT (or have fallen below pTmin).
+ } while (wt < rndmPtr->flat()) ;
+ }
+
+ // QED evolution of photons (so far only for hadron beams).
+ else {
+
+ // Initial values
+ int nFlavour = 3;
+ double kernelPDF = 0.0;
+ double xPDFdaughter = 0.;
+ double xPDFmother[21] = {0.};
+ double xPDFmotherSum = 0.0;
+ double pT2PDF = pT2;
+ double pT2minNow = pT2endDip;
+ bool needNewPDF = true;
+
+ // Begin evolution loop towards smaller pT values.
+ int loopTinyPDFdau = 0;
+ bool hasTinyPDFdau = false;
+ do {
+ wt = 0.;
+
+ // Bad sign if repeated looping with small daughter PDF, so fail.
+ if (hasTinyPDFdau) ++loopTinyPDFdau;
+ if (loopTinyPDFdau > MAXLOOPTINYPDF) {
+ infoPtr->errorMsg("Warning in SpaceShower::pT2nextQED: "
+ "small daughter PDF");
+ return;
+ }
+
+ // Initialize integrals of splitting kernels and evaluate parton
+ // densities at the beginning. Reinitialize after long evolution
+ // in pT2 or when crossing c and b flavour thresholds.
+ if (needNewPDF || pT2 < EVALPDFSTEP * pT2PDF) {
+
+ pT2PDF = pT2;
+ hasTinyPDFdau = false;
+
+ // Determine overestimated z range; switch at c and b masses.
+ if (pT2 > m2b && nQuarkIn >= 5) {
+ nFlavour = 5;
+ pT2minNow = max( m2b, pT2endDip);
+ } else if (pT2 > m2c && nQuarkIn >= 4) {
+ nFlavour = 4;
+ pT2minNow = max( m2c, pT2endDip);
+ } else {
+ nFlavour = 3;
+ pT2minNow = pT2endDip;
+ }
+
+ // Compute upper z limit
+ zMaxAbs = 1. - 0.5 * (pT2minNow / m2Dip) *
+ ( sqrt( 1. + 4. * m2Dip / pT2minNow ) - 1. );
+
+ // Parton density of daughter at current scale.
+ xPDFdaughter = beam.xfISR(iSysNow, idDaughter, xDaughter,
+ factorMultFac * pT2);
+ if (xPDFdaughter < TINYPDF) {
+ xPDFdaughter = TINYPDF;
+ hasTinyPDFdau = true;
+ }
+
+ // Integral over f -> gamma f splitting kernel.
+ // Normalized so: 4/3 aS/2pi P(z) -> eq^2 * aEM/2pi P(z).
+ // (Charge-weighting happens below.)
+ double q2gInt = 4. * (1./sqrt(zMinAbs) - 1./sqrt(zMaxAbs));
+
+ // Charge-weighted Parton density of potential quark mothers.
+ xPDFmotherSum = 0.;
+ for (int i = -nFlavour; i <= nFlavour; ++i) {
+ if (i == 0) {
+ xPDFmother[10] = 0.;
+ } else {
+ xPDFmother[i+10] = pow2((abs(i+1) % 2 + 1)/3.0)
+ * beam.xfISR(iSysNow, i, xDaughter, factorMultFac * pT2);
+ xPDFmotherSum += xPDFmother[i+10];
+ }
+ }
+
+ // Total QED evolution coefficient for a photon.
+ kernelPDF = q2gInt * xPDFmotherSum / xPDFdaughter;
+
+ // End evaluation of splitting kernels and parton densities.
+ needNewPDF = false;
+ }
+ if (kernelPDF < TINYKERNELPDF) return;
+
+ // Select pT2 for next trial branching
+ pT2 *= pow( rndmPtr->flat(), 1. / (alphaEM2pi * kernelPDF));
+
+ // If crossed b threshold, continue evolution from this threshold.
+ if (nFlavour == 5 && pT2 < m2b) {
+ needNewPDF = true;
+ pT2 = m2b;
+ continue;
+ }
+
+ // If crossed c threshold, continue evolution from this threshold.
+ else if (nFlavour == 4 && pT2 < m2c) {
+ needNewPDF = true;
+ pT2 = m2c;
+ continue;
+ }
+
+ // Abort evolution if below cutoff scale, or below another branching.
+ else if (pT2 < pT2endDip) return;
+
+ // Select flavour for trial branching
+ double temp = xPDFmotherSum * rndmPtr->flat();
+ idMother = -nQuarkIn - 1;
+ do {
+ temp -= xPDFmother[(++idMother) + 10];
+ } while (temp > 0. && idMother < nQuarkIn);
+
+ // Sister is same as mother, but can have m2 > 0
+ idSister = idMother;
+ mSister = particleDataPtr->m0(idSister);
+ m2Sister = pow2(mSister);
+
+ // Select z for trial branching
+ z = (zMinAbs * zMaxAbs) / pow2( sqrt(zMinAbs) + rndmPtr->flat()
+ * ( sqrt(zMaxAbs)- sqrt(zMinAbs) ));
+
+ // Trial weight: splitting kernel
+ wt = 0.5 * (1. + pow2(1. - z)) * sqrt(z);
+
+ // Trial weight: running alpha_EM
+ double alphaEMnow = alphaEM.alphaEM(renormMultFac * pT2);
+ wt *= (alphaEMnow / alphaEMmax);
+
+ // Derive Q2 and x of mother from pT2 and z
+ Q2 = pT2 / (1. - z);
+ xMother = xDaughter / z;
+ // Correction to x for massive recoiler from rescattering.
+ if (!dipEndNow->normalRecoil) {
+ if (sideA) xMother += (m2Rec / (x2Now * sCM)) * (1. / z - 1.);
+ else xMother += (m2Rec / (x1Now * sCM)) * (1. / z - 1.);
+ }
+
+ // Compute pT2corr
+ pT2corr = Q2 - z * (m2Dip + Q2) * (Q2 + m2Sister) / m2Dip;
+ if(pT2corr < TINYPT2) { wt = 0.; continue; }
+
+ // If creating heavy quark by Q -> gamma + Q then next g -> Q + Qbar.
+ // So minimum total mass2 is 4 * m2Sister, but use more to be safe.
+ if ( abs(idMother) == 4 || abs(idMother) == 5 ) {
+ double m2QQsister = EXTRASPACEQ * 4. * m2Sister;
+ double pT2QQcorr = Q2 - z * (m2Dip + Q2) * (Q2 + m2QQsister) / m2Dip;
+ if(pT2QQcorr < TINYPT2) { wt = 0.; continue; }
+ }
+
+ // Optional dampening of large pT values in first radiation.
+ if (dopTdamp && iSysNow == 0 && MEtype == 0 && nRad == 0)
+ wt *= pT2damp / (pT2 + pT2damp);
+
+ // Evaluation of new daughter PDF
+ double xPDFdaughterNew = beam.xfISR(iSysNow, idDaughter, xDaughter,
+ factorMultFac * pT2);
+ if (xPDFdaughterNew < TINYPDF) {
+ xPDFdaughterNew = TINYPDF;
+ }
+
+ // Evaluation of new charge-weighted mother PDF
+ double xPDFmotherNew = pow2( (abs(idMother+1) % 2 + 1)/3.0 )
+ * beam.xfISR(iSysNow, idMother, xMother, factorMultFac * pT2);
+
+ // Trial weight: divide out old pdf ratio
+ wt *= xPDFdaughter / xPDFmother[idMother + 10];
+
+ // Trial weight: new pdf ratio
+ wt *= xPDFmotherNew / xPDFdaughterNew;
+
+ // Iterate until acceptable pT (or have fallen below pTmin).
+ } while (wt < rndmPtr->flat()) ;
+ }
+
+ // Save values for (so far) acceptable branching.
+ dipEndNow->store( idDaughter, idMother, idSister, x1Now, x2Now, m2Dip,
+ pT2, z, xMother, Q2, mSister, m2Sister, pT2corr);
+
+}
+
+//--------------------------------------------------------------------------
+
+// Kinematics of branching.
+// Construct mother -> daughter + sister, with recoiler on other side.
+
+bool SpaceShower::branch( Event& event) {
+
+ // Side on which branching occured.
+ int side = abs(dipEndSel->side);
+ double sideSign = (side == 1) ? 1. : -1.;
+
+ // Read in flavour and colour variables.
+ int iDaughter = partonSystemsPtr->getInA(iSysSel);
+ int iRecoiler = partonSystemsPtr->getInB(iSysSel);
+ if (side == 2) swap(iDaughter, iRecoiler);
+ int idDaughterNow = dipEndSel->idDaughter;
+ int idMother = dipEndSel->idMother;
+ int idSister = dipEndSel->idSister;
+ int colDaughter = event[iDaughter].col();
+ int acolDaughter = event[iDaughter].acol();
+
+ // Recoil parton may be rescatterer, requiring special processing.
+ bool normalRecoil = dipEndSel->normalRecoil;
+ int iRecoilMother = event[iRecoiler].mother1();
+
+ // Read in kinematical variables.
+ double x1 = dipEndSel->x1;
+ double x2 = dipEndSel->x2;
+ double xMo = dipEndSel->xMo;
+ double m2 = dipEndSel->m2Dip;
+ double m = sqrt(m2);
+ double pT2 = dipEndSel->pT2;
+ double z = dipEndSel->z;
+ double Q2 = dipEndSel->Q2;
+ double mSister = dipEndSel->mSister;
+ double m2Sister = dipEndSel->m2Sister;
+ double pT2corr = dipEndSel->pT2corr;
+ double x1New = (side == 1) ? xMo : x1;
+ double x2New = (side == 2) ? xMo : x2;
+
+ // Rescatter: kinematics may fail; use the rescatterFail flag to tell
+ // parton level to try again.
+ rescatterFail = false;
+
+ // Construct kinematics of mother, sister and recoiler in old rest frame.
+ // Normally both mother and recoiler are taken massless.
+ double eNewRec, pzNewRec, pTbranch, pzMother;
+ if (normalRecoil) {
+ eNewRec = 0.5 * (m2 + Q2) / m;
+ pzNewRec = -sideSign * eNewRec;
+ pTbranch = sqrt(pT2corr) * m2 / ( z * (m2 + Q2) );
+ pzMother = sideSign * 0.5 * m * ( (m2 - Q2) / ( z * (m2 + Q2) )
+ + (Q2 + m2Sister) / m2 );
+ // More complicated kinematics when recoiler not massless. May fail.
+ } else {
+ m2Rec = event[iRecoiler].m2();
+ double s1Tmp = m2 + Q2 - m2Rec;
+ double s3Tmp = m2 / z - m2Rec;
+ double r1Tmp = sqrt(s1Tmp * s1Tmp + 4. * Q2 * m2Rec);
+ eNewRec = 0.5 * (m2 + m2Rec + Q2) / m;
+ pzNewRec = -sideSign * 0.5 * r1Tmp / m;
+ double pT2br = Q2 * s3Tmp * (m2 / z - m2 - Q2)
+ - m2Sister * s1Tmp * s3Tmp - m2Rec * pow2(Q2 + m2Sister);
+ if (pT2br <= 0.) return false;
+ pTbranch = sqrt(pT2br) / r1Tmp;
+ pzMother = sideSign * (m * s3Tmp
+ - eNewRec * (m2 / z - Q2 - m2Rec - m2Sister)) / r1Tmp;
+ }
+ // Common final kinematics steps for both normal and rescattering.
+ double eMother = sqrt( pow2(pTbranch) + pow2(pzMother) );
+ double pzSister = pzMother + pzNewRec;
+ double eSister = sqrt( pow2(pTbranch) + pow2(pzSister) + m2Sister );
+ Vec4 pMother( pTbranch, 0., pzMother, eMother );
+ Vec4 pSister( pTbranch, 0., pzSister, eSister );
+ Vec4 pNewRec( 0., 0., pzNewRec, eNewRec );
+
+ // Current event and subsystem size.
+ int eventSizeOld = event.size();
+ int systemSizeOld = partonSystemsPtr->sizeAll(iSysSel);
+
+ // Save properties to be restored in case of user-hook veto of emission.
+ int beamOff1 = 1 + beamOffset;
+ int beamOff2 = 2 + beamOffset;
+ int ev1Dau1V = event[beamOff1].daughter1();
+ int ev2Dau1V = event[beamOff2].daughter1();
+ vector<int> statusV, mother1V, mother2V, daughter1V, daughter2V;
+ if (canVetoEmission) {
+ for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
+ int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
+ statusV.push_back( event[iOldCopy].status());
+ mother1V.push_back( event[iOldCopy].mother1());
+ mother2V.push_back( event[iOldCopy].mother2());
+ daughter1V.push_back( event[iOldCopy].daughter1());
+ daughter2V.push_back( event[iOldCopy].daughter2());
+ }
+ }
+
+ // Take copy of existing system, to be given modified kinematics.
+ // Incoming negative status. Rescattered also negative, but after copy.
+ for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
+ int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
+ int statusOld = event[iOldCopy].status();
+ int statusNew = (iOldCopy == iDaughter
+ || iOldCopy == iRecoiler) ? statusOld : 44;
+ int iNewCopy = event.copy(iOldCopy, statusNew);
+ if (statusOld < 0) event[iNewCopy].statusNeg();
+ }
+
+ // Define colour flow in branching.
+ // Default corresponds to f -> f + gamma.
+ int colMother = colDaughter;
+ int acolMother = acolDaughter;
+ int colSister = 0;
+ int acolSister = 0;
+ if (idSister == 22) ;
+ // q -> q + g and 50% of g -> g + g; need new colour.
+ else if (idSister == 21 && ( (idMother > 0 && idMother < 9)
+ || (idMother == 21 && rndmPtr->flat() < 0.5) ) ) {
+ colMother = event.nextColTag();
+ colSister = colMother;
+ acolSister = colDaughter;
+ // qbar -> qbar + g and other 50% of g -> g + g; need new colour.
+ } else if (idSister == 21) {
+ acolMother = event.nextColTag();
+ acolSister = acolMother;
+ colSister = acolDaughter;
+ // q -> g + q.
+ } else if (idDaughterNow == 21 && idMother > 0) {
+ colMother = colDaughter;
+ acolMother = 0;
+ colSister = acolDaughter;
+ // qbar -> g + qbar
+ } else if (idDaughterNow == 21) {
+ acolMother = acolDaughter;
+ colMother = 0;
+ acolSister = colDaughter;
+ // g -> q + qbar.
+ } else if (idDaughterNow > 0 && idDaughterNow < 9) {
+ acolMother = event.nextColTag();
+ acolSister = acolMother;
+ // g -> qbar + q.
+ } else if (idDaughterNow < 0 && idDaughterNow > -9) {
+ colMother = event.nextColTag();
+ colSister = colMother;
+ // q -> gamma + q.
+ } else if (idDaughterNow == 22 && idMother > 0) {
+ colMother = event.nextColTag();
+ colSister = colMother;
+ // qbar -> gamma + qbar.
+ } else if (idDaughterNow == 22) {
+ acolMother = event.nextColTag();
+ acolSister = acolMother;
+ }
+
+ // Indices of partons involved. Add new sister.
+ int iMother = eventSizeOld + side - 1;
+ int iNewRecoiler = eventSizeOld + 2 - side;
+ int iSister = event.append( idSister, 43, iMother, 0, 0, 0,
+ colSister, acolSister, pSister, mSister, sqrt(pT2) );
+
+ // References to the partons involved.
+ Particle& daughter = event[iDaughter];
+ Particle& mother = event[iMother];
+ Particle& newRecoiler = event[iNewRecoiler];
+ Particle& sister = event.back();
+
+ // Replace old by new mother; update new recoiler.
+ mother.id( idMother );
+ mother.status( -41);
+ mother.cols( colMother, acolMother);
+ mother.p( pMother);
+ newRecoiler.status( (normalRecoil) ? -42 : -46 );
+ newRecoiler.p( pNewRec);
+ if (!normalRecoil) newRecoiler.m( event[iRecoiler].m() );
+
+ // Update mother and daughter pointers; also for beams.
+ daughter.mothers( iMother, 0);
+ mother.daughters( iSister, iDaughter);
+ if (iSysSel == 0) {
+ event[beamOff1].daughter1( (side == 1) ? iMother : iNewRecoiler );
+ event[beamOff2].daughter1( (side == 2) ? iMother : iNewRecoiler );
+ }
+
+ // Find boost to old rest frame.
+ RotBstMatrix Mtot;
+ if (normalRecoil) Mtot.bst(0., 0., (x2 - x1) / (x1 + x2) );
+ else if (side == 1)
+ Mtot.toCMframe( event[iDaughter].p(), event[iRecoiler].p() );
+ else Mtot.toCMframe( event[iRecoiler].p(), event[iDaughter].p() );
+
+ // Initially select phi angle of branching at random.
+ double phi = 2. * M_PI * rndmPtr->flat();
+
+ // Evaluate coefficient of azimuthal asymmetry from gluon polarization.
+ findAsymPol( event, dipEndSel);
+ int iFinPol = dipEndSel->iFinPol;
+ double cPol = dipEndSel->asymPol;
+ double phiPol = (iFinPol == 0) ? 0. : event[iFinPol].phi();
+
+ // If interference: try to match sister (anti)colour to final state.
+ int iFinInt = 0;
+ double cInt = 0.;
+ double phiInt = 0.;
+ if (doPhiIntAsym) {
+ for (int i = 0; i < partonSystemsPtr->sizeOut(iSysSel); ++ i) {
+ int iOut = partonSystemsPtr->getOut(iSysSel, i);
+ if ( (acolSister != 0 && event[iOut].col() == acolSister)
+ || (colSister != 0 && event[iOut].acol() == colSister) )
+ iFinInt = iOut;
+ }
+ if (iFinInt != 0) {
+ // Boost final-state parton to current frame of new kinematics.
+ Vec4 pFinTmp = event[iFinInt].p();
+ pFinTmp.rotbst(Mtot);
+ double theFin = pFinTmp.theta();
+ if (side == 2) theFin = M_PI - theFin;
+ double theSis = pSister.theta();
+ if (side == 2) theSis = M_PI - theSis;
+ cInt = strengthIntAsym * 2. * theSis * theFin
+ / (pow2(theSis) + pow2(theFin));
+ phiInt = event[iFinInt].phi();
+ }
+ }
+
+ // Bias phi distribution for polarization and interference.
+ if (iFinPol != 0 || iFinInt != 0) {
+ double cPhiPol, cPhiInt, weight;
+ do {
+ phi = 2. * M_PI * rndmPtr->flat();
+ weight = 1.;
+ if (iFinPol !=0 ) {
+ cPhiPol = cos(phi - phiPol);
+ weight *= ( 1. + cPol * (2. * pow2(cPhiPol) - 1.) )
+ / ( 1. + abs(cPol) );
+ }
+ if (iFinInt !=0 ) {
+ cPhiInt = cos(phi - phiInt);
+ weight *= (1. - cInt) * (1. - cInt * cPhiInt)
+ / (1. + pow2(cInt) - 2. * cInt * cPhiInt);
+ }
+ } while (weight < rndmPtr->flat());
+ }
+
+ // Include rotation -phi on boost to old rest frame.
+ Mtot.rot(0., -phi);
+
+ // Find boost from old rest frame to event cm frame.
+ RotBstMatrix MfromRest;
+ // The boost to the new rest frame.
+ Vec4 sumNew = pMother + pNewRec;
+ double betaX = sumNew.px() / sumNew.e();
+ double betaZ = sumNew.pz() / sumNew.e();
+ MfromRest.bst( -betaX, 0., -betaZ);
+ // Alignment of radiator + recoiler to +- z axis, and rotation +phi.
+ // Note: with spacelike (E < 0) recoiler p'_x_mother < 0 can happen!
+ pMother.rotbst(MfromRest);
+ double theta = pMother.theta();
+ if (pMother.px() < 0.) theta = -theta;
+ if (side == 2) theta += M_PI;
+ MfromRest.rot(-theta, phi);
+ // Boost to radiator + recoiler in event cm frame.
+ if (normalRecoil) {
+ MfromRest.bst( 0., 0., (x1New - x2New) / (x1New + x2New) );
+ } else if (side == 1) {
+ Vec4 pMotherWanted( 0., 0., 0.5 * eCM * x1New, 0.5 * eCM * x1New);
+ MfromRest.fromCMframe( pMotherWanted, event[iRecoiler].p() );
+
+ } else {
+ Vec4 pMotherWanted( 0., 0., -0.5 * eCM * x2New, 0.5 * eCM * x2New);
+ MfromRest.fromCMframe( event[iRecoiler].p(), pMotherWanted );
+ }
+ Mtot.rotbst(MfromRest);
+
+ // Perform cumulative rotation/boost operation.
+ // Mother, recoiler and sister from old rest frame to event cm frame.
+ mother.rotbst(MfromRest);
+ newRecoiler.rotbst(MfromRest);
+ sister.rotbst(MfromRest);
+ // The rest from (and to) event cm frame.
+ for ( int i = eventSizeOld + 2; i < eventSizeOld + systemSizeOld; ++i)
+ event[i].rotbst(Mtot);
+
+ // Allow veto of branching. If so restore event record to before emission.
+ if ( canVetoEmission
+ && userHooksPtr->doVetoISREmission(eventSizeOld, event, iSysSel) ) {
+ event.popBack( event.size() - eventSizeOld);
+ event[beamOff1].daughter1( ev1Dau1V);
+ event[beamOff2].daughter1( ev2Dau1V);
+ for ( int iCopy = 0; iCopy < systemSizeOld; ++iCopy) {
+ int iOldCopy = partonSystemsPtr->getAll(iSysSel, iCopy);
+ event[iOldCopy].status( statusV[iCopy]);
+ event[iOldCopy].mothers( mother1V[iCopy], mother2V[iCopy]);
+ event[iOldCopy].daughters( daughter1V[iCopy], daughter2V[iCopy]);
+ }
+ return false;
+ }
+
+ // Update list of partons in system; adding newly produced one.
+ partonSystemsPtr->setInA(iSysSel, eventSizeOld);
+ partonSystemsPtr->setInB(iSysSel, eventSizeOld + 1);
+ for (int iCopy = 2; iCopy < systemSizeOld; ++iCopy)
+ partonSystemsPtr->setOut(iSysSel, iCopy - 2, eventSizeOld + iCopy);
+ partonSystemsPtr->addOut(iSysSel, eventSizeOld + systemSizeOld);
+ partonSystemsPtr->setSHat(iSysSel, m2 / z);
+
+ // Update info on radiating dipole ends (QCD or QED).
+ for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip)
+ if ( dipEnd[iDip].system == iSysSel) {
+ if (abs(dipEnd[iDip].side) == side) {
+ dipEnd[iDip].iRadiator = iMother;
+ dipEnd[iDip].iRecoiler = iNewRecoiler;
+ if (dipEnd[iDip].side > 0)
+ dipEnd[iDip].colType = mother.colType();
+ else {
+ dipEnd[iDip].chgType = 0;
+ if ( (mother.isQuark() && doQEDshowerByQ)
+ || (mother.isLepton() && doQEDshowerByL) )
+ dipEnd[iDip].chgType = mother.chargeType();
+ }
+ // Kill ME corrections after first emission.
+ dipEnd[iDip].MEtype = 0;
+
+ // Update info on recoiling dipole ends (QCD or QED).
+ } else {
+ dipEnd[iDip].iRadiator = iNewRecoiler;
+ dipEnd[iDip].iRecoiler = iMother;
+ // Optionally also kill recoiler ME corrections after first emission.
+ if (!doMEafterFirst) dipEnd[iDip].MEtype = 0;
+ }
+ }
+
+ // Update info on beam remnants.
+ BeamParticle& beamNow = (side == 1) ? *beamAPtr : *beamBPtr;
+ double xNew = (side == 1) ? x1New : x2New;
+ beamNow[iSysSel].update( iMother, idMother, xNew);
+ // Redo choice of companion kind whenever new flavour.
+ if (idMother != idDaughterNow) {
+ beamNow.xfISR( iSysSel, idMother, xNew, factorMultFac * pT2);
+ beamNow.pickValSeaComp();
+ }
+ BeamParticle& beamRec = (side == 1) ? *beamBPtr : *beamAPtr;
+ beamRec[iSysSel].iPos( iNewRecoiler);
+
+ // Store branching values of current dipole. (For rapidity ordering.)
+ ++dipEndSel->nBranch;
+ dipEndSel->pT2Old = pT2;
+ dipEndSel->zOld = z;
+
+ // Update history if recoiler rescatters.
+ if (!normalRecoil)
+ event[iRecoilMother].daughters( iNewRecoiler, iNewRecoiler);
+
+ // Start list of rescatterers that force changed kinematics.
+ vector<int> iRescatterer;
+ for ( int i = 0; i < systemSizeOld - 2; ++i) {
+ int iOutNew = partonSystemsPtr->getOut( iSysSel, i);
+ if (!event[iOutNew].isFinal()) iRescatterer.push_back(iOutNew);
+ }
+
+ // Start iterate over list of such rescatterers.
+ int iRescNow = -1;
+ while (++iRescNow < int(iRescatterer.size())) {
+
+ // Identify partons that induce or are affected by rescatter shift.
+ // In following Old is before change of kinematics, New after,
+ // Out scatterer in outstate and In in instate of another system.
+ // Daughter sequence is (iOutOld ->) iOutNew -> iInNew -> iInOld.
+ int iOutNew = iRescatterer[iRescNow];
+ int iInOld = event[iOutNew].daughter1();
+ int iSysResc = partonSystemsPtr->getSystemOf(iInOld, true);
+
+ // Copy incoming partons of rescattered system and hook them up.
+ int iOldA = partonSystemsPtr->getInA(iSysResc);
+ int iOldB = partonSystemsPtr->getInB(iSysResc);
+ bool rescSideA = event[iOldA].isRescatteredIncoming();
+ int statusNewA = (rescSideA) ? -45 : -42;
+ int statusNewB = (rescSideA) ? -42 : -45;
+ int iNewA = event.copy(iOldA, statusNewA);
+ int iNewB = event.copy(iOldB, statusNewB);
+
+ // Copy outgoing partons of rescattered system and hook them up.
+ int eventSize = event.size();
+ int sizeOutAB = partonSystemsPtr->sizeOut(iSysResc);
+ int iOldAB, statusOldAB, iNewAB;
+ for (int iOutAB = 0; iOutAB < sizeOutAB; ++iOutAB) {
+ iOldAB = partonSystemsPtr->getOut(iSysResc, iOutAB);
+ statusOldAB = event[iOldAB].status();
+ iNewAB = event.copy(iOldAB, 44);
+ // Status could be negative for parton that rescatters in its turn.
+ if (statusOldAB < 0) {
+ event[iNewAB].statusNeg();
+ iRescatterer.push_back(iNewAB);
+ }
+ }
+
+ // Hook up new outgoing with new incoming parton.
+ int iInNew = (rescSideA) ? iNewA : iNewB;
+ event[iOutNew].daughters( iInNew, iInNew);
+ event[iInNew].mothers( iOutNew, iOutNew);
+
+ // Rescale recoiling incoming parton for correct invariant mass.
+ event[iInNew].p( event[iOutNew].p() );
+ double momFac = (rescSideA)
+ ? event[iInOld].pPos() / event[iInNew].pPos()
+ : event[iInOld].pNeg() / event[iInNew].pNeg();
+ int iInRec = (rescSideA) ? iNewB : iNewA;
+
+ // Rescatter: A previous boost may cause the light cone momentum of a
+ // rescattered parton to change sign. If this happens, tell
+ // parton level to try again.
+ if (momFac < 0.0) {
+ infoPtr->errorMsg("Warning in SpaceShower::branch: "
+ "change in lightcone momentum sign; retrying parton level");
+ rescatterFail = true;
+ return false;
+ }
+
+ event[iInRec].rescale4( momFac);
+
+ // Boost outgoing partons to new frame of incoming.
+ RotBstMatrix MmodResc;
+ MmodResc.toCMframe( event[iOldA].p(), event[iOldB].p());
+ MmodResc.fromCMframe(event[iNewA].p(), event[iNewB].p());
+ for (int iOutAB = 0; iOutAB < sizeOutAB; ++iOutAB)
+ event[eventSize + iOutAB].rotbst(MmodResc);
+
+ // Update list of partons in system.
+ partonSystemsPtr->setInA(iSysResc, iNewA);
+ partonSystemsPtr->setInB(iSysResc, iNewB);
+ for (int iCopy = 0; iCopy < sizeOutAB; ++iCopy)
+ partonSystemsPtr->setOut(iSysResc, iCopy, eventSize + iCopy);
+
+ // Update info on radiating dipole ends (QCD or QED).
+ for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip)
+ if ( dipEnd[iDip].system == iSysResc) {
+ bool sideAnow = (abs(dipEnd[iDip].side) == 1);
+ dipEnd[iDip].iRadiator = (sideAnow) ? iNewA : iNewB;
+ dipEnd[iDip].iRecoiler = (sideAnow) ? iNewB : iNewA;
+ }
+
+ // Update info on beam remnants.
+ BeamParticle& beamResc = (rescSideA) ? *beamAPtr : *beamBPtr;
+ beamResc[iSysResc].iPos( iInNew);
+ beamResc[iSysResc].p( event[iInNew].p() );
+ beamResc[iSysResc].scaleX( 1. / momFac );
+ BeamParticle& beamReco = (rescSideA) ? *beamBPtr : *beamAPtr;
+ beamReco[iSysResc].iPos( iInRec);
+ beamReco[iSysResc].scaleX( momFac);
+
+ // End iterate over list of rescatterers.
+ }
+
+ // Check that beam momentum not used up by rescattered-system boosts.
+ if (beamAPtr->xMax(-1) < 0.0 || beamBPtr->xMax(-1) < 0.0) {
+ infoPtr->errorMsg("Warning in SpaceShower::branch: "
+ "used up beam momentum; retrying parton level");
+ rescatterFail = true;
+ return false;
+ }
+
+ // Done without any errors.
+ return true;
+
+}
+
+//--------------------------------------------------------------------------
+
+// Find class of ME correction.
+
+int SpaceShower::findMEtype( int iSys, Event& event) {
+
+ // Default values and no action.
+ int MEtype = 0;
+ if (!doMEcorrections) ;
+
+ // Identify systems producing a single resonance.
+ else if (partonSystemsPtr->sizeOut( iSys) == 1) {
+ int idIn1 = event[partonSystemsPtr->getInA(iSys)].id();
+ int idIn2 = event[partonSystemsPtr->getInA(iSys)].id();
+ int idRes = event[partonSystemsPtr->getOut(iSys, 0)].id();
+ if (iSys == 0) idResFirst = abs(idRes);
+ if (iSys == 1) idResSecond = abs(idRes);
+
+ // f + fbar -> vector boson.
+ if ( (idRes == 23 || abs(idRes) == 24 || idRes == 32
+ || idRes == 33 || abs(idRes) == 34 || abs(idRes) == 41)
+ && abs(idIn1) < 20 && abs(idIn2) < 20 ) MEtype = 1;
+
+ // g + g, gamma + gamma -> Higgs boson.
+ if ( (idRes == 25 || idRes == 35 || idRes == 36)
+ && ( ( idIn1 == 21 && idIn2 == 21 )
+ || ( idIn1 == 22 && idIn2 == 22 ) ) ) MEtype = 2;
+
+ // f + fbar -> Higgs boson.
+ if ( (idRes == 25 || idRes == 35 || idRes == 36)
+ && abs(idIn1) < 20 && abs(idIn2) < 20 ) MEtype = 3;
+ }
+
+ // Done.
+ return MEtype;
+
+}
+
+//--------------------------------------------------------------------------
+
+// Provide maximum of expected ME weight; for preweighting of evolution.
+
+double SpaceShower::calcMEmax( int MEtype, int idMother, int idDaughterIn) {
+
+ // Currently only one non-unity case: g(gamma) f -> V f'.
+ if (MEtype == 1 && idMother > 20 && idDaughterIn < 20) return 3.;
+ return 1.;
+
+}
+
+//--------------------------------------------------------------------------
+
+// Provide actual ME weight for current branching.
+// Note: currently ME corrections are only allowed for first branching
+// on each side, so idDaughter is essentially known and checks overkill.
+
+double SpaceShower::calcMEcorr(int MEtype, int idMother, int idDaughterIn,
+ double M2, double z, double Q2) {
+
+ // Convert to Mandelstam variables. Sometimes may need to swap later.
+ double sH = M2 / z;
+ double tH = -Q2;
+ double uH = Q2 - M2 * (1. - z) / z;
+ int idMabs = abs(idMother);
+ int idDabs = abs(idDaughterIn);
+
+ // Corrections for f + fbar -> s-channel vector boson.
+ if (MEtype == 1) {
+ if (idMabs < 20 && idDabs < 20) {
+ return (tH*tH + uH*uH + 2. * M2 * sH) / (sH*sH + M2*M2);
+ } else if (idDabs < 20) {
+ // g(gamma) f -> V f': -Q2 = (p_g - p_f')^2 in PS while
+ // tHat = (p_f - p_f')^2 in ME so need to swap tHat <-> uHat.
+ swap( tH, uH);
+ return (sH*sH + uH*uH + 2. * M2 * tH) / (pow2(sH - M2) + M2*M2);
+ }
+
+ // Corrections for g + g -> Higgs boson.
+ } else if (MEtype == 2) {
+ if (idMabs < 20 && idDabs > 20) {
+ return (sH*sH + uH*uH) / (sH*sH + pow2(sH - M2));
+ } else if (idDabs > 20) {
+ return 0.5 * (pow4(sH) + pow4(tH) + pow4(uH) + pow4(M2))
+ / pow2(sH*sH - M2 * (sH - M2));
+ }
+
+ // Corrections for f + fbar -> Higgs boson (f = b mainly).
+ } else if (MEtype == 3) {
+ if (idMabs < 20 && idDabs < 20) {
+ // The PS and ME answers agree for f fbar -> H g/gamma.
+ return 1.;
+ } else if (idDabs < 20) {
+ // Need to swap tHat <-> uHat, cf. vector-boson production above.
+ swap( tH, uH);
+ return (sH*sH + uH*uH + 2. * (M2 - uH) * (M2 - sH))
+ / (pow2(sH - M2) + M2*M2);
+ }
+ }
+
+ return 1.;
+
+}
+
+//--------------------------------------------------------------------------
+
+// Find coefficient of azimuthal asymmetry from gluon polarization.
+
+void SpaceShower::findAsymPol( Event& event, SpaceDipoleEnd* dip) {
+
+ // Default is no asymmetry. Only gluons are studied.
+ dip->iFinPol = 0;
+ dip->asymPol = 0.;
+ int iRad = dip->iRadiator;
+ if (!doPhiPolAsym || dip->idDaughter != 21) return;
+
+ // At least two particles in final state, whereof at least one coloured.
+ int systemSizeOut = partonSystemsPtr->sizeOut( iSysSel);
+ if (systemSizeOut < 2) return;
+ bool foundColOut = false;
+ for (int ii = 0; ii < systemSizeOut; ++ii) {
+ int i = partonSystemsPtr->getOut( iSysSel, ii);
+ if (event[i].col() != 0 || event[i].acol() != 0) foundColOut = true;
+ }
+ if (!foundColOut) return;
+
+ // Check if granddaughter in final state of hard scattering.
+ // (May need to trace across carbon copies to find granddaughters.)
+ // If so, only accept 2 -> 2 scatterings with gg or qq in final state.
+ int iGrandD1 = event[iRad].daughter1();
+ int iGrandD2 = event[iRad].daughter2();
+ bool traceCopy = false;
+ do {
+ traceCopy = false;
+ if (iGrandD1 > 0 && iGrandD2 == iGrandD1) {
+ iGrandD1 = event[iGrandD2].daughter1();
+ iGrandD2 = event[iGrandD2].daughter2();
+ traceCopy = true;
+ }
+ } while (traceCopy);
+ int statusGrandD1 = event[ iGrandD1 ].statusAbs();
+ bool isHardProc = (statusGrandD1 == 23 || statusGrandD1 == 33);
+ if (isHardProc) {
+ if (iGrandD2 != iGrandD1 + 1) return;
+ if (event[iGrandD1].isGluon() && event[iGrandD2].isGluon());
+ else if (event[iGrandD1].isQuark() && event[iGrandD2].isQuark());
+ else return;
+ }
+ dip->iFinPol = iGrandD1;
+
+ // Coefficient from gluon production.
+ if (dip->idMother == 21) dip->asymPol = pow2( (1. - dip->z)
+ / (1. - dip->z * (1. - dip->z) ) );
+ else dip->asymPol = 2. * (1. - dip->z) / (1. + pow2(1. - dip->z) );
+
+ // Coefficients from gluon decay. Put z = 1/2 for hard process.
+ double zDau = (isHardProc) ? 0.5 : dip->zOld;
+ if (event[iGrandD1].isGluon()) dip->asymPol *= pow2( (1. - zDau)
+ / (1. - zDau * (1. - zDau) ) );
+ else dip->asymPol *= -2. * zDau *( 1. - zDau )
+ / (1. - 2. * zDau * (1. - zDau) );
+
+}
+
+//--------------------------------------------------------------------------
+
+// Print the list of dipoles.
+
+void SpaceShower::list(ostream& os) const {
+
+ // Header.
+ os << "\n -------- PYTHIA SpaceShower Dipole Listing -------------- \n"
+ << "\n i syst side rad rec pTmax col chg ME rec \n"
+ << fixed << setprecision(3);
+
+ // Loop over dipole list and print it.
+ for (int i = 0; i < int(dipEnd.size()); ++i)
+ os << setw(5) << i << setw(6) << dipEnd[i].system
+ << setw(6) << dipEnd[i].side << setw(6) << dipEnd[i].iRadiator
+ << setw(6) << dipEnd[i].iRecoiler << setw(12) << dipEnd[i].pTmax
+ << setw(5) << dipEnd[i].colType << setw(5) << dipEnd[i].chgType
+ << setw(5) << dipEnd[i].MEtype << setw(4)
+ << dipEnd[i].normalRecoil << "\n";
+
+ // Done.
+ os << "\n -------- End PYTHIA SpaceShower Dipole Listing ----------"
+ << endl;
+
+
+}
+
+//==========================================================================
+
+} // end namespace Pythia8
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff