9 #include "Pythia8/BeamRemnants.h"
22 BeamDipole(
int colIn = 0,
int iColIn = 0,
int iAcolIn = 0)
23 : col(colIn), iCol(iColIn), iAcol(iAcolIn) {}
41 const bool BeamRemnants::ALLOWCOLOURTWICE =
true;
44 const int BeamRemnants::NTRYCOLMATCH = 10;
45 const int BeamRemnants::NTRYKINMATCH = 10;
50 const bool BeamRemnants::CORRECTMISMATCH =
false;
56 bool BeamRemnants::init( Info* infoPtrIn, Settings& settings, Rndm* rndmPtrIn,
57 BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
58 PartonSystems* partonSystemsPtrIn) {
63 beamAPtr = beamAPtrIn;
64 beamBPtr = beamBPtrIn;
65 partonSystemsPtr = partonSystemsPtrIn;
68 doPrimordialKT = settings.flag(
"BeamRemnants:primordialKT");
69 primordialKTsoft = settings.parm(
"BeamRemnants:primordialKTsoft");
70 primordialKThard = settings.parm(
"BeamRemnants:primordialKThard");
71 primordialKTremnant = settings.parm(
"BeamRemnants:primordialKTremnant");
72 halfScaleForKT = settings.parm(
"BeamRemnants:halfScaleForKT");
73 halfMassForKT = settings.parm(
"BeamRemnants:halfMassForKT");
76 allowRescatter = settings.flag(
"MultipartonInteractions:allowRescatter");
77 doRescatterRestoreY = settings.flag(
"BeamRemnants:rescatterRestoreY");
81 doReconnect = settings.flag(
"BeamRemnants:reconnectColours");
82 reconnectRange = settings.parm(
"BeamRemnants:reconnectRange");
83 pT0Ref = settings.parm(
"MultipartonInteractions:pT0Ref");
84 ecmRef = settings.parm(
"MultipartonInteractions:ecmRef");
85 ecmPow = settings.parm(
"MultipartonInteractions:ecmPow");
92 pT0 = pT0Ref * pow(eCM / ecmRef, ecmPow);
93 pT20Rec = pow2(reconnectRange * pT0);
106 bool BeamRemnants::add(
Event& event) {
109 eCM = infoPtr->eCM();
113 for (
int i = 0; i < beamAPtr->size(); ++i) {
114 int j = (*beamAPtr)[i].iPos();
115 if ((*beamAPtr)[i].
id() != event[j].
id()) {
116 infoPtr->errorMsg(
"Error in BeamRemnants::add: "
117 "event and beam flavours do not match");
121 for (
int i = 0; i < beamBPtr->size(); ++i) {
122 int j = (*beamBPtr)[i].iPos();
123 if ((*beamBPtr)[i].
id() != event[j].
id()) {
124 infoPtr->errorMsg(
"Error in BeamRemnants::add: "
125 "event and beam flavours do not match");
131 nSys = partonSystemsPtr->sizeSys();
132 oldSize =
event.size();
136 if ( !beamAPtr->remnantFlavours(event)
137 || !beamBPtr->remnantFlavours(event) ) {
138 infoPtr->errorMsg(
"Error in BeamRemnants::add:"
139 " remnant flavour setup failed");
144 if (!setKinematics(event))
return false;
147 if (doReconnect) reconnectColours(event);
151 vector<int> acolSave;
152 for (
int i = oldSize; i <
event.size(); ++i) {
153 colSave.push_back( event[i].col() );
154 acolSave.push_back( event[i].acol() );
156 event.saveJunctionSize();
161 bool physical =
true;
162 for (
int iTry = 0; iTry < NTRYCOLMATCH; ++iTry) {
170 if (!beamAPtr->remnantColours(event, colFrom, colTo))
172 if (!beamBPtr->remnantColours(event, colFrom, colTo))
176 if ( physical && !checkColours(event) ) physical =
false;
180 for (
int i = oldSize; i <
event.size(); ++i)
181 event[i].cols( colSave[i - oldSize], acolSave[i - oldSize] );
182 event.restoreJunctionSize();
183 infoPtr->errorMsg(
"Warning in BeamRemnants::add:"
184 " colour tracing failed; will try again");
189 infoPtr->errorMsg(
"Error in BeamRemnants::add:"
190 " colour tracing failed after several attempts");
204 bool BeamRemnants::setKinematics(
Event& event) {
207 BeamParticle& beamA = *beamAPtr;
208 BeamParticle& beamB = *beamBPtr;
211 if ( beamA.isUnresolvedLepton() && beamB.isUnresolvedLepton() )
215 if ( (!beamA.isLepton() && beamA.xMax(-1) <= 0.)
216 || (!beamB.isLepton() && beamB.xMax(-1) <= 0.) ) {
217 infoPtr->errorMsg(
"Error in BeamRemnants::setKinematics:"
218 " no momentum left for beam remnants");
224 for (
int i = 3; i <
event.size(); ++i)
225 if (event[i].statusAbs() < 20) nBeams = i + 1;
226 int nOffset = nBeams - 3;
229 int nMaxBeam = max( beamA.size(), beamB.size() );
230 vector<double> sHatSys(nMaxBeam);
231 vector<double> kTwidth(nMaxBeam);
232 vector<double> kTcomp(nMaxBeam);
233 vector<RotBstMatrix> Msys(nSys);
236 double kTcompSumA = 0.;
237 double kTcompSumB = 0.;
238 for (
int iSys = 0; iSys < nSys; ++iSys) {
239 double kTwidthNow = 0.;
240 double mHatDamp = 1.;
241 int iInA = partonSystemsPtr->getInA(iSys);
242 int iInB = partonSystemsPtr->getInB(iSys);
243 double sHatNow = (
event[iInA].p() +
event[iInB].p()).m2Calc();
247 if (doPrimordialKT) {
248 double mHat = sqrt(sHatNow);
249 mHatDamp = mHat / (mHat + halfMassForKT);
250 double scale = (iSys == 0) ? infoPtr->QRen(iDS)
251 : partonSystemsPtr->getPTHat(iSys);
252 kTwidthNow = ( (halfScaleForKT * primordialKTsoft
253 + scale * primordialKThard) / (halfScaleForKT + scale) ) * mHatDamp;
258 sHatSys[iSys] = sHatNow;
259 kTwidth[iSys] = kTwidthNow ;
260 kTcomp[iSys] = mHatDamp;
261 if (beamA[iSys].isFromBeam()) kTcompSumA += mHatDamp;
262 else beamA[iSys].m( event[iInA].m() );
263 if (beamB[iSys].isFromBeam()) kTcompSumB += mHatDamp;
264 else beamB[iSys].m( event[iInB].m() );
268 double kTwidthNow = (doPrimordialKT) ? primordialKTremnant : 0.;
269 for (
int iRem = nSys; iRem < nMaxBeam; ++iRem) {
271 kTwidth[iRem] = kTwidthNow ;
274 kTcompSumA += beamA.size() - nSys;
275 kTcompSumB += beamB.size() - nSys;
279 double xSum[2], xInvM[2], w2Beam[2], wPosRem, wNegRem, w2Rem;
280 for (
int iTry = 0; iTry < NTRYKINMATCH; ++iTry) {
284 for (
int iBeam = 0; iBeam < 2; ++iBeam) {
285 BeamParticle&
beam = (iBeam == 0) ? beamA : beamB;
286 int nPar = beam.size();
290 if (beam.isHadron() && doPrimordialKT) {
293 for (
int iPar = 0; iPar < nPar; ++iPar) {
294 if ( beam[iPar].isFromBeam() ) {
295 pair<double, double> gauss2 = rndmPtr->gauss2();
296 double px = kTwidth[iPar] * gauss2.first;
297 double py = kTwidth[iPar] * gauss2.second;
303 int iInAB = (iBeam == 0) ? partonSystemsPtr->getInA(iPar)
304 : partonSystemsPtr->getInB(iPar);
305 beam[iPar].p( event[iInAB].p() );
310 double kTcompSum = (iBeam == 0) ? kTcompSumA : kTcompSumB;
311 for (
int iPar = 0; iPar < nPar; ++iPar)
312 if (beam[iPar].isFromBeam() ) {
313 beam[iPar].px( beam[iPar].px() - pxSum * kTcomp[iPar] / kTcompSum );
314 beam[iPar].py( beam[iPar].py() - pySum * kTcomp[iPar] / kTcompSum );
318 }
else if (beam.isHadron()) {
319 for (
int iPar = 0; iPar < nPar; ++iPar)
320 if ( !beam[iPar].isFromBeam() ) {
321 int iInAB = (iBeam == 0) ? partonSystemsPtr->getInA(iPar)
322 : partonSystemsPtr->getInB(iPar);
323 beam[iPar].p( event[iInAB].p() );
330 for (
int iRem = nSys; iRem < nPar; ++iRem) {
331 double xPrel = beam.xRemnant( iRem);
333 xSum[iBeam] += xPrel;
334 xInvM[iBeam] += beam[iRem].mT2()/xPrel;
338 w2Beam[iBeam] = xSum[iBeam] * xInvM[iBeam];
346 for (
int iSys = 0; iSys < nSys; ++iSys) {
347 int iA = beamA[iSys].iPos();
348 int iB = beamB[iSys].iPos();
349 double sHat = sHatSys[iSys];
352 bool normalA = beamA[iSys].isFromBeam();
353 bool normalB = beamB[iSys].isFromBeam();
354 bool normalSys = normalA && normalB;
355 bool doubleRes = !normalA && !normalB;
358 if (allowRescatter && CORRECTMISMATCH) {
359 Vec4 pInitial =
event[iA].p() +
event[iB].p();
361 for (
int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) {
362 int iAB = partonSystemsPtr->getOut(iSys, iMem);
363 if (event[iAB].isFinal()) pFinal +=
event[iAB].p();
367 if (normalA && pFinal.pPos() > pInitial.pPos())
368 beamA[iSys].scalePT( pInitial.pPos() / pFinal.pPos() );
371 if (normalB && pFinal.pNeg() > pInitial.pNeg())
372 beamB[iSys].scalePT( pInitial.pNeg() / pFinal.pNeg() );
379 && (event[iA].pPos() / beamA[iSys].pPos() < 0
380 || event[iB].pNeg() / beamB[iSys].pNeg() < 0) ) {
381 infoPtr->errorMsg(
"Warning in BeamRemnants::setKinematics:"
382 " change in lightcone momentum sign; retrying kinematics");
388 double sHatTAft = sHat + pow2( beamA[iSys].px() + beamB[iSys].px())
389 + pow2( beamA[iSys].py() + beamB[iSys].py());
390 double w2A = beamA[iSys].mT2();
391 double w2B = beamB[iSys].mT2();
392 double w2Diff = sHatTAft - w2A - w2B;
393 double lambda = pow2(w2Diff) - 4. * w2A * w2B;
400 double lamRoot = sqrtpos( lambda );
403 if (allowRescatter && doubleRes && (event[iA].pPos() * event[iB].pNeg()
404 < event[iA].pNeg() * event[iB].pPos()) ) lamRoot = -lamRoot;
409 if (normalSys || doRescatterRestoreY || doubleRes) {
413 double sHatTBef = sHat;
414 double wPosBef, wNegBef, wPosBefRes, wNegBefRes;
417 wPosBef = beamA[iSys].x() * eCM;
418 wNegBef = beamB[iSys].x() * eCM;
422 }
else if (normalB) {
423 sHatTBef +=
event[iA].pT2();
424 wPosBef =
event[iA].pPos();
425 wNegBef =
event[iA].pNeg() + beamB[iSys].x() * eCM;
426 wPosBefRes = beamA[iSys].pPos();
427 wNegBefRes = beamA[iSys].pNeg();
429 }
else if (normalA) {
430 sHatTBef +=
event[iB].pT2();
431 wPosBef = beamA[iSys].x() * eCM +
event[iB].pPos();
432 wNegBef =
event[iB].pNeg();
433 wPosBefRes = beamB[iSys].pPos();
434 wNegBefRes = beamB[iSys].pNeg();
437 sHatTBef += pow2( event[iA].px() + event[iB].px())
438 + pow2( event[iA].py() + event[iB].py());
439 wPosBef =
event[iA].pPos() +
event[iB].pPos();
440 wNegBef =
event[iA].pNeg() +
event[iB].pNeg();
441 wPosBefRes = beamA[iSys].pPos() + beamB[iSys].pPos();
442 wNegBefRes = beamA[iSys].pNeg() + beamB[iSys].pNeg();
447 double rescale = sqrt(sHatTAft / sHatTBef);
448 double wPosAft = rescale * wPosBef;
449 double wNegAft = rescale * wNegBef;
450 wPosRem -= wPosAft - wPosBefRes;
451 wNegRem -= wNegAft - wNegBefRes;
452 double wPosA = 0.5 * (sHatTAft + w2A - w2B + lamRoot) / wNegAft;
453 double wNegB = 0.5 * (sHatTAft + w2B - w2A + lamRoot) / wPosAft;
456 beamA[iSys].e( 0.5 * (wPosA + w2A / wPosA) );
457 beamA[iSys].pz( 0.5 * (wPosA - w2A / wPosA) );
458 beamB[iSys].e( 0.5 * (w2B / wNegB + wNegB) );
459 beamB[iSys].pz( 0.5 * (w2B / wNegB - wNegB) );
467 double wPosA = beamA[iSys].pPos();
468 double wNegB = 0.5 * (w2Diff + lamRoot) / wPosA;
469 beamB[iSys].e( 0.5 * (w2B / wNegB + wNegB) );
470 beamB[iSys].pz( 0.5 * (w2B / wNegB - wNegB) );
471 wPosRem -= w2B / wNegB;
476 }
else if (normalA) {
477 double wNegB = beamB[iSys].pNeg();
478 double wPosA = 0.5 * (w2Diff + lamRoot) / wNegB;
479 beamA[iSys].e( 0.5 * (wPosA + w2A / wPosA) );
480 beamA[iSys].pz( 0.5 * (wPosA - w2A / wPosA) );
482 wNegRem -= w2A / wPosA;
494 Msys[iSys].toCMframe( event[iA].p(), event[iB].p() );
495 Msys[iSys].fromCMframe( beamA[iSys].p(), beamB[iSys].p() );
498 for (
int iSys2 = iSys + 1; iSys2 < nSys; ++iSys2) {
499 if (!beamA[iSys2].isFromBeam()) {
500 int iBefResc =
event[ beamA[iSys2].iPos() ].mother1();
501 for (
int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem)
502 if (partonSystemsPtr->getOut(iSys, iMem) == iBefResc) {
503 Vec4 pTemp =
event[iBefResc].p();
504 pTemp.rotbst( Msys[iSys] );
505 beamA[iSys2].p( pTemp );
510 if (!beamB[iSys2].isFromBeam()) {
511 int iBefResc =
event[ beamB[iSys2].iPos() ].mother1();
512 for (
int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem)
513 if (partonSystemsPtr->getOut(iSys, iMem) == iBefResc) {
514 Vec4 pTemp =
event[iBefResc].p();
515 pTemp.rotbst( Msys[iSys] );
516 beamB[iSys2].p( pTemp );
523 if (wPosRem < 0. || wNegRem < 0.) physical =
false;
524 w2Rem = wPosRem * wNegRem;
525 if (sqrtpos(w2Rem) < sqrt(w2Beam[0]) + sqrt(w2Beam[1]))
534 infoPtr->errorMsg(
"Error in BeamRemnants::setKinematics:"
535 " kinematics construction failed");
541 for (
int iSys = 0; iSys < nSys; ++iSys) {
546 if (beamA[iSys].isFromBeam()) {
547 int iA = beamA[iSys].iPos();
548 int iAcopy =
event.copy(iA, -61);
549 event[iAcopy].rotbst(Msys[iSys]);
550 partonSystemsPtr->setInA(iSys, iAcopy);
551 beamA[iSys].iPos( iAcopy);
553 int mother =
event[iAcopy].mother1();
554 event[mother].daughter1(iAcopy);
557 if (beamB[iSys].isFromBeam()) {
558 int iB = beamB[iSys].iPos();
559 int iBcopy =
event.copy(iB, -61);
560 event[iBcopy].rotbst(Msys[iSys]);
561 partonSystemsPtr->setInB(iSys, iBcopy);
562 beamB[iSys].iPos( iBcopy);
564 int mother =
event[iBcopy].mother1();
565 event[mother].daughter1(iBcopy);
569 for (
int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) {
570 int iAB = partonSystemsPtr->getOut(iSys, iMem);
571 if (event[iAB].isFinal()) {
572 int iABcopy =
event.copy(iAB, 62);
573 event[iABcopy].rotbst(Msys[iSys]);
574 partonSystemsPtr->setOut(iSys, iMem, iABcopy);
575 pSumOut +=
event[iABcopy].p();
583 if (allowRescatter && CORRECTMISMATCH) {
588 for (
int iRem = nSys; iRem < beamA.size(); ++iRem) {
589 pxBeams += beamA[iRem].px();
590 pyBeams += beamA[iRem].py();
592 for (
int iRem = nSys; iRem < beamB.size(); ++iRem) {
593 pxBeams += beamB[iRem].px();
594 pyBeams += beamB[iRem].py();
598 Vec4 pSumTo( -pxBeams, -pyBeams, pSumOut.pz(), sqrt( pow2(pxBeams)
599 + pow2(pyBeams) + pow2(pSumOut.pz()) + pSumOut.m2Calc() ) );
600 RotBstMatrix Mmismatch;
601 Mmismatch.bst( pSumOut, pSumTo);
602 for (
int iSys = 0; iSys < nSys; ++iSys)
603 for (
int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) {
604 int iAB = partonSystemsPtr->getOut(iSys, iMem);
605 if (event[iAB].isFinal())
event[iAB].rotbst(Mmismatch);
607 pSumOut.rotbst(Mmismatch);
610 wPosRem = eCM - (pSumOut.e() + pSumOut.pz());
611 wNegRem = eCM - (pSumOut.e() - pSumOut.pz());
612 w2Rem = wPosRem * wNegRem;
613 if ( wPosRem < 0. || wNegRem < 0.
614 || sqrtpos(w2Rem) < sqrt(w2Beam[0]) + sqrt(w2Beam[1])) {
615 infoPtr->errorMsg(
"Error in BeamRemnants::setKinematics:"
616 " kinematics construction failed owing to recoil mismatch");
622 double lambdaRoot = sqrtpos( pow2(w2Rem - w2Beam[0] - w2Beam[1])
623 - 4. * w2Beam[0] * w2Beam[1] );
624 double rescaleA = (w2Rem + w2Beam[0] - w2Beam[1] + lambdaRoot)
625 / (2. * w2Rem * xSum[0]) ;
626 double rescaleB = (w2Rem + w2Beam[1] - w2Beam[0] + lambdaRoot)
627 / (2. * w2Rem * xSum[1]) ;
630 for (
int iRem = nSys; iRem < beamA.size(); ++iRem) {
631 double pPos = rescaleA * beamA[iRem].x() * wPosRem;
632 double pNeg = beamA[iRem].mT2() / pPos;
633 beamA[iRem].e( 0.5 * (pPos + pNeg) );
634 beamA[iRem].pz( 0.5 * (pPos - pNeg) );
637 int iNew =
event.append( beamA[iRem].
id(), 63, 1 + nOffset, 0, 0, 0,
638 beamA[iRem].col(), beamA[iRem].acol(), beamA[iRem].p(),
640 beamA[iRem].iPos( iNew);
644 for (
int iRem = nSys; iRem < beamB.size(); ++iRem) {
645 double pNeg = rescaleB * beamB[iRem].x() * wNegRem;
646 double pPos = beamB[iRem].mT2() / pNeg;
647 beamB[iRem].e( 0.5 * (pPos + pNeg) );
648 beamB[iRem].pz( 0.5 * (pPos - pNeg) );
651 int iNew =
event.append( beamB[iRem].
id(), 63, 2 + nOffset, 0, 0, 0,
652 beamB[iRem].col(), beamB[iRem].acol(), beamB[iRem].p(),
654 beamB[iRem].iPos( iNew);
671 bool BeamRemnants::reconnectColours(
Event& event) {
674 BeamParticle& beamA = *beamAPtr;
675 BeamParticle& beamB = *beamBPtr;
679 vector<int> iMerge(nSys);
680 vector<bool> hasColour(nSys);
681 for (
int iSys = 0; iSys < nSys; ++iSys) {
684 for (
int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) {
685 int iNow = partonSystemsPtr->getOut( iSys, iMem);
686 if (event[iNow].isFinal() && (event[iNow].col() > 0
687 || event[iNow].acol() > 0) ) {
692 hasColour[iSys] = hasCol;
696 for (
int iRec = nSys - 1; iRec > 0; --iRec) {
699 double pT2Rec = pow2( partonSystemsPtr->getPTHat(iRec) );
700 double probRec = pT20Rec / (pT20Rec + pT2Rec);
704 for (
int iSys = iRec - 1; iSys >= 0; --iSys)
705 if (hasColour[iSys] && probRec > rndmPtr->flat()) {
709 for (
int iRec2 = iRec + 1; iRec2 < nSys; ++iRec2)
710 if (iMerge[iRec2] == iRec) iMerge[iRec2] = iSys;
718 for (
int iSys = 0; iSys < nSys; ++iSys) {
720 for (
int iRec = iSys + 1; iRec < nSys; ++iRec)
721 if (iMerge[iRec] == iSys) ++nMerge;
722 if (nMerge == 0)
continue;
725 int iInASys = partonSystemsPtr->getInA(iSys);
726 bool hasInA = (beamA[iSys].isFromBeam());
727 int iInBSys = partonSystemsPtr->getInB(iSys);
728 bool hasInB = (beamB[iSys].isFromBeam());
731 vector<BeamDipole> dipoles;
732 int sizeOut = partonSystemsPtr->sizeOut(iSys);
733 for (
int iMem = 0; iMem < sizeOut; ++iMem) {
736 int iNow = partonSystemsPtr->getOut( iSys, iMem);
737 if (!event[iNow].isFinal())
continue;
738 int col =
event[iNow].col();
740 if (hasInA && event[iInASys].col() == col)
741 dipoles.push_back( BeamDipole( col, iNow, iInASys ) );
742 else if (hasInB && event[iInBSys].col() == col)
743 dipoles.push_back( BeamDipole( col, iNow, iInBSys ) );
746 else for (
int iMem2 = 0; iMem2 < sizeOut; ++iMem2)
748 int iNow2 = partonSystemsPtr->getOut( iSys, iMem2);
749 if (!event[iNow2].isFinal())
continue;
750 if (event[iNow2].acol() == col) {
751 dipoles.push_back( BeamDipole( col, iNow, iNow2) );
758 int acol =
event[iNow].acol();
760 if (hasInA && event[iInASys].acol() == acol)
761 dipoles.push_back( BeamDipole( acol, iInASys, iNow ) );
762 else if (hasInB && event[iInBSys].acol() == acol)
763 dipoles.push_back( BeamDipole( acol, iInBSys, iNow ) );
768 if (dipoles.size() == 0)
continue;
771 for (
int iDip = 0; iDip < int(dipoles.size()); ++iDip)
772 dipoles[iDip].p1p2 = event[dipoles[iDip].iCol].p()
773 *
event[dipoles[iDip].iAcol].p();
776 for (
int iRec = iSys + 1; iRec < nSys; ++iRec) {
777 if (iMerge[iRec] != iSys)
continue;
780 int sizeRec = partonSystemsPtr->sizeOut(iRec);
781 int iInARec = partonSystemsPtr->getInA(iRec);
782 int iInBRec = partonSystemsPtr->getInB(iRec);
785 vector<double> pT2GluRec;
788 vector<bool> freeAnyRec;
791 for (
int iMem = 0; iMem < sizeRec; ++iMem) {
792 int iNow = partonSystemsPtr->getOut( iRec, iMem);
793 if (!event[iNow].isFinal())
continue;
794 if (event[iNow].isGluon()) {
796 iGluRec.push_back( iNow );
797 pT2GluRec.push_back( event[iNow].pT2() );
798 for (
int i = nGluRec - 1; i > 1; --i) {
799 if (pT2GluRec[i - 1] > pT2GluRec[i])
break;
800 swap( iGluRec[i - 1], iGluRec[i] );
801 swap( pT2GluRec[i - 1], pT2GluRec[i] );
806 iAnyRec.push_back( iNow );
807 freeAnyRec.push_back(
true );
813 for (
int iGRec = 0; iGRec < nGluRec; ++iGRec) {
814 int iGlu = iGluRec[iGRec];
815 Vec4 pGlu =
event[iGlu].p();
817 double pT2DipMin = sCM;
818 for (
int iDip = 0; iDip < int(dipoles.size()); ++iDip) {
819 double pT2Dip = (pGlu *
event[dipoles[iDip].iCol].p())
820 * (pGlu * event[dipoles[iDip].iAcol].p()) / dipoles[iDip].p1p2;
821 if (pT2Dip < pT2DipMin) {
828 int colGlu =
event[iGlu].col();
829 int acolGlu =
event[iGlu].acol();
830 int colDip = dipoles[iDipMin].col;
831 int iColDip = dipoles[iDipMin].iCol;
832 int iAcolDip = dipoles[iDipMin].iAcol;
833 event[iGlu].acol( colDip );
834 if (event[iAcolDip].acol() == colDip)
835 event[iAcolDip].acol( colGlu );
836 else event[iAcolDip].col( colGlu );
837 dipoles[iDipMin].iAcol = iGlu;
838 dipoles[iDipMin].p1p2 =
event[iColDip].p() * pGlu;
839 dipoles.push_back( BeamDipole( colGlu, iGlu, iAcolDip ) );
840 dipoles.back().p1p2 = pGlu *
event[iAcolDip].p();
843 for (
int i = oldSize; i <
event.size(); ++i)
844 if (i != iGlu && i != iAcolDip) {
845 if (event[i].isFinal()) {
846 if (event[i].acol() == colGlu) event[i].acol( acolGlu );
848 if (event[i].col() == colGlu) event[i].col( acolGlu );
853 for (
int iJun = 0; iJun <
event.sizeJunction(); ++iJun) {
856 if (event.kindJunction(iJun) % 2 == 0)
continue;
857 for (
int leg = 0; leg < 3; ++leg) {
858 int col =
event.colJunction(iJun, leg);
860 event.colJunction(iJun, leg, colGlu);
867 for (
int iQRec = 0; iQRec < nAnyRec; ++iQRec) {
868 int iQ = iAnyRec[iQRec];
869 int idQ =
event[iQ].id();
870 if (freeAnyRec[iQRec] && idQ > 0 && idQ < 6)
871 for (
int iQbarRec = 0; iQbarRec < nAnyRec; ++iQbarRec) {
872 int iQbar = iAnyRec[iQbarRec];
873 if (freeAnyRec[iQbarRec] && event[iQbar].
id() == -idQ) {
877 int iTopQ =
event.iTopCopyId(iQ);
878 int iTopQbar =
event.iTopCopyId(iQbar);
879 int iMother =
event[iTopQ].mother1();
880 if (event[iTopQbar].mother1() == iMother
881 && event[iMother].isGluon() && event[iMother].status() != -34
882 && event[iMother + 1].status() != -34 ) {
886 Vec4 pGlu =
event[iQ].p() +
event[iQbar].p();
888 double pT2DipMin = sCM;
889 for (
int iDip = 0; iDip < int(dipoles.size()); ++iDip) {
890 double pT2Dip = (pGlu *
event[dipoles[iDip].iCol].p())
891 * (pGlu * event[dipoles[iDip].iAcol].p())
892 / dipoles[iDip].p1p2;
893 if (pT2Dip < pT2DipMin) {
900 int colGlu =
event[iQ].col();
901 int acolGlu =
event[iQbar].acol();
902 int colDip = dipoles[iDipMin].col;
903 int iColDip = dipoles[iDipMin].iCol;
904 int iAcolDip = dipoles[iDipMin].iAcol;
905 event[iQbar].acol( colDip );
906 if (event[iAcolDip].acol() == colDip)
907 event[iAcolDip].acol( colGlu );
908 else event[iAcolDip].col( colGlu );
909 dipoles[iDipMin].iAcol = iQbar;
910 dipoles[iDipMin].p1p2 =
event[iColDip].p() *
event[iQbar].p();
911 dipoles.push_back( BeamDipole( colGlu, iQ, iAcolDip ) );
912 dipoles.back().p1p2 =
event[iQ].p() *
event[iAcolDip].p();
915 freeAnyRec[iQRec] =
false;
916 freeAnyRec[iQbarRec] =
false;
917 for (
int i = oldSize; i <
event.size(); ++i)
918 if (i != iQRec && i != iQbarRec && i != iColDip
920 if (event[i].isFinal()) {
921 if (event[i].acol() == colGlu) event[i].acol( acolGlu );
923 if (event[i].col() == colGlu) event[i].col( acolGlu );
928 for (
int iJun = 0; iJun <
event.sizeJunction(); ++iJun) {
931 if (event.kindJunction(iJun) % 2 == 0)
continue;
932 for (
int leg = 0; leg < 3; ++leg) {
933 int col =
event.colJunction(iJun, leg);
935 event.colJunction(iJun, leg, colGlu);
947 if ( event[iInARec].isGluon() && !event[iInARec].isRescatteredIncoming()
948 && event[iInBRec].isGluon() && !event[iInBRec].isRescatteredIncoming()
949 && event[iInARec].col() == event[iInBRec].acol()
950 && event[iInARec].acol() == event[iInBRec].col() ) {
951 event[iInARec].acol( event[iInARec].col() );
952 event[iInBRec].acol( event[iInBRec].col() );
968 bool BeamRemnants::checkColours(
Event& event) {
971 if (beamAPtr->isLepton() && beamBPtr->isLepton())
return true;
975 for (
int iCol = 1; iCol < int(colFrom.size()); ++iCol)
976 for (
int iColRef = 0; iColRef < iCol; ++iColRef) {
977 if (colFrom[iCol] == colFrom[iColRef]) {
978 colFrom[iCol] = colTo[iCol];
979 colTo[iCol] = colTo[iColRef];
981 if (colTo[iCol] == colFrom[iColRef]) colTo[iCol] = colTo[iColRef];
985 for (
int i = oldSize; i <
event.size(); ++i) {
986 int col =
event[i].col();
987 int acol =
event[i].acol();
988 for (
int iCol = 0; iCol < int(colFrom.size()); ++iCol) {
989 if (col == colFrom[iCol]) {col = colTo[iCol];
event[i].col(col);}
990 if (acol == colFrom[iCol]) {acol = colTo[iCol];
event[i].acol(acol);}
995 for (
int iJun = 0; iJun <
event.sizeJunction(); ++iJun)
996 for (
int leg = 0; leg < 3; ++leg) {
997 int col =
event.colJunction(iJun, leg);
998 for (
int iCol = 0; iCol < int(colFrom.size()); ++iCol) {
999 if (col == colFrom[iCol]) {
1001 event.colJunction(iJun, leg, col);
1007 vector<int> colList;
1008 vector<int> acolList;
1009 vector<int> iSingletGluon;
1012 for (
int i = oldSize; i <
event.size(); ++i)
1013 if (event[i].isFinal()) {
1014 int id =
event[i].id();
1015 int col =
event[i].col();
1016 int acol =
event[i].acol();
1017 int colType =
event[i].colType();
1020 if ( (
id > 0 &&
id < 9 && (col <= 0 || acol != 0) )
1021 || (id < 0 && id > -9 && (col != 0 || acol <= 0) )
1022 || (
id == 21 && (col <= 0 || acol <= 0) ) ) {
1023 infoPtr->errorMsg(
"Error in BeamRemnants::checkColours: "
1024 "q/qbar/g has wrong colour slots set");
1029 if ( (colType == 3 && (col <= 0 || acol >= 0))
1030 || (colType == -3 && (col >= 0 || acol <= 0)) ) {
1031 infoPtr->errorMsg(
"Error in BeamRemnants::checkColours: "
1032 "sextet has wrong colours");
1036 if ( col > 0) colList.push_back( col );
1037 if (acol > 0) acolList.push_back( acol );
1038 if (col > 0 && acol == col) iSingletGluon.push_back(i);
1040 if ( col < 0) acolList.push_back( -col );
1041 if (acol < 0) colList.push_back( -acol );
1046 for (
int iS = 0; iS < int(iSingletGluon.size()); ++iS) {
1047 int iGlu = iSingletGluon[iS];
1049 double pT2DipMin = sCM;
1050 for (
int iC = oldSize; iC <
event.size(); ++iC)
1051 if (iC != iGlu && event[iC].isFinal()) {
1052 int colDip =
event[iC].col();
1053 if (colDip > 0 && event[iC].acol() !=colDip)
1054 for (
int iA = oldSize; iA <
event.size(); ++iA)
1055 if (iA != iGlu && iA != iC && event[iA].isFinal()
1056 &&
event[iA].acol() == colDip &&
event[iA].col() !=colDip) {
1057 double pT2Dip = (
event[iGlu].p() *
event[iC].p())
1058 * (event[iGlu].p() *
event[iA].p())
1059 / (event[iC].p() *
event[iA].p());
1060 if (pT2Dip < pT2DipMin) {
1068 if (iAcolDip == -1)
return false;
1069 event[iGlu].acol( event[iAcolDip].acol() );
1070 event[iAcolDip].acol( event[iGlu].col() );
1073 for (
int iJun = 0; iJun <
event.sizeJunction(); ++iJun) {
1076 if (event.kindJunction(iJun) % 2 == 0)
continue;
1077 for (
int leg = 0; leg < 3; ++leg) {
1078 int col =
event.colJunction(iJun, leg);
1079 if (col == event[iGlu].acol())
1080 event.colJunction(iJun, leg, event[iGlu].col());
1087 for (
int iCol = 0; iCol < int(colList.size()) - 1; ++iCol) {
1088 int col = colList[iCol];
1089 for (
int iCol2 = iCol + 1; iCol2 < int(colList.size()); ++iCol2)
1090 if (colList[iCol2] == col) {
1091 infoPtr->errorMsg(
"Warning in BeamRemnants::checkColours:"
1092 " colour appears twice");
1093 if (!ALLOWCOLOURTWICE)
return false;
1096 for (
int iAcol = 0; iAcol < int(acolList.size()) - 1; ++iAcol) {
1097 int acol = acolList[iAcol];
1098 for (
int iAcol2 = iAcol + 1; iAcol2 < int(acolList.size()); ++iAcol2)
1099 if (acolList[iAcol2] == acol) {
1100 infoPtr->errorMsg(
"Warning in BeamRemnants::checkColours:"
1101 " anticolour appears twice");
1102 if (!ALLOWCOLOURTWICE)
return false;
1107 bool foundPair =
true;
1108 while (foundPair && colList.size() > 0 && acolList.size() > 0) {
1110 for (
int iCol = 0; iCol < int(colList.size()); ++iCol) {
1111 for (
int iAcol = 0; iAcol < int(acolList.size()); ++iAcol) {
1112 if (acolList[iAcol] == colList[iCol]) {
1113 colList[iCol] = colList.back(); colList.pop_back();
1114 acolList[iAcol] = acolList.back(); acolList.pop_back();
1119 if (foundPair)
break;
1124 for (
int iJun = 0; iJun <
event.sizeJunction(); ++iJun) {
1125 int kindJun =
event.kindJunction(iJun);
1126 for (
int leg = 0; leg < 3; ++leg) {
1127 int colEnd =
event.colJunction(iJun, leg);
1131 for (
int iCol = 0; iCol < int(colList.size()); ++iCol)
1132 if (colList[iCol] == colEnd) {
1134 colList[iCol] = colList.back();
1141 else if (kindJun == 2) {
1142 for (
int iAcol = 0; iAcol < int(acolList.size()); ++iAcol)
1143 if (acolList[iAcol] == colEnd) {
1145 acolList[iAcol] = acolList.back();
1146 acolList.pop_back();
1152 else if ( kindJun == 3 || kindJun == 5) {
1153 for (
int iCol = 0; iCol < int(colList.size()); ++iCol)
1154 if (colList[iCol] == colEnd) {
1156 colList[iCol] = colList.back();
1163 else if ( kindJun == 4 || kindJun == 6) {
1164 for (
int iAcol = 0; iAcol < int(acolList.size()); ++iAcol)
1165 if (acolList[iAcol] == colEnd) {
1167 acolList[iAcol] = acolList.back();
1168 acolList.pop_back();
1179 if (colList.size() > 0 || acolList.size() > 0) {
1180 infoPtr->errorMsg(
"Warning in BeamRemnants::checkColours:"
1181 " need to repair unmatched colours");
1183 while (colList.size() > 0 && acolList.size() > 0) {
1186 int colMatch = colList.back();
1187 int acolMatch = acolList.back();
1188 int colNew =
event.nextColTag();
1190 acolList.pop_back();
1191 for (
int i = oldSize; i <
event.size(); ++i)
1192 if (event[i].isFinal() &&
event[i].col() == colMatch) {
1193 event[i].col( colNew);
1196 for (
int i = oldSize; i <
event.size(); ++i)
1197 if (event[i].isFinal() &&
event[i].acol() == acolMatch) {
1198 event[i].acol( colNew);
1204 return (colList.size() == 0 && acolList.size() == 0);