PYTHIA8/pythia8130/src/BeamRemnants.cxx

   1 // BeamRemnants.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2008 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 // Maximum number of tries to match colours and kinematics in the event.
  41 const int    BeamRemnants::NTRYCOLMATCH   = 10;
  42 const int    BeamRemnants::NTRYKINMATCH   = 10;
  43
  44 //*********
  45
  46 // Initialization.
  47
  48 bool BeamRemnants::init( Info* infoPtrIn, BeamParticle* beamAPtrIn,
  49   BeamParticle* beamBPtrIn) {
  50
  51   // Save pointers.
  52   infoPtr             = infoPtrIn;
  53   beamAPtr            = beamAPtrIn;
  54   beamBPtr            = beamBPtrIn;
  55
  56   // Width of primordial kT distribution.
  57   doPrimordialKT      = Settings::flag("BeamRemnants:primordialKT");
  58   primordialKTsoft    = Settings::parm("BeamRemnants:primordialKTsoft");
  59   primordialKThard    = Settings::parm("BeamRemnants:primordialKThard");
  60   primordialKTremnant = Settings::parm("BeamRemnants:primordialKTremnant");
  61   halfScaleForKT      = Settings::parm("BeamRemnants:halfScaleForKT");
  62   halfMassForKT       = Settings::parm("BeamRemnants:halfMassForKT");
  63
  64   // Parameters for colour reconnection scenario, partly borrowed from
  65   // multiple interactions not to introduce too many new ones.
  66   doReconnect         = Settings::flag("BeamRemnants:reconnectColours");
  67   reconnectRange      = Settings::parm("BeamRemnants:reconnectRange");
  68   pT0Ref              = Settings::parm("MultipleInteractions:pT0Ref");
  69   ecmRef              = Settings::parm("MultipleInteractions:ecmRef");
  70   ecmPow              = Settings::parm("MultipleInteractions:ecmPow");
  71
  72   // Total and squared CM energy at nominal energy.
  73   eCM                 = infoPtr->eCM();
  74   sCM                 = eCM * eCM;
  75
  76   // The MI pT0 smoothening scale and its reconnection-strength combination.
  77   pT0                 = pT0Ref * pow(eCM / ecmRef, ecmPow);
  78   pT20Rec             = pow2(reconnectRange * pT0);
  79
  80   // Done.
  81   return true;
  82
  83 }
  84
  85 //*********
  86
  87 // Select the flavours/kinematics/colours of the two beam remnants.
  88 // Notation: iPar = all partons, iSys = matched systems of two beams,
  89 //           iRem = additional partons in remnants.
  90
  91 bool BeamRemnants::add( Event& event) {
  92
  93   // Update to current CM energy.
  94   eCM     = infoPtr->eCM();
  95   sCM     = eCM * eCM;
  96
  97   // Number of scattering subsystems. Size of event record before treatment.
  98   nSys    = event.sizeSystems();
  99   oldSize = event.size();
 100
 101   // Add required extra remnant flavour content.
 102   // Start all over if fails (in option where junctions not allowed).
 103   if ( !beamAPtr->remnantFlavours(event)
 104     || !beamBPtr->remnantFlavours(event) ) {
 105     infoPtr->errorMsg("Error in BeamRemnants::add:"
 106       " remnant flavour setup failed");
 107     return false;
 108   }
 109
 110   // Do the kinematics of the collision subsystems and two beam remnants
 111   if (!setKinematics(event)) return false;
 112
 113   // Allow colour reconnections.
 114   if (doReconnect) reconnectColours(event);
 115
 116   // Save current modifiable colour configuration for fast restoration.
 117   vector<int> colSave;
 118   vector<int> acolSave;
 119   for (int i = oldSize; i < event.size(); ++i) {
 120     colSave.push_back( event[i].col() );
 121     acolSave.push_back( event[i].acol() );
 122   }
 123   event.saveJunctionSize();
 124
 125   // Allow several tries to match colours of initiators and remnants.
 126   // Frequent "failures" since shortcutting colours separately on
 127   // the two event sides may give "colour singlet gluons" etc.
 128   bool physical = true;
 129   for (int iTry = 0; iTry < NTRYCOLMATCH; ++iTry) {
 130     physical = true;
 131
 132     // Reset list of colour "collapses" (transformations).
 133     colFrom.resize(0);
 134     colTo.resize(0);
 135
 136     // First process each set of beam colours on its own.
 137     if (!beamAPtr->remnantColours(event, colFrom, colTo))
 138       physical = false;
 139     if (!beamBPtr->remnantColours(event, colFrom, colTo))
 140       physical = false;
 141
 142     // Then check that colours and anticolours are matched in whole event.
 143     if ( physical && !checkColours(event) ) physical = false;
 144
 145     // If no problems then done, else restore and loop.
 146     if (physical) break;
 147     for (int i = oldSize; i < event.size(); ++i)
 148       event[i].cols( colSave[i - oldSize], acolSave[i - oldSize] );
 149     event.restoreJunctionSize();
 150   }
 151
 152   // If no solution after several tries then failed.
 153   if (!physical) {
 154     infoPtr->errorMsg("Error in BeamRemnants::add:"
 155       " colour tracing failed");
 156     return false;
 157   }
 158
 159   // Done.
 160   return true;
 161
 162 }
 163
 164 //*********
 165
 166 // Set up trial transverse and longitudinal kinematics for each beam
 167 // separately. Final decisions involve comparing the two beams.
 168
 169 bool BeamRemnants::setKinematics( Event& event) {
 170
 171   // References to beams to simplify indexing.
 172   BeamParticle& beamA = *beamAPtr;
 173   BeamParticle& beamB = *beamBPtr;
 174
 175   // Nothing to do for lepton-lepton scattering with all energy already used.
 176   if ( beamA.isUnresolvedLepton() && beamB.isUnresolvedLepton() )
 177     return true;
 178
 179   // Allow primordial kT reduction for small-mass and small-pT systems
 180   // (for hardest interaction pT -> renormalization scale so also 2 -> 1).
 181   vector<double> kTwidth;
 182   vector<double> kTcomp;
 183   double kTcompSumSys = 0.;
 184   for (int iSys = 0; iSys < nSys; ++iSys) {
 185     double kTwidthNow = 0.;
 186     double mHatDamp   = 1.;
 187     if (doPrimordialKT) {
 188       double mHat     = sqrtpos( beamA[iSys].x() * beamB[iSys].x() * sCM );
 189       mHatDamp        = mHat / (mHat + halfMassForKT);
 190       double scale    = (iSys == 0) ? infoPtr->QRen() : infoPtr->pTMI(iSys);
 191       kTwidthNow      = ( (halfScaleForKT * primordialKTsoft
 192       + scale * primordialKThard) / (halfScaleForKT + scale) ) * mHatDamp;
 193     }
 194     kTwidth.push_back( kTwidthNow );
 195     kTcomp.push_back( mHatDamp );
 196     kTcompSumSys += mHatDamp;
 197   }
 198   double kTwidthNow = (doPrimordialKT) ? primordialKTremnant : 0.;
 199   for (int iRem = nSys; iRem < max( beamA.size(), beamB.size() ); ++iRem) {
 200     kTwidth.push_back( kTwidthNow );
 201     kTcomp.push_back( 1. );
 202   }
 203
 204   // Allow ten tries to construct kinematics (but normally works first).
 205   bool physical;
 206   double xSum[2], xInvM[2], w2Beam[2], wPosRem, wNegRem, w2Rem;
 207   for (int iTry = 0; iTry < NTRYKINMATCH; ++iTry) {
 208     physical = true;
 209
 210     // Loop over the two beams. Sum px and py separately within each.
 211     for (int iBeam = 0; iBeam < 2; ++iBeam) {
 212       BeamParticle& beam = (iBeam == 0) ? beamA : beamB;
 213       int nPar = beam.size();
 214       double pxSum = 0.;
 215       double pySum = 0.;
 216
 217       // Loop over the partons in a beam. Generate Gaussian pT for hadrons.
 218       for (int iPar = 0; iPar < nPar; ++iPar) {
 219         double px = 0.;
 220         double py = 0.;
 221         if (beam.isHadron() && doPrimordialKT) {
 222           px = kTwidth[iPar] * Rndm::gauss();
 223           py = kTwidth[iPar] * Rndm::gauss();
 224         }
 225         beam[iPar].px(px);
 226         beam[iPar].py(py);
 227         pxSum += px;
 228         pySum += py;
 229       }
 230
 231       // Share recoil between all partons.
 232       if (doPrimordialKT) {
 233         double kTcompSum = kTcompSumSys + (nPar - nSys);
 234         for (int iPar = 0; iPar < nPar; ++iPar) {
 235           beam[iPar].px( beam[iPar].px() - pxSum * kTcomp[iPar] / kTcompSum );
 236           beam[iPar].py( beam[iPar].py() - pySum * kTcomp[iPar] / kTcompSum );
 237         }
 238       }
 239
 240       // Pick unrescaled x values for remnants. Sum up (unscaled) p+ and p-.
 241       xSum[iBeam] = 0.;
 242       xInvM[iBeam] = 0.;
 243       for (int iRem = nSys; iRem < nPar; ++iRem) {
 244         double xPrel = beam.xRemnant( iRem);
 245         beam[iRem].x(xPrel);
 246         xSum[iBeam] += xPrel;
 247         xInvM[iBeam] += beam[iRem].mT2()/xPrel;
 248       }
 249
 250       // Squared transverse mass for each beam, using lightcone x.
 251       w2Beam[iBeam] = xSum[iBeam] * xInvM[iBeam];
 252
 253     // End separate treatment of the two beams.
 254     }
 255
 256     // Recalculate kinematics of initiator systems with primordial kT.
 257     wPosRem = eCM;
 258     wNegRem = eCM;
 259     for (int iSys = 0; iSys < nSys; ++iSys) {
 260       double sHat = beamA[iSys].x() * beamB[iSys].x() * sCM;
 261       double sHatT = sHat
 262         + pow2( beamA[iSys].px() + beamB[iSys].px())
 263         + pow2( beamA[iSys].py() + beamB[iSys].py());
 264       double rescale = sqrt( sHatT / sHat);
 265       wPosRem -= rescale * beamA[iSys].x() * eCM;
 266       wNegRem -= rescale * beamB[iSys].x() * eCM;
 267
 268       // Kinematics forbids too large primordial kT.
 269       if (sqrt(sHatT) < beamA[iSys].pT() + beamB[iSys].pT())
 270         physical = false;
 271     }
 272
 273     // Check that remaining momentum is enough for remnants.
 274     if (wPosRem < 0. || wNegRem < 0.) physical = false;
 275     w2Rem = wPosRem * wNegRem;
 276     if (sqrtpos(w2Rem) < sqrt(w2Beam[0]) + sqrt(w2Beam[1]))
 277       physical = false;
 278
 279     // End of loop over ten tries. Do not loop when solution found.
 280     if (physical) break;
 281   }
 282
 283   // If no solution after ten tries then failed.
 284   if (!physical) {
 285     infoPtr->errorMsg("Error in BeamRemnants::add:"
 286       " kinematics construction failed");
 287     return false;
 288   }
 289
 290   // Construct energy and pz for each initiator pair.
 291   for (int iSys = 0; iSys < nSys; ++iSys) {
 292     double sHat = beamA[iSys].x() * beamB[iSys].x() * sCM;
 293     double sHatT = sHat + pow2( beamA[iSys].px() + beamB[iSys].px())
 294       + pow2( beamA[iSys].py() + beamB[iSys].py());
 295     double rescale = sqrt( sHatT / sHat);
 296     double wPos = rescale * beamA[iSys].x() * eCM;
 297     double wNeg = rescale * beamB[iSys].x() * eCM;
 298     double w2A = beamA[iSys].mT2();
 299     double w2B = beamB[iSys].mT2();
 300     double lambdaRoot = sqrtpos( pow2( sHatT - w2A - w2B) - 4. * w2A * w2B );
 301     double pPosA = 0.5 * (sHatT + w2A - w2B + lambdaRoot) / sHatT * wPos;
 302     beamA[iSys].e( 0.5 * (pPosA + w2A / pPosA) );
 303     beamA[iSys].pz( 0.5 * (pPosA - w2A / pPosA) );
 304     double pNegB = 0.5 * (sHatT + w2B - w2A + lambdaRoot) / sHatT * wNeg;
 305     beamB[iSys].e( 0.5 * (pNegB + w2B / pNegB) );
 306     beamB[iSys].pz( 0.5 * (w2B / pNegB - pNegB) );
 307
 308     // Construct rotations and boosts caused by primordial kT.
 309     int iA = beamA[iSys].iPos();
 310     int iB = beamB[iSys].iPos();
 311     RotBstMatrix M;
 312     M.toCMframe( event[iA].p(), event[iB].p() );
 313     M.fromCMframe( beamA[iSys].p(), beamB[iSys].p() );
 314
 315     // Copy initiators and their systems and boost them accordingly.
 316     // Update subsystem and beams info on new positions of partons.
 317     int iAcopy = event.copy(iA, -61);
 318     event[iAcopy].rotbst(M);
 319     event.setInSystem(iSys, 0, iAcopy);
 320     beamA[iSys].iPos( iAcopy);
 321     int iBcopy = event.copy(iB, -61);
 322     event[iBcopy].rotbst(M);
 323     event.setInSystem(iSys, 1, iBcopy);
 324     beamB[iSys].iPos( iBcopy);
 325     for (int iABsys = 2; iABsys < event.sizeSystem(iSys); ++iABsys) {
 326       int iAB = event.getInSystem(iSys, iABsys);
 327       int iABcopy = event.copy(iAB, 62);
 328       event[iABcopy].rotbst(M);
 329       event.setInSystem(iSys, iABsys, iABcopy);
 330     }
 331
 332     // Update daughter info of mothers, i.e. of beams, for hardest interaction.
 333     if (iSys == 0) {
 334       int mother = event[iAcopy].mother1();
 335       event[mother].daughter1(iAcopy);
 336       mother = event[iBcopy].mother1();
 337       event[mother].daughter1(iBcopy);
 338     }
 339   }
 340
 341   // Construct x rescaling factors for the two remants.
 342   double lambdaRoot = sqrtpos( pow2(w2Rem - w2Beam[0] - w2Beam[1])
 343     - 4. * w2Beam[0] * w2Beam[1] );
 344   double rescaleA = (w2Rem + w2Beam[0] - w2Beam[1] + lambdaRoot)
 345     / (2. * w2Rem * xSum[0]) ;
 346   double rescaleB = (w2Rem + w2Beam[1] - w2Beam[0] + lambdaRoot)
 347     / (2. * w2Rem * xSum[1]) ;
 348
 349   // Construct energy and pz for remnants in first beam.
 350   for (int iRem = nSys; iRem < beamA.size(); ++iRem) {
 351     double pPos = rescaleA * beamA[iRem].x() * wPosRem;
 352     double pNeg = beamA[iRem].mT2() / pPos;
 353     beamA[iRem].e( 0.5 * (pPos + pNeg) );
 354     beamA[iRem].pz( 0.5 * (pPos - pNeg) );
 355
 356     // Add these partons to the normal event record.
 357     int iNew = event.append( beamA[iRem].id(), 63, 1, 0, 0, 0,
 358       beamA[iRem].col(), beamA[iRem].acol(), beamA[iRem].p(),
 359       beamA[iRem].m() );
 360     beamA[iRem].iPos( iNew);
 361   }
 362
 363   // Construct energy and pz for remnants in second beam.
 364   for (int iRem = nSys; iRem < beamB.size(); ++iRem) {
 365     double pNeg = rescaleB * beamB[iRem].x() * wNegRem;
 366     double pPos = beamB[iRem].mT2() / pNeg;
 367     beamB[iRem].e( 0.5 * (pPos + pNeg) );
 368     beamB[iRem].pz( 0.5 * (pPos - pNeg) );
 369
 370     // Add these partons to the normal event record.
 371     int iNew = event.append( beamB[iRem].id(), 63, 2, 0, 0, 0,
 372       beamB[iRem].col(), beamB[iRem].acol(), beamB[iRem].p(),
 373       beamB[iRem].m() );
 374     beamB[iRem].iPos( iNew);
 375   }
 376
 377   // Done.
 378   return true;
 379
 380 }
 381
 382 //*********
 383
 384 // Allow colour reconnections by mergings of collision subsystems.
 385 // iRec is system that may be reconnected, by moving its gluons to iSys,
 386 // where minimal pT (or equivalently Lambda) is used to pick location.
 387 // Therefore all dipoles in iSys have to be found, and all gluons in iRec.
 388 // Matching q-qbar pairs are treated by analogy with gluons.
 389
 390 bool BeamRemnants::reconnectColours( Event&  event) {
 391
 392   // Prepare record of which systems should be merged onto another.
 393   // The iSys system must have colour in final state to attach to it.
 394   vector<int>  iMerge;
 395   vector<bool> hasColour;
 396   for (int iSys = 0; iSys < nSys; ++iSys) {
 397     iMerge.push_back( iSys );
 398     bool hasCol = false;
 399     for (int iMem = 2; iMem < event.sizeSystem( iSys); ++iMem) {
 400       int iNow = event.getInSystem( iSys, iMem);
 401       if (event[iNow].col() > 0 || event[iNow].acol() > 0) {
 402         hasCol = true;
 403         break;
 404       }
 405     }
 406     hasColour.push_back( hasCol );
 407   }
 408
 409   // Loop over system that may be reconnected.
 410   for (int iRec = nSys - 1; iRec > 0; --iRec) {
 411
 412     // Determine reconnection strength from pT scale of system.
 413     double pT2Rec  = pow2( infoPtr->pTMI(iRec) );
 414     double probRec = pT20Rec / (pT20Rec + pT2Rec);
 415
 416     // Loop over other systems iSys at higher pT scale and
 417     // decide whether to reconnect the iRec gluons onto one of them.
 418     for (int iSys = iRec - 1; iSys >= 0; --iSys)
 419     if (hasColour[iSys] && probRec > Rndm::flat()) {
 420
 421       // The iRec system and all merged with it to be merged with iSys.
 422       iMerge[iRec] = iSys;
 423       for (int iRec2 = iRec + 1; iRec2 < nSys; ++iRec2)
 424       if (iMerge[iRec2] == iRec) iMerge[iRec2] = iSys;
 425
 426       // Once a system has been merged do not test it anymore.
 427       break;
 428     }
 429   }
 430
 431   // Loop over systems. Check whether other systems to be merged with it.
 432   for (int iSys = 0; iSys < nSys; ++iSys) {
 433
 434     int nMerge = 0;
 435     for (int iRec = iSys + 1; iRec < nSys; ++iRec)
 436     if (iMerge[iRec] == iSys) ++nMerge;
 437     if (nMerge == 0) continue;
 438
 439     // Begin find dipoles in iSys system.
 440     vector<BeamDipole> dipoles;
 441     int sizeSys = event.sizeSystem( iSys);
 442     int iInASys = event.getInSystem( iSys, 0);
 443     int iInBSys = event.getInSystem( iSys, 1);
 444     for (int iMem = 2; iMem < sizeSys; ++iMem) {
 445
 446       // Find colour dipoles to beam remnant.
 447       int iNow = event.getInSystem( iSys, iMem);
 448       int col = event[iNow].col();
 449       if (col > 0) {
 450         if (event[iInASys].col() == col)
 451           dipoles.push_back( BeamDipole( col, iNow, iInASys ) );
 452         else if (event[iInBSys].col() == col)
 453           dipoles.push_back( BeamDipole( col, iNow, iInBSys ) );
 454
 455         // Find colour dipole between final-state partons.
 456         else for (int iMem2 = 2; iMem2 < sizeSys; ++iMem2)
 457         if (iMem2 != iMem) {
 458           int iNow2 = event.getInSystem( iSys, iMem2);
 459           if (event[iNow2].acol() == col) {
 460             dipoles.push_back( BeamDipole( col, iNow, iNow2) );
 461             break;
 462           }
 463         }
 464       }
 465
 466       // Find anticolour dipoles to beam remnant.
 467       int acol = event[iNow].acol();
 468       if (acol > 0) {
 469         if (event[iInASys].acol() == acol)
 470           dipoles.push_back( BeamDipole( acol, iInASys, iNow ) );
 471         else if (event[iInBSys].acol() == acol)
 472           dipoles.push_back( BeamDipole( acol, iInBSys, iNow ) );
 473       }
 474     }
 475
 476     // Find dipole sizes.
 477     for (int iDip = 0; iDip < int(dipoles.size()); ++iDip)
 478       dipoles[iDip].p1p2 = event[dipoles[iDip].iCol].p()
 479                          * event[dipoles[iDip].iAcol].p();
 480
 481     // Loop over systems iRec to be merged with iSys.
 482     for (int iRec = iSys + 1; iRec < nSys; ++iRec) {
 483       if (iMerge[iRec] != iSys) continue;
 484
 485       // Information on iRec. Vectors for gluons and anything else.
 486       int sizeRec = event.sizeSystem( iRec);
 487       int iInARec = event.getInSystem( iRec, 0);
 488       int iInBRec = event.getInSystem( iRec, 1);
 489       int nGluRec = 0;
 490       vector<int>    iGluRec;
 491       vector<double> pT2GluRec;
 492       int nAnyRec = 0;
 493       vector<int>    iAnyRec;
 494       vector<bool>   freeAnyRec;
 495
 496       // Copy of gluon positions in descending order.
 497       for (int iMem = 2; iMem < sizeRec; ++iMem) {
 498         int iNow = event.getInSystem( iRec, iMem);
 499         if (event[iNow].isGluon()) {
 500           ++nGluRec;
 501           iGluRec.push_back( iNow );
 502           pT2GluRec.push_back( event[iNow].pT2() );
 503           for (int i = nGluRec - 1; i > 1; --i) {
 504             if (pT2GluRec[i - 1] > pT2GluRec[i]) break;
 505             swap(   iGluRec[i - 1],   iGluRec[i] );
 506             swap( pT2GluRec[i - 1], pT2GluRec[i] );
 507           }
 508         // Copy of anything else, mainly quarks, in no particular order.
 509         } else {
 510           ++nAnyRec;
 511           iAnyRec.push_back( iNow );
 512           freeAnyRec.push_back( true );
 513         }
 514       }
 515
 516       // For each gluon in iRec now find the dipole that gives the smallest
 517       // (pGlu * pI) (pGlu * pJ) / (pI * pJ), i.e. minimal pT (and Lambda).
 518       for (int iGRec = 0; iGRec < nGluRec; ++iGRec) {
 519         int    iGlu      = iGluRec[iGRec];
 520         Vec4   pGlu      = event[iGlu].p();
 521         int    iDipMin   = 0;
 522         double pT2DipMin = sCM;
 523         for (int iDip = 0; iDip < int(dipoles.size()); ++iDip) {
 524           double pT2Dip = (pGlu * event[dipoles[iDip].iCol].p())
 525             * (pGlu * event[dipoles[iDip].iAcol].p()) / dipoles[iDip].p1p2;
 526           if (pT2Dip < pT2DipMin) {
 527             iDipMin   = iDip;
 528             pT2DipMin = pT2Dip;
 529           }
 530         }
 531
 532         // Attach the gluon to the dipole, i.e. split the dipole in two.
 533         int colGlu   = event[iGlu].col();
 534         int acolGlu  = event[iGlu].acol();
 535         int colDip   = dipoles[iDipMin].col;
 536         int iColDip  = dipoles[iDipMin].iCol;
 537         int iAcolDip = dipoles[iDipMin].iAcol;
 538         event[iGlu].acol( colDip );
 539         if (event[iAcolDip].acol() == colDip)
 540              event[iAcolDip].acol( colGlu );
 541         else event[iAcolDip].col(  colGlu );
 542         dipoles[iDipMin].iAcol = iGlu;
 543         dipoles[iDipMin].p1p2 = event[iColDip].p() * pGlu;
 544         dipoles.push_back( BeamDipole( colGlu, iGlu, iAcolDip ) );
 545         dipoles.back().p1p2 = pGlu * event[iAcolDip].p();
 546
 547         // Remove gluon from old system: reconnect colours.
 548         if      (event[iInARec].col() == colGlu)
 549           event[iInARec].col( acolGlu );
 550         else if (event[iInBRec].col() == colGlu)
 551           event[iInBRec].col( acolGlu );
 552         else for (int i = iGRec + 1; i < nGluRec; ++i)
 553         if (event[iGluRec[i]].acol() == colGlu) {
 554           event[iGluRec[i]].acol( acolGlu );
 555           break;
 556         }
 557         for (int i = 0; i < nAnyRec; ++i)
 558         if (event[iAnyRec[i]].acol() == colGlu) {
 559           event[iAnyRec[i]].acol( acolGlu );
 560           break;
 561         }
 562       }
 563
 564       // See if any matching quark-antiquark pairs among the rest.
 565       for (int iQRec = 0; iQRec < nAnyRec; ++iQRec) {
 566         int iQ  = iAnyRec[iQRec];
 567         int idQ = event[iQ].id();
 568         if (freeAnyRec[iQRec] && idQ > 0 && idQ < 6)
 569         for (int iQbarRec = 0; iQbarRec < nAnyRec; ++iQbarRec) {
 570           int iQbar  = iAnyRec[iQbarRec];
 571           if (freeAnyRec[iQbarRec] && event[iQbar].id() == -idQ) {
 572
 573             // Check that these can be traced back to same gluon splitting.
 574             int iTopQ    = event.iTopCopyId(iQ);
 575             int iTopQbar = event.iTopCopyId(iQbar);
 576             int iMother  = event[iTopQ].mother1();
 577             if (event[iTopQbar].mother1() == iMother
 578               && event[iMother].isGluon()) {
 579
 580               // Now find the dipole that gives the smallest
 581               // ((pQ + pQbar) * pI) ((pQ + pQbar) * pJ) / (pI * pJ).
 582               Vec4   pGlu      = event[iQ].p() + event[iQbar].p();
 583               int    iDipMin   = 0;
 584               double pT2DipMin = sCM;
 585               for (int iDip = 0; iDip < int(dipoles.size()); ++iDip) {
 586                 double pT2Dip = (pGlu * event[dipoles[iDip].iCol].p())
 587                   * (pGlu * event[dipoles[iDip].iAcol].p())
 588                   / dipoles[iDip].p1p2;
 589                 if (pT2Dip < pT2DipMin) {
 590                   iDipMin   = iDip;
 591                   pT2DipMin = pT2Dip;
 592                 }
 593               }
 594
 595               // Attach the q-qbar pair to the dipole, i.e. split the dipole.
 596               int colGlu   = event[iQ].col();
 597               int acolGlu  = event[iQbar].acol();
 598               int colDip   = dipoles[iDipMin].col;
 599               int iColDip  = dipoles[iDipMin].iCol;
 600               int iAcolDip = dipoles[iDipMin].iAcol;
 601               event[iQbar].acol( colDip );
 602               if (event[iAcolDip].acol() == colDip)
 603                    event[iAcolDip].acol( colGlu );
 604               else event[iAcolDip].col(  colGlu );
 605               dipoles[iDipMin].iAcol = iQbar;
 606               dipoles[iDipMin].p1p2 = event[iColDip].p() * event[iQbar].p();
 607               dipoles.push_back( BeamDipole( colGlu, iQ, iAcolDip ) );
 608               dipoles.back().p1p2 = event[iQ].p() * event[iAcolDip].p();
 609
 610               // Remove q-qbar pair from old system: reconnect colours.
 611               freeAnyRec[iQRec]    = false;
 612               freeAnyRec[iQbarRec] = false;
 613               if      (event[iInARec].col() == colGlu)
 614                 event[iInARec].col( acolGlu );
 615               else if (event[iInBRec].col() == colGlu)
 616                 event[iInBRec].col( acolGlu );
 617               else for (int i = 0; i < nAnyRec; ++i)
 618               if (freeAnyRec[i] && event[iAnyRec[i]].acol() == colGlu) {
 619                 event[iAnyRec[i]].acol( acolGlu );
 620                 break;
 621               }
 622
 623             // Done with processing of q-qbar pairs.
 624             }
 625           }
 626         }
 627       }
 628
 629       // If only two beam gluons left of system, set their colour = anticolour.
 630       // Used by BeamParticle::remnantColours to skip irrelevant gluons.
 631       if ( event[iInARec].isGluon() && event[iInBRec].isGluon()
 632         && event[iInARec].col() == event[iInBRec].acol()
 633         && event[iInARec].acol() == event[iInBRec].col() ) {
 634           event[iInARec].acol( event[iInARec].col() );
 635           event[iInBRec].acol( event[iInBRec].col() );
 636       }
 637
 638     // End of loops over iRec and iSys systems.
 639     }
 640   }
 641
 642   // Done.
 643   return true;
 644
 645 }
 646
 647 //*********
 648
 649 // Collapse colours and check that they are consistent.
 650
 651 bool BeamRemnants::checkColours( Event& event) {
 652
 653   // No colours in lepton beams so no need to do anything.
 654   if (beamAPtr->isLepton() && beamBPtr->isLepton()) return true;
 655
 656   // Remove ambiguities when one colour collapses two ways.
 657   // Resolve chains where one colour is mapped to another.
 658   for (int iCol = 1; iCol < int(colFrom.size()); ++iCol)
 659   for (int iColRef = 0; iColRef < iCol; ++iColRef) {
 660     if (colFrom[iCol] == colFrom[iColRef]) {
 661       colFrom[iCol] = colTo[iCol];
 662       colTo[iCol] = colTo[iColRef];
 663     }
 664     if (colTo[iCol] == colFrom[iColRef]) colTo[iCol] = colTo[iColRef];
 665   }
 666
 667   // Transform event record colours from beam remnant colour collapses.
 668   for (int i = oldSize; i < event.size(); ++i) {
 669     int col = event[i].col();
 670     int acol = event[i].acol();
 671     for (int iCol = 0; iCol < int(colFrom.size()); ++iCol) {
 672       if (col == colFrom[iCol]) {col = colTo[iCol]; event[i].col(col);}
 673       if (acol == colFrom[iCol]) {acol = colTo[iCol]; event[i].acol(acol);}
 674     }
 675   }
 676
 677   // Transform junction colours from beam remnant colour collapses.
 678   for (int iJun = 0; iJun < event.sizeJunction(); ++iJun)
 679   for (int leg = 0; leg < 3; ++leg) {
 680     int col = event.colJunction(iJun, leg);
 681     for (int iCol = 0; iCol < int(colFrom.size()); ++iCol)
 682     if (col == colFrom[iCol]) {
 683       col = colTo[iCol];
 684       event.colJunction(iJun, leg, col);
 685     }
 686   }
 687
 688   // Arrays for current colours and anticolours, and for singlet gluons.
 689   vector<int> colList;
 690   vector<int> acolList;
 691   vector<int> iSingletGluon;
 692
 693   // Find current colours and anticolours in the event record.
 694   for (int i = oldSize; i < event.size(); ++i)
 695   if (event[i].status() > 0) {
 696     int id   = event[i].id();
 697     int col  = event[i].col();
 698     int acol = event[i].acol();
 699
 700     // Quarks must have colour set, antiquarks anticolour, gluons both.
 701     if ( (id > 0 && id < 9 && (col <= 0 || acol != 0) )
 702       || (id < 0 && id > -9 && (col != 0 || acol <= 0) )
 703       || (id == 21 && (col <= 0 || acol <= 0) ) ) {
 704       infoPtr->errorMsg("Error in BeamRemnants::checkColours: "
 705         "q/qbar/g has wrong colour slots set");
 706       return false;
 707     }
 708
 709     // Save colours/anticolours, and position of colour singlet gluons.
 710     if ( col > 0)  colList.push_back(  col );
 711     if (acol > 0) acolList.push_back( acol );
 712     if (col > 0 && acol == col) iSingletGluon.push_back(i);
 713   }
 714
 715   // Run though list of singlet gluons and put them on final-state dipole
 716   // (i,j) that offers smallest (p_g p_i) * (p_g p_j) / (p_i p_j).
 717   for (int iS = 0; iS < int(iSingletGluon.size()); ++iS) {
 718     int    iGlu      = iSingletGluon[iS];
 719     int    iAcolDip  = 0;
 720     double pT2DipMin = sCM;
 721     for (int iC = oldSize; iC < event.size(); ++iC)
 722     if (iC != iGlu && event[iC].status() > 0) {
 723       int colDip = event[iC].col();
 724       if (colDip > 0 && event[iC].acol() !=colDip)
 725       for (int iA = oldSize; iA < event.size(); ++iA)
 726       if (iA != iGlu && iA != iC && event[iA].status() > 0
 727         && event[iA].acol() == colDip && event[iA].col() !=colDip) {
 728         double pT2Dip = (event[iGlu].p() * event[iC].p())
 729           * (event[iGlu].p() * event[iA].p())
 730           / (event[iC].p() * event[iA].p());
 731         if (pT2Dip < pT2DipMin) {
 732           iAcolDip  = iA;
 733           pT2DipMin = pT2Dip;
 734         }
 735       }
 736     }
 737     event[iGlu].acol( event[iAcolDip].acol() );
 738     event[iAcolDip].acol( event[iGlu].col() );
 739   }
 740
 741   // Check that not the same colour or anticolour appears twice.
 742   for (int iCol = 0; iCol < int(colList.size()) - 1; ++iCol) {
 743     int col = colList[iCol];
 744     for (int iCol2 = iCol + 1; iCol2 < int(colList.size()); ++iCol2)
 745       if (colList[iCol2] == col) return false;
 746   }
 747   for (int iAcol = 0; iAcol < int(acolList.size()) - 1; ++iAcol) {
 748     int acol = acolList[iAcol];
 749     for (int iAcol2 = iAcol + 1; iAcol2 < int(acolList.size()); ++iAcol2)
 750       if (acolList[iAcol2] == acol) return false;
 751   }
 752
 753   // Remove all matching colour-anticolour pairs.
 754   bool foundPair = true;
 755   while (foundPair && colList.size() > 0 && acolList.size() > 0) {
 756     foundPair = false;
 757     for (int iCol = 0; iCol < int(colList.size()); ++iCol) {
 758       for (int iAcol = 0; iAcol < int(acolList.size()); ++iAcol) {
 759         if (acolList[iAcol] == colList[iCol]) {
 760           colList[iCol] = colList.back(); colList.pop_back();
 761           acolList[iAcol] = acolList.back(); acolList.pop_back();
 762           foundPair = true; break;
 763         }
 764       } if (foundPair) break;
 765     }
 766   }
 767
 768   // Check that remaining (anti)colours are accounted for by junctions.
 769   for (int iJun = 0; iJun < event.sizeJunction(); ++iJun) {
 770     int kindJun = event.kindJunction(iJun);
 771     for (int leg = 0; leg < 3; ++leg) {
 772       int colEnd = event.colJunction(iJun, leg);
 773
 774       // Junction connected to three colours.
 775       if (kindJun == 1) {
 776         bool foundCol = false;
 777         for (int iCol = 0; iCol < int(colList.size()); ++iCol)
 778         if (colList[iCol] == colEnd) {
 779           colList[iCol] = colList.back();
 780           colList.pop_back();
 781           foundCol = true;
 782           break;
 783         }
 784       }
 785
 786       // Junction connected to three anticolours.
 787       else if (kindJun == 2) {
 788         bool foundCol = false;
 789         for (int iAcol = 0; iAcol < int(acolList.size()); ++iAcol)
 790         if (acolList[iAcol] == colEnd) {
 791           acolList[iAcol] = acolList.back();
 792           acolList.pop_back();
 793           foundCol = true;
 794           break;
 795         }
 796       }
 797
 798     // End junction check. More junction cases to come??
 799     }
 800   }
 801
 802   // Done.
 803   return (colList.size() == 0 && acolList.size() == 0);
 804
 805 }
 806
 807 //**************************************************************************
 808
 809 } // end namespace Pythia8