// TimeShower.cc is a part of the PYTHIA event generator. // Copyright (C) 2010 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 TimeShower class. #include "TimeShower.h" namespace Pythia8 { //========================================================================== // The TimeShower class. //-------------------------------------------------------------------------- // Constants: could be changed here if desired, but normally should not. // These are of technical nature, as described for each. // For small x approximate 1 - sqrt(1 - x) by x/2. const double TimeShower::SIMPLIFYROOT = 1e-8; // Do not allow x too close to 0 or 1 in matrix element expressions. // Warning: cuts into phase space for E_CM > 2 * pTmin * sqrt(1/XMARGIN), // i.e. will become problem roughly for E_CM > 10^6 GeV. const double TimeShower::XMARGIN = 1e-12; const double TimeShower::XMARGINCOMB = 1e-4; // Lower limit on PDF value in order to avoid division by zero. const double TimeShower::TINYPDF = 1e-10; // Big starting value in search for smallest invariant-mass pair. const double TimeShower::LARGEM2 = 1e20; // In g -> q qbar or gamma -> f fbar require m2_pair > this * m2_q/f. const double TimeShower::THRESHM2 = 4.004; // Never pick pT so low that alphaS is evaluated too close to Lambda_3. const double TimeShower::LAMBDA3MARGIN = 1.1; // Rescatter: rescattering + ISR + FSR + primordial kT can lead to // systems not locally conserving momentum. // Fix up momentum in intermediate systems with rescattering const bool TimeShower::FIXRESCATTER = true; // Veto negative energies when using FIXRESCATTER option. const bool TimeShower::VETONEGENERGY = false; // Do not allow too large time- or spacelike virtualities in fixing-up. const double TimeShower::MAXVIRTUALITYFRACTION = 0.5; // Do not allow too large negative spacelike energy in system rest frame. const double TimeShower::MAXNEGENERGYFRACTION = 0.7; //-------------------------------------------------------------------------- // Initialize alphaStrong, alphaEM and related pTmin parameters. void TimeShower::init( BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn) { // Store input pointers for future use. beamAPtr = beamAPtrIn; beamBPtr = beamBPtrIn; // Main flags. doQCDshower = settingsPtr->flag("TimeShower:QCDshower"); doQEDshowerByQ = settingsPtr->flag("TimeShower:QEDshowerByQ"); doQEDshowerByL = settingsPtr->flag("TimeShower:QEDshowerByL"); doQEDshowerByGamma = settingsPtr->flag("TimeShower:QEDshowerByGamma"); doMEcorrections = settingsPtr->flag("TimeShower:MEcorrections"); doPhiPolAsym = settingsPtr->flag("TimeShower:phiPolAsym"); doInterleave = settingsPtr->flag("TimeShower:interleave"); allowBeamRecoil = settingsPtr->flag("TimeShower:allowBeamRecoil"); dampenBeamRecoil = settingsPtr->flag("TimeShower:dampenBeamRecoil"); // Matching in pT of hard interaction or MI to shower evolution. pTmaxMatch = settingsPtr->mode("TimeShower:pTmaxMatch"); pTdampMatch = settingsPtr->mode("TimeShower:pTdampMatch"); pTmaxFudge = settingsPtr->parm("TimeShower:pTmaxFudge"); pTmaxFudgeMI = settingsPtr->parm("TimeShower:pTmaxFudgeMI"); pTdampFudge = settingsPtr->parm("TimeShower:pTdampFudge"); // Charm and bottom mass thresholds. mc = particleDataPtr->m0(4); mb = particleDataPtr->m0(5); m2c = mc * mc; m2b = mb * mb; // Parameters of alphaStrong generation . alphaSvalue = settingsPtr->parm("TimeShower:alphaSvalue"); alphaSorder = settingsPtr->mode("TimeShower:alphaSorder"); alphaS2pi = 0.5 * alphaSvalue / M_PI; // Initialize alphaStrong generation. alphaS.init( alphaSvalue, alphaSorder); // Lambda for 5, 4 and 3 flavours. Lambda3flav = alphaS.Lambda3(); Lambda4flav = alphaS.Lambda4(); Lambda5flav = alphaS.Lambda5(); Lambda5flav2 = pow2(Lambda5flav); Lambda4flav2 = pow2(Lambda4flav); Lambda3flav2 = pow2(Lambda3flav); // Parameters of QCD evolution. Warn if pTmin must be raised. nGluonToQuark = settingsPtr->mode("TimeShower:nGluonToQuark"); pTcolCutMin = settingsPtr->parm("TimeShower:pTmin"); if (pTcolCutMin > LAMBDA3MARGIN * Lambda3flav) pTcolCut = pTcolCutMin; else { pTcolCut = LAMBDA3MARGIN * Lambda3flav; ostringstream newPTcolCut; newPTcolCut << fixed << setprecision(3) << pTcolCut; infoPtr->errorMsg("Warning in TimeShower::init: pTmin too low", ", raised to " + newPTcolCut.str() ); infoPtr->setTooLowPTmin(true); } pT2colCut = pow2(pTcolCut); // Parameters of alphaEM generation . alphaEMorder = settingsPtr->mode("TimeShower:alphaEMorder"); // Initialize alphaEM generation. alphaEM.init( alphaEMorder, settingsPtr); // Parameters of QED evolution. nGammaToQuark = settingsPtr->mode("TimeShower:nGammaToQuark"); nGammaToLepton = settingsPtr->mode("TimeShower:nGammaToLepton"); pTchgQCut = settingsPtr->parm("TimeShower:pTminChgQ"); pT2chgQCut = pow2(pTchgQCut); pTchgLCut = settingsPtr->parm("TimeShower:pTminChgL"); pT2chgLCut = pow2(pTchgLCut); mMaxGamma = settingsPtr->parm("TimeShower:mMaxGamma"); m2MaxGamma = pow2(mMaxGamma); // Consisteny check for gamma -> f fbar variables. if (nGammaToQuark <= 0 && nGammaToLepton <= 0) doQEDshowerByGamma = false; // Fraction and colour factor of gluon emission off onium octat state. octetOniumFraction = settingsPtr->parm("TimeShower:octetOniumFraction"); octetOniumColFac = settingsPtr->parm("TimeShower:octetOniumColFac"); // Z0 properties needed for gamma/Z0 mixing. mZ = particleDataPtr->m0(23); gammaZ = particleDataPtr->mWidth(23); thetaWRat = 1. / (16. * coupSMPtr->sin2thetaW() * coupSMPtr->cos2thetaW()); // May have to fix up recoils related to rescattering. allowRescatter = settingsPtr->flag("PartonLevel:MI") && settingsPtr->flag("MultipleInteractions:allowRescatter"); // Hidden Valley scenario with further shower activity. doHVshower = settingsPtr->flag("HiddenValley:FSR"); nCHV = settingsPtr->mode("HiddenValley:Ngauge"); alphaHVfix = settingsPtr->parm("HiddenValley:alphaFSR"); pThvCut = settingsPtr->parm("HiddenValley:pTminFSR"); pT2hvCut = pThvCut * pThvCut; CFHV = (nCHV == 1) ? 1. : (nCHV * nCHV - 1.)/(2. * nCHV); idHV = (nCHV == 1) ? 4900022 : 4900021; mHV = particleDataPtr->m0(idHV); brokenHVsym = (nCHV == 1 && mHV > 0.); // Possibility to allow user veto of emission step. canVetoEmission = (userHooksPtr > 0) ? userHooksPtr->canVetoFSREmission() : false; } //-------------------------------------------------------------------------- // Find whether to limit maximum scale of emissions. // Also allow for dampening at factorization or renormalization scale. bool TimeShower::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; } //-------------------------------------------------------------------------- // Top-level routine to do a full time-like shower in resonance decay. int TimeShower::shower( int iBeg, int iEnd, Event& event, double pTmax, int nBranchMax) { // Add new system, automatically with two empty beam slots. int iSys = partonSystemsPtr->addSys(); // Loop over allowed range to find all final-state particles. Vec4 pSum; for (int i = iBeg; i <= iEnd; ++i) if (event[i].isFinal()) { partonSystemsPtr->addOut( iSys, i); pSum += event[i].p(); } partonSystemsPtr->setSHat( iSys, pSum.m2Calc() ); // Let prepare routine do the setup. prepare( iSys, event, true); // Begin evolution down in pT from hard pT scale. int nBranch = 0; pTLastBranch = 0.; do { double pTtimes = pTnext( event, pTmax, 0.); // Do a final-state emission (if allowed). if (pTtimes > 0.) { if (branch( event)) { ++nBranch; pTLastBranch = pTtimes; } pTmax = pTtimes; } // Keep on evolving until nothing is left to be done. else pTmax = 0.; } while (pTmax > 0. && (nBranchMax <= 0 || nBranch < nBranchMax)); // Return number of emissions that were performed. return nBranch; } //-------------------------------------------------------------------------- // Prepare system for evolution; identify ME. void TimeShower::prepare( int iSys, Event& event, bool limitPTmaxIn) { // Reset dipole-ends list for first interaction and for resonance decays. int iInA = partonSystemsPtr->getInA(iSys); int iInB = partonSystemsPtr->getInB(iSys); if (iSys == 0 || iInA == 0) dipEnd.resize(0); int dipEndSizeBeg = dipEnd.size(); // Loop through final state of system to find possible dipole ends. for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) { int iRad = partonSystemsPtr->getOut( iSys, i); if (event[iRad].isFinal() && event[iRad].scale() > 0.) { // Identify colour octet onium state. Check whether QCD shower allowed. int idRad = event[iRad].id(); int idRadAbs = abs(idRad); bool isOctetOnium = ( idRad == 9900441 || idRad == 9900443 || idRad == 9910441 || idRad == 9900551 || idRad == 9900553 || idRad == 9910551 ); bool doQCD = doQCDshower; if (doQCD && isOctetOnium) doQCD = (rndmPtr->flat() < octetOniumFraction); // Find dipole end formed by colour index. int colTag = event[iRad].col(); if (doQCD && colTag > 0) setupQCDdip( iSys, i, colTag, 1, event, isOctetOnium, limitPTmaxIn); // Find dipole end formed by anticolour index. int acolTag = event[iRad].acol(); if (doQCD && acolTag > 0) setupQCDdip( iSys, i, acolTag, -1, event, isOctetOnium, limitPTmaxIn); // Find "charge-dipole" and "photon-dipole" ends. int chgType = event[iRad].chargeType(); bool doChgDip = (chgType != 0) && ( ( doQEDshowerByQ && event[iRad].isQuark() ) || ( doQEDshowerByL && event[iRad].isLepton() ) ); int gamType = (idRad == 22) ? 1 : 0; bool doGamDip = (gamType == 1) && doQEDshowerByGamma; if (doChgDip || doGamDip) setupQEDdip( iSys, i, chgType, gamType, event, limitPTmaxIn); // Find Hidden Valley dipole ends. bool isHVrad = (idRadAbs > 4900000 && idRadAbs < 4900007) || (idRadAbs > 4900010 && idRadAbs < 4900017) || idRadAbs == 4900101; if (doHVshower && isHVrad) setupHVdip( iSys, i, event, limitPTmaxIn); // End loop over system final state. Have now found the dipole ends. } } // Loop through dipole ends to find matrix element corrections. for (int iDip = dipEndSizeBeg; iDip < int(dipEnd.size()); ++iDip) findMEtype( event, dipEnd[iDip]); // Update dipole list after a multiple interactions rescattering. if (iSys > 0 && ( (iInA > 0 && event[iInA].status() == -34) || (iInB > 0 && event[iInB].status() == -34) ) ) rescatterUpdate( iSys, event); } //-------------------------------------------------------------------------- // Update dipole list after a multiple interactions rescattering. void TimeShower::rescatterUpdate( int iSys, Event& event) { // Loop over two incoming partons in system; find their rescattering mother. // (iOut is outgoing from old system, = incoming iIn of rescattering system.) for (int iResc = 0; iResc < 2; ++iResc) { int iIn = (iResc == 0) ? partonSystemsPtr->getInA(iSys) : partonSystemsPtr->getInB(iSys); if (iIn == 0 || event[iIn].status() != -34) continue; int iOut = event[iIn].mother1(); // Loop over all dipoles. int dipEndSize = dipEnd.size(); for (int iDip = 0; iDip < dipEndSize; ++iDip) { TimeDipoleEnd& dipNow = dipEnd[iDip]; // Kill dipoles where rescattered parton is radiator. if (dipNow.iRadiator == iOut) { dipNow.colType = 0; dipNow.chgType = 0; dipNow.gamType = 0; continue; } // No matrix element for dipoles between scatterings. if (dipNow.iMEpartner == iOut) { dipNow.MEtype = 0; dipNow.iMEpartner = -1; } // Update dipoles where outgoing rescattered parton is recoiler. if (dipNow.iRecoiler == iOut) { int iRad = dipNow.iRadiator; // Colour dipole: recoil in final state, initial state or new if (dipNow.colType > 0) { int colTag = event[iRad].col(); bool done = false; for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) { int iRecNow = partonSystemsPtr->getOut( iSys, i); if (event[iRecNow].acol() == colTag) { dipNow.iRecoiler = iRecNow; dipNow.systemRec = iSys; dipNow.MEtype = 0; done = true; break; } } if (!done) { int iIn2 = (iResc == 0) ? partonSystemsPtr->getInB(iSys) : partonSystemsPtr->getInA(iSys); if (event[iIn2].col() == colTag) { dipNow.iRecoiler = iIn2; dipNow.systemRec = iSys; dipNow.MEtype = 0; int isrType = event[iIn2].mother1(); // This line in case mother is a rescattered parton. while (isrType > 2 + beamOffset) isrType = event[isrType].mother1(); if (isrType > 2) isrType -= beamOffset; dipNow.isrType = isrType; done = true; } } // If above options failed, then create new dipole if (!done) { int iRadNow = partonSystemsPtr->getIndexOfOut(dipNow.system, iRad); if (iRadNow != -1) setupQCDdip(dipNow.system, iRadNow, event[iRad].col(), 1, event, dipNow.isOctetOnium, true); else infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: " "failed to locate radiator in system"); dipNow.colType = 0; dipNow.chgType = 0; dipNow.gamType = 0; infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: " "failed to locate new recoiling colour partner"); } // Anticolour dipole: recoil in final state, initial state or new } else if (dipNow.colType < 0) { int acolTag = event[iRad].acol(); bool done = false; for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) { int iRecNow = partonSystemsPtr->getOut( iSys, i); if (event[iRecNow].col() == acolTag) { dipNow.iRecoiler = iRecNow; dipNow.systemRec = iSys; dipNow.MEtype = 0; done = true; break; } } if (!done) { int iIn2 = (iResc == 0) ? partonSystemsPtr->getInB(iSys) : partonSystemsPtr->getInA(iSys); if (event[iIn2].acol() == acolTag) { dipNow.iRecoiler = iIn2; dipNow.systemRec = iSys; dipNow.MEtype = 0; int isrType = event[iIn2].mother1(); // This line in case mother is a rescattered parton. while (isrType > 2 + beamOffset) isrType = event[isrType].mother1(); if (isrType > 2) isrType -= beamOffset; dipNow.isrType = isrType; done = true; } } // If above options failed, then create new dipole if (!done) { int iRadNow = partonSystemsPtr->getIndexOfOut(dipNow.system, iRad); if (iRadNow != -1) setupQCDdip(dipNow.system, iRadNow, event[iRad].acol(), -1, event, dipNow.isOctetOnium, true); else infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: " "failed to locate radiator in system"); dipNow.colType = 0; dipNow.chgType = 0; dipNow.gamType = 0; infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: " "failed to locate new recoiling colour partner"); } // Charge or photon dipoles: same flavour in final or initial state. } else if (dipNow.chgType != 0 || dipNow.gamType != 0) { int idTag = event[dipNow.iRecoiler].id(); bool done = false; for (int i = 0; i < partonSystemsPtr->sizeOut(iSys); ++i) { int iRecNow = partonSystemsPtr->getOut( iSys, i); if (event[iRecNow].id() == idTag) { dipNow.iRecoiler = iRecNow; dipNow.systemRec = iSys; dipNow.MEtype = 0; done = true; break; } } if (!done) { int iIn2 = (iResc == 0) ? partonSystemsPtr->getInB(iSys) : partonSystemsPtr->getInA(iSys); if (event[iIn2].id() == -idTag) { dipNow.iRecoiler = iIn2; dipNow.systemRec = iSys; dipNow.MEtype = 0; int isrType = event[iIn2].mother1(); // This line in case mother is a rescattered parton. while (isrType > 2 + beamOffset) isrType = event[isrType].mother1(); if (isrType > 2) isrType -= beamOffset; dipNow.isrType = isrType; done = true; } } // If above options failed, then create new dipole if (!done) { int iRadNow = partonSystemsPtr->getIndexOfOut(dipNow.system, iRad); if (iRadNow != -1) setupQEDdip(dipNow.system, iRadNow, dipNow.chgType, dipNow.gamType, event, true); else infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: " "failed to locate radiator in system"); dipNow.colType = 0; dipNow.chgType = 0; dipNow.gamType = 0; infoPtr->errorMsg("Warning in TimeShower::rescatterUpdate: " "failed to locate new recoiling charge partner"); } } } // End of loop over dipoles and two incoming sides. } } } //-------------------------------------------------------------------------- // Update dipole list after each ISR emission (so not used for resonances). void TimeShower::update( int iSys, Event& event) { // Start list of rescatterers that gave further changed systems in ISR. vector iRescatterer; // Find new and old positions of incoming partons in the system. vector iNew, iOld; iNew.push_back( partonSystemsPtr->getInA(iSys) ); iOld.push_back( event[iNew[0]].daughter2() ); iNew.push_back( partonSystemsPtr->getInB(iSys) ); iOld.push_back( event[iNew[1]].daughter2() ); // Ditto for outgoing partons, except the newly created one. int sizeOut = partonSystemsPtr->sizeOut(iSys) - 1; for (int i = 0; i < sizeOut; ++i) { int iNow = partonSystemsPtr->getOut(iSys, i); iNew.push_back( iNow ); iOld.push_back( event[iNow].mother1() ); // Add non-final to list of rescatterers. if (!event[iNow].isFinal()) iRescatterer.push_back( iNow ); } int iNewNew = partonSystemsPtr->getOut(iSys, sizeOut); // Swap beams to let 0 be side on which branching occured. if (event[iNew[0]].status() != -41) { swap( iNew[0], iNew[1]); swap( iOld[0], iOld[1]); } // Loop over all dipole ends belonging to the system // or to the recoil system, if different. for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip) if (dipEnd[iDip].system == iSys || dipEnd[iDip].systemRec == iSys) { TimeDipoleEnd& dipNow = dipEnd[iDip]; // Replace radiator (always in final state so simple). for (int i = 2; i < 2 + sizeOut; ++i) if (dipNow.iRadiator == iOld[i]) { dipNow.iRadiator = iNew[i]; break; } // Replace ME partner (always in final state, if exists, so simple). for (int i = 2; i < 2 + sizeOut; ++i) if (dipNow.iMEpartner == iOld[i]) { dipNow.iMEpartner = iNew[i]; break; } // Recoiler: by default pick old one, only moved. Note excluded beam. int iRec = 0; if (dipNow.systemRec == iSys) { for (int i = 1; i < 2 + sizeOut; ++i) if (dipNow.iRecoiler == iOld[i]) { iRec = iNew[i]; break; } // QCD recoiler: check if colour hooks up with new final parton. if ( dipNow.colType > 0 && event[dipNow.iRadiator].col() == event[iNewNew].acol() ) { iRec = iNewNew; dipNow.isrType = 0; } if ( dipNow.colType < 0 && event[dipNow.iRadiator].acol() == event[iNewNew].col() ) { iRec = iNewNew; dipNow.isrType = 0; } // QCD recoiler: check if colour hooks up with new beam parton. if ( iRec == 0 && dipNow.colType > 0 && event[dipNow.iRadiator].col() == event[iNew[0]].col() ) iRec = iNew[0]; if ( iRec == 0 && dipNow.colType < 0 && event[dipNow.iRadiator].acol() == event[iNew[0]].acol() ) iRec = iNew[0]; // QED/photon recoiler: either to new particle or remains to beam. if ( iRec == 0 && (dipNow.chgType != 0 || dipNow.gamType != 0) ) { if ( event[iNew[0]].chargeType() == 0 ) { iRec = iNewNew; dipNow.isrType = 0; } else { iRec = iNew[0]; } } // Recoiler in another system: keep it as is. } else iRec = dipNow.iRecoiler; // Done. Kill dipole if failed to find new recoiler. dipNow.iRecoiler = iRec; if ( iRec == 0 && (dipNow.colType != 0 || dipNow.chgType != 0 || dipNow.gamType != 0) ) { dipNow.colType = 0; dipNow.chgType = 0; dipNow.gamType = 0; infoPtr->errorMsg("Error in TimeShower::update: " "failed to locate new recoiling partner"); } } // Find new dipole end formed by colour index. int colTag = event[iNewNew].col(); if (doQCDshower && colTag > 0) setupQCDdip( iSys, sizeOut, colTag, 1, event, false, true); // Find new dipole end formed by anticolour index. int acolTag = event[iNewNew].acol(); if (doQCDshower && acolTag > 0) setupQCDdip( iSys, sizeOut, acolTag, -1, event, false, true); // Find new "charge-dipole" and "photon-dipole" ends. int chgType = event[iNewNew].chargeType(); bool doChgDip = (chgType != 0) && ( ( doQEDshowerByQ && event[iNewNew].isQuark() ) || ( doQEDshowerByL && event[iNewNew].isLepton() ) ); int gamType = (event[iNewNew].id() == 22) ? 1 : 0; bool doGamDip = (gamType == 1) && doQEDshowerByGamma; if (doChgDip || doGamDip) setupQEDdip( iSys, sizeOut, chgType, gamType, event, true); // Start iterate over list of rescatterers - may be empty. int iRescNow = -1; while (++iRescNow < int(iRescatterer.size())) { // Identify systems that rescatterers belong to. int iOutNew = iRescatterer[iRescNow]; int iInNew = event[iOutNew].daughter1(); int iSysResc = partonSystemsPtr->getSystemOf(iInNew, true); // Find new and old positions of incoming partons in the system. iNew.resize(0); iOld.resize(0); iNew.push_back( partonSystemsPtr->getInA(iSysResc) ); iOld.push_back( event[iNew[0]].daughter1() ); iNew.push_back( partonSystemsPtr->getInB(iSysResc) ); iOld.push_back( event[iNew[1]].daughter1() ); // Ditto for outgoing partons. sizeOut = partonSystemsPtr->sizeOut(iSysResc); for (int i = 0; i < sizeOut; ++i) { int iNow = partonSystemsPtr->getOut(iSysResc, i); iNew.push_back( iNow ); iOld.push_back( event[iNow].mother1() ); // Add non-final to list of rescatterers. if (!event[iNow].isFinal()) iRescatterer.push_back( iNow ); } // Loop over all dipole ends belonging to the system // or to the recoil system, if different. for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip) if (dipEnd[iDip].system == iSysResc || dipEnd[iDip].systemRec == iSysResc) { TimeDipoleEnd& dipNow = dipEnd[iDip]; // Replace radiator (always in final state so simple). for (int i = 2; i < 2 + sizeOut; ++i) if (dipNow.iRadiator == iOld[i]) { dipNow.iRadiator = iNew[i]; break; } // Replace ME partner (always in final state, if exists, so simple). for (int i = 2; i < 2 + sizeOut; ++i) if (dipNow.iMEpartner == iOld[i]) { dipNow.iMEpartner = iNew[i]; break; } // Replace recoiler. for (int i = 0; i < 2 + sizeOut; ++i) if (dipNow.iRecoiler == iOld[i]) { dipNow.iRecoiler = iNew[i]; break; } } // End iterate over list of rescatterers. } } //-------------------------------------------------------------------------- // Setup a dipole end for a QCD colour charge. void TimeShower::setupQCDdip( int iSys, int i, int colTag, int colSign, Event& event, bool isOctetOnium, bool limitPTmaxIn) { // Initial values. Find if allowed to hook up beams. int iRad = partonSystemsPtr->getOut(iSys, i); int iRec = 0; int sizeAllA = partonSystemsPtr->sizeAll(iSys); int sizeOut = partonSystemsPtr->sizeOut(iSys); int sizeAll = ( allowBeamRecoil ) ? sizeAllA : sizeOut; int sizeIn = sizeAll - sizeOut; int sizeInA = sizeAllA - sizeIn - sizeOut; int iOffset = i + sizeAllA - sizeOut; bool otherSystemRec = false; bool allowInitial = (partonSystemsPtr->hasInAB(iSys)) ? true : false; // Colour: other end by same index in beam or opposite in final state. // Exclude rescattered incoming and not final outgoing. if (colSign > 0) for (int j = 0; j < sizeAll; ++j) if (j + sizeInA != iOffset) { int iRecNow = partonSystemsPtr->getAll(iSys, j + sizeInA); if ( ( j < sizeIn && event[iRecNow].col() == colTag && !event[iRecNow].isRescatteredIncoming() ) || ( j >= sizeIn && event[iRecNow].acol() == colTag && event[iRecNow].isFinal() ) ) { iRec = iRecNow; break; } } // Anticolour: other end by same index in beam or opposite in final state. // Exclude rescattered incoming and not final outgoing. if (colSign < 0) for (int j = 0; j < sizeAll; ++j) if (j + sizeInA != iOffset) { int iRecNow = partonSystemsPtr->getAll(iSys, j + sizeInA); if ( ( j < sizeIn && event[iRecNow].acol() == colTag && !event[iRecNow].isRescatteredIncoming() ) || ( j >= sizeIn && event[iRecNow].col() == colTag && event[iRecNow].isFinal() ) ) { iRec = iRecNow; break; } } // Resonance decays (= no instate), without a junction colour: // other end to nearest recoiler in same system final state, // by (p_i + p_j)^2 - (m_i + m_j)^2 = 2 (p_i p_j - m_i m_j). if (iRec == 0 && !allowInitial) { bool hasJunction = false; for (int iJun = 0; iJun < event.sizeJunction(); ++ iJun) { // For types 1&2, all legs in final state // For types 3&4, two legs in final state // For types 5&6, one leg in final state int iBeg = (event.kindJunction(iJun)-1)/2; for (int iLeg = iBeg; iLeg < 3; ++iLeg) { if (event.endColJunction( iJun, iLeg) == colTag) hasJunction = true; } } if (!hasJunction) { double ppMin = LARGEM2; for (int j = 0; j < sizeOut; ++j) if (j != i) { int iRecNow = partonSystemsPtr->getOut(iSys, j); if (!event[iRecNow].isFinal()) continue; double ppNow = event[iRecNow].p() * event[iRad].p() - event[iRecNow].m() * event[iRad].m(); if (ppNow < ppMin) { iRec = iRecNow; ppMin = ppNow; } } } } // If no success then look for matching (anti)colour anywhere in final state. if ( iRec == 0 || (!doInterleave && !event[iRec].isFinal()) ) { iRec = 0; for (int j = 0; j < event.size(); ++j) if (event[j].isFinal()) if ( (colSign > 0 && event[j].acol() == colTag) || (colSign < 0 && event[j].col() == colTag) ) { iRec = j; otherSystemRec = true; break; } // If no success then look for match to non-rescattered in initial state. if (iRec == 0 && allowInitial) { for (int iSysR = 0; iSysR < partonSystemsPtr->sizeSys(); ++iSysR) if (iSysR != iSys) { int j = partonSystemsPtr->getInA(iSysR); if (j > 0 && event[j].isRescatteredIncoming()) j = 0; if (j > 0 && ( (colSign > 0 && event[j].col() == colTag) || (colSign < 0 && event[j].acol() == colTag) ) ) { iRec = j; otherSystemRec = true; break; } j = partonSystemsPtr->getInB(iSysR); if (j > 0 && event[j].isRescatteredIncoming()) j = 0; if (j > 0 && ( (colSign > 0 && event[j].col() == colTag) || (colSign < 0 && event[j].acol() == colTag) ) ) { iRec = j; otherSystemRec = true; break; } } } // For junction pick any of other leg as recoiler. Note colour sign. if (iRec == 0) { for (int iJun = 0; iJun < event.sizeJunction(); ++ iJun) for (int iLeg = 0; iLeg < 3; ++iLeg) if (event.endColJunction( iJun, iLeg) == colTag) { int iLegRec = iLeg + 1 + int(2. * rndmPtr->flat()); if (iLegRec >= 3) iLegRec -= 3; int colTagRec = event.endColJunction( iJun, iLegRec); for (int j = 0; j < event.size(); ++j) if (event[j].isFinal()) if ( (colSign > 0 && event[j].col() == colTagRec) || (colSign < 0 && event[j].acol() == colTagRec) ) { iRec = j; otherSystemRec = true; break; } } } } // If fail, then other end to nearest recoiler in same system final state, // by (p_i + p_j)^2 - (m_i + m_j)^2 = 2 (p_i p_j - m_i m_j). if (iRec == 0) { double ppMin = LARGEM2; for (int j = 0; j < sizeOut; ++j) if (j != i) { int iRecNow = partonSystemsPtr->getOut(iSys, j); if (!event[iRecNow].isFinal()) continue; double ppNow = event[iRecNow].p() * event[iRad].p() - event[iRecNow].m() * event[iRad].m(); if (ppNow < ppMin) { iRec = iRecNow; ppMin = ppNow; } } } // If fail, then other end to nearest recoiler in any system final state, // by (p_i + p_j)^2 - (m_i + m_j)^2 = 2 (p_i p_j - m_i m_j). if (iRec == 0) { double ppMin = LARGEM2; for (int iRecNow = 0; iRecNow < event.size(); ++iRecNow) if (iRecNow != iRad && event[iRecNow].isFinal()) { double ppNow = event[iRecNow].p() * event[iRad].p() - event[iRecNow].m() * event[iRad].m(); if (ppNow < ppMin) { iRec = iRecNow; otherSystemRec = true; ppMin = ppNow; } } } // Store dipole colour end. if (iRec > 0) { // Max scale either by parton scale or by half dipole mass. double pTmax = event[iRad].scale(); if (limitPTmaxIn) { if (iSys == 0) pTmax *= pTmaxFudge; if (iSys > 0 && sizeIn > 0) pTmax *= pTmaxFudgeMI; } else pTmax = 0.5 * m( event[iRad], event[iRec]); int colType = (event[iRad].id() == 21) ? 2 * colSign : colSign; int isrType = (event[iRec].isFinal()) ? 0 : event[iRec].mother1(); // This line in case mother is a rescattered parton. while (isrType > 2 + beamOffset) isrType = event[isrType].mother1(); if (isrType > 2) isrType -= beamOffset; dipEnd.push_back( TimeDipoleEnd( iRad, iRec, pTmax, colType, 0, 0, isrType, iSys, -1, -1, isOctetOnium) ); // If hooked up with other system then find which. if (otherSystemRec) { int systemRec = partonSystemsPtr->getSystemOf(iRec, true); if (systemRec >= 0) dipEnd.back().systemRec = systemRec; dipEnd.back().MEtype = 0; } } else { infoPtr->errorMsg("Error in TimeShower::setupQCDdip: " "failed to locate any recoiling partner"); } } //-------------------------------------------------------------------------- // Setup a dipole end for a QED colour charge or a photon. // No failsafe choice of recoiler, so gradually widen search. void TimeShower::setupQEDdip( int iSys, int i, int chgType, int gamType, Event& event, bool limitPTmaxIn) { // Initial values. Find if allowed to hook up beams. int iRad = partonSystemsPtr->getOut(iSys, i); int idRad = event[iRad].id(); int iRec = 0; int sizeAllA = partonSystemsPtr->sizeAll(iSys); int sizeOut = partonSystemsPtr->sizeOut(iSys); int sizeAll = ( allowBeamRecoil ) ? sizeAllA : sizeOut; int sizeIn = sizeAll - sizeOut; int sizeInA = sizeAllA - sizeIn - sizeOut; int iOffset = i + sizeAllA - sizeOut; double ppMin = LARGEM2; bool hasRescattered = false; bool otherSystemRec = false; // Find nearest same- (opposide-) flavour recoiler in initial (final) // state of same system, excluding rescattered (in or out) partons. // Also find if system is involved in rescattering. // Note: (p_i + p_j)2 - (m_i + m_j)2 = 2 (p_i p_j - m_i m_j). for (int j = 0; j < sizeAll; ++j) if (j + sizeInA != iOffset) { int iRecNow = partonSystemsPtr->getAll(iSys, j + sizeInA); if ( (j < sizeIn && !event[iRecNow].isRescatteredIncoming()) || (j >= sizeIn && event[iRecNow].isFinal()) ) { if ( (j < sizeIn && event[iRecNow].id() == idRad) || (j >= sizeIn && event[iRecNow].id() == -idRad) ) { double ppNow = event[iRecNow].p() * event[iRad].p() - event[iRecNow].m() * event[iRad].m(); if (ppNow < ppMin) { iRec = iRecNow; ppMin = ppNow; } } } else hasRescattered = true; } // If rescattering then find nearest opposite-flavour recoiler // anywhere in final state. if (iRec == 0 && hasRescattered) { for (int iRecNow = 0; iRecNow < event.size(); ++iRecNow) if (event[iRecNow].id() == -idRad && event[iRecNow].isFinal()) { double ppNow = event[iRecNow].p() * event[iRad].p() - event[iRecNow].m() * event[iRad].m(); if (ppNow < ppMin) { iRec = iRecNow; ppMin = ppNow; otherSystemRec = true; } } } // Find nearest recoiler in same system, charge-squared-weighted, // including initial state, but excluding rescatterer. if (iRec == 0) for (int j = 0; j < sizeAll; ++j) if (j + sizeInA != iOffset) { int iRecNow = partonSystemsPtr->getAll(iSys, j + sizeInA); int chgTypeRecNow = event[iRecNow].chargeType(); if (chgTypeRecNow == 0) continue; if ( (j < sizeIn && !event[iRecNow].isRescatteredIncoming()) || (j >= sizeIn && event[iRecNow].isFinal()) ) { double ppNow = (event[iRecNow].p() * event[iRad].p() - event[iRecNow].m() * event[iRad].m()) / pow2(chgTypeRecNow); if (ppNow < ppMin) { iRec = iRecNow; ppMin = ppNow; } } } // If rescattering then find nearest recoiler in the final state, // charge-squared-weighted. if (iRec == 0 && hasRescattered) { for (int iRecNow = 0; iRecNow < event.size(); ++iRecNow) if (iRecNow != iRad && event[iRecNow].isFinal()) { int chgTypeRecNow = event[iRecNow].chargeType(); if (chgTypeRecNow != 0 && event[iRecNow].isFinal()) { double ppNow = (event[iRecNow].p() * event[iRad].p() - event[iRecNow].m() * event[iRad].m()) / pow2(chgTypeRecNow); if (ppNow < ppMin) { iRec = iRecNow; ppMin = ppNow; otherSystemRec = true; } } } } // Find any nearest recoiler in final state of same system. if (iRec == 0) for (int j = 0; j < sizeOut; ++j) if (j != i) { int iRecNow = partonSystemsPtr->getOut(iSys, j); double ppNow = event[iRecNow].p() * event[iRad].p() - event[iRecNow].m() * event[iRad].m(); if (ppNow < ppMin) { iRec = iRecNow; ppMin = ppNow; } } // Find any nearest recoiler in final state. if (iRec == 0) for (int iRecNow = 0; iRecNow < event.size(); ++iRecNow) if (iRecNow != iRad && event[iRecNow].isFinal()) { double ppNow = event[iRecNow].p() * event[iRad].p() - event[iRecNow].m() * event[iRad].m(); if (ppNow < ppMin) { iRec = iRecNow; ppMin = ppNow; otherSystemRec = true; } } // Fill charge-dipole or photon-dipole end. if (iRec > 0) { // Max scale either by parton scale or by half dipole mass. double pTmax = event[iRad].scale(); if (limitPTmaxIn) { if (iSys == 0) pTmax *= pTmaxFudge; if (iSys > 0 && sizeIn > 0) pTmax *= pTmaxFudgeMI; } else pTmax = 0.5 * m( event[iRad], event[iRec]); int isrType = (event[iRec].isFinal()) ? 0 : event[iRec].mother1(); // This line in case mother is a rescattered parton. while (isrType > 2 + beamOffset) isrType = event[isrType].mother1(); if (isrType > 2) isrType -= beamOffset; dipEnd.push_back( TimeDipoleEnd(iRad, iRec, pTmax, 0, chgType, gamType, isrType, iSys, -1) ); // If hooked up with other system then find which. if (otherSystemRec) { int systemRec = partonSystemsPtr->getSystemOf(iRec); if (systemRec >= 0) dipEnd.back().systemRec = systemRec; dipEnd.back().MEtype = 0; } // Failure to find other end of dipole. } else { infoPtr->errorMsg("Error in TimeShower::setupQEDdip: " "failed to locate any recoiling partner"); } } //-------------------------------------------------------------------------- // Setup a dipole end for a Hidden Valley colour charge. void TimeShower::setupHVdip( int iSys, int i, Event& event, bool limitPTmaxIn) { // Initial values. int iRad = partonSystemsPtr->getOut(iSys, i); int iRec = 0; int idRad = event[iRad].id(); int sizeOut = partonSystemsPtr->sizeOut(iSys); // Hidden Valley colour positive for positive id, and vice versa. // Find opposte HV colour in final state of same system. for (int j = 0; j < sizeOut; ++j) if (j != i) { int iRecNow = partonSystemsPtr->getOut(iSys, j); int idRec = event[iRecNow].id(); if ( (abs(idRec) > 4900000 && abs(idRec) < 4900017) && idRad * idRec < 0) { iRec = iRecNow; break; } } // Else find heaviest other final-state in same system. // (Intended for decays; should mainly be two-body so unique.) double mMax = -sqrt(LARGEM2); if (iRec == 0) for (int j = 0; j < sizeOut; ++j) if (j != i) { int iRecNow = partonSystemsPtr->getOut(iSys, j); if (event[iRecNow].m() > mMax) { iRec = iRecNow; mMax = event[iRecNow].m(); } } // Set up dipole end, or report failure. if (iRec > 0) { // Max scale either by parton scale or by half dipole mass. double pTmax = event[iRad].scale(); if (limitPTmaxIn) { if (iSys == 0) pTmax *= pTmaxFudge; } else pTmax = 0.5 * m( event[iRad], event[iRec]); int colvType = (event[iRad].id() > 0) ? 1 : -1; dipEnd.push_back( TimeDipoleEnd( iRad, iRec, pTmax, 0, 0, 0, 0, iSys, -1, -1, false, true, colvType) ); } else infoPtr->errorMsg("Error in TimeShower::setupHVdip: " "failed to locate any recoiling partner"); } //-------------------------------------------------------------------------- // Select next pT in downwards evolution of the existing dipoles. double TimeShower::pTnext( Event& event, double pTbegAll, double pTendAll) { // Begin loop over all possible radiating dipole ends. dipSel = 0; iDipSel = -1; double pT2sel = pTendAll * pTendAll; for (int iDip = 0; iDip < int(dipEnd.size()); ++iDip) { TimeDipoleEnd& dip = dipEnd[iDip]; // Dipole properties. (Could partly be moved up to prepare??) dip.mRad = event[dip.iRadiator].m(); dip.m2Rad = pow2(dip.mRad); dip.mRec = event[dip.iRecoiler].m(); dip.m2Rec = pow2(dip.mRec); dip.mDip = m( event[dip.iRadiator], event[dip.iRecoiler] ); dip.m2Dip = pow2(dip.mDip); // Find maximum evolution scale for dipole. dip.m2DipCorr = pow2(dip.mDip - dip.mRec) - dip.m2Rad; double pTbegDip = min( pTbegAll, dip.pTmax ); double pT2begDip = min( pow2(pTbegDip), 0.25 * dip.m2DipCorr); // Do QCD, QED or HV evolution if it makes sense. dip.pT2 = 0.; if (pT2begDip > pT2sel) { if (dip.colType != 0) pT2nextQCD(pT2begDip, pT2sel, dip, event); else if (dip.chgType != 0 || dip.gamType != 0) pT2nextQED(pT2begDip, pT2sel, dip, event); else if (dip.colvType != 0) pT2nextHV(pT2begDip, pT2sel, dip, event); // Update if found larger pT than current maximum. if (dip.pT2 > pT2sel) { pT2sel = dip.pT2; dipSel = &dip; iDipSel = iDip; } } } // Return nonvanishing value if found pT bigger than already found. return (dipSel == 0) ? 0. : sqrt(pT2sel); } //-------------------------------------------------------------------------- // Evolve a QCD dipole end. void TimeShower::pT2nextQCD(double pT2begDip, double pT2sel, TimeDipoleEnd& dip, Event& event) { // Lower cut for evolution. Return if no evolution range. double pT2endDip = max( pT2sel, pT2colCut ); if (pT2begDip < pT2endDip) return; // Upper estimate for matrix element weighting and colour factor. // Special cases for triplet recoiling against gluino and octet onia. // Note that g -> g g and g -> q qbar are split on two sides. int colTypeAbs = abs(dip.colType); double wtPSglue = 2.; double colFac = (colTypeAbs == 1) ? 4./3. : 3./2.; if (dip.MEgluinoRec) colFac = 3.; if (dip.isOctetOnium) colFac *= 0.5 * octetOniumColFac; double wtPSqqbar = (colTypeAbs == 2) ? 0.25 * nGluonToQuark : 0.; // Variables used inside evolution loop. (Mainly dummy start values.) dip.pT2 = pT2begDip; int nFlavour = 3; double zMinAbs = 0.5; double pT2min = pT2endDip; double b0 = 4.5; double Lambda2 = Lambda3flav2; double emitCoefGlue = 0.; double emitCoefQqbar = 0.; double emitCoefTot = 0.; double wt = 0.; bool mustFindRange = true; // Begin evolution loop towards smaller pT values. do { // Initialize evolution coefficients at the beginning and // reinitialize when crossing c and b flavour thresholds. if (mustFindRange) { // Determine overestimated z range; switch at c and b masses. if (dip.pT2 > m2b) { nFlavour = 5; pT2min = max( m2b, pT2endDip); b0 = 23./6.; Lambda2 = Lambda5flav2; } else if (dip.pT2 > m2c) { nFlavour = 4; pT2min = max( m2c, pT2endDip); b0 = 25./6.; Lambda2 = Lambda4flav2; } else { nFlavour = 3; pT2min = pT2endDip; b0 = 27./6.; Lambda2 = Lambda3flav2; } zMinAbs = 0.5 - sqrtpos( 0.25 - pT2min / dip.m2DipCorr ); if (zMinAbs < SIMPLIFYROOT) zMinAbs = pT2min / dip.m2DipCorr; // Find emission coefficients for X -> X g and g -> q qbar. emitCoefGlue = wtPSglue * colFac * log(1. / zMinAbs - 1.); emitCoefTot = emitCoefGlue; if (colTypeAbs == 2 && event[dip.iRadiator].id() == 21) { emitCoefQqbar = wtPSqqbar * (1. - 2. * zMinAbs); emitCoefTot += emitCoefQqbar; } // Initialization done for current range. mustFindRange = false; } // Pick pT2 (in overestimated z range) for fixed alpha_strong. if (alphaSorder == 0) { dip.pT2 = dip.pT2 * pow( rndmPtr->flat(), 1. / (alphaS2pi * emitCoefTot) ); // Ditto for first-order alpha_strong. } else if (alphaSorder == 1) { dip.pT2 = Lambda2 * pow( dip.pT2 / Lambda2, pow( rndmPtr->flat(), b0 / emitCoefTot) ); // For second order reject by second term in alpha_strong expression. } else { do dip.pT2 = Lambda2 * pow( dip.pT2 / Lambda2, pow( rndmPtr->flat(), b0 / emitCoefTot) ); while (alphaS.alphaS2OrdCorr(dip.pT2) < rndmPtr->flat() && dip.pT2 > pT2min); } wt = 0.; // If crossed c or b thresholds: continue evolution from threshold. if (nFlavour == 5 && dip.pT2 < m2b) { mustFindRange = true; dip.pT2 = m2b; } else if ( nFlavour == 4 && dip.pT2 < m2c) { mustFindRange = true; dip.pT2 = m2c; // Abort evolution if below cutoff scale, or below another branching. } else { if ( dip.pT2 < pT2endDip) { dip.pT2 = 0.; return; } // Pick kind of branching: X -> X g or g -> q qbar. dip.flavour = 21; dip.mFlavour = 0.; if (colTypeAbs == 2 && emitCoefQqbar > rndmPtr->flat() * emitCoefTot) dip.flavour = 0; // Pick z: either dz/(1-z) or flat dz. if (dip.flavour == 21) { dip.z = 1. - zMinAbs * pow( 1. / zMinAbs - 1., rndmPtr->flat() ); } else { dip.z = zMinAbs + (1. - 2. * zMinAbs) * rndmPtr->flat(); } // Do not accept branching if outside allowed z range. double zMin = 0.5 - sqrtpos( 0.25 - dip.pT2 / dip.m2DipCorr ); if (zMin < SIMPLIFYROOT) zMin = dip.pT2 / dip.m2DipCorr; dip.m2 = dip.m2Rad + dip.pT2 / (dip.z * (1. - dip.z)); if (dip.z > zMin && dip.z < 1. - zMin && dip.m2 * dip.m2Dip < dip.z * (1. - dip.z) * pow2(dip.m2Dip + dip.m2 - dip.m2Rec) ) { // Flavour choice for g -> q qbar. if (dip.flavour == 0) { dip.flavour = min(5, 1 + int(nGluonToQuark * rndmPtr->flat())); dip.mFlavour = particleDataPtr->m0(dip.flavour); } // No z weight, except threshold, if to do ME corrections later on. if (dip.MEtype > 0) { wt = 1.; if (dip.flavour < 10 && dip.m2 < THRESHM2 * pow2(dip.mFlavour)) wt = 0.; // z weight for X -> X g. } else if (dip.flavour == 21 && colTypeAbs == 1) { wt = (1. + pow2(dip.z)) / wtPSglue; } else if (dip.flavour == 21) { wt = (1. + pow3(dip.z)) / wtPSglue; // z weight for g -> q qbar. } else { double beta = sqrtpos( 1. - 4. * pow2(dip.mFlavour) / dip.m2 ); wt = beta * ( pow2(dip.z) + pow2(1. - dip.z) ); } // Suppression factors for dipole to beam remnant. if (dip.isrType != 0) { BeamParticle& beam = (dip.isrType == 1) ? *beamAPtr : *beamBPtr; int iSysRec = dip.systemRec; double xOld = beam[iSysRec].x(); double xNew = xOld * (1. + (dip.m2 - dip.m2Rad) / (dip.m2Dip - dip.m2Rad)); double xMaxAbs = beam.xMax(iSysRec); if (xMaxAbs < 0.) { infoPtr->errorMsg("Warning in TimeShower::pT2nextQCD: " "xMaxAbs negative"); return; } // Firstly reduce by PDF ratio. if (xNew > xMaxAbs) wt = 0.; else { int idRec = event[dip.iRecoiler].id(); double pdfOld = max ( TINYPDF, beam.xfISR( iSysRec, idRec, xOld, dip.pT2) ); double pdfNew = beam.xfISR( iSysRec, idRec, xNew, dip.pT2); wt *= min( 1., pdfNew / pdfOld); } // Secondly optionally reduce by 4 pT2_hard / (4 pT2_hard + m2). if (dampenBeamRecoil) { double pTpT = sqrt(event[dip.iRadiator].pT2() * dip.pT2); wt *= pTpT / (pTpT + dip.m2); } } } } // Iterate until acceptable pT (or have fallen below pTmin). } while (wt < rndmPtr->flat()); } //-------------------------------------------------------------------------- // Evolve a QED dipole end, either charged or photon. void TimeShower::pT2nextQED(double pT2begDip, double pT2sel, TimeDipoleEnd& dip, Event& event) { // Lower cut for evolution. Return if no evolution range. double pT2chgCut = (dip.chgType != 0 && abs(dip.chgType) != 3) ? pT2chgQCut : pT2chgLCut; double pT2endDip = max( pT2sel, pT2chgCut ); if (pT2begDip < pT2endDip) return; // Emission of photon or photon branching. bool hasCharge = (dip.chgType != 0); // Default values. double wtPSgam = 0.; double chg2Sum = 0.; double chg2SumL = 0.; double chg2SumQ = 0.; double zMinAbs = 0.; double emitCoefTot = 0.; // alpha_em at maximum scale provides upper estimate. double alphaEMmax = alphaEM.alphaEM(pT2begDip); double alphaEM2pi = alphaEMmax / (2. * M_PI); // Emission: upper estimate for matrix element weighting; charge factor. if (hasCharge) { wtPSgam = 2.; double chg2 = pow2(dip.chgType / 3.); // Determine overestimated z range. Find evolution coefficient. zMinAbs = 0.5 - sqrtpos( 0.25 - pT2endDip / dip.m2DipCorr ); if (zMinAbs < SIMPLIFYROOT) zMinAbs = pT2endDip / dip.m2DipCorr; emitCoefTot = alphaEM2pi * chg2 * wtPSgam * log(1. / zMinAbs - 1.); // Branching: sum of squared charge factors for lepton and quark daughters. } else { chg2SumL = max(0, min(3, nGammaToLepton)); if (nGammaToQuark > 4) chg2SumQ = 11. / 9.; else if (nGammaToQuark > 3) chg2SumQ = 10. / 9.; else if (nGammaToQuark > 2) chg2SumQ = 6. / 9.; else if (nGammaToQuark > 1) chg2SumQ = 5. / 9.; else if (nGammaToQuark > 0) chg2SumQ = 1. / 9.; // Total sum of squared charge factors. Find evolution coefficient. chg2Sum = chg2SumL + 3. * chg2SumQ; emitCoefTot = alphaEM2pi * chg2Sum; } // Variables used inside evolution loop. dip.pT2 = pT2begDip; double wt; // Begin evolution loop towards smaller pT values. do { // Pick pT2 (in overestimated z range). dip.pT2 = dip.pT2 * pow(rndmPtr->flat(), 1. / emitCoefTot); wt = 0.; // Abort evolution if below cutoff scale, or below another branching. if ( dip.pT2 < pT2endDip) { dip.pT2 = 0.; return; } // Pick z according to dz/(1-z) or flat. if (hasCharge) dip.z = 1. - zMinAbs * pow( 1. / zMinAbs - 1., rndmPtr->flat() ); else dip.z = rndmPtr->flat(); // Do not accept branching if outside allowed z range. double zMin = 0.5 - sqrtpos( 0.25 - dip.pT2 / dip.m2DipCorr ); if (zMin < SIMPLIFYROOT) zMin = dip.pT2 / dip.m2DipCorr; dip.m2 = dip.m2Rad + dip.pT2 / (dip.z * (1. - dip.z)); if (dip.z > zMin && dip.z < 1. - zMin && dip.m2 * dip.m2Dip < dip.z * (1. - dip.z) * pow2(dip.m2Dip + dip.m2 - dip.m2Rec) // For gamma -> f fbar also impose maximum mass. && (hasCharge || dip.m2 < m2MaxGamma) ) { // Photon emission: unique flavour choice. if (hasCharge) { dip.flavour = 22; dip.mFlavour = 0.; // Photon branching: either lepton or quark flavour choice. } else { if (rndmPtr->flat() * chg2Sum < chg2SumL) dip.flavour = 9 + 2 * min(3, 1 + int(chg2SumL * rndmPtr->flat())); else { double rndmQ = 9. * chg2SumQ * rndmPtr->flat(); if (rndmQ < 1.) dip.flavour = 1; else if (rndmQ < 5.) dip.flavour = 2; else if (rndmQ < 6.) dip.flavour = 3; else if (rndmQ < 10.) dip.flavour = 4; else dip.flavour = 5; } dip.mFlavour = particleDataPtr->m0(dip.flavour); } // No z weight, except threshold, if to do ME corrections later on. if (dip.MEtype > 0) { wt = 1.; if (dip.flavour < 20 && dip.m2 < THRESHM2 * pow2(dip.mFlavour)) wt = 0.; // z weight for X -> X gamma. } else if (hasCharge) { wt = (1. + pow2(dip.z)) / wtPSgam; // z weight for gamma -> f fbar. } else { double beta = sqrtpos( 1. - 4. * pow2(dip.mFlavour) / dip.m2 ); wt = beta * ( pow2(dip.z) + pow2(1. - dip.z) ); } // Correct to current value of alpha_EM. double alphaEMnow = alphaEM.alphaEM(dip.pT2); wt *= (alphaEMnow / alphaEMmax); // Suppression factors for dipole to beam remnant. if (dip.isrType != 0) { BeamParticle& beam = (dip.isrType == 1) ? *beamAPtr : *beamBPtr; int iSys = dip.system; double xOld = beam[iSys].x(); double xNew = xOld * (1. + (dip.m2 - dip.m2Rad) / (dip.m2Dip - dip.m2Rad)); double xMaxAbs = beam.xMax(iSys); if (xMaxAbs < 0.) { infoPtr->errorMsg("Warning in TimeShower::pT2nextQED: " "xMaxAbs negative"); return; } // Firstly reduce by PDF ratio. if (xNew > xMaxAbs) wt = 0.; else { int idRec = event[dip.iRecoiler].id(); double pdfOld = max ( TINYPDF, beam.xfISR( iSys, idRec, xOld, dip.pT2) ); double pdfNew = beam.xfISR( iSys, idRec, xNew, dip.pT2); wt *= min( 1., pdfNew / pdfOld); } // Secondly optionally reduce by 4 pT2_hard / (4 pT2_hard + m2). if (dampenBeamRecoil) { double pT24 = 4. * event[dip.iRadiator].pT2(); wt *= pT24 / (pT24 + dip.m2); } } } // Iterate until acceptable pT (or have fallen below pTmin). } while (wt < rndmPtr->flat()); } //-------------------------------------------------------------------------- // Evolve a Hidden Valley dipole end. void TimeShower::pT2nextHV(double pT2begDip, double pT2sel, TimeDipoleEnd& dip, Event& ) { // Lower cut for evolution. Return if no evolution range. double pT2endDip = max( pT2sel, pT2hvCut ); if (pT2begDip < pT2endDip) return; // C_F * alpha_HV/2 pi. int colvTypeAbs = abs(dip.colvType); double colvFac = (colvTypeAbs == 1) ? CFHV : 0.5 * nCHV; double alphaHV2pi = colvFac * (alphaHVfix / (2. * M_PI)); // Determine overestimated z range. Find evolution coefficient. double zMinAbs = 0.5 - sqrtpos( 0.25 - pT2endDip / dip.m2DipCorr ); if (zMinAbs < SIMPLIFYROOT) zMinAbs = pT2endDip / dip.m2DipCorr; double emitCoefTot = alphaHV2pi * 2. * log(1. / zMinAbs - 1.); // Variables used inside evolution loop. dip.pT2 = pT2begDip; double wt; // Begin evolution loop towards smaller pT values. do { // Pick pT2 (in overestimated z range). dip.pT2 = dip.pT2 * pow(rndmPtr->flat(), 1. / emitCoefTot); wt = 0.; // Abort evolution if below cutoff scale, or below another branching. if ( dip.pT2 < pT2endDip) { dip.pT2 = 0.; return; } // Pick z according to dz/(1-z). dip.z = 1. - zMinAbs * pow( 1. / zMinAbs - 1., rndmPtr->flat() ); // Do not accept branching if outside allowed z range. double zMin = 0.5 - sqrtpos( 0.25 - dip.pT2 / dip.m2DipCorr ); if (zMin < SIMPLIFYROOT) zMin = dip.pT2 / dip.m2DipCorr; dip.m2 = dip.m2Rad + dip.pT2 / (dip.z * (1. - dip.z)); if (dip.z > zMin && dip.z < 1. - zMin && dip.m2 * dip.m2Dip < dip.z * (1. - dip.z) * pow2(dip.m2Dip + dip.m2 - dip.m2Rec) ) { // HV gamma or gluon emission: unique flavour choice. dip.flavour = idHV; dip.mFlavour = mHV; // No z weight, except threshold, if to do ME corrections later on. if (dip.MEtype > 0) wt = 1.; // z weight for X -> X g_HV. else if (colvTypeAbs == 1) wt = (1. + pow2(dip.z)) / 2.; else wt = (1. + pow3(dip.z)) / 2.; } // Iterate until acceptable pT (or have fallen below pTmin). } while (wt < rndmPtr->flat()); } //-------------------------------------------------------------------------- // ME corrections and kinematics that may give failure. // Notation: radBef, recBef = radiator, recoiler before emission, // rad, rec, emt = radiator, recoiler, emitted efter emission. // (rad, emt distinguished by colour flow for g -> q qbar.) bool TimeShower::branch( Event& event, bool isInterleaved) { // Find initial particles in dipole branching. int iRadBef = dipSel->iRadiator; int iRecBef = dipSel->iRecoiler; Particle& radBef = event[iRadBef]; Particle& recBef = event[iRecBef]; // Default flavours and colour tags for new particles in dipole branching. int idRad = radBef.id(); int idEmt = dipSel->flavour; int colRad = radBef.col(); int acolRad = radBef.acol(); int colEmt = 0; int acolEmt = 0; iSysSel = dipSel->system; int iSysSelRec = dipSel->systemRec; // Default OK for photon, photon_HV or gluon_HV emission. if (dipSel->flavour == 22 || dipSel->flavour == idHV) { // New colour tag required for gluon emission. } else if (dipSel->flavour == 21 && dipSel->colType > 0) { colEmt = colRad; colRad = event.nextColTag(); acolEmt = colRad; } else if (dipSel->flavour == 21) { acolEmt = acolRad; acolRad = event.nextColTag(); colEmt = acolRad; // New flavours for g -> q qbar; split colours. } else if (dipSel->colType > 0) { idEmt = dipSel->flavour ; idRad = -idEmt; colEmt = colRad; colRad = 0; } else if (dipSel->colType < 0) { idEmt = -dipSel->flavour ; idRad = -idEmt; acolEmt = acolRad; acolRad = 0; // New flavours for gamma -> f fbar, and maybe also colours. } else if (dipSel->gamType == 1 && rndmPtr->flat() > 0.5) { idEmt = -dipSel->flavour ; idRad = -idEmt; if (idRad < 10) colRad = event.nextColTag(); acolEmt = colRad; } else if (dipSel->gamType == 1) { idEmt = dipSel->flavour ; idRad = -idEmt; if (idEmt < 10) colEmt = event.nextColTag(); acolRad = colEmt; } // Construct kinematics in dipole rest frame: // begin simple (like g -> g g). double eRadPlusEmt = 0.5 * (dipSel->m2Dip + dipSel->m2 - dipSel->m2Rec) / dipSel->mDip; double e2RadPlusEmt = pow2(eRadPlusEmt); double pzRadPlusEmt = 0.5 * sqrtpos( pow2(dipSel->m2Dip - dipSel->m2 - dipSel->m2Rec) - 4. * dipSel->m2 * dipSel->m2Rec ) / dipSel->mDip; double pT2corr = dipSel->m2 * (e2RadPlusEmt * dipSel->z * (1. - dipSel->z) - 0.25 * dipSel->m2) / pow2(pzRadPlusEmt); double pTcorr = sqrtpos( pT2corr ); double pzRad = (e2RadPlusEmt * dipSel->z - 0.5 * dipSel->m2) / pzRadPlusEmt; double pzEmt = (e2RadPlusEmt * (1. - dipSel->z) - 0.5 * dipSel->m2) / pzRadPlusEmt; double mRad = dipSel->mRad; double mEmt = 0.; // Kinematics reduction for q -> q gamma_v when m_q > 0 and m_gamma_v > 0. if ( abs(dipSel->colvType) == 1 && dipSel->mFlavour > 0.) { mEmt = dipSel->mFlavour; if (pow2(mRad + mEmt) > dipSel->m2) return false; double m2Emt = pow2(mEmt); double lambda = sqrtpos( pow2(dipSel->m2 - dipSel->m2Rad - m2Emt) - 4. * dipSel->m2Rad * m2Emt ); kRad = 0.5 * (dipSel->m2 - lambda + m2Emt - dipSel->m2Rad) / dipSel->m2; kEmt = 0.5 * (dipSel->m2 - lambda + dipSel->m2Rad - m2Emt) / dipSel->m2; pTcorr *= 1. - kRad - kEmt; double pzMove = kRad * pzRad - kEmt * pzEmt; pzRad -= pzMove; pzEmt += pzMove; // Kinematics reduction for q -> q g/gamma/g_HV when m_q > 0. } else if (abs(dipSel->colType) == 1 || dipSel->chgType != 0 || abs(dipSel->colvType) == 1) { pTcorr *= 1. - dipSel->m2Rad / dipSel->m2; pzRad += pzEmt * dipSel->m2Rad / dipSel->m2; pzEmt *= 1. - dipSel->m2Rad / dipSel->m2; // Kinematics reduction for g -> q qbar or gamma -> f fbar when m_f > 0; } else if (abs(dipSel->flavour) < 20) { mEmt = dipSel->mFlavour; mRad = mEmt; double beta = sqrtpos( 1. - 4. * pow2(mEmt) / dipSel->m2 ); pTcorr *= beta; pzRad = 0.5 * ( (1. + beta) * pzRad + (1. - beta) * pzEmt ); pzEmt = pzRadPlusEmt - pzRad; } // Find rest frame and angles of original dipole. RotBstMatrix M; M.fromCMframe(radBef.p(), recBef.p()); // Evaluate coefficient of azimuthal asymmetry from gluon polarization. findAsymPol( event, dipSel); // Begin construction of new dipole kinematics: pick azimuthal angle. Vec4 pRad, pEmt, pRec; double wtPhi = 1.; do { double phi = 2. * M_PI * rndmPtr->flat(); // Define kinematics of branching in dipole rest frame. pRad = Vec4( pTcorr * cos(phi), pTcorr * sin(phi), pzRad, sqrt( pow2(pTcorr) + pow2(pzRad) + pow2(mRad) ) ); pEmt = Vec4( -pRad.px(), -pRad.py(), pzEmt, sqrt( pow2(pTcorr) + pow2(pzEmt) + pow2(mEmt) ) ); pRec = Vec4( 0., 0., -pzRadPlusEmt, sqrt( pow2(pzRadPlusEmt) + dipSel->m2Rec ) ); // Rotate and boost dipole products to the event frame. pRad.rotbst(M); pEmt.rotbst(M); pRec.rotbst(M); // Azimuthal phi weighting: loop to new phi value if required. if (dipSel->asymPol != 0.) { Vec4 pRadBef = event[iRadBef].p(); Vec4 pAunt = event[dipSel->iAunt].p(); double cosPhi = cosphi( pRad, pAunt, pRadBef ); wtPhi = ( 1. + dipSel->asymPol * (2. * pow2(cosPhi) - 1.) ) / ( 1. + abs(dipSel->asymPol) ); } } while (wtPhi < rndmPtr->flat()) ; // Kinematics when recoiler is initial-state parton. int isrTypeNow = dipSel->isrType; if (isrTypeNow != 0) pRec = 2. * recBef.p() - pRec; // Define new particles from dipole branching. double pTsel = sqrt(dipSel->pT2); Particle rad = Particle(idRad, 51, iRadBef, 0, 0, 0, colRad, acolRad, pRad, mRad, pTsel); Particle emt = Particle(idEmt, 51, iRadBef, 0, 0, 0, colEmt, acolEmt, pEmt, mEmt, pTsel); // Recoiler either in final or in initial state Particle rec = (isrTypeNow == 0) ? Particle(recBef.id(), 52, iRecBef, iRecBef, 0, 0, recBef.col(), recBef.acol(), pRec, dipSel->mRec, pTsel) : Particle(recBef.id(), -53, 0, 0, iRecBef, iRecBef, recBef.col(), recBef.acol(), pRec, 0., 0.); // ME corrections can lead to branching being rejected. if (dipSel->MEtype > 0) { Particle& partner = (dipSel->iMEpartner == iRecBef) ? rec : event[dipSel->iMEpartner]; if ( findMEcorr( dipSel, rad, partner, emt) < rndmPtr->flat() ) return false; } // Rescatter: if the recoiling partner is not in the same system // as the radiator, fix up intermediate systems (can lead // to emissions being vetoed) if (allowRescatter && FIXRESCATTER && isInterleaved && iSysSel != iSysSelRec) { Vec4 pNew = rad.p() + emt.p(); if (!rescatterPropagateRecoil(event, pNew)) return false; } // Save properties to be restored in case of user-hook veto of emission. int eventSizeOld = event.size(); int iRadStatusV = event[iRadBef].status(); int iRadDau1V = event[iRadBef].daughter1(); int iRadDau2V = event[iRadBef].daughter2(); int iRecStatusV = event[iRecBef].status(); int iRecMot1V = event[iRecBef].mother1(); int iRecMot2V = event[iRecBef].mother2(); int iRecDau1V = event[iRecBef].daughter1(); int iRecDau2V = event[iRecBef].daughter2(); int beamOff1 = 1 + beamOffset; int beamOff2 = 2 + beamOffset; int ev1Dau1V = event[beamOff1].daughter1(); int ev2Dau1V = event[beamOff2].daughter1(); // Put new particles into the event record. int iRad = event.append(rad); int iEmt = event.append(emt); int iRec = event.append(rec); // Mark original dipole partons as branched and set daughters/mothers. event[iRadBef].statusNeg(); event[iRadBef].daughters( iRad, iEmt); if (isrTypeNow == 0) { event[iRecBef].statusNeg(); event[iRecBef].daughters( iRec, iRec); } else { event[iRecBef].mothers( iRec, iRec); event[iRec].mothers( iRecMot1V, iRecMot2V); if (iRecMot1V == beamOff1) event[beamOff1].daughter1( iRec); if (iRecMot1V == beamOff2) event[beamOff2].daughter1( iRec); } // Allow veto of branching. If so restore event record to before emission. if ( canVetoEmission && userHooksPtr->doVetoFSREmission(eventSizeOld, event) ) { event.popBack( event.size() - eventSizeOld); event[iRadBef].status( iRadStatusV); event[iRadBef].daughters( iRadDau1V, iRadDau2V); if (isrTypeNow == 0) { event[iRecBef].status( iRecStatusV); event[iRecBef].daughters( iRecDau1V, iRecDau2V); } else { event[iRecBef].mothers( iRecMot1V, iRecMot2V); if (iRecMot1V == beamOff1) event[beamOff1].daughter1( ev1Dau1V); if (iRecMot1V == beamOff2) event[beamOff2].daughter1( ev2Dau1V); } return false; } // For initial-state recoiler also update beam and sHat info. if (isrTypeNow != 0) { BeamParticle& beamRec = (isrTypeNow == 1) ? *beamAPtr : *beamBPtr; double xOld = beamRec[iSysSelRec].x(); double xRec = 2. * pRec.e() / (beamAPtr->e() + beamBPtr->e()); beamRec[iSysSelRec].iPos( iRec); beamRec[iSysSelRec].x( xRec); partonSystemsPtr->setSHat( iSysSelRec, partonSystemsPtr->getSHat(iSysSelRec) * xRec / xOld); } // Photon emission: update to new dipole ends; add new photon "dipole". if (dipSel->flavour == 22) { dipSel->iRadiator = iRad; dipSel->iRecoiler = iRec; dipSel->pTmax = pTsel; if (doQEDshowerByGamma) dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel, 0, 0, 1, 0, iSysSel, 0)); // Gluon emission: update both dipole ends and add two new ones. } else if (dipSel->flavour == 21) { dipSel->iRadiator = iRad; dipSel->iRecoiler = iEmt; dipSel->systemRec = iSysSel; dipSel->isrType = 0; dipSel->pTmax = pTsel; for (int i = 0; i < int(dipEnd.size()); ++i) { if (dipEnd[i].iRadiator == iRecBef && dipEnd[i].iRecoiler == iRadBef && dipEnd[i].colType != 0) { dipEnd[i].iRadiator = iRec; dipEnd[i].iRecoiler = iEmt; // Strive to match colour to anticolour inside closed system. if (dipEnd[i].colType * dipSel->colType > 0) dipEnd[i].iRecoiler = iRad; dipEnd[i].pTmax = pTsel; } } int colType = (dipSel->colType > 0) ? 2 : -2 ; dipEnd.push_back( TimeDipoleEnd(iEmt, iRec, pTsel, colType, 0, 0, isrTypeNow, iSysSel, 0)); dipEnd.back().systemRec = iSysSelRec; dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel, -colType, 0, 0, 0, iSysSel, 0)); // Gluon branching to q qbar: update current dipole and other of gluon. } else if (dipSel->colType != 0) { for (int i = 0; i < int(dipEnd.size()); ++i) { // Strive to match colour to anticolour inside closed system. if ( dipEnd[i].iRecoiler == iRadBef && dipEnd[i].colType * dipSel->colType < 0 ) dipEnd[i].iRecoiler = iEmt; if (dipEnd[i].iRadiator == iRadBef && abs(dipEnd[i].colType) == 2) { dipEnd[i].colType /= 2; if (dipEnd[i].system != dipEnd[i].systemRec) continue; // Note: gluino -> quark + squark gives a deeper radiation dip than // the more obvious alternative photon decay, so is more realistic. dipEnd[i].MEtype = 66; if (&dipEnd[i] == dipSel) dipEnd[i].iMEpartner = iRad; else dipEnd[i].iMEpartner = iEmt; } } dipSel->iRadiator = iEmt; dipSel->iRecoiler = iRec; dipSel->pTmax = pTsel; // Gluon branching to q qbar: also add two charge dipole ends. // Note: gluino -> quark + squark gives a deeper radiation dip than // the more obvious alternative photon decay, so is more realistic. if (doQEDshowerByQ) { int chgType = event[iRad].chargeType(); dipEnd.push_back( TimeDipoleEnd(iRad, iEmt, pTsel, 0, chgType, 0, 0, iSysSel, 66, iEmt)); dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel, 0, -chgType, 0, 0, iSysSel, 66, iRad)); } // Photon branching to f fbar: inactivate photon "dipole"; // optionally add new charge and colour dipole ends. } else if (dipSel->gamType != 0) { dipSel->gamType = 0; int chgType = event[iRad].chargeType(); int colType = event[iRad].colType(); // MEtype = 102 for charge in vector decay. if ( chgType != 0 && ( ( doQEDshowerByQ && colType != 0 ) || ( doQEDshowerByL && colType == 0 ) ) ) { dipEnd.push_back( TimeDipoleEnd(iRad, iEmt, pTsel, 0, chgType, 0, 0, iSysSel, 102, iEmt)); dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel, 0, -chgType, 0, 0, iSysSel, 102, iRad)); } // MEtype = 11 for colour in vector decay. if (colType != 0 && doQCDshower) { dipEnd.push_back( TimeDipoleEnd(iRad, iEmt, pTsel, colType, 0, 0, 0, iSysSel, 11, iEmt)); dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel, -colType, 0, 0, 0, iSysSel, 11, iRad)); } // Photon_HV emission: update to new dipole ends. } else if (dipSel->flavour == 4900022) { dipSel->iRadiator = iRad; dipSel->iRecoiler = iRec; dipSel->pTmax = pTsel; // Gluon_HV emission: update to new dipole ends. } else if (dipSel->flavour == 4900021) { dipSel->iRadiator = iRad; dipSel->iRecoiler = iEmt; dipSel->pTmax = pTsel; for (int i = 0; i < int(dipEnd.size()); ++i) if (dipEnd[i].iRadiator == iRecBef && dipEnd[i].iRecoiler == iRadBef && dipEnd[i].isHiddenValley) { dipEnd[i].iRadiator = iRec; dipEnd[i].iRecoiler = iEmt; dipEnd[i].pTmax = pTsel; } int colvType = (dipSel->colvType > 0) ? 2 : -2 ; dipEnd.push_back( TimeDipoleEnd(iEmt, iRec, pTsel, 0, 0, 0, isrTypeNow, iSysSel, 0, -1, false, true, colvType) ); dipEnd.back().systemRec = iSysSelRec; dipEnd.push_back( TimeDipoleEnd(iEmt, iRad, pTsel, 0, 0, 0, 0, iSysSel, 0, -1, false, true, -colvType) ); } // Copy or set lifetime for new final state. if (event[iRad].id() == event[iRadBef].id()) { event[iRad].tau( event[iRadBef].tau() ); } else { event[iRad].tau( event[iRad].tau0() * rndmPtr->exp() ); event[iEmt].tau( event[iEmt].tau0() * rndmPtr->exp() ); } event[iRec].tau( event[iRecBef].tau() ); // Now update other dipoles that also involved the radiator or recoiler. for (int i = 0; i < int(dipEnd.size()); ++i) { if (dipEnd[i].iRadiator == iRadBef) dipEnd[i].iRadiator = iRad; if (dipEnd[i].iRadiator == iRecBef) dipEnd[i].iRadiator = iRec; if (dipEnd[i].iRecoiler == iRadBef) dipEnd[i].iRecoiler = iRad; if (dipEnd[i].iRecoiler == iRecBef) dipEnd[i].iRecoiler = iRec; if (dipEnd[i].iMEpartner == iRadBef) dipEnd[i].iMEpartner = iRad; if (dipEnd[i].iMEpartner == iRecBef) dipEnd[i].iMEpartner = iRec; } // PS Sept 2010 // Update any junctions downstream of this branching (if necessary) // (This happens, e.g., via LHEF, when adding showers to intermediate // coloured resonances whose decays involved junctions) for (int iJun = 0; iJun < event.sizeJunction(); iJun++) { // Number of incoming colour lines for this junction. int nIncoming = (event.kindJunction(iJun)-1)/2; // Check radiator colour or anticolour, depending on junction kind // (if junction, incoming = anticolours, and vice versa) int colChk = 0; colChk = ( event.kindJunction(iJun) % 2 == 0 ) ? event[iRadBef].col() : event[iRadBef].acol(); // Loop over incoming junction ends for (int iCol = 0; iCol < nIncoming; iCol++) { int colJun = event.endColJunction( iJun, iCol); // If match, update junction end with new upstream (anti)colour if (colJun == colChk) { int colNew = 0; if ( event.kindJunction(iJun) % 2 == 0 ) colNew = colRad; else colNew = acolRad; event.endColJunction( iJun, iCol, colNew ); } } } // Finally update the list of all partons in all systems. partonSystemsPtr->replace(iSysSel, iRadBef, iRad); partonSystemsPtr->addOut(iSysSel, iEmt); partonSystemsPtr->replace(iSysSelRec, iRecBef, iRec); // Done. return true; } //-------------------------------------------------------------------------- // Rescatter: If a dipole stretches between two different systems, those // systems will no longer locally conserve momentum. These // imbalances become problematic when ISR or primordial kT // is switched on as these steps involve Lorentz boosts. // // 'rescatterPropagateRecoil' tries to fix momentum in all // systems by propogating recoil momentum through all // intermediate systems. As the momentum transfer is already // defined, this can lead to internal lines gaining a // virtuality. // Useful definitions for a pair of integers and a vector of pairs typedef pair < int, int > pairInt; typedef vector < pairInt > vectorPairInt; //-------------------------------------------------------------------------- // findParentSystems // Utility routine to find all parent systems of a given system // Returns a vector of pairs of integers with: // a) The system index, including the starting system (negative // if (b) points to a parent system, positive if (b) points // to a daughter system // b) The event record index that is the path out of the system // (if forwards == false, this is an incoming parton to the // system, and is +ve if side A or -ve if side B, // if forwards == true, this is an outgoing parton from the // system). // Returns as empty vector on failure // Note: this assumes single rescattering only and therefore only // one possible parent system inline vectorPairInt findParentSystems(const int sys, Event& event, PartonSystems* partonSystemsPtr, bool forwards) { vectorPairInt parentSystems; parentSystems.reserve(10); int iSysCur = sys; while (true) { // Get two incoming partons int iInA = partonSystemsPtr->getInA(iSysCur); int iInB = partonSystemsPtr->getInB(iSysCur); // Check if either of these links to another system int iIn = 0; if (event[iInA].isRescatteredIncoming()) iIn = iInA; if (event[iInB].isRescatteredIncoming()) iIn = -iInB; // Save the current system to the vector parentSystems.push_back( pairInt(-iSysCur, iIn) ); if (iIn == 0) break; int iInAbs = abs(iIn); int iMother = event[iInAbs].mother1(); iSysCur = partonSystemsPtr->getSystemOf(iMother); if (iSysCur == -1) { parentSystems.clear(); break; } } // while (true) // If forwards is set, change all event record indices to go to daughter // systems rather than parent systems if (forwards) { vectorPairInt::reverse_iterator rit; for (rit = parentSystems.rbegin(); rit < (parentSystems.rend() - 1); ++rit) { pairInt &cur = *rit; pairInt &next = *(rit + 1); cur.first = -cur.first; cur.second = (next.second < 0) ? -event[abs(next.second)].mother1() : event[abs(next.second)].mother1(); } } return parentSystems; } //-------------------------------------------------------------------------- // rescatterPropagateRecoil // Fix up momentum in all intermediate systems when radiator and recoiler // systems are different. The strategy is to look at all parent systems // from the radiator system and the recoiler system and find where they // intersect. bool TimeShower::rescatterPropagateRecoil( Event& event, Vec4& pNew) { // Some useful variables for later int iRadBef = dipSel->iRadiator; iSysSel = dipSel->system; int iSysSelRec = dipSel->systemRec; Vec4 pImbal = pNew - event[iRadBef].p(); // Store changes locally at first in case we veto the branching // eventMod stores index into the event record and new 4-vector vector < pair < int, Vec4 > > eventMod; eventMod.reserve(10); // systemMod stores system index (iSys) and system-parton index (iMem) // iMem >= 0 - index into outgoing partons (iOut) // iMem == -1 - incoming A // iMem == -2 - incoming B vectorPairInt systemMod; systemMod.reserve(10); // Find all parent systems from radiating and recoiling systems vectorPairInt radParent = findParentSystems(iSysSel, event, partonSystemsPtr, false); vectorPairInt recParent = findParentSystems(iSysSelRec, event, partonSystemsPtr, true); if (radParent.size() == 0 || recParent.size() == 0) { // This should never happen infoPtr->errorMsg("Error in TimeShower::rescatterPropagateRecoil: " "couldn't find parent system; branching vetoed"); return false; } // Find the system that connects radiating and recoiling system bool foundPath = false; unsigned int iRadP = 0; unsigned int iRecP = 0; for (iRadP = 0; iRadP < radParent.size(); iRadP++) { for (iRecP = 0; iRecP < recParent.size(); iRecP++) if (abs(radParent[iRadP].first) == abs(recParent[iRecP].first)) { foundPath = true; break; } if (foundPath) break; } if (!foundPath) { // Can fail e.g. for QED dipoles where there is no connection // between radiator and recoiler systems infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " "couldn't find recoil path; branching vetoed"); return false; } // Join together to form complete path from radiating system // to recoiling system vectorPairInt path; if (radParent.size() > 1) path.assign(radParent.begin(), radParent.begin() + iRadP); if (recParent.size() > 1) path.insert(path.end(), recParent.rend() - iRecP - 1, recParent.rend() - 1); // Follow the path fixing up momenta as we go for (unsigned int i = 0; i < path.size(); i++) { // Line out of the current system bool isIncoming = (path[i].first < 0) ? true : false; int iSysCur = abs(path[i].first); bool isIncomingA = (path[i].second > 0) ? true : false; int iLink = abs(path[i].second); int iMemCur; if (isIncoming) iMemCur = (isIncomingA) ? -1 : -2; else { iMemCur = -1; for (int j = 0; j < partonSystemsPtr->sizeOut(iSysCur); j++) if (partonSystemsPtr->getOut(iSysCur, j) == iLink) { iMemCur = j; break; } if (iMemCur == -1) { // This should never happen infoPtr->errorMsg("Error in TimeShower::rescatterPropagateRecoil: " "couldn't find parton system; branching vetoed"); return false; } } Vec4 pMod = (isIncoming) ? event[iLink].p() + pImbal : event[iLink].p() - pImbal; eventMod.push_back(pair (iLink, pMod)); systemMod.push_back(pairInt(iSysCur, iMemCur)); // Calculate sHat of iSysCur int iInCurA = partonSystemsPtr->getInA(iSysCur); int iInCurB = partonSystemsPtr->getInB(iSysCur); Vec4 pTotCur = event[iInCurA].p() + event[iInCurB].p(); // If iMemCur is -1 or -2, then we must have changed the sHat of iSysCur if (iMemCur < 0) pTotCur += (isIncoming) ? pImbal : -pImbal; double sHatCur = pTotCur.m2Calc(); // The fixed-up incoming and outgoing partons should not have // too large a virtuality in relation to the system mass-square. if (abs(pMod.m2Calc()) > MAXVIRTUALITYFRACTION * sHatCur) { infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " "virtuality much larger than sHat; branching vetoed"); return false; } // Outgoing ones should also not have too large negative energy // in the rest frame of the system. if (!isIncoming && pMod * pTotCur < -MAXNEGENERGYFRACTION * sHatCur) { infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " "rest frame energy too negative; branching vetoed"); return false; } // Veto negative sHat. if (sHatCur < 0.0) { infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " "sHat became negative; branching vetoed"); return false; } // Line into the new current system iLink = (isIncoming) ? event[iLink].mother1() : event[iLink].daughter1(); iSysCur = partonSystemsPtr->getSystemOf(iLink, true); if (!isIncoming) iMemCur = (isIncomingA) ? -1 : -2; else { iMemCur = -1; for (int j = 0; j < partonSystemsPtr->sizeOut(iSysCur); j++) if (partonSystemsPtr->getOut(iSysCur, j) == iLink) { iMemCur = j; break; } if (iMemCur == -1) { // This should never happen infoPtr->errorMsg("Error in TimeShower::rescatterPropagateRecoil: " "couldn't find parton system; branching vetoed"); return false; } } pMod = (isIncoming) ? event[iLink].p() + pImbal : event[iLink].p() - pImbal; eventMod.push_back(pair (iLink, pMod)); systemMod.push_back(pairInt(iSysCur, iMemCur)); // Calculate sHat of iSysCur iInCurA = partonSystemsPtr->getInA(iSysCur); iInCurB = partonSystemsPtr->getInB(iSysCur); pTotCur = event[iInCurA].p() + event[iInCurB].p(); // If iMemCur is -1 or -2, then we must have changed the sHat of iSysCur if (iMemCur < 0) pTotCur += (isIncoming) ? pImbal : -pImbal; sHatCur = pTotCur.m2Calc(); // The fixed-up incoming and outgoing partons should not have // too large a virtuality in relation to the system mass-square. if (abs(pMod.m2Calc()) > MAXVIRTUALITYFRACTION * sHatCur) { infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " "virtuality much larger than sHat; branching vetoed"); return false; } // Outgoing ones should also not have too large negative energy // in the rest frame of the system. if (!isIncoming && pMod * pTotCur < -MAXNEGENERGYFRACTION * sHatCur) { infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " "rest frame energy too negative; branching vetoed"); return false; } // Veto negative sHat if (sHatCur < 0.0) { infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " "sHat became negative; branching vetoed"); return false; } // Do negative energy veto if (VETONEGENERGY && pMod.e() < 0.0) { infoPtr->errorMsg("Warning in TimeShower::rescatterPropagateRecoil: " "energy became negative; branching vetoed"); return false; } } // for (unsigned int i = 0; i < path.size(); i++) // If no vetos by this point, apply the changes to the event record // An incoming particle with changed momentum is given status code -54, // an outgoing particle with changed momentum is given status code -55 for (unsigned int i = 0; i < eventMod.size(); i++) { int idx = eventMod[i].first; Vec4 &pMod = eventMod[i].second; int iSys = systemMod[i].first; int iMem = systemMod[i].second; // If incoming to a process then set the copy to be the mother if (event[idx].isRescatteredIncoming()) { int mother1 = event[idx].mother1(); idx = event.copy(idx, -54); event[mother1].daughters(idx, idx); // Update beam information if necessary double eCM = sqrt(m2( beamAPtr->p(), beamBPtr->p())); if (iMem == -1) { partonSystemsPtr->setInA(iSys, idx); (*beamAPtr)[iSys].x((pMod.e() + pMod.pz()) / eCM); (*beamAPtr)[iSys].m(pMod.mCalc()); (*beamAPtr)[iSys].p(pMod); (*beamAPtr)[iSys].iPos(idx); } else if (iMem == -2) { partonSystemsPtr->setInB(iSys, idx); (*beamBPtr)[iSys].x((pMod.e() - pMod.pz()) / eCM); (*beamBPtr)[iSys].m(pMod.mCalc()); (*beamBPtr)[iSys].p(pMod); (*beamBPtr)[iSys].iPos(idx); } else { // This should never happen infoPtr->errorMsg("Error in TimeShower::rescatterPropagateRecoil: " "internal bookeeping error"); } // Otherwise set the new event record entry to be the daughter } else { int daughter1 = event[idx].daughter1(); idx = event.copy(idx, 55); event[idx].statusNeg(); event[daughter1].mothers(idx, idx); partonSystemsPtr->setOut(iSys, iMem, idx); } event[idx].p( eventMod[i].second ); event[idx].m( event[idx].mCalc() ); } return true; } //-------------------------------------------------------------------------- // Find class of QCD ME correction. // MEtype classification follow codes in Norrbin article, // additionally -1 = try to find type, 0 = no ME corrections. // Warning: not yet tried out to do a correct assignment in // arbitrary multiparton configurations! ?? void TimeShower::findMEtype( Event& event, TimeDipoleEnd& dip) { // Initial value. Mark if no ME corrections to be applied. bool setME = true; if (!doMEcorrections) setME = false; int iMother = event[dip.iRadiator].mother1(); int iMother2 = event[dip.iRadiator].mother2(); // Allow ME corrections for Hidden Valley pair in 2 -> 2. if (dip.isHiddenValley && event[dip.iRecoiler].id() == -event[dip.iRadiator].id()); // Else no ME corrections in 2 -> n processes. else { if (iMother2 != iMother && iMother2 != 0) setME = false; if (event[dip.iRecoiler].mother1() != iMother) setME = false; if (event[dip.iRecoiler].mother2() != iMother2) setME = false; } // No ME corrections for recoiler in initial state. if (event[dip.iRecoiler].status() < 0) setME = false; // No ME corrections for recoiler not in same system if (dip.system != dip.systemRec) setME = false; // Done if no ME to be set. if (!setME) { dip.MEtype = 0; return; } // If no ME partner set, assume it is the recoiler. if (dip.iMEpartner < 0) dip.iMEpartner = dip.iRecoiler; // Now begin processing of colour dipole, including Hidden Valley. if (dip.colType != 0 || dip.colvType != 0) { bool isHiddenColour = (dip.colvType != 0); // Find daughter types (may or may not be used later on). int idDau1 = event[dip.iRadiator].id(); int idDau2 = event[dip.iMEpartner].id(); int dau1Type = findMEparticle(idDau1, isHiddenColour); int dau2Type = findMEparticle(idDau2, isHiddenColour); int minDauType = min(dau1Type, dau2Type); int maxDauType = max(dau1Type, dau2Type); // Reorder dipole ends in kinematics. Split ME expression in two sides. dip.MEorder = (dau2Type >= dau1Type); dip.MEsplit = (maxDauType <= 6); dip.MEgluinoRec = false; // If type already set (or set not to have) then done. if (minDauType == 0 && dip.MEtype < 0) dip.MEtype = 0; if (dip.MEtype >= 0) return; dip.MEtype = 0; // For H -> gg -> ggg we found that DGLAP kernels do better than eikonal. if (dau1Type == 4 && dau2Type == 4) return; // Find mother type. int idMother = 0; if ( event[dip.iRecoiler].mother1() == iMother && iMother >= 0) idMother = event[iMother].id(); int motherType = (idMother != 0) ? findMEparticle(idMother, isHiddenColour) : 0; // When a mother if not known then use colour and spin content to guess. if (motherType == 0) { int col1 = event[dip.iRadiator].col(); int acol1 = event[dip.iRadiator].acol(); int col2 = event[dip.iMEpartner].col(); int acol2 = event[dip.iMEpartner].acol(); // spinT = 0/1 = integer or half-integer. int spinT = ( event[dip.iRadiator].spinType() + event[dip.iMEpartner].spinType() )%2; // Colour singlet mother. if ( col1 == acol2 && acol1 == col2 ) motherType = (spinT == 0) ? 7 : 9; // Colour octet mother. else if ( (col1 == acol2 && acol1 != 0 && col2 != 0) || (acol1 == col2 && col1 != 0 && acol2 != 0) ) motherType = (spinT == 0) ? 4 : 5; // Colour triplet mother. else if ( (col1 == acol2 && acol1 != col2) || (acol1 == col2 && col1 != acol2) ) motherType = (spinT == 0) ? 2 : 1; // If no colours are matched then cannot have common mother, so done. else return; } // Now start from default, which is eikonal ME corrections, // and try to find matching ME cases below. int MEkind = 0; int MEcombi = 4; dip.MEmix = 0.5; // Hidden Valley with massive gamma_v covered by two special cases. if (isHiddenColour && brokenHVsym) { MEkind = (dau2Type == 0 || dau2Type > 6) ? 30 : 31; dip.MEtype = 5 * MEkind + 1; return; } // Triplet recoiling against gluino needs enhanced radiation // to match to matrix elements. dip.MEgluinoRec = (dau1Type >= 1 && dau1Type <= 3 && dau2Type == 5); // Vector/axial vector -> q + qbar. if (minDauType == 1 && maxDauType == 1 && (motherType == 4 || motherType == 7) ) { MEkind = 2; if (idMother == 21 || idMother == 22) MEcombi = 1; else if (idMother == 23 || idDau1 + idDau2 == 0) { MEcombi = 3; dip.MEmix = gammaZmix( event, iMother, dip.iRadiator, dip.iRecoiler ); } else if (idMother == 24) MEcombi = 4; } // For chi -> chi q qbar, use V/A -> q qbar as first approximation. else if (minDauType == 1 && maxDauType == 1 && motherType == 9) MEkind = 2; // q -> q + V. else if (minDauType == 1 && maxDauType == 7 && motherType == 1) MEkind = 3; if (idDau1 == 22 || idDau2 == 22) MEcombi = 1; // Scalar/pseudoscalar -> q + qbar; q -> q + S. else if (minDauType == 1 && maxDauType == 1 && motherType == 8) { MEkind = 4; if (idMother == 25 || idMother == 35 || idMother == 37) MEcombi = 1; else if (idMother == 36) MEcombi = 2; } else if (minDauType == 1 && maxDauType == 8 && motherType == 1) MEkind = 5; // V -> ~q + ~qbar; ~q -> ~q + V; S -> ~q + ~qbar; ~q -> ~q + S. else if (minDauType == 2 && maxDauType == 2 && (motherType == 4 || motherType == 7) ) MEkind = 6; else if (minDauType == 2 && (maxDauType == 4 || maxDauType == 7) && motherType == 2) MEkind = 7; else if (minDauType == 2 && maxDauType == 2 && motherType == 8) MEkind = 8; else if (minDauType == 2 && maxDauType == 8 && motherType == 2) MEkind = 9; // chi -> q + ~qbar; ~q -> q + chi; q -> ~q + chi. else if (minDauType == 1 && maxDauType == 2 && motherType == 9) MEkind = 10; else if (minDauType == 1 && maxDauType == 9 && motherType == 2) MEkind = 11; else if (minDauType == 2 && maxDauType == 9 && motherType == 1) MEkind = 12; // ~g -> q + ~qbar; ~q -> q + ~g; q -> ~q + ~g. else if (minDauType == 1 && maxDauType == 2 && motherType == 5) MEkind = 13; else if (minDauType == 1 && maxDauType == 5 && motherType == 2) MEkind = 14; else if (minDauType == 2 && maxDauType == 5 && motherType == 1) MEkind = 15; // In cases where coloured spin 1 particle involved use spin 0. // V_coloured -> q + l. else if (minDauType == 1 && maxDauType == 9 && motherType == 3) MEkind = 11; // q -> V_coloured + l; else if (minDauType == 3 && maxDauType == 9 && motherType == 1) MEkind = 12; // g (+V, S) -> ~g + ~g (eikonal approximation). else if (minDauType == 5 && maxDauType == 5) MEkind = 16; // Save ME type and gamma_5 admixture. dip.MEtype = 5 * MEkind + MEcombi; // Now begin processing of charge dipole - still primitive. } else if (dip.chgType != 0) { // Set defaults for QED case; then possibly done. dip.MEorder = true; dip.MEsplit = true; if (dip.MEtype >= 0) return; // So far only ME corrections for q qbar or l lbar. int idDau1 = event[dip.iRadiator].id(); int idDau2 = event[dip.iMEpartner].id(); if (abs(idDau1) < 9 && abs(idDau2) < 9 && idDau1 * idDau2 < 0) ; else if (abs(idDau1) > 10 && abs(idDau1) < 19 && abs(idDau2) > 10 && abs(idDau2) < 19 && idDau1 * idDau2 < 0) ; else { dip.MEtype = 0; return; } // Distinguish charge sum != 0 or = 0; in latter assume vector source. dip.MEtype = 101; if (idDau1 + idDau2 == 0) dip.MEtype = 102; dip.MEmix = 1.; } } //-------------------------------------------------------------------------- // Find type of particle for ME type: 0 = unknown, 1 = quark, 2 = squark, // 3 = spare triplet, 4 = gluon, 5 = gluino, 6 = spare octet, // 7 = vector boson, 8 = colourless scalar, 9 = colourless spin 1/2. int TimeShower::findMEparticle( int id, bool isHiddenColour) { // find colour and spin of particle. int type = 0; int colType = abs(particleDataPtr->colType(id)); int spinType = particleDataPtr->spinType(id); // For hidden valley particle treat HV colour as normal one. // Note: no need to assign gv/gammav since not in ME. if (isHiddenColour) { colType = 0; int idAbs = abs(id); if ( (idAbs > 4900000 && idAbs < 4900007) || (idAbs > 4900010 && idAbs < 4900017) || idAbs == 4900101) colType = 1; } // Find particle type from colour and spin. if (colType == 1 && spinType == 2) type = 1; else if (colType == 1 && spinType == 1) type = 2; else if (colType == 1) type = 3; else if (colType == 2 && spinType == 3) type = 4; else if (colType == 2 && spinType == 2) type = 5; else if (colType == 2) type = 6; else if (colType == 0 && spinType == 3) type = 7; else if (colType == 0 && spinType == 1) type = 8; else if (colType == 0 && spinType == 2) type = 9; // Done. return type; } //-------------------------------------------------------------------------- // Find mixture of V and A in gamma/Z: energy- and flavour-dependent. double TimeShower::gammaZmix( Event& event, int iRes, int iDau1, int iDau2) { // Try to identify initial flavours; use e+e- as default. int idIn1 = -11; int idIn2 = 11; int iIn1 = (iRes >= 0) ? event[iRes].mother1() : -1; int iIn2 = (iRes >= 0) ? event[iRes].mother2() : -1; if (iIn1 >=0) idIn1 = event[iIn1].id(); if (iIn2 >=0) idIn2 = event[iIn1].id(); // In processes f + g/gamma -> f + Z only need find one fermion. if (idIn1 == 21 || idIn1 == 22) idIn1 = -idIn2; if (idIn2 == 21 || idIn2 == 22) idIn2 = -idIn1; // Initial flavours and couplings; return if don't make sense. if (idIn1 + idIn2 != 0 ) return 0.5; int idInAbs = abs(idIn1); if (idInAbs == 0 || idInAbs > 18 ) return 0.5; double ei = coupSMPtr->ef(idInAbs); double vi = coupSMPtr->vf(idInAbs); double ai = coupSMPtr->af(idInAbs); // Final flavours and couplings; return if don't make sense. if (event[iDau1].id() + event[iDau2].id() != 0) return 0.5; int idOutAbs = abs(event[iDau1].id()); if (idOutAbs == 0 || idOutAbs >18 ) return 0.5; double ef = coupSMPtr->ef(idOutAbs); double vf = coupSMPtr->vf(idOutAbs); double af = coupSMPtr->af(idOutAbs); // Calculate prefactors for interference and resonance part. Vec4 psum = event[iDau1].p() + event[iDau2].p(); double sH = psum.m2Calc(); double intNorm = 2. * thetaWRat * sH * (sH - mZ*mZ) / ( pow2(sH - mZ*mZ) + pow2(sH * gammaZ / mZ) ); double resNorm = pow2(thetaWRat * sH) / ( pow2(sH - mZ*mZ) + pow2(sH * gammaZ / mZ) ); // Calculate vector and axial expressions and find mix. double vect = ei*ei * ef*ef + ei*vi * intNorm * ef*vf + (vi*vi + ai*ai) * resNorm * vf*vf; double axiv = (vi*vi + ai*ai) * resNorm * af*af; return vect / (vect + axiv); } //-------------------------------------------------------------------------- // Set up to calculate QCD ME correction with calcMEcorr. // Normally for primary particles, but also from g/gamma -> f fbar. double TimeShower::findMEcorr(TimeDipoleEnd* dip, Particle& rad, Particle& partner, Particle& emt) { // Initial values and matrix element kind. double wtME = 1.; double wtPS = 1.; int MEkind = dip->MEtype / 5; int MEcombi = dip->MEtype % 5; // Construct ME variables. Vec4 sum = rad.p() + partner.p() + emt.p(); double eCMME = sum.mCalc(); double x1 = 2. * (sum * rad.p()) / pow2(eCMME); double x2 = 2. * (sum * partner.p()) / pow2(eCMME); double r1 = rad.m() / eCMME; double r2 = partner.m() / eCMME; double r3 = 0.; // Evaluate kinematics for Hidden Valley with massive gamma_v. double gammavCorr = 1.; if (dip->colvType != 0 && brokenHVsym) { r3 = emt.m() / eCMME; double x3Tmp = 2. - x1 - x2; gammavCorr = x3Tmp / (x3Tmp - kRad * (x1 + x3Tmp)); // For Q_v Qbar_v pair correct kinematics to common average mass. if (MEkind == 31) { double m2Pair = (rad.p() + partner.p()).m2Calc(); double m2Avg = 0.5 * (rad.m2() + partner.m2()) - 0.25 * pow2(rad.m2() - partner.m2()) / m2Pair; r1 = sqrt(m2Avg) / eCMME; r2 = r1; double xShift = 0.5 * (x1 + x2) * (partner.m2() - rad.m2()) / m2Pair; x1 += xShift; x2 -= xShift; } } // Derived ME variables, suitably protected. double x1minus = max(XMARGIN, 1. + r1*r1 - r2*r2 - x1); double x2minus = max(XMARGIN, 1. + r2*r2 - r1*r1 - x2) ; double x3 = max(XMARGIN, 2. - x1 - x2); // Begin processing of QCD dipoles. if (dip->colType !=0 || dip->colvType != 0) { // Evaluate normal ME, for proper order of particles. if (dip->MEorder) wtME = calcMEcorr(MEkind, MEcombi, dip->MEmix, x1, x2, r1, r2, r3); else wtME = calcMEcorr(MEkind, MEcombi, dip->MEmix, x2, x1, r2, r1, r3); // Split up total ME when two radiating particles. if (dip->MEsplit) wtME = wtME * x1minus / x3; // Evaluate shower rate to be compared with. wtPS = 2. / ( x3 * x2minus ); if (dip->MEgluinoRec) wtPS *= 9./4.; if (dip->colvType != 0 && brokenHVsym) wtPS *= gammavCorr; // For generic charge combination currently only massless expression. // (Masses included only to respect phase space boundaries.) } else if (dip->chgType !=0 && dip->MEtype == 101) { double chg1 = particleDataPtr->charge(rad.id()); double chg2 = particleDataPtr->charge(partner.id()); wtME = (x1*x1 + x2*x2) * pow2( chg1 * x1minus / x3 - chg2 * x2minus / x3 ); wtPS = 2. * ( chg1*chg1 * x1minus / x3 + chg2*chg2 * x2minus / x3 ); // For flavour neutral system assume vector source and include masses. } else if (dip->chgType !=0 && dip->MEtype == 102) { wtME = calcMEcorr(2, 1, dip->MEmix, x1, x2, r1, r2) * x1minus / x3; wtPS = 2. / ( x3 * x2minus ); } if (wtME > wtPS) infoPtr->errorMsg("Warning in TimeShower::findMEcorr: " "ME weight above PS one"); // Return ratio of actual ME to assumed PS rate of emission. return wtME / wtPS; } //-------------------------------------------------------------------------- // Matrix elements for gluon (or photon) emission from // a two-body state; to be used by the parton shower routine. // Here x_i = 2 E_i/E_cm, r_i = m_i/E_cm and // 1/sigma_0 d(sigma)/d(x_1)d(x_2) = (alpha-strong/2 pi) * C_F * (this), // i.e. normalization is such that one recovers the familiar // (x_1^2 + x_2^2)/((1-x_1)*(1-x_2)) for the massless case. // Coupling structure: // kind = 1 : eikonal soft-gluon expression (spin-independent) // = 2 : V -> q qbar (V = vector/axial vector colour singlet) // = 3 : q -> q V // = 4 : S -> q qbar (S = scalar/pseudoscalar colour singlet) // = 5 : q -> q S // = 6 : V -> ~q ~qbar (~q = squark) // = 7 : ~q -> ~q V // = 8 : S -> ~q ~qbar // = 9 : ~q -> ~q S // = 10 : chi -> q ~qbar (chi = neutralino/chargino) // = 11 : ~q -> q chi // = 12 : q -> ~q chi // = 13 : ~g -> q ~qbar // = 14 : ~q -> q ~g // = 15 : q -> ~q ~g // = 16 : (9/4)*(eikonal) for gg -> ~g ~g // = 30 : Dv -> d qv (Dv= hidden valley fermion, qv= valley scalar) // = 31 : S -> Dv Dvbar (S=scalar color singlet) // Note that the order of the decay products is important. // combi = 1 : pure non-gamma5, i.e. vector/scalar/... // = 2 : pure gamma5, i.e. axial vector/pseudoscalar/.... // = 3 : mixture mix*(combi=1) + (1-mix)*(combi=2) // = 4 : mixture (combi=1) +- (combi=2) double TimeShower::calcMEcorr( int kind, int combiIn, double mixIn, double x1, double x2, double r1, double r2, double r3) { // Frequent variable combinations. double x3 = 2. - x1 - x2; double x1s = x1 * x1; double x2s = x2 * x2; double x3s = x3 * x3; double x1c = x1 * x1s; double x2c = x2 * x2s; double x3c = x3 * x3s; double r1s = r1 * r1; double r2s = r2 * r2; double r1c = r1 * r1s; double r2c = r2 * r2s; double r1q = r1s * r1s; double r2q = r2s * r2s; double prop1 = 1. + r1s - r2s - x1; double prop2 = 1. + r2s - r1s - x2; double prop1s = prop1 * prop1; double prop2s = prop2 * prop2; double prop12 = prop1 * prop2; double prop13 = prop1 * x3; double prop23 = prop2 * x3; // Special case: Hidden-Valley massive photon. double r3s = r3 * r3; double prop3 = r3s - x3; double prop3s = prop3 * prop3; if (kind == 30) prop23 = prop2 * prop3; // Check input values. Return zero outside allowed phase space. if (x1 - 2.*r1 < XMARGIN || prop1 < XMARGIN) return 0.; if (x2 - 2.*r2 < XMARGIN || prop2 < XMARGIN) return 0.; // Limits not worked out for r3 > 0. if (kind != 30 && kind != 31) { if (x1 + x2 - 1. - pow2(r1+r2) < XMARGIN) return 0.; // Note: equivalent rewritten form 4. * ( (1. - x1) * (1. - x2) // * (1. - r1s - r2s - x3) + r1s * (1. - x2s - x3) + r2s // * (1. - x1s - x3) - pow2(r1s - r2s) ) gives about same result. if ( (x1s - 4.*r1s) * (x2s - 4.*r2s) - pow2( 2. * (1. - x1 - x2 + r1s + r2s) + x1*x2 ) < XMARGIN * (XMARGINCOMB + r1 + r2) ) return 0.; } // Initial values; phase space. int combi = max(1, min(4, combiIn) ); double mix = max(0., min(1., mixIn) ); bool isSet1 = false; bool isSet2 = false; bool isSet4 = false; double ps = sqrtpos( pow2(1. - r1*r1 - r2*r2) - pow2(2. * r1 * r2) ); double rLO = 0., rFO = 0., rLO1 = 0., rFO1 = 0., rLO2 = 0., rFO2 = 0., rLO4 = 0., rFO4 = 0.; double offset = 0; // Select which kind of ME to use. switch (kind) { // case 1 is equal to default, i.e. eikonal expression. // V -> q qbar (V = gamma*/Z0/W+-/...). case 2: if (combi == 1 || combi == 3) { rLO1 = ps*(2.-r1s-r1q+6.*r1*r2-r2s+2.*r1s*r2s-r2q)/2.; rFO1 = -(3.+6.*r1s+r1q-6.*r1*r2+6.*r1c*r2-2.*r2s-6.*r1s*r2s +6.*r1*r2c+r2q-3.*x1+6.*r1*r2*x1+2.*r2s*x1+x1s-2.*r1s*x1s +3.*r1s*x3+6.*r1*r2*x3-r2s*x3-2.*x1*x3-5.*r1s*x1*x3 +r2s*x1*x3+x1s*x3-3.*x3s-3.*r1s*x3s+r2s*x3s +2.*x1*x3s+x3c-x2) /prop2s -2.*(-3.+r1s-6.*r1*r2+6.*r1c*r2+3.*r2s-4.*r1s*r2s +6.*r1*r2c+2.*x1+3.*r1s*x1+r2s*x1-x1s-r1s*x1s -r2s*x1s+4.*x3+2.*r1s*x3+3.*r1*r2*x3-r2s*x3-3.*x1*x3 -2.*r1s*x1*x3+x1s*x3-x3s-r1s*x3s+r1*r2*x3s+x1*x3s) /prop12 -(-1.+2.*r1s+r1q+6.*r1*r2+6.*r1c*r2-2.*r2s-6.*r1s*r2s +6.*r1*r2c+r2q-x1-2.*r1s*x1-6.*r1*r2*x1+8.*r2s*x1+x1s -2.*r2s*x1s-r1s*x3+r2s*x3-r1s*x1*x3+r2s*x1*x3+x1s*x3+x2) /prop1s; rFO1 = rFO1/2.; isSet1 = true; } if (combi == 2 || combi == 3) { rLO2 = ps*(2.-r1s-r1q-6.*r1*r2-r2s+2.*r1s*r2s-r2q)/2.; rFO2 = -(3.+6.*r1s+r1q+6.*r1*r2-6.*r1c*r2-2.*r2s-6.*r1s*r2s -6.*r1*r2c+r2q-3.*x1-6.*r1*r2*x1+2.*r2s*x1+x1s-2.*r1s*x1s +3.*r1s*x3-6.*r1*r2*x3-r2s*x3-2.*x1*x3-5.*r1s*x1*x3 +r2s*x1*x3+x1s*x3-3.*x3s-3.*r1s*x3s+r2s*x3s+2.*x1*x3s+x3c-x2) /prop2s -2.*(-3+r1s+6.*r1*r2-6.*r1c*r2+3.*r2s-4.*r1s*r2s-6.*r1*r2c +2.*x1+3.*r1s*x1+r2s*x1-x1s-r1s*x1s-r2s*x1s+4.*x3+2.*r1s*x3 -3.*r1*r2*x3-r2s*x3-3.*x1*x3-2.*r1s*x1*x3+x1s*x3-x3s-r1s*x3s -r1*r2*x3s+x1*x3s) /prop12 -(-1.+2.*r1s+r1q-6.*r1*r2-6.*r1c*r2-2.*r2s-6.*r1s*r2s -6.*r1*r2c+r2q-x1-2.*r1s*x1+6.*r1*r2*x1+8.*r2s*x1+x1s -2.*r2s*x1s-r1s*x3+r2s*x3-r1s*x1*x3+r2s*x1*x3+x1s*x3+x2) /prop1s; rFO2 = rFO2/2.; isSet2 = true; } if (combi == 4) { rLO4 = ps*(2.-r1s-r1q-r2s+2.*r1s*r2s-r2q)/2.; rFO4 = (1.-r1q+6.*r1s*r2s-r2q+x1+3.*r1s*x1-9.*r2s*x1-3.*x1s -r1s*x1s+3.*r2s*x1s+x1c-x2-r1s*x2+r2s*x2-r1s*x1*x2+r2s*x1*x2 +x1s*x2) /prop1s -2.*(1.+r1s+r2s-4.*r1s*r2s+r1s*x1+2.*r2s*x1-x1s-r2s*x1s +2.*r1s*x2+r2s*x2-3.*x1*x2+x1s*x2-x2s-r1s*x2s+x1*x2s) /prop12 +(1.-r1q+6.*r1s*r2s-r2q-x1+r1s*x1-r2s*x1+x2-9.*r1s*x2 +3.*r2s*x2+r1s*x1*x2-r2s*x1*x2-3.*x2s+3.*r1s*x2s-r2s*x2s +x1*x2s+x2c) /prop2s; rFO4 = rFO4/2.; isSet4 = true; } break; // q -> q V. case 3: if (combi == 1 || combi == 3) { rLO1 = ps*(1.-2.*r1s+r1q+r2s-6.*r1*r2s+r1s*r2s-2.*r2q); rFO1 = -2.*(-1.+r1-2.*r1s+2.*r1c-r1q+pow5(r1)-r2s+r1*r2s -5.*r1s*r2s+r1c*r2s-2.*r1*r2q+2.*x1-2.*r1*x1+2.*r1s*x1 -2.*r1c*x1+2.*r2s*x1+5.*r1*r2s*x1+r1s*r2s*x1+2.*r2q*x1 -x1s+r1*x1s-r2s*x1s+3.*x2+4.*r1s*x2+r1q*x2+2.*r2s*x2 +2.*r1s*r2s*x2-4.*x1*x2-2.*r1s*x1*x2-r2s*x1*x2+x1s*x2 -2.*x2s-2.*r1s*x2s+x1*x2s) /prop23 +(2.*r2s+6.*r1*r2s-6.*r1s*r2s+6.*r1c*r2s+2.*r2q+6.*r1*r2q -r2s*x1+r1s*r2s*x1-r2q*x1+x2-r1q*x2-3.*r2s*x2-6.*r1*r2s*x2 +9.*r1s*r2s*x2-2.*r2q*x2-x1*x2+r1s*x1*x2-x2s-3.*r1s*x2s +2.*r2s*x2s+x1*x2s) /prop2s +(-4.-8.*r1s-4.*r1q+4.*r2s-4.*r1s*r2s+8.*r2q+9.*x1+10.*r1s*x1 +r1q*x1-3.*r2s*x1+6.*r1*r2s*x1+r1s*r2s*x1-2.*r2q*x1-6.*x1s- 2.*r1s*x1s+x1c+7.*x2+8.*r1s*x2+r1q*x2-7.*r2s*x2+6.*r1*r2s*x2 +r1s*r2s*x2-2.*r2q*x2-9.*x1*x2-3.*r1s*x1*x2+2.*r2s*x1*x2 +2.*x1s*x2-3.*x2s-r1s*x2s+2.*r2s*x2s+x1*x2s) /x3s; isSet1 = true; } if (combi == 2 || combi == 3) { rLO2 = ps*(1.-2.*r1s+r1q+r2s+6.*r1*r2s+r1s*r2s-2.*r2q); rFO2 = 2*(1.+r1+2.*r1s+2.*r1c+r1q+pow5(r1)+r2s+r1*r2s +5.*r1s*r2s+r1c*r2s-2.*r1*r2q-2.*x1-2.*r1*x1-2.*r1s*x1 -2.*r1c*x1-2.*r2s*x1+5.*r1*r2s*x1-r1s*r2s*x1-2.*r2q*x1+x1s +r1*x1s+r2s*x1s-3.*x2-4.*r1s*x2-r1q*x2-2.*r2s*x2 -2.*r1s*r2s*x2+4.*x1*x2+2.*r1s*x1*x2+r2s*x1*x2-x1s*x2 +2.*x2s+2.*r1s*x2s-x1*x2s) /prop23 +(2.*r2s-6.*r1*r2s-6.*r1s*r2s-6.*r1c*r2s+2.*r2q-6.*r1*r2q -r2s*x1+r1s*r2s*x1-r2q*x1+x2-r1q*x2-3.*r2s*x2+6.*r1*r2s*x2 +9.*r1s*r2s*x2-2.*r2q*x2-x1*x2+r1s*x1*x2-x2s-3.*r1s*x2s +2.*r2s*x2s+x1*x2s) /prop2s +(-4.-8.*r1s-4.*r1q+4.*r2s-4.*r1s*r2s+8.*r2q+9.*x1+10.*r1s*x1 +r1q*x1-3.*r2s*x1-6.*r1*r2s*x1+r1s*r2s*x1-2.*r2q*x1-6.*x1s -2.*r1s*x1s+x1c+7.*x2+8.*r1s*x2+r1q*x2-7.*r2s*x2-6.*r1*r2s*x2 +r1s*r2s*x2-2.*r2q*x2-9.*x1*x2-3.*r1s*x1*x2+2.*r2s*x1*x2 +2.*x1s*x2-3.*x2s-r1s*x2s+2.*r2s*x2s+x1*x2s) /x3s; isSet2 = true; } if (combi == 4) { rLO4 = ps*(1.-2.*r1s+r1q+r2s+r1s*r2s-2.*r2q); rFO4 = 2*(1.+2.*r1s+r1q+r2s+5.*r1s*r2s-2.*x1-2.*r1s*x1 -2.*r2s*x1-r1s*r2s*x1-2.*r2q*x1+x1s+r2s*x1s-3.*x2-4.*r1s*x2 -r1q*x2-2.*r2s*x2-2.*r1s*r2s*x2+4.*x1*x2+2.*r1s*x1*x2+r2s*x1*x2 -x1s*x2+2.*x2s+2.*r1s*x2s-x1*x2s) /prop23 +(2.*r2s-6.*r1s*r2s+2.*r2q-r2s*x1+r1s*r2s*x1-r2q*x1+x2-r1q*x2 -3.*r2s*x2+9.*r1s*r2s*x2-2.*r2q*x2-x1*x2+r1s*x1*x2-x2s-3.*r1s*x2s +2.*r2s*x2s+x1*x2s) /prop2s +(-4.-8.*r1s-4.*r1q+4.*r2s-4.*r1s*r2s+8.*r2q+9.*x1+10.*r1s*x1 +r1q*x1-3.*r2s*x1+r1s*r2s*x1-2.*r2q*x1-6.*x1s-2.*r1s*x1s+x1c +7.*x2+8.*r1s*x2+r1q*x2-7.*r2s*x2+r1s*r2s*x2-2.*r2q*x2-9.*x1*x2 -3.*r1s*x1*x2+2.*r2s*x1*x2+2.*x1s*x2-3.*x2s-r1s*x2s+2.*r2s*x2s +x1*x2s) /x3s; isSet4 = true; } break; // S -> q qbar (S = h0/H0/A0/H+-/...). case 4: if (combi == 1 || combi == 3) { rLO1 = ps*(1.-r1s-r2s-2.*r1*r2); rFO1 = -(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q+x1 -r1s*x1+2.*r1*r2*x1+3.*r2s*x1+x2+r1s*x2-r2s*x2-x1*x2) /prop1s -2.*(r1s+r1q-2.*r1c*r2+r2s-6.*r1s*r2s-2.*r1*r2c+r2q-r1s*x1 +r1*r2*x1+2.*r2s*x1+2.*r1s*x2+r1*r2*x2-r2s*x2-x1*x2) /prop12 -(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q+x1-r1s*x1 +r2s*x1+x2+3.*r1s*x2+2.*r1*r2*x2-r2s*x2-x1*x2) /prop2s; isSet1 = true; } if (combi == 2 || combi == 3) { rLO2 = ps*(1.-r1s-r2s+2.*r1*r2); rFO2 = -(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1 -r1s*x1-2.*r1*r2*x1+3.*r2s*x1+x2+r1s*x2-r2s*x2-x1*x2) /prop1s -(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1 -r1s*x1+r2s*x1+x2+3.*r1s*x2-2.*r1*r2*x2-r2s*x2-x1*x2) /prop2s +2.*(-r1s-r1q-2.*r1c*r2-r2s+6.*r1s*r2s-2.*r1*r2c-r2q+r1s*x1 +r1*r2*x1-2.*r2s*x1-2.*r1s*x2+r1*r2*x2+r2s*x2+x1*x2) /prop12; isSet2 = true; } if (combi == 4) { rLO4 = ps*(1.-r1s-r2s); rFO4 = -(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+3.*r2s*x1+x2 +r1s*x2-r2s*x2-x1*x2) /prop1s -2.*(r1s+r1q+r2s-6.*r1s*r2s+r2q-r1s*x1 +2.*r2s*x1+2.*r1s*x2-r2s*x2-x1*x2) /prop12 -(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+r2s*x1 +x2+3.*r1s*x2-r2s*x2-x1*x2) /prop2s; isSet4 = true; } break; // q -> q S. case 5: if (combi == 1 || combi == 3) { rLO1 = ps*(1.+r1s-r2s+2.*r1); rFO1 = (4.-4.*r1s+4.*r2s-3.*x1-2.*r1*x1+r1s*x1-r2s*x1-5.*x2 -2.*r1*x2+r1s*x2-r2s*x2+x1*x2+x2s) /x3s -2.*(3.-r1-5.*r1s-r1c+3.*r2s+r1*r2s-2.*x1-r1*x1 +r1s*x1-4.*x2+2.*r1s*x2-r2s*x2+x1*x2+x2s) /prop23 +(2.-2.*r1-6.*r1s-2.*r1c+2.*r2s-2.*r1*r2s-x1+r1s*x1 -r2s*x1-3.*x2+2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s) /prop2s; isSet1 = true; } if (combi == 2 || combi == 3) { rLO2 = ps*(1.+r1s-r2s-2.*r1); rFO2 = (4.-4.*r1s+4.*r2s-3.*x1+2.*r1*x1+r1s*x1-r2s*x1-5.*x2 +2.*r1*x2+r1s*x2-r2s*x2+x1*x2+x2s) /x3s -2.*(3.+r1-5.*r1s+r1c+3.*r2s-r1*r2s-2.*x1+r1*x1 +r1s*x1-4.*x2+2.*r1s*x2-r2s*x2+x1*x2+x2s) /prop23 +(2.+2.*r1-6.*r1s+2.*r1c+2.*r2s+2.*r1*r2s-x1+r1s*x1 -r2s*x1-3.*x2-2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s) /prop2s; isSet2 = true; } if (combi == 4) { rLO4 = ps*(1.+r1s-r2s); rFO4 = (4.-4.*r1s+4.*r2s-3.*x1+r1s*x1-r2s*x1-5.*x2+r1s*x2 -r2s*x2+x1*x2+x2s) /x3s -2.*(3.-5.*r1s+3.*r2s-2.*x1+r1s*x1-4.*x2+2.*r1s*x2 -r2s*x2+x1*x2+x2s) /prop23 +(2.-6.*r1s+2.*r2s-x1+r1s*x1-r2s*x1-3.*x2+3.*r1s*x2 -r2s*x2+x1*x2+x2s) /prop2s; isSet4 = true; } break; // V -> ~q ~qbar (~q = squark). case 6: rLO1 = ps*(1.-2.*r1s+r1q-2.*r2s-2.*r1s*r2s+r2q); rFO1 = 2.*3.+(1.+r1s+r2s-x1)*(4.*r1s-x1s) /prop1s +2.*(-1.-3.*r1s-r2s+x1+x1s*0.5+x2-x1*x2*0.5) /prop1 +(1.+r1s+r2s-x2)*(4.*r2s-x2s) /prop2s +2.*(-1.-r1s-3.*r2s+x1+x2-x1*x2*0.5+x2s*0.5) /prop2 -(-4.*r1s-4.*r1q-4.*r2s-8.*r1s*r2s-4.*r2q+2.*x1+6.*r1s*x1 +6.*r2s*x1-2.*x1s+2.*x2+6.*r1s*x2+6.*r2s*x2-4.*x1*x2 -2.*r1s*x1*x2-2.*r2s*x1*x2+x1s*x2-2.*x2s+x1*x2s) /prop12; isSet1 = true; break; // ~q -> ~q V. case 7: rLO1 = ps*(1.-2.*r1s+r1q-2.*r2s-2.*r1s*r2s+r2q); rFO1 = 16.*r2s-8.*(4.*r2s+2.*r2s*x1+x2+r1s*x2+r2s*x2-x1*x2 -2.*x2s) /(3.*prop2) +8.*(1.+r1s+r2s-x2)*(4.*r2s-x2s) /(3.*prop2s) +8.*(x1+x2)*(-1.-2.*r1s-r1q-2.*r2s+2.*r1s*r2s-r2q+2.*x1 +2.*r1s*x1+2.*r2s*x1-x1s+2.*x2+2.*r1s*x2+2.*r2s*x2-2.*x1*x2-x2s) /(3.*x3s) +8.*(-1.-r1s+r2s-x1)*(2.*r2s*x1+x2+r1s*x2+r2s*x2-x1*x2-x2s) /(3.*prop2*x3) -8.*(1.+2.*r1s+r1q+2.*r2s-2.*r1s*r2s+r2q-2.*x1-2.*r1s*x1 -4.*r2s*x1+x1s-3.*x2-3.*r1s*x2-3.*r2s*x2+3.*x1*x2+2.*x2s) /(3.*x3); rFO1 = 3.*rFO1/8.; isSet1 = true; break; // S -> ~q ~qbar. case 8: rLO1 = ps; rFO1 = (-1.-2.*r1s-r1q-2.*r2s+2.*r1s*r2s-r2q+2.*x1+2.*r1s*x1 +2.*r2s*x1-x1s-r2s*x1s+2.*x2+2.*r1s*x2+2.*r2s*x2-3.*x1*x2 -r1s*x1*x2-r2s*x1*x2+x1s*x2-x2s-r1s*x2s+x1*x2s) /(prop1s*prop2s); rFO1 = 2.*rFO1; isSet1 = true; break; // ~q -> ~q S. case 9: rLO1 = ps; rFO1 = (-1.-r1s-r2s+x2) /prop2s +(1.+r1s-r2s+x1) /prop23 -(x1+x2) /x3s; isSet1 = true; break; // chi -> q ~qbar (chi = neutralino/chargino). case 10: if (combi == 1 || combi == 3) { rLO1 = ps*(1.+r1s-r2s+2.*r1); rFO1 = (2.*r1+x1)*(-1.-r1s-r2s+x1) /prop1s +2.*(-1.-r1s-2.*r1c-r2s-2.*r1*r2s+3.*x1*0.5+r1*x1 -r1s*x1*0.5-r2s*x1*0.5+x2+r1*x2+r1s*x2-x1*x2*0.5) /prop12 +(2.-2.*r1-6.*r1s-2.*r1c+2.*r2s-2.*r1*r2s-x1+r1s*x1 -r2s*x1-3.*x2+2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s) /prop2s; isSet1 = true; } if (combi == 2 || combi == 3) { rLO2 = ps*(1.-2.*r1+r1s-r2s); rFO2 = (2.*r1-x1)*(1.+r1s+r2s-x1) /prop1s +2.*(-1.-r1s+2.*r1c-r2s+2.*r1*r2s+3.*x1*0.5-r1*x1 -r1s*x1*0.5-r2s*x1*0.5+x2-r1*x2+r1s*x2-x1*x2*0.5) /prop12 +(2.+2.*r1-6.*r1s+2.*r1c+2.*r2s+2.*r1*r2s-x1+r1s*x1 -r2s*x1-3.*x2-2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s)/ prop2s; isSet2 = true; } if (combi == 4) { rLO4 = ps*(1.+r1s-r2s); rFO4 = x1*(-1.-r1s-r2s+x1) /prop1s +2.*(-1.-r1s-r2s+3.*x1*0.5-r1s*x1*0.5-r2s*x1*0.5 +x2+r1s*x2-x1*x2*0.5) /prop12 +(2.-6.*r1s+2.*r2s-x1+r1s*x1-r2s*x1-3.*x2+3.*r1s*x2 -r2s*x2+x1*x2+x2s) /prop2s; isSet4 = true; } break; // ~q -> q chi. case 11: if (combi == 1 || combi == 3) { rLO1 = ps*(1.-pow2(r1+r2)); rFO1 = (1.+r1s+2.*r1*r2+r2s-x1-x2)*(x1+x2) /x3s -(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q+x1 -r1s*x1+r2s*x1+x2+3.*r1s*x2+2.*r1*r2*x2-r2s*x2-x1*x2) /prop2s +(-1.-2.*r1s-r1q-2.*r1*r2-2.*r1c*r2+2.*r1*r2c+r2q+x1+r1s*x1 -2.*r1*r2*x1-3.*r2s*x1+2.*r1s*x2-2.*r2s*x2+x1*x2+x2s) /prop23; isSet1 = true; } if (combi == 2 || combi == 3) { rLO2 = ps*(1.-pow2(r1-r2)); rFO2 = (1.+r1s-2.*r1*r2+r2s-x1-x2)*(x1+x2) /x3s -(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1 -r1s*x1+r2s*x1+x2+3.*r1s*x2-2.*r1*r2*x2-r2s*x2-x1*x2) /prop2s +(-1.-2.*r1s-r1q+2.*r1*r2+2.*r1c*r2-2.*r1*r2c+r2q+x1+r1s*x1 +2.*r1*r2*x1-3.*r2s*x1+2.*r1s*x2-2.*r2s*x2+x1*x2+x2s) /prop23; isSet2 = true; } if (combi == 4) { rLO4 = ps*(1.-r1s-r2s); rFO4 = (1.+r1s+r2s-x1-x2)*(x1+x2) /x3s -(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+r2s*x1+x2 +3.*r1s*x2-r2s*x2-x1*x2) /prop2s +(-1.-2.*r1s-r1q+r2q+x1+r1s*x1-3.*r2s*x1 +2.*r1s*x2-2.*r2s*x2+x1*x2+x2s) /prop23; isSet4 = true; } break; // q -> ~q chi. case 12: if (combi == 1 || combi == 3) { rLO1 = ps*(1.-r1s+r2s+2.*r2); rFO1 = (2.*r2+x2)*(-1.-r1s-r2s+x2) /prop2s +(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1-2.*r2*x1+r2s*x1+x1s -3.*x2-r1s*x2-2.*r2*x2+r2s*x2+x1*x2) /x3s +2.*(-1.-r1s+r2+r1s*r2-r2s-r2c+x1+r2*x1+r2s*x1+2.*x2 +r1s*x2-x1*x2*0.5-x2s*0.5) /prop23; isSet1 = true; } if (combi == 2 || combi == 3) { rLO2 = ps*(1.-r1s+r2s-2.*r2); rFO2 = (2.*r2-x2)*(1.+r1s+r2s-x2) /prop2s +(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1+2.*r2*x1+r2s*x1+x1s -3.*x2-r1s*x2+2.*r2*x2+r2s*x2+x1*x2) /x3s +2.*(-1.-r1s-r2-r1s*r2-r2s+r2c+x1-r2*x1+r2s*x1+2.*x2 +r1s*x2-x1*x2*0.5-x2s*0.5) /prop23; isSet2 = true; } if (combi == 4) { rLO4 = ps*(1.-r1s+r2s); rFO4 = x2*(-1.-r1s-r2s+x2) /prop2s +(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1+r2s*x1+x1s -3.*x2-r1s*x2+r2s*x2+x1*x2) /x3s +2.*(-1.-r1s-r2s+x1+r2s*x1+2.*x2 +r1s*x2-x1*x2*0.5-x2s*0.5) /prop23; isSet4 = true; } break; // ~g -> q ~qbar. case 13: if (combi == 1 || combi == 3) { rLO1 = ps*(1.+r1s-r2s+2.*r1); rFO1 = 4.*(2.*r1+x1)*(-1.-r1s-r2s+x1) /(3.*prop1s) -(-1.-r1s-2.*r1c-r2s-2.*r1*r2s+3.*x1*0.5+r1*x1-r1s*x1*0.5 -r2s*x1*0.5+x2+r1*x2+r1s*x2-x1*x2*0.5) /(3.*prop12) +3.*(-1.+r1-r1s-r1c-r2s+r1*r2s+2.*x1+r2s*x1-x1s*0.5+x2+r1*x2 +r1s*x2-x1*x2*0.5) /prop13 +3.*(4.-4.*r1s+4.*r2s-3.*x1-2.*r1*x1+r1s*x1-r2s*x1-5.*x2 -2.*r1*x2+r1s*x2-r2s*x2+x1*x2+x2s) /x3s -3.*(3.-r1-5.*r1s-r1c+3.*r2s+r1*r2s-2.*x1-r1*x1+r1s*x1 -4.*x2+2.*r1s*x2-r2s*x2+x1*x2+x2s) /prop23 +4.*(2.-2.*r1-6.*r1s-2.*r1c+2.*r2s-2.*r1*r2s-x1+r1s*x1-r2s*x1 -3.*x2+2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s) /(3.*prop2s); rFO1 = 3.*rFO1/4.; isSet1 = true; } if (combi == 2 || combi == 3) { rLO2 = ps*(1.+r1s-r2s-2.*r1); rFO2 = 4.*(2.*r1-x1)*(1.+r1s+r2s-x1) /(3.*prop1s) +3.*(-1.-r1-r1s+r1c-r2s-r1*r2s+2.*x1+r2s*x1-x1s*0.5 +x2-r1*x2+r1s*x2-x1*x2*0.5) /prop13 +(2.+2.*r1s-4.*r1c+2.*r2s-4.*r1*r2s-3.*x1+2.*r1*x1 +r1s*x1+r2s*x1-2.*x2+2.*r1*x2-2.*r1s*x2+x1*x2) /(6.*prop12) +3.*(4.-4.*r1s+4.*r2s-3.*x1+2.*r1*x1+r1s*x1-r2s*x1-5.*x2 +2.*r1*x2+r1s*x2-r2s*x2+x1*x2+x2s) /x3s -3.*(3.+r1-5.*r1s+r1c+3.*r2s-r1*r2s-2.*x1+r1*x1+r1s*x1-4.*x2 +2.*r1s*x2-r2s*x2+x1*x2+x2s) /prop23 +4.*(2.+2.*r1-6.*r1s+2.*r1c+2.*r2s+2.*r1*r2s-x1+r1s*x1-r2s*x1 -3.*x2-2.*r1*x2+3.*r1s*x2-r2s*x2+x1*x2+x2s) /(3.*prop2s); rFO2 = 3.*rFO2/4.; isSet2 = true; } if (combi == 4) { rLO4 = ps*(1.+r1s-r2s); rFO4 = 8.*x1*(-1.-r1s-r2s+x1) /(3.*prop1s) +6.*(-1-r1s-r2s+2.*x1+r2s*x1-x1s*0.5+x2+r1s*x2-x1*x2*0.5) /prop13 +(2.+2.*r1s+2.*r2s-3.*x1+r1s*x1+r2s*x1-2.*x2-2.*r1s*x2+x1*x2) /(3.*prop12) +6.*(4.-4.*r1s+4.*r2s-3.*x1+r1s*x1-r2s*x1-5.*x2+r1s*x2-r2s*x2 +x1*x2+x2s) /x3s -6.*(3.-5.*r1s+3.*r2s-2.*x1+r1s*x1-4.*x2+2.*r1s*x2-r2s*x2+x1*x2+x2s) /prop23 +8.*(2.-6.*r1s+2.*r2s-x1+r1s*x1-r2s*x1-3.*x2+3.*r1s*x2-r2s*x2 +x1*x2+x2s) /(3.*prop2s); rFO4 = 3.*rFO4/8.; isSet4 = true; } break; // ~q -> q ~g. case 14: if (combi == 1 || combi == 3) { rLO1 = ps*(1.-r1s-r2s-2.*r1*r2); rFO1 = 64.*(1.+r1s+2.*r1*r2+r2s-x1-x2)*(x1+x2) /(9.*x3s) -16.*(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q +x1-r1s*x1+2.*r1*r2*x1+3.*r2s*x1+x2+r1s*x2-r2s*x2-x1*x2) /prop1s -16.*(r1s+r1q-2.*r1c*r2+r2s-6.*r1s*r2s-2.*r1*r2c+r2q-r1s*x1 +r1*r2*x1+2.*r2s*x1+2.*r1s*x2+r1*r2*x2-r2s*x2-x1*x2) /prop12 -64.*(-1.+r1q-2.*r1*r2-2.*r1c*r2-6.*r1s*r2s-2.*r1*r2c+r2q+x1 -r1s*x1+r2s*x1+x2+3.*r1s*x2+2.*r1*r2*x2-r2s*x2-x1*x2) /(9.*prop2s) +8.*(-1.+r1q-2.*r1*r2+2.*r1c*r2-2.*r2s-2.*r1*r2c-r2q-2.*r1s*x1 +2.*r2s*x1+x1s+x2-3.*r1s*x2-2.*r1*r2*x2+r2s*x2+x1*x2) /prop13 -8.*(-1.-2.*r1s-r1q-2.*r1*r2-2.*r1c*r2+2.*r1*r2c+r2q+x1+r1s*x1 -2.*r1*r2*x1-3.*r2s*x1+2.*r1s*x2-2.*r2s*x2+x1*x2+x2s) /(9.*prop23); rFO1 = 9.*rFO1/64.; isSet1 = true; } if (combi == 2 || combi == 3) { rLO2 = ps*(1.-r1s-r2s+2.*r1*r2); rFO2 = 64.*(1.+r1s-2.*r1*r2+r2s-x1-x2)*(x1+x2) /(9.*x3s) -16.*(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1 -r1s*x1-2.*r1*r2*x1+3.*r2s*x1+x2+r1s*x2-r2s*x2-x1*x2) /prop1s -64.*(-1.+r1q+2.*r1*r2+2.*r1c*r2-6.*r1s*r2s+2.*r1*r2c+r2q+x1 -r1s*x1+r2s*x1+x2+3.*r1s*x2-2.*r1*r2*x2-r2s*x2-x1*x2) /(9.*prop2s) +16.*(-r1s-r1q-2.*r1c*r2-r2s+6.*r1s*r2s-2.*r1*r2c-r2q+r1s*x1 +r1*r2*x1-2.*r2s*x1-2.*r1s*x2+r1*r2*x2+r2s*x2+x1*x2) /prop12 +8.*(-1.+r1q+2.*r1*r2-2.*r1c*r2-2.*r2s+2.*r1*r2c-r2q-2.*r1s*x1 +2.*r2s*x1+x1s+x2-3.*r1s*x2+2.*r1*r2*x2+r2s*x2+x1*x2) /prop13 -8.*(-1.-2.*r1s-r1q+2.*r1*r2+2.*r1c*r2-2.*r1*r2c+r2q+x1+r1s*x1+ 2.*r1*r2*x1-3.*r2s*x1+2.*r1s*x2-2.*r2s*x2+x1*x2+x2s) /(9.*prop23); rFO2 = 9.*rFO2/64.; isSet2 = true; } if (combi == 4) { rLO4 = ps*(1.-r1s-r2s); rFO4 = 128.*(1.+r1s+r2s-x1-x2)*(x1+x2) /(9.*x3s) -32*(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+3.*r2s*x1+x2 +r1s*x2-r2s*x2-x1*x2) /prop1s -32.*(r1s+r1q+r2s-6.*r1s*r2s+r2q-r1s*x1+2.*r2s*x1+2.*r1s*x2 -r2s*x2-x1*x2) /prop12 -128.*(-1.+r1q-6.*r1s*r2s+r2q+x1-r1s*x1+r2s*x1+x2+3.*r1s*x2 -r2s*x2-x1*x2) /(9.*prop2s) +16.*(-1.+r1q-2.*r2s-r2q-2.*r1s*x1+2.*r2s*x1+x1s +x2-3.*r1s*x2+r2s*x2+x1*x2) /prop13 -16.*(-1.-2.*r1s-r1q+r2q+x1+r1s*x1-3.*r2s*x1 +2.*r1s*x2-2.*r2s*x2+x1*x2+x2s) /(9.*prop23); rFO4 = 9.*rFO4/128.; isSet4 = true; } break; // q -> ~q ~g. case 15: if (combi == 1 || combi == 3) { rLO1 = ps*(1.-r1s+r2s+2.*r2); rFO1 = 32*(2.*r2+x2)*(-1.-r1s-r2s+x2) /(9.*prop2s) +8.*(-1.-r1s-2.*r1s*r2-r2s-2.*r2c+x1+r2*x1+r2s*x1 +3.*x2*0.5-r1s*x2*0.5+r2*x2-r2s*x2*0.5-x1*x2*0.5) /prop12 +8.*(2.+2.*r1s-2.*r2-2.*r1s*r2-6.*r2s-2.*r2c-3.*x1-r1s*x1 +2.*r2*x1+3.*r2s*x1+x1s-x2-r1s*x2+r2s*x2+x1*x2) /prop1s +32.*(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1-2.*r2*x1+r2s*x1+x1s -3.*x2-r1s*x2-2.*r2*x2+r2s*x2+x1*x2) /(9.*x3s) -8.*(3.+3.*r1s-r2+r1s*r2-5.*r2s-r2c-4.*x1-r1s*x1 +2.*r2s*x1+x1s-2.*x2-r2*x2+r2s*x2+x1*x2) /prop13 -8.*(-1.-r1s+r2+r1s*r2-r2s-r2c+x1+r2*x1+r2s*x1+2.*x2+r1s*x2 -x1*x2*0.5-x2s*0.5) /(9.*prop23); rFO1 = 9.*rFO1/32.; isSet1 = true; } if (combi == 2 || combi == 3) { rLO2 = ps*(1.-r1s+r2s-2.*r2); rFO2 = 32*(2.*r2-x2)*(1.+r1s+r2s-x2) /(9.*prop2s) +8.*(-1.-r1s+2.*r1s*r2-r2s+2.*r2c+x1-r2*x1+r2s*x1 +3.*x2*0.5-r1s*x2*0.5-r2*x2-r2s*x2*0.5-x1*x2*0.5) /prop12 +8.*(2.+2.*r1s+2.*r2+2.*r1s*r2-6.*r2s+2.*r2c-3.*x1-r1s*x1 -2.*r2*x1+3.*r2s*x1+x1s-x2-r1s*x2+r2s*x2+x1*x2) /prop1s -8.*(3.+3.*r1s+r2-r1s*r2-5.*r2s+r2c-4.*x1-r1s*x1+2.*r2s*x1+x1s -2.*x2+r2*x2+r2s*x2+x1*x2) /prop13 +32*(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1+2.*r2*x1+r2s*x1 +x1s-3.*x2-r1s*x2+2.*r2*x2+r2s*x2+x1*x2) /(9.*x3s) -8.*(-1.-r1s-r2-r1s*r2-r2s+r2c+x1-r2*x1+r2s*x1+2.*x2+r1s*x2 -x1*x2*0.5-x2s*0.5) /(9.*prop23); rFO2 = 9.*rFO2/32.; isSet2 = true; } if (combi == 4) { rLO4 = ps*(1.-r1s+r2s); rFO4 = 64.*x2*(-1.-r1s-r2s+x2) /(9.*prop2s) +16.*(-1.-r1s-r2s+x1+r2s*x1+3.*x2*0.5-r1s*x2*0.5 -r2s*x2*0.5-x1*x2*0.5) /prop12 -16.*(3.+3.*r1s-5.*r2s-4.*x1-r1s*x1+2.*r2s*x1+x1s-2.*x2+r2s*x2 +x1*x2) /prop13 +64.*(4.+4.*r1s-4.*r2s-5.*x1-r1s*x1+r2s*x1+x1s-3.*x2 -r1s*x2+r2s*x2+x1*x2) /(9.*x3s) +16.*(2.+2.*r1s-6.*r2s-3.*x1-r1s*x1+3.*r2s*x1+x1s -x2-r1s*x2+r2s*x2+x1*x2) /prop1s -16.*(-1.-r1s-r2s+x1+r2s*x1+2.*x2+r1s*x2-x1*x2*0.5-x2s*0.5) /(9.*prop23); rFO4 = 9.*rFO4/64.; isSet4 = true; } break; // g -> ~g ~g. Use (9/4)*eikonal. May be changed in the future. case 16: rLO = ps; if (combi == 2) offset = x3s; else if (combi == 3) offset = mix * x3s; else if (combi == 4) offset = 0.5 * x3s; rFO = ps * 4.5 * ( (x1+x2-1.+offset-r1s-r2s)/prop12 - r1s/prop2s - r2s/prop1s ); break; // Dv -> qv d. case 30: rLO = ps*(1.-r1s+r2s+2.*r2); rFO = ( 0.5*r3s + 2.*r1q + 0.5*r2s*r3s + r2*r3s - 2.*r1s - 0.5*r1s*r3s - 2.*r1s*r2s - 4.*r1s*r2 ) / prop2s + ( -2. + 2.*r2q + 2.*r1q + 2.*r2s*r3s - 4.*r2 + 2.*r2*r3s + 4.*r2*r2s - 4.*r1s*r2s - 4.*r1s*r2 ) /prop23 + ( -2. - 0.5*r3s - 2.*r2s - 4.*r2 + 2.*r1s ) / prop2 + ( -2. - r3s - 2.*r2s - r2s*r3s - 4.*r2 - 2.*r2*r3s + 2.*r1s + r1s*r3s ) / prop3s + ( -1. - r3s - r2s - 4.*r2 + r1s - x2 ) / prop3 + 1.; break; // S -> Dv Dvbar case 31: rLO = ps*(1.-4.*r1s); rFO = (r3s + 2.*r1s) * (-1. + 4.*r1s) * (1./prop1s + 1./prop2s) + (-1. + 8.*r1s - x2) / prop1 + (-1. + 8.*r1s - x1) / prop2 + 2. * (1. - 6.*r1s + 8.*r1q + 4.*r3s*r1s) / prop12 + 2.; break; // Eikonal expression for kind == 1; also acts as default. default: rLO = ps; if (combi == 2) offset = x3s; else if (combi == 3) offset = mix * x3s; else if (combi == 4) offset = 0.5 * x3s; rFO = ps * 2. * ( (x1+x2-1.+offset-r1s-r2s)/prop12 - r1s/prop2s - r2s/prop1s ); break; // End of ME cases. } // Find relevant leading and first order expressions. if (combi == 1 && isSet1) { rLO = rLO1; rFO = rFO1; } else if (combi == 2 && isSet2) { rLO = rLO2; rFO = rFO2; } else if (combi == 3 && isSet1 && isSet2) { rLO = mix * rLO1 + (1.-mix) * rLO2; rFO = mix * rFO1 + (1.-mix) * rFO2; } else if (isSet4) { rLO = rLO4; rFO = rFO4; } else if (combi == 4 && isSet1 && isSet2) { rLO = 0.5 * (rLO1 + rLO2); rFO = 0.5 * (rFO1 + rFO2); } else if (isSet1) { rLO = rLO1; rFO = rFO1; } // Return ratio of first to leading order cross section. return rFO / rLO; } //-------------------------------------------------------------------------- // Find coefficient of azimuthal asymmetry from gluon polarization. void TimeShower::findAsymPol( Event& event, TimeDipoleEnd* dip) { // Default is no asymmetry. Only gluons are studied. dip->asymPol = 0.; dip->iAunt = 0; int iRad = dip->iRadiator; if (!doPhiPolAsym || event[iRad].id() != 21) return; // Trace grandmother via possibly intermediate recoil copies. int iMother = event.iTopCopy(iRad); int iGrandM = event[iMother].mother1(); // If grandmother in initial state of hard scattering, // then only keep gg and qq initial states. int statusGrandM = event[iGrandM].status(); bool isHardProc = (statusGrandM == -21 || statusGrandM == -31); if (isHardProc) { if (event[iGrandM + 1].status() != statusGrandM) return; if (event[iGrandM].isGluon() && event[iGrandM + 1].isGluon()); else if (event[iGrandM].isQuark() && event[iGrandM + 1].isQuark()); else return; } // Set aunt by history or, for hard scattering, by colour flow. if (isHardProc) dip->iAunt = dip->iRecoiler; else dip->iAunt = (event[iGrandM].daughter1() == iMother) ? event[iGrandM].daughter2() : event[iGrandM].daughter1(); // Coefficient from gluon production (approximate z by energy). // For hard process arbitrarily put z = 1/2. double zProd = (isHardProc) ? 0.5 : event[iRad].e() / (event[iRad].e() + event[dip->iAunt].e()); if (event[iGrandM].isGluon()) dip->asymPol = pow2( (1. - zProd) / (1. - zProd * (1. - zProd) ) ); else dip->asymPol = 2. * (1. - zProd) / (1. + pow2(1. - zProd) ); // Coefficients from gluon decay. if (dip->flavour == 21) dip->asymPol *= pow2( (1. - dip->z) / (1. - dip->z * (1. - dip->z) ) ); else dip->asymPol *= -2. * dip->z *( 1. - dip->z ) / (1. - 2. * dip->z * (1. - dip->z) ); } //-------------------------------------------------------------------------- // Print the list of dipoles. void TimeShower::list(ostream& os) const { // Header. os << "\n -------- PYTHIA TimeShower Dipole Listing ----------------" << "--------------------------------------------- \n \n i rad" << " rec pTmax col chg gam oni hv isr sys sysR typ" << "e MErec mix ord spl ~gR \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(7) << dipEnd[i].iRadiator << setw(7) << dipEnd[i].iRecoiler << setw(12) << dipEnd[i].pTmax << setw(5) << dipEnd[i].colType << setw(5) << dipEnd[i].chgType << setw(5) << dipEnd[i].gamType << setw(5) << dipEnd[i].isOctetOnium << setw(5) << dipEnd[i].isHiddenValley << setw(5) << dipEnd[i].isrType << setw(5) << dipEnd[i].system << setw(5) << dipEnd[i].systemRec << setw(5) << dipEnd[i].MEtype << setw(7) << dipEnd[i].iMEpartner << setw(8) << dipEnd[i].MEmix << setw(5) << dipEnd[i].MEorder << setw(5) << dipEnd[i].MEsplit << setw(5) << dipEnd[i].MEgluinoRec << "\n"; // Done. os << "\n -------- End PYTHIA TimeShower Dipole Listing ------------" << "---------------------------------------------" << endl; } //========================================================================== } // end namespace Pythia8