// BeamRemnants.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 // BeamRemnants class. #include "BeamRemnants.h" namespace Pythia8 { //========================================================================== // The BeamDipole class is purely internal to reconnectColours. class BeamDipole { public: // Constructor. BeamDipole( int colIn = 0, int iColIn = 0, int iAcolIn = 0) : col(colIn), iCol(iColIn), iAcol(iAcolIn) {} // Members. int col, iCol, iAcol; double p1p2; }; //========================================================================== // The BeamRemnants class. //-------------------------------------------------------------------------- // Constants: could be changed here if desired, but normally should not. // These are of technical nature, as described for each. // If same (anti)colour appears twice in final state, repair or reject. const bool BeamRemnants::ALLOWCOLOURTWICE = true; // Maximum number of tries to match colours and kinematics in the event. const int BeamRemnants::NTRYCOLMATCH = 10; const int BeamRemnants::NTRYKINMATCH = 10; // Overall correction step for energy-momentum conservation; only // becomes relevant in rescattering scenarios when FSR dipole emissions // and primordial kT is added. Should hopefully not be needed. const bool BeamRemnants::CORRECTMISMATCH = false; //-------------------------------------------------------------------------- // Initialization. bool BeamRemnants::init( Info* infoPtrIn, Settings& settings, Rndm* rndmPtrIn, BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn, PartonSystems* partonSystemsPtrIn) { // Save pointers. infoPtr = infoPtrIn; rndmPtr = rndmPtrIn; beamAPtr = beamAPtrIn; beamBPtr = beamBPtrIn; partonSystemsPtr = partonSystemsPtrIn; // Width of primordial kT distribution. doPrimordialKT = settings.flag("BeamRemnants:primordialKT"); primordialKTsoft = settings.parm("BeamRemnants:primordialKTsoft"); primordialKThard = settings.parm("BeamRemnants:primordialKThard"); primordialKTremnant = settings.parm("BeamRemnants:primordialKTremnant"); halfScaleForKT = settings.parm("BeamRemnants:halfScaleForKT"); halfMassForKT = settings.parm("BeamRemnants:halfMassForKT"); // Handling of rescattering kinematics uncertainties from primodial kT. allowRescatter = settings.flag("MultipleInteractions:allowRescatter"); doRescatterRestoreY = settings.flag("BeamRemnants:rescatterRestoreY"); // Parameters for colour reconnection scenario, partly borrowed from // multiple interactions not to introduce too many new ones. doReconnect = settings.flag("BeamRemnants:reconnectColours"); reconnectRange = settings.parm("BeamRemnants:reconnectRange"); pT0Ref = settings.parm("MultipleInteractions:pT0Ref"); ecmRef = settings.parm("MultipleInteractions:ecmRef"); ecmPow = settings.parm("MultipleInteractions:ecmPow"); // Total and squared CM energy at nominal energy. eCM = infoPtr->eCM(); sCM = eCM * eCM; // The MI pT0 smoothening scale and its reconnection-strength combination. pT0 = pT0Ref * pow(eCM / ecmRef, ecmPow); pT20Rec = pow2(reconnectRange * pT0); // Done. return true; } //-------------------------------------------------------------------------- // Select the flavours/kinematics/colours of the two beam remnants. // Notation: iPar = all partons, iSys = matched systems of two beams, // iRem = additional partons in remnants. bool BeamRemnants::add( Event& event) { // Update to current CM energy. eCM = infoPtr->eCM(); sCM = eCM * eCM; // Check that flavour bookkept in event and in beam particles agree. for (int i = 0; i < beamAPtr->size(); ++i) { int j = (*beamAPtr)[i].iPos(); if ((*beamAPtr)[i].id() != event[j].id()) { infoPtr->errorMsg("Error in BeamRemnants::add: " "event and beam flavours do not match"); return false; } } for (int i = 0; i < beamBPtr->size(); ++i) { int j = (*beamBPtr)[i].iPos(); if ((*beamBPtr)[i].id() != event[j].id()) { infoPtr->errorMsg("Error in BeamRemnants::add: " "event and beam flavours do not match"); return false; } } // Number of scattering subsystems. Size of event record before treatment. nSys = partonSystemsPtr->sizeSys(); oldSize = event.size(); // Add required extra remnant flavour content. // Start all over if fails (in option where junctions not allowed). if ( !beamAPtr->remnantFlavours(event) || !beamBPtr->remnantFlavours(event) ) { infoPtr->errorMsg("Error in BeamRemnants::add:" " remnant flavour setup failed"); return false; } // Do the kinematics of the collision subsystems and two beam remnants. if (!setKinematics(event)) return false; // Allow colour reconnections. if (doReconnect) reconnectColours(event); // Save current modifiable colour configuration for fast restoration. vector colSave; vector acolSave; for (int i = oldSize; i < event.size(); ++i) { colSave.push_back( event[i].col() ); acolSave.push_back( event[i].acol() ); } event.saveJunctionSize(); // Allow several tries to match colours of initiators and remnants. // Frequent "failures" since shortcutting colours separately on // the two event sides may give "colour singlet gluons" etc. bool physical = true; for (int iTry = 0; iTry < NTRYCOLMATCH; ++iTry) { physical = true; // Reset list of colour "collapses" (transformations). colFrom.resize(0); colTo.resize(0); // First process each set of beam colours on its own. if (!beamAPtr->remnantColours(event, colFrom, colTo)) physical = false; if (!beamBPtr->remnantColours(event, colFrom, colTo)) physical = false; // Then check that colours and anticolours are matched in whole event. if ( physical && !checkColours(event) ) physical = false; // If no problems then done, else restore and loop. if (physical) break; for (int i = oldSize; i < event.size(); ++i) event[i].cols( colSave[i - oldSize], acolSave[i - oldSize] ); event.restoreJunctionSize(); infoPtr->errorMsg("Warning in BeamRemnants::add:" " colour tracing failed; will try again"); } // If no solution after several tries then failed. if (!physical) { infoPtr->errorMsg("Error in BeamRemnants::add:" " colour tracing failed after several attempts"); return false; } // Done. return true; } //-------------------------------------------------------------------------- // Set up trial transverse and longitudinal kinematics for each beam // separately. Final decisions involve comparing the two beams. bool BeamRemnants::setKinematics( Event& event) { // References to beams to simplify indexing. BeamParticle& beamA = *beamAPtr; BeamParticle& beamB = *beamBPtr; // Nothing to do for lepton-lepton scattering with all energy already used. if ( beamA.isUnresolvedLepton() && beamB.isUnresolvedLepton() ) return true; // Check that has not already used up beams. if ( (!beamA.isLepton() && beamA.xMax(-1) <= 0.) || (!beamB.isLepton() && beamB.xMax(-1) <= 0.) ) { infoPtr->errorMsg("Error in BeamRemnants::setKinematics:" " no momentum left for beam remnants"); return false; } // Last beam-status particles. Offset relative to normal beam locations. int nBeams = 3; for (int i = 3; i < event.size(); ++i) if (event[i].statusAbs() < 20) nBeams = i + 1; int nOffset = nBeams - 3; // Reserve space for extra information on the systems and beams. int nMaxBeam = max( beamA.size(), beamB.size() ); vector sHatSys(nMaxBeam); vector kTwidth(nMaxBeam); vector kTcomp(nMaxBeam); vector Msys(nSys); // Loop over all subsystems. Default values. Find invariant mass. double kTcompSumA = 0.; double kTcompSumB = 0.; for (int iSys = 0; iSys < nSys; ++iSys) { double kTwidthNow = 0.; double mHatDamp = 1.; int iInA = partonSystemsPtr->getInA(iSys); int iInB = partonSystemsPtr->getInB(iSys); double sHatNow = (event[iInA].p() + event[iInB].p()).m2Calc(); // Allow primordial kT reduction for small-mass and small-pT systems // (for hardest interaction pT -> renormalization scale so also 2 -> 1). if (doPrimordialKT) { double mHat = sqrt(sHatNow); mHatDamp = mHat / (mHat + halfMassForKT); double scale = (iSys == 0) ? infoPtr->QRen() : infoPtr->pTMI(iSys); kTwidthNow = ( (halfScaleForKT * primordialKTsoft + scale * primordialKThard) / (halfScaleForKT + scale) ) * mHatDamp; } // Store properties of compensation systems and total compensation power. // Rescattered partons do not compensate, but may be massive. sHatSys[iSys] = sHatNow; kTwidth[iSys] = kTwidthNow ; kTcomp[iSys] = mHatDamp; if (beamA[iSys].isFromBeam()) kTcompSumA += mHatDamp; else beamA[iSys].m( event[iInA].m() ); if (beamB[iSys].isFromBeam()) kTcompSumB += mHatDamp; else beamB[iSys].m( event[iInB].m() ); } // Primordial kT and compensation power among remnants. double kTwidthNow = (doPrimordialKT) ? primordialKTremnant : 0.; for (int iRem = nSys; iRem < nMaxBeam; ++iRem) { sHatSys[iRem] = 0.; kTwidth[iRem] = kTwidthNow ; kTcomp[iRem] = 1.; } kTcompSumA += beamA.size() - nSys; kTcompSumB += beamB.size() - nSys; // Allow ten tries to construct kinematics (but normally works first). bool physical; double xSum[2], xInvM[2], w2Beam[2], wPosRem, wNegRem, w2Rem; for (int iTry = 0; iTry < NTRYKINMATCH; ++iTry) { physical = true; // Loop over the two beams. for (int iBeam = 0; iBeam < 2; ++iBeam) { BeamParticle& beam = (iBeam == 0) ? beamA : beamB; int nPar = beam.size(); // Generate Gaussian pT for initiator partons inside hadrons. // Store/restore rescattered parton momenta before primordial kT. if (beam.isHadron() && doPrimordialKT) { double pxSum = 0.; double pySum = 0.; for (int iPar = 0; iPar < nPar; ++iPar) { if ( beam[iPar].isFromBeam() ) { pair gauss2 = rndmPtr->gauss2(); double px = kTwidth[iPar] * gauss2.first; double py = kTwidth[iPar] * gauss2.second; beam[iPar].px(px); beam[iPar].py(py); pxSum += px; pySum += py; } else { int iInAB = (iBeam == 0) ? partonSystemsPtr->getInA(iPar) : partonSystemsPtr->getInB(iPar); beam[iPar].p( event[iInAB].p() ); } } // Share recoil between all initiator partons, rescatterers excluded. double kTcompSum = (iBeam == 0) ? kTcompSumA : kTcompSumB; for (int iPar = 0; iPar < nPar; ++iPar) if (beam[iPar].isFromBeam() ) { beam[iPar].px( beam[iPar].px() - pxSum * kTcomp[iPar] / kTcompSum ); beam[iPar].py( beam[iPar].py() - pySum * kTcomp[iPar] / kTcompSum ); } // Without primordial kT: still need to store rescattered partons. } else if (beam.isHadron()) { for (int iPar = 0; iPar < nPar; ++iPar) if ( !beam[iPar].isFromBeam() ) { int iInAB = (iBeam == 0) ? partonSystemsPtr->getInA(iPar) : partonSystemsPtr->getInB(iPar); beam[iPar].p( event[iInAB].p() ); } } // Pick unrescaled x values for remnants. Sum up (unscaled) p+ and p-. xSum[iBeam] = 0.; xInvM[iBeam] = 0.; for (int iRem = nSys; iRem < nPar; ++iRem) { double xPrel = beam.xRemnant( iRem); beam[iRem].x(xPrel); xSum[iBeam] += xPrel; xInvM[iBeam] += beam[iRem].mT2()/xPrel; } // Squared transverse mass for each beam, using lightcone x. w2Beam[iBeam] = xSum[iBeam] * xInvM[iBeam]; // End separate treatment of the two beams. } // Recalculate kinematics of initiator systems with primordial kT. wPosRem = eCM; wNegRem = eCM; for (int iSys = 0; iSys < nSys; ++iSys) { int iA = beamA[iSys].iPos(); int iB = beamB[iSys].iPos(); double sHat = sHatSys[iSys]; // Classify system: rescattering on either or both sides? bool normalA = beamA[iSys].isFromBeam(); bool normalB = beamB[iSys].isFromBeam(); bool normalSys = normalA && normalB; bool doubleRes = !normalA && !normalB; // Check whether final-state system momentum matches initial-state one. if (allowRescatter && CORRECTMISMATCH) { Vec4 pInitial = event[iA].p() + event[iB].p(); Vec4 pFinal; for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) { int iAB = partonSystemsPtr->getOut(iSys, iMem); if (event[iAB].isFinal()) pFinal += event[iAB].p(); } // Scale down primordial kT from side A if p+ increased. if (normalA && pFinal.pPos() > pInitial.pPos()) beamA[iSys].scalePT( pInitial.pPos() / pFinal.pPos() ); // Scale down primordial kT from side B if p- increased. if (normalB && pFinal.pNeg() > pInitial.pNeg()) beamB[iSys].scalePT( pInitial.pNeg() / pFinal.pNeg() ); } // Rescatter: possible change in sign of lightcone momentum of a // rescattered parton. If this happens, try to pick // new primordial kT values if (allowRescatter && (event[iA].pPos() / beamA[iSys].pPos() < 0 || event[iB].pNeg() / beamB[iSys].pNeg() < 0) ) { infoPtr->errorMsg("Warning in BeamRemnants::setKinematics:" " change in lightcone momentum sign; retrying kinematics"); physical = false; break; } // Begin kinematics of partons after primordial kT has been added. double sHatTAft = sHat + pow2( beamA[iSys].px() + beamB[iSys].px()) + pow2( beamA[iSys].py() + beamB[iSys].py()); double w2A = beamA[iSys].mT2(); double w2B = beamB[iSys].mT2(); double w2Diff = sHatTAft - w2A - w2B; double lambda = pow2(w2Diff) - 4. * w2A * w2B; // Too large transverse momenta means that kinematics will not work. if (lambda <= 0.) { physical = false; break; } double lamRoot = sqrtpos( lambda ); // Mirror solution if the two incoming have reverse rapidity ordering. if (allowRescatter && doubleRes && (event[iA].pPos() * event[iB].pNeg() < event[iA].pNeg() * event[iB].pPos()) ) lamRoot = -lamRoot; // Two procedures, which agree for normal scattering, separate here. // First option keeps rapidity (and mass) of system unchanged by // primordial kT, by modification of rescattered parton. if (normalSys || doRescatterRestoreY || doubleRes) { // Find kinematics of initial system: transverse mass, and // longitudinal momentum carried by all or rescattered partons. double sHatTBef = sHat; double wPosBef, wNegBef, wPosBefRes, wNegBefRes; // Normal scattering. if (normalSys) { wPosBef = beamA[iSys].x() * eCM; wNegBef = beamB[iSys].x() * eCM; wPosBefRes = 0.; wNegBefRes = 0.; // Rescattering on side A. } else if (normalB) { sHatTBef += event[iA].pT2(); wPosBef = event[iA].pPos(); wNegBef = event[iA].pNeg() + beamB[iSys].x() * eCM; wPosBefRes = beamA[iSys].pPos(); wNegBefRes = beamA[iSys].pNeg(); // Rescattering on side B. } else if (normalA) { sHatTBef += event[iB].pT2(); wPosBef = beamA[iSys].x() * eCM + event[iB].pPos(); wNegBef = event[iB].pNeg(); wPosBefRes = beamB[iSys].pPos(); wNegBefRes = beamB[iSys].pNeg(); // Rescattering on both sides. } else { sHatTBef += pow2( event[iA].px() + event[iB].px()) + pow2( event[iA].py() + event[iB].py()); wPosBef = event[iA].pPos() + event[iB].pPos(); wNegBef = event[iA].pNeg() + event[iB].pNeg(); wPosBefRes = beamA[iSys].pPos() + beamB[iSys].pPos(); wNegBefRes = beamA[iSys].pNeg() + beamB[iSys].pNeg(); } // Rescale outgoing momenta to keep same mass and rapidity of system // as originally, and solve for kinematics. double rescale = sqrt(sHatTAft / sHatTBef); double wPosAft = rescale * wPosBef; double wNegAft = rescale * wNegBef; wPosRem -= wPosAft - wPosBefRes; wNegRem -= wNegAft - wNegBefRes; double wPosA = 0.5 * (sHatTAft + w2A - w2B + lamRoot) / wNegAft; double wNegB = 0.5 * (sHatTAft + w2B - w2A + lamRoot) / wPosAft; // Store modified beam parton momenta. beamA[iSys].e( 0.5 * (wPosA + w2A / wPosA) ); beamA[iSys].pz( 0.5 * (wPosA - w2A / wPosA) ); beamB[iSys].e( 0.5 * (w2B / wNegB + wNegB) ); beamB[iSys].pz( 0.5 * (w2B / wNegB - wNegB) ); // Second option keeps rescattered parton (and system mass) unchanged // by primordial kT, by modification of system rapidity. } else { // Rescattering on side A: preserve already scattered parton. if (normalB) { double wPosA = beamA[iSys].pPos(); double wNegB = 0.5 * (w2Diff + lamRoot) / wPosA; beamB[iSys].e( 0.5 * (w2B / wNegB + wNegB) ); beamB[iSys].pz( 0.5 * (w2B / wNegB - wNegB) ); wPosRem -= w2B / wNegB; wNegRem -= wNegB; // Rescattering on side B: preserve already scattered parton. } else if (normalA) { double wNegB = beamB[iSys].pNeg(); double wPosA = 0.5 * (w2Diff + lamRoot) / wNegB; beamA[iSys].e( 0.5 * (wPosA + w2A / wPosA) ); beamA[iSys].pz( 0.5 * (wPosA - w2A / wPosA) ); wPosRem -= wPosA; wNegRem -= w2A / wPosA; // Primordial kT in double rescattering does change the mass of // the system without any possible fix in this framework, which // leads to incorrect boosts. Current choice is therefore always // to handle it with the first procedure, where mass is retained. } else { } } // Construct system rotation and boost caused by primordial kT. Msys[iSys].reset(); Msys[iSys].toCMframe( event[iA].p(), event[iB].p() ); Msys[iSys].fromCMframe( beamA[iSys].p(), beamB[iSys].p() ); // Boost rescattered partons in subsequent beam A list. for (int iSys2 = iSys + 1; iSys2 < nSys; ++iSys2) { if (!beamA[iSys2].isFromBeam()) { int iBefResc = event[ beamA[iSys2].iPos() ].mother1(); for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) if (partonSystemsPtr->getOut(iSys, iMem) == iBefResc) { Vec4 pTemp = event[iBefResc].p(); pTemp.rotbst( Msys[iSys] ); beamA[iSys2].p( pTemp ); } } // Boost rescattered partons in subsequent beam B list. if (!beamB[iSys2].isFromBeam()) { int iBefResc = event[ beamB[iSys2].iPos() ].mother1(); for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) if (partonSystemsPtr->getOut(iSys, iMem) == iBefResc) { Vec4 pTemp = event[iBefResc].p(); pTemp.rotbst( Msys[iSys] ); beamB[iSys2].p( pTemp ); } } } } // Check that remaining momentum is enough for remnants. if (wPosRem < 0. || wNegRem < 0.) physical = false; w2Rem = wPosRem * wNegRem; if (sqrtpos(w2Rem) < sqrt(w2Beam[0]) + sqrt(w2Beam[1])) physical = false; // End of loop over ten tries. Do not loop when solution found. if (physical) break; } // If no solution after ten tries then failed. if (!physical) { infoPtr->errorMsg("Error in BeamRemnants::setKinematics:" " kinematics construction failed"); return false; } // For successful initiator kinematics process whole systems. Vec4 pSumOut; for (int iSys = 0; iSys < nSys; ++iSys) { // Copy initiators and their systems and boost them accordingly. // Update subsystem and beams info on new positions of partons. // Update daughter info of mothers, i.e. of beams, for hardest interaction. if (beamA[iSys].isFromBeam()) { int iA = beamA[iSys].iPos(); int iAcopy = event.copy(iA, -61); event[iAcopy].rotbst(Msys[iSys]); partonSystemsPtr->setInA(iSys, iAcopy); beamA[iSys].iPos( iAcopy); if (iSys == 0) { int mother = event[iAcopy].mother1(); event[mother].daughter1(iAcopy); } } if (beamB[iSys].isFromBeam()) { int iB = beamB[iSys].iPos(); int iBcopy = event.copy(iB, -61); event[iBcopy].rotbst(Msys[iSys]); partonSystemsPtr->setInB(iSys, iBcopy); beamB[iSys].iPos( iBcopy); if (iSys == 0) { int mother = event[iBcopy].mother1(); event[mother].daughter1(iBcopy); } } for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) { int iAB = partonSystemsPtr->getOut(iSys, iMem); if (event[iAB].isFinal()) { int iABcopy = event.copy(iAB, 62); event[iABcopy].rotbst(Msys[iSys]); partonSystemsPtr->setOut(iSys, iMem, iABcopy); pSumOut += event[iABcopy].p(); } } } // Colour dipoles spanning systems gives mismatch between FSR recoils // and primordial kT boosts. if (allowRescatter && CORRECTMISMATCH) { // Find summed pT of beam remnants = - wanted pT of systems. double pxBeams = 0.; double pyBeams = 0.; for (int iRem = nSys; iRem < beamA.size(); ++iRem) { pxBeams += beamA[iRem].px(); pyBeams += beamA[iRem].py(); } for (int iRem = nSys; iRem < beamB.size(); ++iRem) { pxBeams += beamB[iRem].px(); pyBeams += beamB[iRem].py(); } // Boost all final partons in systems transversely, and also their sum. Vec4 pSumTo( -pxBeams, -pyBeams, pSumOut.pz(), sqrt( pow2(pxBeams) + pow2(pyBeams) + pow2(pSumOut.pz()) + pSumOut.m2Calc() ) ); RotBstMatrix Mmismatch; Mmismatch.bst( pSumOut, pSumTo); for (int iSys = 0; iSys < nSys; ++iSys) for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) { int iAB = partonSystemsPtr->getOut(iSys, iMem); if (event[iAB].isFinal()) event[iAB].rotbst(Mmismatch); } pSumOut.rotbst(Mmismatch); // Reset energy and momentum sum, to be compensated by beam remnants. wPosRem = eCM - (pSumOut.e() + pSumOut.pz()); wNegRem = eCM - (pSumOut.e() - pSumOut.pz()); w2Rem = wPosRem * wNegRem; if ( wPosRem < 0. || wNegRem < 0. || sqrtpos(w2Rem) < sqrt(w2Beam[0]) + sqrt(w2Beam[1])) { infoPtr->errorMsg("Error in BeamRemnants::setKinematics:" " kinematics construction failed owing to recoil mismatch"); return false; } } // Construct x rescaling factors for the two remants. double lambdaRoot = sqrtpos( pow2(w2Rem - w2Beam[0] - w2Beam[1]) - 4. * w2Beam[0] * w2Beam[1] ); double rescaleA = (w2Rem + w2Beam[0] - w2Beam[1] + lambdaRoot) / (2. * w2Rem * xSum[0]) ; double rescaleB = (w2Rem + w2Beam[1] - w2Beam[0] + lambdaRoot) / (2. * w2Rem * xSum[1]) ; // Construct energy and pz for remnants in first beam. for (int iRem = nSys; iRem < beamA.size(); ++iRem) { double pPos = rescaleA * beamA[iRem].x() * wPosRem; double pNeg = beamA[iRem].mT2() / pPos; beamA[iRem].e( 0.5 * (pPos + pNeg) ); beamA[iRem].pz( 0.5 * (pPos - pNeg) ); // Add these partons to the normal event record. int iNew = event.append( beamA[iRem].id(), 63, 1 + nOffset, 0, 0, 0, beamA[iRem].col(), beamA[iRem].acol(), beamA[iRem].p(), beamA[iRem].m() ); beamA[iRem].iPos( iNew); } // Construct energy and pz for remnants in second beam. for (int iRem = nSys; iRem < beamB.size(); ++iRem) { double pNeg = rescaleB * beamB[iRem].x() * wNegRem; double pPos = beamB[iRem].mT2() / pNeg; beamB[iRem].e( 0.5 * (pPos + pNeg) ); beamB[iRem].pz( 0.5 * (pPos - pNeg) ); // Add these partons to the normal event record. int iNew = event.append( beamB[iRem].id(), 63, 2 + nOffset, 0, 0, 0, beamB[iRem].col(), beamB[iRem].acol(), beamB[iRem].p(), beamB[iRem].m() ); beamB[iRem].iPos( iNew); } // Done. return true; } //-------------------------------------------------------------------------- // Allow colour reconnections by mergings of collision subsystems. // iRec is system that may be reconnected, by moving its gluons to iSys, // where minimal pT (or equivalently Lambda) is used to pick location. // Therefore all dipoles in iSys have to be found, and all gluons in iRec. // Matching q-qbar pairs are treated by analogy with gluons. // Note: owing to rescatterings some outgoing partons must be skipped. bool BeamRemnants::reconnectColours( Event& event) { // References to beams to simplify indexing. BeamParticle& beamA = *beamAPtr; BeamParticle& beamB = *beamBPtr; // Prepare record of which systems should be merged onto another. // The iSys system must have colour in final state to attach to it. vector iMerge(nSys); vector hasColour(nSys); for (int iSys = 0; iSys < nSys; ++iSys) { iMerge[iSys] = iSys; bool hasCol = false; for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) { int iNow = partonSystemsPtr->getOut( iSys, iMem); if (event[iNow].isFinal() && (event[iNow].col() > 0 || event[iNow].acol() > 0) ) { hasCol = true; break; } } hasColour[iSys] = hasCol; } // Loop over systems to decide which should be reconnected. for (int iRec = nSys - 1; iRec > 0; --iRec) { // Determine reconnection strength from pT scale of system. double pT2Rec = pow2( infoPtr->pTMI(iRec) ); double probRec = pT20Rec / (pT20Rec + pT2Rec); // Loop over other systems iSys at higher pT scale and // decide whether to reconnect the iRec gluons onto one of them. for (int iSys = iRec - 1; iSys >= 0; --iSys) if (hasColour[iSys] && probRec > rndmPtr->flat()) { // The iRec system and all merged with it to be merged with iSys. iMerge[iRec] = iSys; for (int iRec2 = iRec + 1; iRec2 < nSys; ++iRec2) if (iMerge[iRec2] == iRec) iMerge[iRec2] = iSys; // Once a system has been merged do not test it anymore. break; } } // Loop over systems. Check whether other systems to be merged with it. for (int iSys = 0; iSys < nSys; ++iSys) { int nMerge = 0; for (int iRec = iSys + 1; iRec < nSys; ++iRec) if (iMerge[iRec] == iSys) ++nMerge; if (nMerge == 0) continue; // Incoming partons not counted if rescattered. int iInASys = partonSystemsPtr->getInA(iSys); bool hasInA = (beamA[iSys].isFromBeam()); int iInBSys = partonSystemsPtr->getInB(iSys); bool hasInB = (beamB[iSys].isFromBeam()); // Begin find dipoles in iSys system. vector dipoles; int sizeOut = partonSystemsPtr->sizeOut(iSys); for (int iMem = 0; iMem < sizeOut; ++iMem) { // Find colour dipoles to beam remnant. int iNow = partonSystemsPtr->getOut( iSys, iMem); if (!event[iNow].isFinal()) continue; int col = event[iNow].col(); if (col > 0) { if (hasInA && event[iInASys].col() == col) dipoles.push_back( BeamDipole( col, iNow, iInASys ) ); else if (hasInB && event[iInBSys].col() == col) dipoles.push_back( BeamDipole( col, iNow, iInBSys ) ); // Find colour dipole between final-state partons. else for (int iMem2 = 0; iMem2 < sizeOut; ++iMem2) if (iMem2 != iMem) { int iNow2 = partonSystemsPtr->getOut( iSys, iMem2); if (!event[iNow2].isFinal()) continue; if (event[iNow2].acol() == col) { dipoles.push_back( BeamDipole( col, iNow, iNow2) ); break; } } } // Find anticolour dipoles to beam remnant. int acol = event[iNow].acol(); if (acol > 0) { if (hasInA && event[iInASys].acol() == acol) dipoles.push_back( BeamDipole( acol, iInASys, iNow ) ); else if (hasInB && event[iInBSys].acol() == acol) dipoles.push_back( BeamDipole( acol, iInBSys, iNow ) ); } } // Skip mergings if no dipoles found. if (dipoles.size() == 0) continue; // Find dipole sizes. for (int iDip = 0; iDip < int(dipoles.size()); ++iDip) dipoles[iDip].p1p2 = event[dipoles[iDip].iCol].p() * event[dipoles[iDip].iAcol].p(); // Loop over systems iRec to be merged with iSys. for (int iRec = iSys + 1; iRec < nSys; ++iRec) { if (iMerge[iRec] != iSys) continue; // Information on iRec. Vectors for gluons and anything else. int sizeRec = partonSystemsPtr->sizeOut(iRec); int iInARec = partonSystemsPtr->getInA(iRec); int iInBRec = partonSystemsPtr->getInB(iRec); int nGluRec = 0; vector iGluRec; vector pT2GluRec; int nAnyRec = 0; vector iAnyRec; vector freeAnyRec; // Copy of gluon positions in descending order. for (int iMem = 0; iMem < sizeRec; ++iMem) { int iNow = partonSystemsPtr->getOut( iRec, iMem); if (!event[iNow].isFinal()) continue; if (event[iNow].isGluon()) { ++nGluRec; iGluRec.push_back( iNow ); pT2GluRec.push_back( event[iNow].pT2() ); for (int i = nGluRec - 1; i > 1; --i) { if (pT2GluRec[i - 1] > pT2GluRec[i]) break; swap( iGluRec[i - 1], iGluRec[i] ); swap( pT2GluRec[i - 1], pT2GluRec[i] ); } // Copy of anything else, mainly quarks, in no particular order. } else { ++nAnyRec; iAnyRec.push_back( iNow ); freeAnyRec.push_back( true ); } } // For each gluon in iRec now find the dipole that gives the smallest // (pGlu * pI) (pGlu * pJ) / (pI * pJ), i.e. minimal pT (and Lambda). for (int iGRec = 0; iGRec < nGluRec; ++iGRec) { int iGlu = iGluRec[iGRec]; Vec4 pGlu = event[iGlu].p(); int iDipMin = 0; double pT2DipMin = sCM; for (int iDip = 0; iDip < int(dipoles.size()); ++iDip) { double pT2Dip = (pGlu * event[dipoles[iDip].iCol].p()) * (pGlu * event[dipoles[iDip].iAcol].p()) / dipoles[iDip].p1p2; if (pT2Dip < pT2DipMin) { iDipMin = iDip; pT2DipMin = pT2Dip; } } // Attach the gluon to the dipole, i.e. split the dipole in two. int colGlu = event[iGlu].col(); int acolGlu = event[iGlu].acol(); int colDip = dipoles[iDipMin].col; int iColDip = dipoles[iDipMin].iCol; int iAcolDip = dipoles[iDipMin].iAcol; event[iGlu].acol( colDip ); if (event[iAcolDip].acol() == colDip) event[iAcolDip].acol( colGlu ); else event[iAcolDip].col( colGlu ); dipoles[iDipMin].iAcol = iGlu; dipoles[iDipMin].p1p2 = event[iColDip].p() * pGlu; dipoles.push_back( BeamDipole( colGlu, iGlu, iAcolDip ) ); dipoles.back().p1p2 = pGlu * event[iAcolDip].p(); // Remove gluon from old system: reconnect colours. for (int i = oldSize; i < event.size(); ++i) if (i != iGlu && i != iAcolDip) { if (event[i].isFinal()) { if (event[i].acol() == colGlu) event[i].acol( acolGlu ); } else { if (event[i].col() == colGlu) event[i].col( acolGlu ); } } } // See if any matching quark-antiquark pairs among the rest. for (int iQRec = 0; iQRec < nAnyRec; ++iQRec) { int iQ = iAnyRec[iQRec]; int idQ = event[iQ].id(); if (freeAnyRec[iQRec] && idQ > 0 && idQ < 6) for (int iQbarRec = 0; iQbarRec < nAnyRec; ++iQbarRec) { int iQbar = iAnyRec[iQbarRec]; if (freeAnyRec[iQbarRec] && event[iQbar].id() == -idQ) { // Check that these can be traced back to same gluon splitting. // For now also avoid qqbar pairs produced in rescatterings.?? int iTopQ = event.iTopCopyId(iQ); int iTopQbar = event.iTopCopyId(iQbar); int iMother = event[iTopQ].mother1(); if (event[iTopQbar].mother1() == iMother && event[iMother].isGluon() && event[iMother].status() != -34 && event[iMother + 1].status() != -34 ) { // Now find the dipole that gives the smallest // ((pQ + pQbar) * pI) ((pQ + pQbar) * pJ) / (pI * pJ). Vec4 pGlu = event[iQ].p() + event[iQbar].p(); int iDipMin = 0; double pT2DipMin = sCM; for (int iDip = 0; iDip < int(dipoles.size()); ++iDip) { double pT2Dip = (pGlu * event[dipoles[iDip].iCol].p()) * (pGlu * event[dipoles[iDip].iAcol].p()) / dipoles[iDip].p1p2; if (pT2Dip < pT2DipMin) { iDipMin = iDip; pT2DipMin = pT2Dip; } } // Attach the q-qbar pair to the dipole, i.e. split the dipole. int colGlu = event[iQ].col(); int acolGlu = event[iQbar].acol(); int colDip = dipoles[iDipMin].col; int iColDip = dipoles[iDipMin].iCol; int iAcolDip = dipoles[iDipMin].iAcol; event[iQbar].acol( colDip ); if (event[iAcolDip].acol() == colDip) event[iAcolDip].acol( colGlu ); else event[iAcolDip].col( colGlu ); dipoles[iDipMin].iAcol = iQbar; dipoles[iDipMin].p1p2 = event[iColDip].p() * event[iQbar].p(); dipoles.push_back( BeamDipole( colGlu, iQ, iAcolDip ) ); dipoles.back().p1p2 = event[iQ].p() * event[iAcolDip].p(); // Remove q-qbar pair from old system: reconnect colours. freeAnyRec[iQRec] = false; freeAnyRec[iQbarRec] = false; for (int i = oldSize; i < event.size(); ++i) if (i != iQRec && i != iQbarRec && i != iColDip && i != iAcolDip) { if (event[i].isFinal()) { if (event[i].acol() == colGlu) event[i].acol( acolGlu ); } else { if (event[i].col() == colGlu) event[i].col( acolGlu ); } } // Done with processing of q-qbar pairs. } } } } // If only two beam gluons left of system, set their colour = anticolour. // Used by BeamParticle::remnantColours to skip irrelevant gluons. if ( event[iInARec].isGluon() && !event[iInARec].isRescatteredIncoming() && event[iInBRec].isGluon() && !event[iInBRec].isRescatteredIncoming() && event[iInARec].col() == event[iInBRec].acol() && event[iInARec].acol() == event[iInBRec].col() ) { event[iInARec].acol( event[iInARec].col() ); event[iInBRec].acol( event[iInBRec].col() ); } // End of loops over iRec and iSys systems. } } // Done. return true; } //-------------------------------------------------------------------------- // Collapse colours and check that they are consistent. bool BeamRemnants::checkColours( Event& event) { // No colours in lepton beams so no need to do anything. if (beamAPtr->isLepton() && beamBPtr->isLepton()) return true; // Remove ambiguities when one colour collapses two ways. // Resolve chains where one colour is mapped to another. for (int iCol = 1; iCol < int(colFrom.size()); ++iCol) for (int iColRef = 0; iColRef < iCol; ++iColRef) { if (colFrom[iCol] == colFrom[iColRef]) { colFrom[iCol] = colTo[iCol]; colTo[iCol] = colTo[iColRef]; } if (colTo[iCol] == colFrom[iColRef]) colTo[iCol] = colTo[iColRef]; } // Transform event record colours from beam remnant colour collapses. for (int i = oldSize; i < event.size(); ++i) { int col = event[i].col(); int acol = event[i].acol(); for (int iCol = 0; iCol < int(colFrom.size()); ++iCol) { if (col == colFrom[iCol]) {col = colTo[iCol]; event[i].col(col);} if (acol == colFrom[iCol]) {acol = colTo[iCol]; event[i].acol(acol);} } } // Transform junction colours from beam remnant colour collapses. for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) for (int leg = 0; leg < 3; ++leg) { int col = event.colJunction(iJun, leg); //cout<< " BeamRemnants iJun = "< 0) colList.push_back( col ); if (acol > 0) acolList.push_back( acol ); if (col > 0 && acol == col) iSingletGluon.push_back(i); } // Run though list of singlet gluons and put them on final-state dipole // (i,j) that offers smallest (p_g p_i) * (p_g p_j) / (p_i p_j). for (int iS = 0; iS < int(iSingletGluon.size()); ++iS) { int iGlu = iSingletGluon[iS]; int iAcolDip = -1; double pT2DipMin = sCM; for (int iC = oldSize; iC < event.size(); ++iC) if (iC != iGlu && event[iC].isFinal()) { int colDip = event[iC].col(); if (colDip > 0 && event[iC].acol() !=colDip) for (int iA = oldSize; iA < event.size(); ++iA) if (iA != iGlu && iA != iC && event[iA].isFinal() && event[iA].acol() == colDip && event[iA].col() !=colDip) { double pT2Dip = (event[iGlu].p() * event[iC].p()) * (event[iGlu].p() * event[iA].p()) / (event[iC].p() * event[iA].p()); if (pT2Dip < pT2DipMin) { iAcolDip = iA; pT2DipMin = pT2Dip; } } } // Fail if no dipole. Else insert singlet gluon onto relevant dipole. if (iAcolDip == -1) return false; event[iGlu].acol( event[iAcolDip].acol() ); event[iAcolDip].acol( event[iGlu].col() ); } // Check that not the same colour or anticolour appears twice. for (int iCol = 0; iCol < int(colList.size()) - 1; ++iCol) { int col = colList[iCol]; for (int iCol2 = iCol + 1; iCol2 < int(colList.size()); ++iCol2) if (colList[iCol2] == col) { infoPtr->errorMsg("Warning in BeamRemnants::checkColours:" " colour appears twice"); if (!ALLOWCOLOURTWICE) return false; } } for (int iAcol = 0; iAcol < int(acolList.size()) - 1; ++iAcol) { int acol = acolList[iAcol]; for (int iAcol2 = iAcol + 1; iAcol2 < int(acolList.size()); ++iAcol2) if (acolList[iAcol2] == acol) { infoPtr->errorMsg("Warning in BeamRemnants::checkColours:" " anticolour appears twice"); if (!ALLOWCOLOURTWICE) return false; } } // Remove all matching colour-anticolour pairs. bool foundPair = true; while (foundPair && colList.size() > 0 && acolList.size() > 0) { foundPair = false; for (int iCol = 0; iCol < int(colList.size()); ++iCol) { for (int iAcol = 0; iAcol < int(acolList.size()); ++iAcol) { if (acolList[iAcol] == colList[iCol]) { colList[iCol] = colList.back(); colList.pop_back(); acolList[iAcol] = acolList.back(); acolList.pop_back(); foundPair = true; break; } } if (foundPair) break; } } // Check that remaining (anti)colours are accounted for by junctions. for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) { int kindJun = event.kindJunction(iJun); for (int leg = 0; leg < 3; ++leg) { int colEnd = event.colJunction(iJun, leg); // Junction connected to three colours. if (kindJun == 1) { bool foundCol = false; for (int iCol = 0; iCol < int(colList.size()); ++iCol) if (colList[iCol] == colEnd) { colList[iCol] = colList.back(); colList.pop_back(); foundCol = true; break; } } // Junction connected to three anticolours. else if (kindJun == 2) { bool foundCol = false; for (int iAcol = 0; iAcol < int(acolList.size()); ++iAcol) if (acolList[iAcol] == colEnd) { acolList[iAcol] = acolList.back(); acolList.pop_back(); foundCol = true; break; } } // Other junction types else if ( kindJun == 3 || kindJun == 5) { bool foundCol = false; for (int iCol = 0; iCol < int(colList.size()); ++iCol) if (colList[iCol] == colEnd) { colList[iCol] = colList.back(); colList.pop_back(); foundCol = true; break; } } // Other antijunction types else if ( kindJun == 4 || kindJun == 6) { bool foundCol = false; for (int iAcol = 0; iAcol < int(acolList.size()); ++iAcol) if (acolList[iAcol] == colEnd) { acolList[iAcol] = acolList.back(); acolList.pop_back(); foundCol = true; break; } } // End junction check. } } // Repair step - sometimes needed when rescattering allowed. if (colList.size() > 0 || acolList.size() > 0) infoPtr->errorMsg("Warning in BeamRemnants::checkColours:" " need to repair unmatched colours"); while (colList.size() > 0 && acolList.size() > 0) { // Replace one colour and one anticolour index by a new common one. int colMatch = colList.back(); int acolMatch = acolList.back(); int colNew = event.nextColTag(); colList.pop_back(); acolList.pop_back(); for (int i = oldSize; i < event.size(); ++i) if (event[i].isFinal() && event[i].col() == colMatch) { event[i].col( colNew); break; } for (int i = oldSize; i < event.size(); ++i) if (event[i].isFinal() && event[i].acol() == acolMatch) { event[i].acol( colNew); break; } } // Done. return (colList.size() == 0 && acolList.size() == 0); } //========================================================================== } // end namespace Pythia8