PYTHIA8/pythia8145/src/BeamRemnants.cxx

   1 // BeamRemnants.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 Torbjorn Sjostrand.
   3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
   4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
   5
   6 // Function definitions (not found in the header) for the
   7 // BeamRemnants class.
   8
   9 #include "BeamRemnants.h"
  10
  11 namespace Pythia8 {
  12
  13 //==========================================================================
  14
  15 // The BeamDipole class is purely internal to reconnectColours.
  16
  17 class BeamDipole {
  18
  19 public:
  20
  21   // Constructor.
  22   BeamDipole( int colIn = 0, int iColIn = 0, int iAcolIn = 0)
  23     : col(colIn), iCol(iColIn), iAcol(iAcolIn) {}
  24
  25   // Members.
  26   int    col, iCol, iAcol;
  27   double p1p2;
  28
  29 };
  30
  31 //==========================================================================
  32
  33 // The BeamRemnants class.
  34
  35 //--------------------------------------------------------------------------
  36
  37 // Constants: could be changed here if desired, but normally should not.
  38 // These are of technical nature, as described for each.
  39
  40 // If same (anti)colour appears twice in final state, repair or reject.
  41 const bool   BeamRemnants::ALLOWCOLOURTWICE = true;
  42
  43 // Maximum number of tries to match colours and kinematics in the event.
  44 const int    BeamRemnants::NTRYCOLMATCH     = 10;
  45 const int    BeamRemnants::NTRYKINMATCH     = 10;
  46
  47 // Overall correction step for energy-momentum conservation; only
  48 // becomes relevant in rescattering scenarios when FSR dipole emissions
  49 // and primordial kT is added. Should hopefully not be needed.
  50 const bool   BeamRemnants::CORRECTMISMATCH  = false;
  51
  52 //--------------------------------------------------------------------------
  53
  54 // Initialization.
  55
  56 bool BeamRemnants::init( Info* infoPtrIn, Settings& settings, Rndm* rndmPtrIn,
  57   BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
  58   PartonSystems* partonSystemsPtrIn) {
  59
  60   // Save pointers.
  61   infoPtr             = infoPtrIn;
  62   rndmPtr             = rndmPtrIn;
  63   beamAPtr            = beamAPtrIn;
  64   beamBPtr            = beamBPtrIn;
  65   partonSystemsPtr    = partonSystemsPtrIn;
  66
  67   // Width of primordial kT distribution.
  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");
  74
  75   // Handling of rescattering kinematics uncertainties from primodial kT.
  76   allowRescatter      = settings.flag("MultipleInteractions:allowRescatter");
  77   doRescatterRestoreY = settings.flag("BeamRemnants:rescatterRestoreY");
  78
  79   // Parameters for colour reconnection scenario, partly borrowed from
  80   // multiple interactions not to introduce too many new ones.
  81   doReconnect         = settings.flag("BeamRemnants:reconnectColours");
  82   reconnectRange      = settings.parm("BeamRemnants:reconnectRange");
  83   pT0Ref              = settings.parm("MultipleInteractions:pT0Ref");
  84   ecmRef              = settings.parm("MultipleInteractions:ecmRef");
  85   ecmPow              = settings.parm("MultipleInteractions:ecmPow");
  86
  87   // Total and squared CM energy at nominal energy.
  88   eCM                 = infoPtr->eCM();
  89   sCM                 = eCM * eCM;
  90
  91   // The MI pT0 smoothening scale and its reconnection-strength combination.
  92   pT0                 = pT0Ref * pow(eCM / ecmRef, ecmPow);
  93   pT20Rec             = pow2(reconnectRange * pT0);
  94
  95   // Done.
  96   return true;
  97
  98 }
  99
 100 //--------------------------------------------------------------------------
 101
 102 // Select the flavours/kinematics/colours of the two beam remnants.
 103 // Notation: iPar = all partons, iSys = matched systems of two beams,
 104 //           iRem = additional partons in remnants.
 105
 106 bool BeamRemnants::add( Event& event) {
 107
 108   // Update to current CM energy.
 109   eCM     = infoPtr->eCM();
 110   sCM     = eCM * eCM;
 111
 112   // Check that flavour bookkept in event and in beam particles agree.
 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");
 118       return false;
 119     }
 120   }
 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");
 126       return false;
 127     }
 128   }
 129
 130   // Number of scattering subsystems. Size of event record before treatment.
 131   nSys    = partonSystemsPtr->sizeSys();
 132   oldSize = event.size();
 133
 134   // Add required extra remnant flavour content.
 135   // Start all over if fails (in option where junctions not allowed).
 136   if ( !beamAPtr->remnantFlavours(event)
 137     || !beamBPtr->remnantFlavours(event) ) {
 138     infoPtr->errorMsg("Error in BeamRemnants::add:"
 139       " remnant flavour setup failed");
 140     return false;
 141   }
 142
 143   // Do the kinematics of the collision subsystems and two beam remnants.
 144   if (!setKinematics(event)) return false;
 145
 146   // Allow colour reconnections.
 147   if (doReconnect) reconnectColours(event);
 148
 149   // Save current modifiable colour configuration for fast restoration.
 150   vector<int> colSave;
 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() );
 155   }
 156   event.saveJunctionSize();
 157
 158   // Allow several tries to match colours of initiators and remnants.
 159   // Frequent "failures" since shortcutting colours separately on
 160   // the two event sides may give "colour singlet gluons" etc.
 161   bool physical = true;
 162   for (int iTry = 0; iTry < NTRYCOLMATCH; ++iTry) {
 163     physical = true;
 164
 165     // Reset list of colour "collapses" (transformations).
 166     colFrom.resize(0);
 167     colTo.resize(0);
 168
 169     // First process each set of beam colours on its own.
 170     if (!beamAPtr->remnantColours(event, colFrom, colTo))
 171       physical = false;
 172     if (!beamBPtr->remnantColours(event, colFrom, colTo))
 173       physical = false;
 174
 175     // Then check that colours and anticolours are matched in whole event.
 176     if ( physical && !checkColours(event) ) physical = false;
 177
 178     // If no problems then done, else restore and loop.
 179     if (physical) break;
 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");
 185   }
 186
 187   // If no solution after several tries then failed.
 188   if (!physical) {
 189     infoPtr->errorMsg("Error in BeamRemnants::add:"
 190       " colour tracing failed after several attempts");
 191     return false;
 192   }
 193
 194   // Done.
 195   return true;
 196
 197 }
 198
 199 //--------------------------------------------------------------------------
 200
 201 // Set up trial transverse and longitudinal kinematics for each beam
 202 // separately. Final decisions involve comparing the two beams.
 203
 204 bool BeamRemnants::setKinematics( Event& event) {
 205
 206   // References to beams to simplify indexing.
 207   BeamParticle& beamA = *beamAPtr;
 208   BeamParticle& beamB = *beamBPtr;
 209
 210   // Nothing to do for lepton-lepton scattering with all energy already used.
 211   if ( beamA.isUnresolvedLepton() && beamB.isUnresolvedLepton() )
 212     return true;
 213
 214   // Check that has not already used up beams.
 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");
 219     return false;
 220   }
 221
 222   // Last beam-status particles. Offset relative to normal beam locations.
 223   int nBeams   = 3;
 224   for (int i = 3; i < event.size(); ++i)
 225     if (event[i].statusAbs() < 20) nBeams = i + 1;
 226   int nOffset  = nBeams - 3;
 227
 228   // Reserve space for extra information on the systems and beams.
 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);
 234
 235   // Loop over all subsystems. Default values. Find invariant mass.
 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();
 244
 245     // Allow primordial kT reduction for small-mass and small-pT systems
 246     // (for hardest interaction pT -> renormalization scale so also 2 -> 1).
 247     if (doPrimordialKT) {
 248       double mHat     = sqrt(sHatNow);
 249       mHatDamp        = mHat / (mHat + halfMassForKT);
 250       double scale    = (iSys == 0) ? infoPtr->QRen() : infoPtr->pTMI(iSys);
 251       kTwidthNow      = ( (halfScaleForKT * primordialKTsoft
 252       + scale * primordialKThard) / (halfScaleForKT + scale) ) * mHatDamp;
 253     }
 254
 255     // Store properties of compensation systems and total compensation power.
 256     // Rescattered partons do not compensate, but may be massive.
 257     sHatSys[iSys] = sHatNow;
 258     kTwidth[iSys] = kTwidthNow ;
 259     kTcomp[iSys]  = mHatDamp;
 260     if (beamA[iSys].isFromBeam()) kTcompSumA += mHatDamp;
 261     else beamA[iSys].m( event[iInA].m() );
 262     if (beamB[iSys].isFromBeam()) kTcompSumB += mHatDamp;
 263     else beamB[iSys].m( event[iInB].m() );
 264   }
 265
 266   // Primordial kT and compensation power among remnants.
 267   double kTwidthNow = (doPrimordialKT) ? primordialKTremnant : 0.;
 268   for (int iRem = nSys; iRem < nMaxBeam; ++iRem) {
 269     sHatSys[iRem] = 0.;
 270     kTwidth[iRem] = kTwidthNow ;
 271     kTcomp[iRem]  = 1.;
 272   }
 273   kTcompSumA += beamA.size() - nSys;
 274   kTcompSumB += beamB.size() - nSys;
 275
 276   // Allow ten tries to construct kinematics (but normally works first).
 277   bool physical;
 278   double xSum[2], xInvM[2], w2Beam[2], wPosRem, wNegRem, w2Rem;
 279   for (int iTry = 0; iTry < NTRYKINMATCH; ++iTry) {
 280     physical = true;
 281
 282     // Loop over the two beams.
 283     for (int iBeam = 0; iBeam < 2; ++iBeam) {
 284       BeamParticle& beam = (iBeam == 0) ? beamA : beamB;
 285       int nPar = beam.size();
 286
 287       // Generate Gaussian pT for initiator partons inside hadrons.
 288       // Store/restore rescattered parton momenta before primordial kT.
 289       if (beam.isHadron() && doPrimordialKT) {
 290         double pxSum = 0.;
 291         double pySum = 0.;
 292         for (int iPar = 0; iPar < nPar; ++iPar) {
 293           if ( beam[iPar].isFromBeam() ) {
 294             pair<double, double> gauss2 = rndmPtr->gauss2();
 295             double px = kTwidth[iPar] * gauss2.first;
 296             double py = kTwidth[iPar] * gauss2.second;
 297             beam[iPar].px(px);
 298             beam[iPar].py(py);
 299             pxSum += px;
 300             pySum += py;
 301           } else {
 302             int iInAB = (iBeam == 0) ? partonSystemsPtr->getInA(iPar)
 303                                      : partonSystemsPtr->getInB(iPar);
 304             beam[iPar].p( event[iInAB].p() );
 305           }
 306         }
 307
 308         // Share recoil between all initiator partons, rescatterers excluded.
 309         double kTcompSum = (iBeam == 0) ? kTcompSumA : kTcompSumB;
 310         for (int iPar = 0; iPar < nPar; ++iPar)
 311         if (beam[iPar].isFromBeam() ) {
 312           beam[iPar].px( beam[iPar].px() - pxSum * kTcomp[iPar] / kTcompSum );
 313           beam[iPar].py( beam[iPar].py() - pySum * kTcomp[iPar] / kTcompSum );
 314         }
 315
 316       // Without primordial kT: still need to store rescattered partons.
 317       } else if (beam.isHadron()) {
 318         for (int iPar = 0; iPar < nPar; ++iPar)
 319         if ( !beam[iPar].isFromBeam() ) {
 320           int iInAB = (iBeam == 0) ? partonSystemsPtr->getInA(iPar)
 321                                    : partonSystemsPtr->getInB(iPar);
 322           beam[iPar].p( event[iInAB].p() );
 323         }
 324       }
 325
 326       // Pick unrescaled x values for remnants. Sum up (unscaled) p+ and p-.
 327       xSum[iBeam]  = 0.;
 328       xInvM[iBeam] = 0.;
 329       for (int iRem = nSys; iRem < nPar; ++iRem) {
 330         double xPrel = beam.xRemnant( iRem);
 331         beam[iRem].x(xPrel);
 332         xSum[iBeam]  += xPrel;
 333         xInvM[iBeam] += beam[iRem].mT2()/xPrel;
 334       }
 335
 336       // Squared transverse mass for each beam, using lightcone x.
 337       w2Beam[iBeam] = xSum[iBeam] * xInvM[iBeam];
 338
 339     // End separate treatment of the two beams.
 340     }
 341
 342     // Recalculate kinematics of initiator systems with primordial kT.
 343     wPosRem = eCM;
 344     wNegRem = eCM;
 345     for (int iSys = 0; iSys < nSys; ++iSys) {
 346       int iA          = beamA[iSys].iPos();
 347       int iB          = beamB[iSys].iPos();
 348       double sHat     = sHatSys[iSys];
 349
 350       // Classify system: rescattering on either or both sides?
 351       bool normalA    = beamA[iSys].isFromBeam();
 352       bool normalB    = beamB[iSys].isFromBeam();
 353       bool normalSys  = normalA && normalB;
 354       bool doubleRes  = !normalA && !normalB;
 355
 356       // Check whether final-state system momentum matches initial-state one.
 357       if (allowRescatter && CORRECTMISMATCH) {
 358         Vec4 pInitial = event[iA].p() + event[iB].p();
 359         Vec4 pFinal;
 360         for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) {
 361           int iAB      = partonSystemsPtr->getOut(iSys, iMem);
 362           if (event[iAB].isFinal()) pFinal += event[iAB].p();
 363         }
 364
 365         // Scale down primordial kT from side A if p+ increased.
 366         if (normalA && pFinal.pPos() > pInitial.pPos())
 367           beamA[iSys].scalePT( pInitial.pPos() / pFinal.pPos() );
 368
 369         // Scale down primordial kT from side B if p- increased.
 370         if (normalB && pFinal.pNeg() > pInitial.pNeg())
 371           beamB[iSys].scalePT( pInitial.pNeg() / pFinal.pNeg() );
 372       }
 373
 374       // Rescatter: possible change in sign of lightcone momentum of a
 375       //            rescattered parton. If this happens, try to pick
 376       //            new primordial kT values
 377       if (allowRescatter
 378          && (event[iA].pPos() / beamA[iSys].pPos() < 0
 379          ||  event[iB].pNeg() / beamB[iSys].pNeg() < 0) ) {
 380         infoPtr->errorMsg("Warning in BeamRemnants::setKinematics:"
 381           " change in lightcone momentum sign; retrying kinematics");
 382         physical = false;
 383         break;
 384       }
 385
 386       // Begin kinematics of partons after primordial kT has been added.
 387       double sHatTAft = sHat + pow2( beamA[iSys].px() + beamB[iSys].px())
 388                              + pow2( beamA[iSys].py() + beamB[iSys].py());
 389       double w2A      = beamA[iSys].mT2();
 390       double w2B      = beamB[iSys].mT2();
 391       double w2Diff   = sHatTAft - w2A - w2B;
 392       double lambda   = pow2(w2Diff) - 4. * w2A * w2B;
 393
 394       // Too large transverse momenta means that kinematics will not work.
 395       if (lambda <= 0.) {
 396         physical      = false;
 397         break;
 398       }
 399       double lamRoot  = sqrtpos( lambda );
 400
 401       // Mirror solution if the two incoming have reverse rapidity ordering.
 402       if (allowRescatter && doubleRes && (event[iA].pPos() * event[iB].pNeg()
 403         < event[iA].pNeg() * event[iB].pPos()) ) lamRoot = -lamRoot;
 404
 405       // Two procedures, which agree for normal scattering, separate here.
 406       // First option keeps rapidity (and mass) of system unchanged by
 407       // primordial kT, by modification of rescattered parton.
 408       if (normalSys || doRescatterRestoreY || doubleRes) {
 409
 410         // Find kinematics of initial system: transverse mass, and
 411         // longitudinal momentum carried by all or rescattered partons.
 412         double sHatTBef = sHat;
 413         double wPosBef, wNegBef, wPosBefRes, wNegBefRes;
 414         // Normal scattering.
 415         if (normalSys) {
 416           wPosBef       = beamA[iSys].x() * eCM;
 417           wNegBef       = beamB[iSys].x() * eCM;
 418           wPosBefRes    = 0.;
 419           wNegBefRes    = 0.;
 420         // Rescattering on side A.
 421         } else if (normalB) {
 422           sHatTBef     += event[iA].pT2();
 423           wPosBef       = event[iA].pPos();
 424           wNegBef       = event[iA].pNeg() + beamB[iSys].x() * eCM;
 425           wPosBefRes    = beamA[iSys].pPos();
 426           wNegBefRes    = beamA[iSys].pNeg();
 427         // Rescattering on side B.
 428         } else if (normalA) {
 429           sHatTBef     += event[iB].pT2();
 430           wPosBef       = beamA[iSys].x() * eCM + event[iB].pPos();
 431           wNegBef       = event[iB].pNeg();
 432           wPosBefRes    = beamB[iSys].pPos();
 433           wNegBefRes    = beamB[iSys].pNeg();
 434         // Rescattering on both sides.
 435         } else {
 436           sHatTBef     += pow2( event[iA].px() + event[iB].px())
 437                         + pow2( event[iA].py() + event[iB].py());
 438           wPosBef       = event[iA].pPos() + event[iB].pPos();
 439           wNegBef       = event[iA].pNeg() + event[iB].pNeg();
 440           wPosBefRes    = beamA[iSys].pPos() + beamB[iSys].pPos();
 441           wNegBefRes    = beamA[iSys].pNeg() + beamB[iSys].pNeg();
 442         }
 443
 444         // Rescale outgoing momenta to keep same mass and rapidity of system
 445         // as originally, and solve for kinematics.
 446         double rescale  = sqrt(sHatTAft / sHatTBef);
 447         double wPosAft  = rescale * wPosBef;
 448         double wNegAft  = rescale * wNegBef;
 449         wPosRem        -= wPosAft - wPosBefRes;
 450         wNegRem        -= wNegAft - wNegBefRes;
 451         double wPosA    = 0.5 * (sHatTAft + w2A - w2B + lamRoot) / wNegAft;
 452         double wNegB    = 0.5 * (sHatTAft + w2B - w2A + lamRoot) / wPosAft;
 453
 454         // Store modified beam parton momenta.
 455         beamA[iSys].e(  0.5 * (wPosA + w2A / wPosA) );
 456         beamA[iSys].pz( 0.5 * (wPosA - w2A / wPosA) );
 457         beamB[iSys].e(  0.5 * (w2B / wNegB + wNegB) );
 458         beamB[iSys].pz( 0.5 * (w2B / wNegB - wNegB) );
 459
 460       // Second option keeps rescattered parton (and system mass) unchanged
 461       // by primordial kT, by modification of system rapidity.
 462       } else {
 463
 464         // Rescattering on side A: preserve already scattered parton.
 465         if (normalB) {
 466           double wPosA  = beamA[iSys].pPos();
 467           double wNegB  = 0.5 * (w2Diff + lamRoot) / wPosA;
 468           beamB[iSys].e(  0.5 * (w2B / wNegB + wNegB) );
 469           beamB[iSys].pz( 0.5 * (w2B / wNegB - wNegB) );
 470           wPosRem      -= w2B / wNegB;
 471           wNegRem      -= wNegB;
 472
 473
 474         // Rescattering on side B: preserve already scattered parton.
 475         } else if (normalA) {
 476           double wNegB  = beamB[iSys].pNeg();
 477           double wPosA  = 0.5 * (w2Diff + lamRoot) / wNegB;
 478           beamA[iSys].e(  0.5 * (wPosA + w2A / wPosA) );
 479           beamA[iSys].pz( 0.5 * (wPosA - w2A / wPosA) );
 480           wPosRem      -= wPosA;
 481           wNegRem      -= w2A / wPosA;
 482
 483         // Primordial kT in double rescattering does change the mass of
 484         // the system without any possible fix in this framework, which
 485         // leads to incorrect boosts. Current choice is therefore always
 486         // to handle it with the first procedure, where mass is retained.
 487         } else {
 488         }
 489       }
 490
 491       // Construct system rotation and boost caused by primordial kT.
 492       Msys[iSys].reset();
 493       Msys[iSys].toCMframe( event[iA].p(), event[iB].p() );
 494       Msys[iSys].fromCMframe( beamA[iSys].p(), beamB[iSys].p() );
 495
 496       // Boost rescattered partons in subsequent beam A list.
 497       for (int iSys2 = iSys + 1; iSys2 < nSys; ++iSys2) {
 498         if (!beamA[iSys2].isFromBeam()) {
 499           int iBefResc = event[ beamA[iSys2].iPos() ].mother1();
 500           for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem)
 501           if (partonSystemsPtr->getOut(iSys, iMem) == iBefResc) {
 502             Vec4 pTemp = event[iBefResc].p();
 503             pTemp.rotbst( Msys[iSys] );
 504             beamA[iSys2].p( pTemp );
 505           }
 506         }
 507
 508         // Boost rescattered partons in subsequent beam B list.
 509         if (!beamB[iSys2].isFromBeam()) {
 510           int iBefResc = event[ beamB[iSys2].iPos() ].mother1();
 511           for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem)
 512           if (partonSystemsPtr->getOut(iSys, iMem) == iBefResc) {
 513             Vec4 pTemp = event[iBefResc].p();
 514             pTemp.rotbst( Msys[iSys] );
 515             beamB[iSys2].p( pTemp );
 516           }
 517         }
 518       }
 519     }
 520
 521     // Check that remaining momentum is enough for remnants.
 522     if (wPosRem < 0. || wNegRem < 0.) physical = false;
 523     w2Rem = wPosRem * wNegRem;
 524     if (sqrtpos(w2Rem) < sqrt(w2Beam[0]) + sqrt(w2Beam[1]))
 525       physical = false;
 526
 527     // End of loop over ten tries. Do not loop when solution found.
 528     if (physical) break;
 529   }
 530
 531   // If no solution after ten tries then failed.
 532   if (!physical) {
 533     infoPtr->errorMsg("Error in BeamRemnants::setKinematics:"
 534       " kinematics construction failed");
 535     return false;
 536   }
 537
 538   // For successful initiator kinematics process whole systems.
 539   Vec4 pSumOut;
 540   for (int iSys = 0; iSys < nSys; ++iSys) {
 541
 542     // Copy initiators and their systems and boost them accordingly.
 543     // Update subsystem and beams info on new positions of partons.
 544     // Update daughter info of mothers, i.e. of beams, for hardest interaction.
 545     if (beamA[iSys].isFromBeam()) {
 546       int iA       = beamA[iSys].iPos();
 547       int iAcopy   = event.copy(iA, -61);
 548       event[iAcopy].rotbst(Msys[iSys]);
 549       partonSystemsPtr->setInA(iSys, iAcopy);
 550       beamA[iSys].iPos( iAcopy);
 551       if (iSys == 0) {
 552         int mother = event[iAcopy].mother1();
 553         event[mother].daughter1(iAcopy);
 554       }
 555     }
 556     if (beamB[iSys].isFromBeam()) {
 557       int iB       = beamB[iSys].iPos();
 558       int iBcopy   = event.copy(iB, -61);
 559       event[iBcopy].rotbst(Msys[iSys]);
 560       partonSystemsPtr->setInB(iSys, iBcopy);
 561       beamB[iSys].iPos( iBcopy);
 562       if (iSys == 0) {
 563         int mother = event[iBcopy].mother1();
 564         event[mother].daughter1(iBcopy);
 565       }
 566     }
 567
 568     for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) {
 569       int iAB      = partonSystemsPtr->getOut(iSys, iMem);
 570       if (event[iAB].isFinal()) {
 571         int iABcopy = event.copy(iAB, 62);
 572         event[iABcopy].rotbst(Msys[iSys]);
 573         partonSystemsPtr->setOut(iSys, iMem, iABcopy);
 574         pSumOut   += event[iABcopy].p();
 575       }
 576     }
 577
 578   }
 579
 580   // Colour dipoles spanning systems gives mismatch between FSR recoils
 581   // and primordial kT boosts.
 582   if (allowRescatter && CORRECTMISMATCH) {
 583
 584     // Find summed pT of beam remnants = - wanted pT of systems.
 585     double pxBeams = 0.;
 586     double pyBeams = 0.;
 587     for (int iRem = nSys; iRem < beamA.size(); ++iRem) {
 588       pxBeams     += beamA[iRem].px();
 589       pyBeams     += beamA[iRem].py();
 590     }
 591     for (int iRem = nSys; iRem < beamB.size(); ++iRem) {
 592       pxBeams     += beamB[iRem].px();
 593       pyBeams     += beamB[iRem].py();
 594     }
 595
 596     // Boost all final partons in systems transversely, and also their sum.
 597     Vec4 pSumTo( -pxBeams, -pyBeams, pSumOut.pz(), sqrt( pow2(pxBeams)
 598       + pow2(pyBeams) + pow2(pSumOut.pz()) + pSumOut.m2Calc() ) );
 599     RotBstMatrix Mmismatch;
 600     Mmismatch.bst( pSumOut, pSumTo);
 601     for (int iSys = 0; iSys < nSys; ++iSys)
 602     for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) {
 603       int iAB = partonSystemsPtr->getOut(iSys, iMem);
 604       if (event[iAB].isFinal()) event[iAB].rotbst(Mmismatch);
 605     }
 606     pSumOut.rotbst(Mmismatch);
 607
 608    // Reset energy and momentum sum, to be compensated by beam remnants.
 609     wPosRem = eCM - (pSumOut.e() + pSumOut.pz());
 610     wNegRem = eCM - (pSumOut.e() - pSumOut.pz());
 611     w2Rem    = wPosRem * wNegRem;
 612     if ( wPosRem < 0. || wNegRem < 0.
 613       || sqrtpos(w2Rem) < sqrt(w2Beam[0]) + sqrt(w2Beam[1])) {
 614       infoPtr->errorMsg("Error in BeamRemnants::setKinematics:"
 615         " kinematics construction failed owing to recoil mismatch");
 616       return false;
 617     }
 618   }
 619
 620   // Construct x rescaling factors for the two remants.
 621   double lambdaRoot = sqrtpos( pow2(w2Rem - w2Beam[0] - w2Beam[1])
 622     - 4. * w2Beam[0] * w2Beam[1] );
 623   double rescaleA   = (w2Rem + w2Beam[0] - w2Beam[1] + lambdaRoot)
 624     / (2. * w2Rem * xSum[0]) ;
 625   double rescaleB   = (w2Rem + w2Beam[1] - w2Beam[0] + lambdaRoot)
 626     / (2. * w2Rem * xSum[1]) ;
 627
 628   // Construct energy and pz for remnants in first beam.
 629   for (int iRem = nSys; iRem < beamA.size(); ++iRem) {
 630     double pPos = rescaleA * beamA[iRem].x() * wPosRem;
 631     double pNeg = beamA[iRem].mT2() / pPos;
 632     beamA[iRem].e( 0.5 * (pPos + pNeg) );
 633     beamA[iRem].pz( 0.5 * (pPos - pNeg) );
 634
 635     // Add these partons to the normal event record.
 636     int iNew = event.append( beamA[iRem].id(), 63, 1 + nOffset, 0, 0, 0,
 637       beamA[iRem].col(), beamA[iRem].acol(), beamA[iRem].p(),
 638       beamA[iRem].m() );
 639     beamA[iRem].iPos( iNew);
 640   }
 641
 642   // Construct energy and pz for remnants in second beam.
 643   for (int iRem = nSys; iRem < beamB.size(); ++iRem) {
 644     double pNeg = rescaleB * beamB[iRem].x() * wNegRem;
 645     double pPos = beamB[iRem].mT2() / pNeg;
 646     beamB[iRem].e( 0.5 * (pPos + pNeg) );
 647     beamB[iRem].pz( 0.5 * (pPos - pNeg) );
 648
 649     // Add these partons to the normal event record.
 650     int iNew = event.append( beamB[iRem].id(), 63, 2 + nOffset, 0, 0, 0,
 651       beamB[iRem].col(), beamB[iRem].acol(), beamB[iRem].p(),
 652       beamB[iRem].m() );
 653     beamB[iRem].iPos( iNew);
 654   }
 655
 656   // Done.
 657   return true;
 658
 659 }
 660
 661 //--------------------------------------------------------------------------
 662
 663 // Allow colour reconnections by mergings of collision subsystems.
 664 // iRec is system that may be reconnected, by moving its gluons to iSys,
 665 // where minimal pT (or equivalently Lambda) is used to pick location.
 666 // Therefore all dipoles in iSys have to be found, and all gluons in iRec.
 667 // Matching q-qbar pairs are treated by analogy with gluons.
 668 // Note: owing to rescatterings some outgoing partons must be skipped.
 669
 670 bool BeamRemnants::reconnectColours( Event&  event) {
 671
 672   // References to beams to simplify indexing.
 673   BeamParticle& beamA = *beamAPtr;
 674   BeamParticle& beamB = *beamBPtr;
 675
 676   // Prepare record of which systems should be merged onto another.
 677   // The iSys system must have colour in final state to attach to it.
 678   vector<int>  iMerge(nSys);
 679   vector<bool> hasColour(nSys);
 680   for (int iSys = 0; iSys < nSys; ++iSys) {
 681     iMerge[iSys] = iSys;
 682     bool hasCol = false;
 683     for (int iMem = 0; iMem < partonSystemsPtr->sizeOut(iSys); ++iMem) {
 684       int iNow = partonSystemsPtr->getOut( iSys, iMem);
 685       if (event[iNow].isFinal() && (event[iNow].col() > 0
 686         || event[iNow].acol() > 0) ) {
 687         hasCol = true;
 688         break;
 689       }
 690     }
 691     hasColour[iSys] = hasCol;
 692   }
 693
 694   // Loop over systems to decide which should be reconnected.
 695   for (int iRec = nSys - 1; iRec > 0; --iRec) {
 696
 697     // Determine reconnection strength from pT scale of system.
 698     double pT2Rec  = pow2( infoPtr->pTMI(iRec) );
 699     double probRec = pT20Rec / (pT20Rec + pT2Rec);
 700
 701     // Loop over other systems iSys at higher pT scale and
 702     // decide whether to reconnect the iRec gluons onto one of them.
 703     for (int iSys = iRec - 1; iSys >= 0; --iSys)
 704     if (hasColour[iSys] && probRec > rndmPtr->flat()) {
 705
 706       // The iRec system and all merged with it to be merged with iSys.
 707       iMerge[iRec] = iSys;
 708       for (int iRec2 = iRec + 1; iRec2 < nSys; ++iRec2)
 709       if (iMerge[iRec2] == iRec) iMerge[iRec2] = iSys;
 710
 711       // Once a system has been merged do not test it anymore.
 712       break;
 713     }
 714   }
 715
 716   // Loop over systems. Check whether other systems to be merged with it.
 717   for (int iSys = 0; iSys < nSys; ++iSys) {
 718     int nMerge = 0;
 719     for (int iRec = iSys + 1; iRec < nSys; ++iRec)
 720     if (iMerge[iRec] == iSys) ++nMerge;
 721     if (nMerge == 0) continue;
 722
 723     // Incoming partons not counted if rescattered.
 724     int  iInASys = partonSystemsPtr->getInA(iSys);
 725     bool hasInA  = (beamA[iSys].isFromBeam());
 726     int  iInBSys = partonSystemsPtr->getInB(iSys);
 727     bool hasInB  = (beamB[iSys].isFromBeam());
 728
 729     // Begin find dipoles in iSys system.
 730     vector<BeamDipole> dipoles;
 731     int sizeOut = partonSystemsPtr->sizeOut(iSys);
 732     for (int iMem = 0; iMem < sizeOut; ++iMem) {
 733
 734       // Find colour dipoles to beam remnant.
 735       int iNow = partonSystemsPtr->getOut( iSys, iMem);
 736       if (!event[iNow].isFinal()) continue;
 737       int col = event[iNow].col();
 738       if (col > 0) {
 739         if      (hasInA && event[iInASys].col() == col)
 740           dipoles.push_back( BeamDipole( col, iNow, iInASys ) );
 741         else if (hasInB && event[iInBSys].col() == col)
 742           dipoles.push_back( BeamDipole( col, iNow, iInBSys ) );
 743
 744         // Find colour dipole between final-state partons.
 745         else for (int iMem2 = 0; iMem2 < sizeOut; ++iMem2)
 746         if (iMem2 != iMem) {
 747           int iNow2 = partonSystemsPtr->getOut( iSys, iMem2);
 748           if (!event[iNow2].isFinal()) continue;
 749           if (event[iNow2].acol() == col) {
 750             dipoles.push_back( BeamDipole( col, iNow, iNow2) );
 751             break;
 752           }
 753         }
 754       }
 755
 756       // Find anticolour dipoles to beam remnant.
 757       int acol = event[iNow].acol();
 758       if (acol > 0) {
 759         if      (hasInA && event[iInASys].acol() == acol)
 760           dipoles.push_back( BeamDipole( acol, iInASys, iNow ) );
 761         else if (hasInB && event[iInBSys].acol() == acol)
 762           dipoles.push_back( BeamDipole( acol, iInBSys, iNow ) );
 763       }
 764     }
 765
 766     // Skip mergings if no dipoles found.
 767     if (dipoles.size() == 0) continue;
 768
 769     // Find dipole sizes.
 770     for (int iDip = 0; iDip < int(dipoles.size()); ++iDip)
 771       dipoles[iDip].p1p2 = event[dipoles[iDip].iCol].p()
 772                          * event[dipoles[iDip].iAcol].p();
 773
 774     // Loop over systems iRec to be merged with iSys.
 775     for (int iRec = iSys + 1; iRec < nSys; ++iRec) {
 776       if (iMerge[iRec] != iSys) continue;
 777
 778       // Information on iRec. Vectors for gluons and anything else.
 779       int sizeRec = partonSystemsPtr->sizeOut(iRec);
 780       int iInARec = partonSystemsPtr->getInA(iRec);
 781       int iInBRec = partonSystemsPtr->getInB(iRec);
 782       int nGluRec = 0;
 783       vector<int>    iGluRec;
 784       vector<double> pT2GluRec;
 785       int nAnyRec = 0;
 786       vector<int>    iAnyRec;
 787       vector<bool>   freeAnyRec;
 788
 789       // Copy of gluon positions in descending order.
 790       for (int iMem = 0; iMem < sizeRec; ++iMem) {
 791         int iNow = partonSystemsPtr->getOut( iRec, iMem);
 792         if (!event[iNow].isFinal()) continue;
 793         if (event[iNow].isGluon()) {
 794           ++nGluRec;
 795           iGluRec.push_back( iNow );
 796           pT2GluRec.push_back( event[iNow].pT2() );
 797           for (int i = nGluRec - 1; i > 1; --i) {
 798             if (pT2GluRec[i - 1] > pT2GluRec[i]) break;
 799             swap(   iGluRec[i - 1],   iGluRec[i] );
 800             swap( pT2GluRec[i - 1], pT2GluRec[i] );
 801           }
 802         // Copy of anything else, mainly quarks, in no particular order.
 803         } else {
 804           ++nAnyRec;
 805           iAnyRec.push_back( iNow );
 806           freeAnyRec.push_back( true );
 807         }
 808       }
 809
 810       // For each gluon in iRec now find the dipole that gives the smallest
 811       // (pGlu * pI) (pGlu * pJ) / (pI * pJ), i.e. minimal pT (and Lambda).
 812       for (int iGRec = 0; iGRec < nGluRec; ++iGRec) {
 813         int    iGlu      = iGluRec[iGRec];
 814         Vec4   pGlu      = event[iGlu].p();
 815         int    iDipMin   = 0;
 816         double pT2DipMin = sCM;
 817         for (int iDip = 0; iDip < int(dipoles.size()); ++iDip) {
 818           double pT2Dip = (pGlu * event[dipoles[iDip].iCol].p())
 819             * (pGlu * event[dipoles[iDip].iAcol].p()) / dipoles[iDip].p1p2;
 820           if (pT2Dip < pT2DipMin) {
 821             iDipMin   = iDip;
 822             pT2DipMin = pT2Dip;
 823           }
 824         }
 825
 826         // Attach the gluon to the dipole, i.e. split the dipole in two.
 827         int colGlu   = event[iGlu].col();
 828         int acolGlu  = event[iGlu].acol();
 829         int colDip   = dipoles[iDipMin].col;
 830         int iColDip  = dipoles[iDipMin].iCol;
 831         int iAcolDip = dipoles[iDipMin].iAcol;
 832         event[iGlu].acol( colDip );
 833         if (event[iAcolDip].acol() == colDip)
 834              event[iAcolDip].acol( colGlu );
 835         else event[iAcolDip].col(  colGlu );
 836         dipoles[iDipMin].iAcol = iGlu;
 837         dipoles[iDipMin].p1p2 = event[iColDip].p() * pGlu;
 838         dipoles.push_back( BeamDipole( colGlu, iGlu, iAcolDip ) );
 839         dipoles.back().p1p2 = pGlu * event[iAcolDip].p();
 840
 841         // Remove gluon from old system: reconnect colours.
 842         for (int i = oldSize; i < event.size(); ++i)
 843         if (i != iGlu && i != iAcolDip) {
 844           if (event[i].isFinal()) {
 845             if (event[i].acol() == colGlu) event[i].acol( acolGlu );
 846           } else {
 847               if (event[i].col()  == colGlu) event[i].col( acolGlu );
 848           }
 849         }
 850       }
 851
 852       // See if any matching quark-antiquark pairs among the rest.
 853       for (int iQRec = 0; iQRec < nAnyRec; ++iQRec) {
 854         int iQ  = iAnyRec[iQRec];
 855         int idQ = event[iQ].id();
 856         if (freeAnyRec[iQRec] && idQ > 0 && idQ < 6)
 857         for (int iQbarRec = 0; iQbarRec < nAnyRec; ++iQbarRec) {
 858           int iQbar  = iAnyRec[iQbarRec];
 859           if (freeAnyRec[iQbarRec] && event[iQbar].id() == -idQ) {
 860
 861             // Check that these can be traced back to same gluon splitting.
 862             // For now also avoid qqbar pairs produced in rescatterings.??
 863             int iTopQ    = event.iTopCopyId(iQ);
 864             int iTopQbar = event.iTopCopyId(iQbar);
 865             int iMother  = event[iTopQ].mother1();
 866             if (event[iTopQbar].mother1() == iMother
 867               && event[iMother].isGluon() && event[iMother].status() != -34
 868               && event[iMother + 1].status() != -34 ) {
 869
 870               // Now find the dipole that gives the smallest
 871               // ((pQ + pQbar) * pI) ((pQ + pQbar) * pJ) / (pI * pJ).
 872               Vec4   pGlu      = event[iQ].p() + event[iQbar].p();
 873               int    iDipMin   = 0;
 874               double pT2DipMin = sCM;
 875               for (int iDip = 0; iDip < int(dipoles.size()); ++iDip) {
 876                 double pT2Dip = (pGlu * event[dipoles[iDip].iCol].p())
 877                   * (pGlu * event[dipoles[iDip].iAcol].p())
 878                   / dipoles[iDip].p1p2;
 879                 if (pT2Dip < pT2DipMin) {
 880                   iDipMin   = iDip;
 881                   pT2DipMin = pT2Dip;
 882                 }
 883               }
 884
 885               // Attach the q-qbar pair to the dipole, i.e. split the dipole.
 886               int colGlu   = event[iQ].col();
 887               int acolGlu  = event[iQbar].acol();
 888               int colDip   = dipoles[iDipMin].col;
 889               int iColDip  = dipoles[iDipMin].iCol;
 890               int iAcolDip = dipoles[iDipMin].iAcol;
 891               event[iQbar].acol( colDip );
 892               if (event[iAcolDip].acol() == colDip)
 893                    event[iAcolDip].acol( colGlu );
 894               else event[iAcolDip].col(  colGlu );
 895               dipoles[iDipMin].iAcol = iQbar;
 896               dipoles[iDipMin].p1p2 = event[iColDip].p() * event[iQbar].p();
 897               dipoles.push_back( BeamDipole( colGlu, iQ, iAcolDip ) );
 898               dipoles.back().p1p2 = event[iQ].p() * event[iAcolDip].p();
 899
 900               // Remove q-qbar pair from old system: reconnect colours.
 901               freeAnyRec[iQRec]    = false;
 902               freeAnyRec[iQbarRec] = false;
 903               for (int i = oldSize; i < event.size(); ++i)
 904               if (i != iQRec && i != iQbarRec && i != iColDip && i != iAcolDip) {
 905                  if (event[i].isFinal()) {
 906                   if (event[i].acol() == colGlu) event[i].acol( acolGlu );
 907                 } else {
 908                     if (event[i].col()  == colGlu) event[i].col( acolGlu );
 909                 }
 910               }
 911
 912             // Done with processing of q-qbar pairs.
 913             }
 914           }
 915         }
 916       }
 917
 918       // If only two beam gluons left of system, set their colour = anticolour.
 919       // Used by BeamParticle::remnantColours to skip irrelevant gluons.
 920       if ( event[iInARec].isGluon() && !event[iInARec].isRescatteredIncoming()
 921         && event[iInBRec].isGluon() && !event[iInBRec].isRescatteredIncoming()
 922         && event[iInARec].col() == event[iInBRec].acol()
 923         && event[iInARec].acol() == event[iInBRec].col() ) {
 924           event[iInARec].acol( event[iInARec].col() );
 925           event[iInBRec].acol( event[iInBRec].col() );
 926       }
 927
 928     // End of loops over iRec and iSys systems.
 929     }
 930   }
 931
 932   // Done.
 933   return true;
 934
 935 }
 936
 937 //--------------------------------------------------------------------------
 938
 939 // Collapse colours and check that they are consistent.
 940
 941 bool BeamRemnants::checkColours( Event& event) {
 942
 943   // No colours in lepton beams so no need to do anything.
 944   if (beamAPtr->isLepton() && beamBPtr->isLepton()) return true;
 945
 946   // Remove ambiguities when one colour collapses two ways.
 947   // Resolve chains where one colour is mapped to another.
 948   for (int iCol = 1; iCol < int(colFrom.size()); ++iCol)
 949   for (int iColRef = 0; iColRef < iCol; ++iColRef) {
 950     if (colFrom[iCol] == colFrom[iColRef]) {
 951       colFrom[iCol] = colTo[iCol];
 952       colTo[iCol] = colTo[iColRef];
 953     }
 954     if (colTo[iCol] == colFrom[iColRef]) colTo[iCol] = colTo[iColRef];
 955   }
 956
 957   // Transform event record colours from beam remnant colour collapses.
 958   for (int i = oldSize; i < event.size(); ++i) {
 959     int col = event[i].col();
 960     int acol = event[i].acol();
 961     for (int iCol = 0; iCol < int(colFrom.size()); ++iCol) {
 962       if (col == colFrom[iCol]) {col = colTo[iCol]; event[i].col(col);}
 963       if (acol == colFrom[iCol]) {acol = colTo[iCol]; event[i].acol(acol);}
 964     }
 965   }
 966
 967   // Transform junction colours from beam remnant colour collapses.
 968   for (int iJun = 0; iJun < event.sizeJunction(); ++iJun)
 969   for (int leg = 0; leg < 3; ++leg) {
 970     int col = event.colJunction(iJun, leg);
 971     //cout<< " BeamRemnants iJun = "<<iJun<<" leg = "<<leg<<" col = "<<col<<endl;
 972     for (int iCol = 0; iCol < int(colFrom.size()); ++iCol) {
 973       //cout << " iCol = "<<iCol<<" colFrom = "<<colFrom[iCol]<<endl;
 974       if (col == colFrom[iCol]) {
 975         col = colTo[iCol];
 976         event.colJunction(iJun, leg, col);
 977       }
 978     }
 979   }
 980
 981   // Arrays for current colours and anticolours, and for singlet gluons.
 982   vector<int> colList;
 983   vector<int> acolList;
 984   vector<int> iSingletGluon;
 985
 986   // Find current colours and anticolours in the event record.
 987   for (int i = oldSize; i < event.size(); ++i)
 988   if (event[i].isFinal()) {
 989     int id   = event[i].id();
 990     int col  = event[i].col();
 991     int acol = event[i].acol();
 992
 993     // Quarks must have colour set, antiquarks anticolour, gluons both.
 994     if ( (id > 0 && id < 9 && (col <= 0 || acol != 0) )
 995       || (id < 0 && id > -9 && (col != 0 || acol <= 0) )
 996       || (id == 21 && (col <= 0 || acol <= 0) ) ) {
 997       infoPtr->errorMsg("Error in BeamRemnants::checkColours: "
 998         "q/qbar/g has wrong colour slots set");
 999       return false;
1000     }
1001
1002     // Save colours/anticolours, and position of colour singlet gluons.
1003     if ( col > 0)  colList.push_back(  col );
1004     if (acol > 0) acolList.push_back( acol );
1005     if (col > 0 && acol == col) iSingletGluon.push_back(i);
1006   }
1007
1008   // Run though list of singlet gluons and put them on final-state dipole
1009   // (i,j) that offers smallest (p_g p_i) * (p_g p_j) / (p_i p_j).
1010   for (int iS = 0; iS < int(iSingletGluon.size()); ++iS) {
1011     int    iGlu      = iSingletGluon[iS];
1012     int    iAcolDip  = -1;
1013     double pT2DipMin = sCM;
1014     for (int iC = oldSize; iC < event.size(); ++iC)
1015       if (iC != iGlu && event[iC].isFinal()) {
1016       int colDip = event[iC].col();
1017       if (colDip > 0 && event[iC].acol() !=colDip)
1018       for (int iA = oldSize; iA < event.size(); ++iA)
1019         if (iA != iGlu && iA != iC && event[iA].isFinal()
1020         && event[iA].acol() == colDip && event[iA].col() !=colDip) {
1021         double pT2Dip = (event[iGlu].p() * event[iC].p())
1022           * (event[iGlu].p() * event[iA].p())
1023           / (event[iC].p() * event[iA].p());
1024         if (pT2Dip < pT2DipMin) {
1025           iAcolDip  = iA;
1026           pT2DipMin = pT2Dip;
1027         }
1028       }
1029     }
1030
1031     // Fail if no dipole. Else insert singlet gluon onto relevant dipole.
1032     if (iAcolDip == -1)  return false;
1033     event[iGlu].acol( event[iAcolDip].acol() );
1034     event[iAcolDip].acol( event[iGlu].col() );
1035   }
1036
1037   // Check that not the same colour or anticolour appears twice.
1038   for (int iCol = 0; iCol < int(colList.size()) - 1; ++iCol) {
1039     int col = colList[iCol];
1040     for (int iCol2 = iCol + 1; iCol2 < int(colList.size()); ++iCol2)
1041     if (colList[iCol2] == col) {
1042       infoPtr->errorMsg("Warning in BeamRemnants::checkColours:"
1043         " colour appears twice");
1044       if (!ALLOWCOLOURTWICE) return false;
1045     }
1046   }
1047   for (int iAcol = 0; iAcol < int(acolList.size()) - 1; ++iAcol) {
1048     int acol = acolList[iAcol];
1049     for (int iAcol2 = iAcol + 1; iAcol2 < int(acolList.size()); ++iAcol2)
1050     if (acolList[iAcol2] == acol) {
1051       infoPtr->errorMsg("Warning in BeamRemnants::checkColours:"
1052         " anticolour appears twice");
1053       if (!ALLOWCOLOURTWICE) return false;
1054     }
1055   }
1056
1057   // Remove all matching colour-anticolour pairs.
1058   bool foundPair = true;
1059   while (foundPair && colList.size() > 0 && acolList.size() > 0) {
1060     foundPair = false;
1061     for (int iCol = 0; iCol < int(colList.size()); ++iCol) {
1062       for (int iAcol = 0; iAcol < int(acolList.size()); ++iAcol) {
1063         if (acolList[iAcol] == colList[iCol]) {
1064           colList[iCol] = colList.back(); colList.pop_back();
1065           acolList[iAcol] = acolList.back(); acolList.pop_back();
1066           foundPair = true; break;
1067         }
1068       } if (foundPair) break;
1069     }
1070   }
1071
1072   // Check that remaining (anti)colours are accounted for by junctions.
1073   for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) {
1074     int kindJun = event.kindJunction(iJun);
1075     for (int leg = 0; leg < 3; ++leg) {
1076       int colEnd = event.colJunction(iJun, leg);
1077
1078       // Junction connected to three colours.
1079       if (kindJun == 1) {
1080         bool foundCol = false;
1081         for (int iCol = 0; iCol < int(colList.size()); ++iCol)
1082         if (colList[iCol] == colEnd) {
1083           colList[iCol] = colList.back();
1084           colList.pop_back();
1085           foundCol = true;
1086           break;
1087         }
1088       }
1089
1090       // Junction connected to three anticolours.
1091       else if (kindJun == 2) {
1092         bool foundCol = false;
1093         for (int iAcol = 0; iAcol < int(acolList.size()); ++iAcol)
1094         if (acolList[iAcol] == colEnd) {
1095           acolList[iAcol] = acolList.back();
1096           acolList.pop_back();
1097           foundCol = true;
1098           break;
1099         }
1100       }
1101
1102       // Other junction types
1103       else if ( kindJun == 3 || kindJun == 5) {
1104         bool foundCol = false;
1105         for (int iCol = 0; iCol < int(colList.size()); ++iCol)
1106         if (colList[iCol] == colEnd) {
1107           colList[iCol] = colList.back();
1108           colList.pop_back();
1109           foundCol = true;
1110           break;
1111         }
1112       }
1113
1114       // Other antijunction types
1115       else if ( kindJun == 4 || kindJun == 6) {
1116         bool foundCol = false;
1117         for (int iAcol = 0; iAcol < int(acolList.size()); ++iAcol)
1118         if (acolList[iAcol] == colEnd) {
1119           acolList[iAcol] = acolList.back();
1120           acolList.pop_back();
1121           foundCol = true;
1122           break;
1123         }
1124       }
1125
1126       // End junction check.
1127     }
1128   }
1129
1130
1131   // Repair step - sometimes needed when rescattering allowed.
1132   if (colList.size() > 0 || acolList.size() > 0)
1133     infoPtr->errorMsg("Warning in BeamRemnants::checkColours:"
1134     " need to repair unmatched colours");
1135   while (colList.size() > 0 && acolList.size() > 0) {
1136
1137     // Replace one colour and one anticolour index by a new common one.
1138     int  colMatch =  colList.back();
1139     int acolMatch = acolList.back();
1140     int  colNew   = event.nextColTag();
1141     colList.pop_back();
1142     acolList.pop_back();
1143     for (int i = oldSize; i < event.size(); ++i)
1144     if (event[i].isFinal() && event[i].col() == colMatch) {
1145       event[i].col( colNew);
1146       break;
1147     }
1148     for (int i = oldSize; i < event.size(); ++i)
1149     if (event[i].isFinal() && event[i].acol() == acolMatch) {
1150       event[i].acol( colNew);
1151       break;
1152     }
1153   }
1154
1155   // Done.
1156   return (colList.size() == 0 && acolList.size() == 0);
1157
1158 }
1159
1160 //==========================================================================
1161
1162 } // end namespace Pythia8