PYTHIA8/pythia8130/src/BeamParticle.cxx

   1 // BeamParticle.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 // BeamParticle class.
   8
   9 #include "BeamParticle.h"
  10
  11 namespace Pythia8 {
  12
  13 //**************************************************************************
  14
  15 // The BeamParticle class.
  16
  17 //*********
  18
  19 // Constants: could be changed here if desired, but normally should not.
  20 // These are of technical nature, as described for each.
  21
  22 // A lepton that takes (almost) the full beam energy does not leave a remnant.
  23 const double BeamParticle::XMINUNRESOLVED = 1. - 1e-10;
  24
  25 //*********
  26
  27 // Initialize data on a beam particle and save pointers.
  28
  29 void BeamParticle::init( int idIn, double pzIn, double eIn, double mIn,
  30   Info* infoPtrIn, PDF* pdfInPtr, PDF* pdfHardInPtr,
  31   bool isUnresolvedIn, StringFlav* flavSelPtrIn) {
  32
  33   // Store input pointers (and one bool) for future use.
  34   infoPtr           = infoPtrIn;
  35   pdfBeamPtr        = pdfInPtr;
  36   pdfHardBeamPtr    = pdfHardInPtr;
  37   isUnresolvedBeam  = isUnresolvedIn;
  38   flavSelPtr        = flavSelPtrIn;
  39
  40   // Maximum quark kind in allowed incoming beam hadrons.
  41   maxValQuark       = Settings::mode("BeamRemnants:maxValQuark");
  42
  43   // Power of (1-x)^power/sqrt(x) for remnant valence quark distribution.
  44   valencePowerMeson = Settings::parm("BeamRemnants:valencePowerMeson");
  45   valencePowerUinP  = Settings::parm("BeamRemnants:valencePowerUinP");
  46   valencePowerDinP  = Settings::parm("BeamRemnants:valencePowerDinP");
  47
  48   // Enhancement factor of x of diquark.
  49   valenceDiqEnhance = Settings::parm("BeamRemnants:valenceDiqEnhance");
  50
  51   // Assume g(x) ~ (1-x)^power/x to constrain companion to sea quark.
  52   companionPower    = Settings::mode("BeamRemnants:companionPower");
  53
  54   // Assume g(x) ~ (1-x)^power/x to constrain companion to sea quark.
  55   companionPower    = Settings::mode("BeamRemnants:companionPower");
  56
  57   // Allow or not more than one valence quark to be kicked out.
  58   allowJunction     = Settings::flag("BeamRemnants:allowJunction");
  59
  60   // For diffractive system kick out q/g = norm / mass^power.
  61   pickQuarkNorm     = Settings::parm("BeamRemnants:pickQuarkNorm");
  62   pickQuarkPower    = Settings::parm("BeamRemnants:pickQuarkPower");
  63
  64   // Width of primordial kT distribution in diffractive systems.
  65   diffPrimKTwidth   = Settings::parm("BeamRemnants:diffPrimKTwidth");
  66
  67   // Suppress large masses of beam remnant in diffractive systems.
  68   diffLargeMassSuppress = Settings::parm("BeamRemnants:diffLargeMassSuppress");
  69
  70   // Store info on the incoming beam.
  71   idBeam            = idIn;
  72   initBeamKind();
  73   pBeam             = Vec4( 0., 0., pzIn, eIn);
  74   mBeam             = mIn;
  75
  76 }
  77
  78 //*********
  79
  80 // Initialize kind and valence flavour content of incoming beam.
  81 // For recognized hadrons one can generate multiple interactions.
  82 // So far we do not handle diagonal mesons or K0S/K0L (or photons),
  83 // for which flavour content is only known after first interaction.
  84
  85 void BeamParticle::initBeamKind() {
  86
  87   // Reset.
  88   idBeamAbs    = abs(idBeam);
  89   isLeptonBeam = false;
  90   isHadronBeam = false;
  91   isMesonBeam  = false;
  92   isBaryonBeam = false;
  93   nValKinds    = 0;
  94
  95   // Check for leptons.
  96   if (idBeamAbs > 10 && idBeamAbs < 17) {
  97     nValKinds = 1;
  98     nVal[0]   = 1;
  99     idVal[0]  = idBeam;
 100     isLeptonBeam = true;
 101   }
 102
 103   //  Done if cannot be lowest-lying hadron state.
 104   if (idBeamAbs < 101 || idBeamAbs > 9999) return;
 105
 106   // Resolve valence content for assumed meson. Flunk unallowed codes.
 107   if (idBeamAbs < 1000) {
 108     int id1 = idBeamAbs/100;
 109     int id2 = (idBeamAbs/10)%10;
 110     if ( id1 < 1 || id1 > maxValQuark
 111       || id2 < 1 || id2 > maxValQuark ) return;
 112     if (id2 == id1 || idBeamAbs == 130 || idBeamAbs == 310) return;
 113     isMesonBeam = true;
 114
 115     // Store valence content of a confirmed meson.
 116     nValKinds = 2;
 117     nVal[0]   = 1 ;
 118     nVal[1]   = 1;
 119     if (id1%2 == 0) {
 120       idVal[0] = id1;
 121       idVal[1] = -id2;
 122     } else {
 123       idVal[0] = id2;
 124       idVal[1] = -id1;
 125     }
 126
 127   // Resolve valence content for assumed baryon. Flunk unallowed codes.
 128   } else {
 129     int id1 = idBeamAbs/1000;
 130     int id2 = (idBeamAbs/100)%10;
 131     int id3 = (idBeamAbs/10)%10;
 132     if ( id1 < 1 || id1 > maxValQuark || id2 < 1 || id2 > maxValQuark
 133       || id3 < 1 || id3 > maxValQuark) return;
 134     if (id2 > id1 || id3 > id1) return;
 135     isBaryonBeam = true;
 136
 137     // Store valence content of a confirmed baryon.
 138     nValKinds = 1; idVal[0] = id1; nVal[0] = 1;
 139     if (id2 == id1) ++nVal[0];
 140     else {
 141       nValKinds = 2;
 142       idVal[1]  = id2;
 143       nVal[1]   = 1;
 144     }
 145     if (id3 == id1) ++nVal[0];
 146     else if (id3 == id2) ++nVal[1];
 147     else {
 148       idVal[nValKinds] = id3;
 149       nVal[nValKinds]  = 1;
 150       ++nValKinds;
 151     }
 152   }
 153
 154   // Flip flavours for antimeson or antibaryon, and then done.
 155   if (idBeam < 0) for (int i = 0; i < nValKinds; ++i) idVal[i] = -idVal[i];
 156   isHadronBeam = true;
 157   Q2ValFracSav = -1.;
 158
 159 }
 160
 161 //*********
 162
 163 double BeamParticle::xMax(int iSkip) {
 164
 165   // Minimum requirement on remaining energy > nominal mass for hadron.
 166   double xLeft = 1.;
 167   if (isHadron()) xLeft -= m() / e();
 168   if (size() == 0) return xLeft;
 169
 170   // Subtract what was carried away by initiators (to date).
 171   for (int i = 0; i < size(); ++i)
 172     if (i != iSkip) xLeft -= resolved[i].x();
 173   return xLeft;
 174
 175 }
 176
 177 //*********
 178
 179 // Parton distributions, reshaped to take into account previous
 180 // multiple interactions. By picking a non-negative iSkip value,
 181 // one particular interaction is skipped, as needed for ISR
 182
 183 double BeamParticle::xfModified(int iSkip, int idIn, double x, double Q2) {
 184
 185   // Initial values.
 186   idSave    = idIn;
 187   iSkipSave = iSkip;
 188   xqVal     = 0.;
 189   xqgSea    = 0.;
 190   xqCompSum = 0.;
 191
 192   // Fast procedure for first interaction.
 193   if (size() == 0) {
 194     if (x >= 1.) return 0.;
 195     bool canBeVal = false;
 196     for (int i = 0; i < nValKinds; ++i)
 197       if (idIn == idVal[i]) canBeVal = true;
 198     if (canBeVal) {
 199       xqVal     = xfVal( idIn, x, Q2);
 200       xqgSea    = xfSea( idIn, x, Q2);
 201     }
 202     else xqgSea = xf( idIn, x, Q2);
 203
 204   // More complicated procedure for non-first interaction.
 205   } else {
 206
 207     // Sum up the x already removed, and check that remaining x is enough.
 208     double xUsed = 0.;
 209     for (int i = 0; i < size(); ++i)
 210       if (i != iSkip) xUsed += resolved[i].x();
 211     double xLeft = 1. - xUsed;
 212     if (x >= xLeft) return 0.;
 213     double xRescaled = x / xLeft;
 214
 215     // Calculate total and remaining amount of x carried by valence quarks.
 216     double xValTot = 0.;
 217     double xValLeft = 0.;
 218     for (int i = 0; i < nValKinds; ++i) {
 219       nValLeft[i] = nVal[i];
 220       for (int j = 0; j < size(); ++j)
 221       if (j != iSkip && resolved[j].isValence()
 222         && resolved[j].id() == idVal[i]) --nValLeft[i];
 223       double xValNow =  xValFrac(i, Q2);
 224       xValTot += nVal[i] * xValNow;
 225       xValLeft += nValLeft[i] * xValNow;
 226     }
 227
 228     // Calculate total amount of x carried by unmatched companion quarks.
 229     double xCompAdded = 0.;
 230     for (int i = 0; i < size(); ++i)
 231     if (i != iSkip && resolved[i].isUnmatched()) xCompAdded
 232       += xCompFrac( resolved[i].x() / (xLeft + resolved[i].x()) )
 233       // Typo warning: extrafactor missing in Skands&Sjostrand article;
 234       // <x> for companion refers to fraction of x left INCLUDING sea quark.
 235       // To be modified further??
 236       * (1. + resolved[i].x() / xLeft);
 237
 238     // Calculate total rescaling factor and pdf for sea and gluon.
 239     double rescaleGS = max( 0., (1. - xValLeft - xCompAdded)
 240       / (1. - xValTot) );
 241     xqgSea = rescaleGS * xfSea( idIn, xRescaled, Q2);
 242
 243     // Find valence part and rescale it to remaining number of quarks.
 244     for (int i = 0; i < nValKinds; ++i)
 245     if (idIn == idVal[i] && nValLeft[i] > 0)
 246       xqVal = xfVal( idIn, xRescaled, Q2)
 247       * double(nValLeft[i]) / double(nVal[i]);
 248
 249     // Find companion part, by adding all companion contributions.
 250     for (int i = 0; i < size(); ++i)
 251     if (i != iSkip && resolved[i].id() == -idIn
 252       && resolved[i].isUnmatched()) {
 253       double xsRescaled = resolved[i].x() / (xLeft + resolved[i].x());
 254       double xcRescaled = x / (xLeft + resolved[i].x());
 255       double xqCompNow = xCompDist( xcRescaled, xsRescaled);
 256       resolved[i].xqCompanion( xqCompNow);
 257       xqCompSum += xqCompNow;
 258     }
 259   }
 260
 261   // Add total, but only return relevant part for ISR. More cases??
 262   // Watch out, e.g. g can come from either kind of quark.??
 263   xqgTot = xqVal + xqgSea + xqCompSum;
 264   if (iSkip >= 0) {
 265     if (resolved[iSkip].isValence()) return xqVal;
 266     if (resolved[iSkip].isUnmatched()) return xqgSea + xqCompSum;
 267   }
 268   return xqgTot;
 269
 270 }
 271
 272 //*********
 273
 274 // Decide whether a quark extracted from the beam is of valence, sea or
 275 // companion kind; in the latter case also pick its companion.
 276 // Assumes xfModified has already been called.
 277
 278   int BeamParticle::pickValSeaComp() {
 279
 280   // If parton already has a companion than reset code for this.
 281   int oldCompanion = resolved[iSkipSave].companion();
 282   if (oldCompanion >= 0) resolved[oldCompanion].companion(-2);
 283
 284   // Default assignment is sea.
 285   int vsc = -2;
 286
 287   // For gluons or photons no sense of valence or sea.
 288   if (idSave == 21 || idSave == 22) vsc = -1;
 289
 290   // For lepton beam assume same-kind lepton inside is valence.
 291   else if (isLeptonBeam && idSave == idBeam) vsc = -3;
 292
 293   // Decide if valence or sea quark.
 294   else {
 295     double xqRndm = xqgTot * Rndm::flat();
 296     if (xqRndm < xqVal) vsc = -3;
 297     else if (xqRndm < xqVal + xqgSea) vsc = -2;
 298
 299     // If not either, loop over all possible companion quarks.
 300     else {
 301       xqRndm -= xqVal + xqgSea;
 302       for (int i = 0; i < size(); ++i)
 303       if (i != iSkipSave && resolved[i].id() == -idSave
 304         && resolved[i].isUnmatched()) {
 305         xqRndm -= resolved[i].xqCompanion();
 306         if (xqRndm < 0.) vsc = i;
 307         break;
 308       }
 309     }
 310   }
 311
 312   // Bookkeep assignment; for sea--companion pair both ways.
 313   resolved[iSkipSave].companion(vsc);
 314   if (vsc >= 0) resolved[vsc].companion(iSkipSave);
 315
 316   // Done; return code for choice (to distinguish valence/sea in Info).
 317   return vsc;
 318
 319 }
 320
 321 //*********
 322
 323 // Fraction of hadron momentum sitting in a valence quark distribution.
 324 // Based on hardcoded parametrizations of CTEQ 5L numbers.
 325
 326 double BeamParticle::xValFrac(int j, double Q2) {
 327
 328   // Only recalculate when required.
 329   if (Q2 != Q2ValFracSav) {
 330     Q2ValFracSav = Q2;
 331
 332     // Q2-dependence of log-log form; assume fixed Lambda = 0.2.
 333     double llQ2 = log( log( max( 1., Q2) / 0.04 ));
 334
 335     // Fractions carried by u and d in proton.
 336     uValInt =  0.48 / (1. + 1.56 * llQ2);
 337     dValInt = 0.385 / (1. + 1.60 * llQ2);
 338   }
 339
 340   // Baryon with three different quark kinds: (2 * u + d) / 3 of proton.
 341   if (isBaryonBeam && nValKinds == 3) return (2. * uValInt + dValInt) / 3.;
 342
 343   // Baryon with one or two identical: like d or u of proton.
 344   if (isBaryonBeam && nVal[j] == 1) return dValInt;
 345   if (isBaryonBeam && nVal[j] == 2) return uValInt;
 346
 347   // Meson: (2 * u + d) / 2 of proton so same total valence quark fraction.
 348     return 0.5 * (2. * uValInt + dValInt);
 349
 350 }
 351
 352 //*********
 353
 354 // The momentum integral of a companion quark, with its partner at x_s,
 355 // using an approximate gluon density like (1 - x_g)^power / x_g.
 356 // The value corresponds to an unrescaled range between 0 and 1 - x_s.
 357
 358 double BeamParticle::xCompFrac(double xs) {
 359
 360   // Select case by power of gluon (1-x_g) shape.
 361   switch (companionPower) {
 362
 363     case 0:
 364        return xs * ( 5. + xs * (-9. - 2. * xs * (-3. + xs)) + 3. * log(xs) )
 365          / ( (-1. + xs) * (2. + xs * (-1. + 2. * xs)) );
 366
 367     case 1:
 368        return -1. -3. * xs + ( 2. * pow2(-1. + xs) * (1. + xs + xs*xs))
 369          / ( 2. + xs*xs * (xs - 3.) + 3. * xs * log(xs) );
 370
 371     case 2:
 372        return xs * ( (1. - xs) * (19. + xs * (43. + 4. * xs))
 373          + 6. * log(xs) * (1. + 6. * xs + 4.*xs*xs) ) /
 374         ( 4. * ( (xs - 1.) * (1. + xs * (4. + xs) )
 375         - 3. * xs * log(xs) * (1 + xs) ) );
 376
 377     case 3:
 378       return 3. * xs * ( (xs - 1.) * (7. + xs * (28. + 13. * xs))
 379         - 2. * log(xs) * (1. + xs * (9. + 2. * xs * (6. + xs))) )
 380         / ( 4. + 27. * xs - 31. * pow3(xs)
 381         + 6. * xs * log(xs) * (3. + 2. * xs * (3.+xs)) );
 382
 383     default:
 384       return ( -9. * xs * (xs*xs - 1.) * (5. + xs * (24. + xs)) + 12. * xs
 385         * log(xs) * (1. + 2. * xs) * (1. + 2. * xs * (5. + 2. * xs)) )
 386         / ( 8. * (1. + 2. * xs) * ((xs - 1.) * (1. + xs * (10. + xs))
 387         - 6. * xs * log(xs) * (1. + xs)) );
 388
 389   }
 390 }
 391
 392 //*********
 393
 394 // The x*f pdf of a companion quark at x_c, with its sea partner at x_s,
 395 // using an approximate gluon density like (1 - x_g)^power / x_g.
 396 // The value corresponds to an unrescaled range between 0 and 1 - x_s.
 397
 398 double BeamParticle::xCompDist(double xc, double xs) {
 399
 400   // Mother gluon momentum fraction. Check physical limit.
 401   double xg = xc + xs;
 402   if (xg > 1.) return 0.;
 403
 404   // Common factor, including splitting kernel and part of gluon density
 405   // (and that it is x_c * f that is coded).
 406   double fac = 3. * xc * xs * (xc*xc + xs*xs) / pow4(xg);
 407
 408   // Select case by power of gluon (1-x_g) shape.
 409   switch (companionPower) {
 410
 411     case 0:
 412       return fac / ( 2. - xs * (3. - xs * (3. - 2. * xs)) );
 413
 414     case 1:
 415       return fac * (1. - xg) / ( 2. + xs*xs * (-3. + xs) + 3. * xs * log(xs) );
 416
 417     case 2:
 418       return fac * pow2(1. - xg) / ( 2. * ((1. - xs) * (1. + xs * (4. + xs))
 419         + 3. * xs * (1. + xs) * log(xs)) );
 420
 421     case 3:
 422       return fac * pow3(1. - xg) * 2. / ( 4. + 27. * xs - 31. * pow3(xs)
 423         + 6. * xs * log(xs) * (3. + 2. * xs * (3. + xs)) );
 424
 425     default:
 426        return fac * pow4(1. - xg) / ( 2. * (1. + 2. * xs) * ((1. - xs)
 427          * (1. + xs * (10. + xs)) + 6. * xs * log(xs) * (1. + xs)) );
 428
 429   }
 430 }
 431
 432 //*********
 433
 434 // Add required extra remnant flavour content. Also initial colours.
 435
 436 bool BeamParticle::remnantFlavours(Event& event) {
 437
 438   // A baryon will have a junction, unless a diquark is formed later.
 439   hasJunctionBeam = (isBaryon());
 440
 441   // Store how many hard-scattering partons were removed from beam.
 442   nInit = size();
 443
 444   // Find remaining valence quarks.
 445   for (int i = 0; i < nValKinds; ++i) {
 446     nValLeft[i] = nVal[i];
 447     for (int j = 0; j < nInit; ++j) if (resolved[j].isValence()
 448       && resolved[j].id() == idVal[i]) --nValLeft[i];
 449     // Add remaining valence quarks to record. Partly temporary values.
 450     for (int k = 0; k < nValLeft[i]; ++k) append(0, idVal[i], 0., -3);
 451   }
 452
 453   // If at least two valence quarks left in baryon remnant then form diquark.
 454   int nInitPlusVal = size();
 455   if (isBaryon() && nInitPlusVal - nInit >= 2) {
 456
 457     // If three, pick two at random to form diquark, else trivial.
 458     int iQ1 = nInit;
 459     int iQ2 = nInit + 1;
 460     if (nInitPlusVal - nInit == 3) {
 461       double pickDq = 3. * Rndm::flat();
 462       if (pickDq > 1.) iQ2 = nInit + 2;
 463       if (pickDq > 2.) iQ1 = nInit + 1;
 464     }
 465
 466     // Pick spin 0 or 1 according to SU(6) wave function factors.
 467     int idDq = flavSelPtr->makeDiquark( resolved[iQ1].id(),
 468       resolved[iQ2].id(), idBeam);
 469
 470     // Overwrite with diquark flavour and remove one slot. No more junction.
 471     resolved[iQ1].id(idDq);
 472     if (nInitPlusVal - nInit == 3 && iQ2 == nInit + 1)
 473       resolved[nInit + 1].id( resolved[nInit + 2].id() );
 474     resolved.pop_back();
 475     hasJunctionBeam = false;
 476   }
 477
 478   // Find companion quarks to unmatched sea quarks.
 479   for (int i = 0; i < nInit; ++i)
 480   if (resolved[i].isUnmatched()) {
 481
 482     // Add companion quark to record; and bookkeep both ways.
 483     append(0, -resolved[i].id(), 0., i);
 484     resolved[i].companion(size() - 1);
 485   }
 486
 487   // If no other remnants found, add a gluon or photon to carry momentum.
 488   if (size() == nInit) {
 489     int    idRemnant = (isHadronBeam) ? 21 : 22;
 490     append(0, idRemnant, 1., -1);
 491   }
 492
 493   // Set initiator and remnant masses.
 494   for (int i = 0; i < size(); ++i) {
 495     if (i < nInit) resolved[i].m(0.);
 496     else resolved[i].m( ParticleDataTable::m0( resolved[i].id() ) );
 497   }
 498
 499   // For debug purposes: reject beams with resolved junction topology.
 500   if (hasJunctionBeam && !allowJunction) return false;
 501
 502   // Pick initial colours for remnants.
 503   for (int i = nInit; i < size(); ++i) {
 504     int colType = ParticleDataTable::colType( resolved[i].id() );
 505     int col = (colType == 1 || colType == 2) ? event.nextColTag() : 0;
 506     int acol = (colType == -1 || colType == 2) ? event.nextColTag() : 0;
 507     resolved[i].cols( col, acol);
 508   }
 509
 510   // Done.
 511   return true;
 512
 513 }
 514
 515 //*********
 516
 517 // Correlate all initiators and remnants to make a colour singlet.
 518
 519 bool BeamParticle::remnantColours(Event& event, vector<int>& colFrom,
 520   vector<int>& colTo) {
 521
 522   // No colours in lepton beams so no need to do anything.
 523   if (isLeptonBeam) return true;
 524
 525   // Copy initiator colour info from the event record to the beam.
 526   for (int i = 0; i < size(); ++i) {
 527     int j =  resolved[i].iPos();
 528     resolved[i].cols( event[j].col(), event[j].acol());
 529   }
 530
 531   // Find number and position of valence quarks, of gluons, and
 532   // of sea-companion pairs (counted as gluons) in the beam remnants.
 533   // Skip gluons with same colour as anticolour.
 534   vector<int> iVal;
 535   vector<int> iGlu;
 536   for (int i = 0; i < size(); ++i) {
 537     if ( resolved[i].isValence() ) iVal.push_back(i);
 538     else if ( resolved[i].isCompanion() && resolved[i].companion() > i )
 539       iGlu.push_back(i);
 540     else if ( resolved[i].id() == 21
 541       && resolved[i].col() != resolved[i].acol() ) iGlu.push_back(i);
 542   }
 543
 544   // Pick a valence quark to which gluons are attached.
 545   // Do not resolve quarks in diquark. (More sophisticated??)
 546   int iValSel= iVal[0];
 547   if (iVal.size() == 2) {
 548     if ( abs(resolved[iValSel].id()) > 10 ) iValSel = iVal[1];
 549   } else {
 550     double rndmValSel = 3. * Rndm::flat();
 551     if (rndmValSel > 1.) iValSel= iVal[1];
 552     if (rndmValSel > 2.) iValSel= iVal[2];
 553   }
 554
 555   // This valence quark defines initial (anti)colour.
 556   int iBeg = iValSel;
 557   bool hasCol = (resolved[iBeg].col() > 0);
 558   int begCol = (hasCol) ? resolved[iBeg].col() : resolved[iBeg].acol();
 559
 560   // Do random stepping through gluon/(sea+companion) list.
 561   vector<int> iGluRndm;
 562   for (int i = 0; i < int(iGlu.size()); ++i)
 563     iGluRndm.push_back( iGlu[i] );
 564   for (int iOrder = 0; iOrder < int(iGlu.size()); ++iOrder) {
 565     int iRndm = int( double(iGluRndm.size()) * Rndm::flat());
 566     int iGluSel = iGluRndm[iRndm];
 567     iGluRndm[iRndm] = iGluRndm[iGluRndm.size() - 1];
 568     iGluRndm.pop_back();
 569
 570     // Find matching anticolour/colour to current colour/anticolour.
 571     int iEnd = iGluSel;
 572     int endCol = (hasCol) ? resolved[iEnd].acol() : resolved[iEnd].col();
 573     // Not gluon but sea+companion pair: go to other.
 574     if (endCol == 0) {
 575       iEnd = resolved[iEnd].companion();
 576       endCol = (hasCol) ? resolved[iEnd].acol() : resolved[iEnd].col();
 577     }
 578
 579     // Collapse this colour-anticolour pair to the lowest one.
 580     if (begCol < endCol) {
 581       if (hasCol) resolved[iEnd].acol(begCol);
 582       else resolved[iEnd].col(begCol);
 583       colFrom.push_back(endCol);
 584       colTo.push_back(begCol);
 585     } else {
 586       if (hasCol) resolved[iBeg].col(endCol);
 587       else resolved[iBeg].acol(endCol);
 588       colFrom.push_back(begCol);
 589       colTo.push_back(endCol);
 590     }
 591
 592     // Pick up the other colour of the recent gluon and repeat.
 593     iBeg = iEnd;
 594     begCol = (hasCol) ? resolved[iBeg].col() : resolved[iBeg].acol();
 595     // Not gluon but sea+companion pair: go to other.
 596     if (begCol == 0) {
 597       iBeg = resolved[iBeg].companion();
 598       begCol = (hasCol) ? resolved[iBeg].col() : resolved[iBeg].acol();
 599     }
 600
 601   // At end of gluon/(sea+companion) list.
 602   }
 603
 604   // Now begin checks, and also finding junction information.
 605   // Loop through remnant partons; isolate all colours and anticolours.
 606   vector<int> colList;
 607   vector<int> acolList;
 608   for (int i = 0; i < size(); ++i)
 609   if ( resolved[i].col() != resolved[i].acol() ) {
 610     if (resolved[i].col() > 0) colList.push_back( resolved[i].col() );
 611     if (resolved[i].acol() > 0) acolList.push_back( resolved[i].acol() );
 612   }
 613
 614   // Remove all matching colour-anticolour pairs.
 615   bool foundPair = true;
 616   while (foundPair && colList.size() > 0 && acolList.size() > 0) {
 617     foundPair = false;
 618     for (int iCol = 0; iCol < int(colList.size()); ++iCol) {
 619       for (int iAcol = 0; iAcol < int(acolList.size()); ++iAcol) {
 620         if (acolList[iAcol] == colList[iCol]) {
 621           colList[iCol] = colList.back(); colList.pop_back();
 622           acolList[iAcol] = acolList.back(); acolList.pop_back();
 623           foundPair = true; break;
 624         }
 625       } if (foundPair) break;
 626     }
 627   }
 628
 629   // Usually one unmatched pair left to collapse.
 630   if (colList.size() == 1 && acolList.size() == 1) {
 631     int finalFrom = max( colList[0], acolList[0]);
 632     int finalTo   = min( colList[0], acolList[0]);
 633     for (int i = 0; i < size(); ++i) {
 634       if (resolved[i].col()  == finalFrom) resolved[i].col(finalTo);
 635       if (resolved[i].acol() == finalFrom) resolved[i].acol(finalTo);
 636     }
 637     colFrom.push_back(finalFrom);
 638     colTo.push_back(finalTo);
 639
 640   // Store an (anti)junction when three (anti)coloured daughters.
 641   } else if (hasJunctionBeam && colList.size() == 3
 642     && acolList.size() == 0) {
 643     event.appendJunction( 1, colList[0], colList[1], colList[2]);
 644     junCol[0] = colList[0];
 645     junCol[1] = colList[1];
 646     junCol[2] = colList[2];
 647   } else if (hasJunctionBeam && acolList.size() == 3
 648     && colList.size() == 0) {
 649     event.appendJunction( 2, acolList[0], acolList[1], acolList[2]);
 650     junCol[0] = acolList[0];
 651     junCol[1] = acolList[1];
 652     junCol[2] = acolList[2];
 653
 654   // Any other nonvanishing values indicate failure.
 655   } else if (colList.size() > 0 || acolList.size() > 0) {
 656     infoPtr->errorMsg("Error in BeamParticle::remnantColours: "
 657       "leftover unmatched colours");
 658     return false;
 659   }
 660
 661   // Store colour assignment of beam particles.
 662   for (int i = nInit; i < size(); ++i)
 663     event[resolved[i].iPos()].cols( resolved[i].col(), resolved[i].acol() );
 664
 665   // Done.
 666   return true;
 667 }
 668
 669
 670 //*********
 671
 672 // Pick unrescaled x values for beam remnant sharing.
 673
 674 double BeamParticle::xRemnant( int i) {
 675
 676   double x = 0.;
 677
 678   // Calculation of x of valence quark or diquark, for latter as sum.
 679   if (resolved[i].isValence()) {
 680
 681     // Resolve diquark into sum of two quarks.
 682     int id1 = resolved[i].id();
 683     int id2 = 0;
 684     if (abs(id1) > 10) {
 685       id2 = (id1 > 0) ? (id1/100)%10 : -(((-id1)/100)%10);
 686       id1 = (id1 > 0) ? id1/1000 : -((-id1)/1000);
 687     }
 688
 689     // Loop over (up to) two quarks; add their contributions.
 690     for (int iId = 0; iId < 2; ++iId) {
 691       int idNow = (iId == 0) ? id1 : id2;
 692       if (idNow == 0) break;
 693       double xPart = 0.;
 694
 695       // Assume form (1-x)^a / sqrt(x).
 696       double xPow = valencePowerMeson;
 697       if (isBaryonBeam) {
 698         if (nValKinds == 3 || nValKinds == 1)
 699           xPow = (3. * Rndm::flat() < 2.)
 700             ? valencePowerUinP : valencePowerDinP ;
 701         else if (nValence(idNow) == 2) xPow = valencePowerUinP;
 702         else xPow = valencePowerDinP;
 703       }
 704       do xPart = pow2( Rndm::flat() );
 705       while ( pow(1. - xPart, xPow) < Rndm::flat() );
 706
 707       // End loop over (up to) two quarks. Possibly enhancement for diquarks.
 708       x += xPart;
 709     }
 710    if (id2 != 0) x *= valenceDiqEnhance;
 711
 712   // Calculation of x of sea quark, based on companion association.
 713   } else if (resolved[i].isCompanion()) {
 714
 715     // Find rescaled x value of companion.
 716     double xLeft = 1.;
 717     for (int iInit = 0; iInit < nInit; ++iInit)
 718       xLeft -= resolved[iInit].x();
 719     double xCompanion = resolved[ resolved[i].companion() ].x();
 720     xCompanion /= (xLeft + xCompanion);
 721
 722     // Now use ansatz q(x; x_c) < N/(x +x_c) to pick x.
 723     do x = pow( xCompanion, Rndm::flat()) - xCompanion;
 724     while ( pow( (1. - x - xCompanion) / (1. - xCompanion), companionPower)
 725       * (pow2(x) + pow2(xCompanion)) / pow2(x + xCompanion) < Rndm::flat() );
 726
 727   // Else, rarely, a single gluon remnant, so value does not matter.
 728   } else x = 1.;
 729   return x;
 730
 731 }
 732
 733 //*********
 734
 735 // Print the list of resolved partons in a beam.
 736
 737 void BeamParticle::list(ostream& os) {
 738
 739   // Header.
 740   os << "\n --------  PYTHIA Partons resolved in beam  ------------"
 741      << "--------------------------------------------------------\n"
 742      << "\n    i  iPos      id       x    comp   xqcomp      colours"
 743      << "      p_x        p_y        p_z         e          m \n";
 744
 745   // Loop over list of removed partons and print it.
 746   double xSum = 0.;
 747   Vec4 pSum;
 748   for (int i = 0; i < size(); ++i) {
 749     ResolvedParton res = resolved[i];
 750     os << fixed << setprecision(6) << setw(5) << i << setw(6) << res.iPos()
 751        << setw(8) << res.id() << setw(10) << res.x() << setw(6)
 752        << res.companion() << setw(10) << res.xqCompanion()
 753        << setprecision(3) << setw(6) << res.col() << setw(6) << res.acol()
 754        << setw(11) << res.px() << setw(11) << res.py() << setw(11)
 755        << res.pz() << setw(11) << res.e() << setw(11) << res.m() << "\n";
 756
 757     // Also find and print sum of x and p values. Endline.
 758     xSum += res.x();
 759     pSum += res.p();
 760   }
 761   os << setprecision(6) << "             x sum:" << setw(10) << xSum
 762      << setprecision(3) << "                      p sum:" << setw(11)
 763      << pSum.px() << setw(11) << pSum.py() << setw(11) << pSum.pz()
 764      << setw(11) << pSum.e()
 765      << "\n\n --------  End PYTHIA Partons resolved in beam  ------"
 766      << "----------------------------------------------------------"
 767      << endl;
 768 }
 769
 770 //*********
 771
 772 // Test whether a lepton is to be considered as unresolved.
 773
 774 bool BeamParticle::isUnresolvedLepton() {
 775
 776   // Require record to consist of lepton with full energy plus a photon.
 777   if (!isLeptonBeam || resolved.size() > 2 || resolved[1].id() != 22
 778     || resolved[0].x() < XMINUNRESOLVED) return false;
 779   return true;
 780
 781 }
 782
 783 //*********
 784
 785 // For a diffractive system, decide whether to kick out gluon or quark.
 786
 787 bool BeamParticle::pickGluon(double mDiff) {
 788
 789   // Relative weight to pick a quark, assumed falling with energy.
 790   double probPickQuark = pickQuarkNorm / pow( mDiff, pickQuarkPower);
 791   return  ( (1. + probPickQuark) * Rndm::flat() < 1. );
 792
 793 }
 794
 795 //*********
 796
 797 // Pick a valence quark at random. (Used for diffractive systems.)
 798
 799 int BeamParticle::pickValence() {
 800
 801   // Pick one valence quark at random.
 802   int nTotVal = (isBaryonBeam) ? 3 : 2;
 803   double rnVal = Rndm::flat() * nTotVal;
 804   int iVal = (rnVal < 1.) ? 1 : ( (rnVal < 2.) ? 2 : 3 );
 805
 806   // This valence in slot 1, the rest thereafter.
 807   idVal1 = 0;
 808   idVal2 = 0;
 809   idVal3 = 0;
 810   int iNow = 0;
 811   for (int i = 0; i < nValKinds; ++i)
 812   for (int j = 0; j < nVal[i]; ++j) {
 813     ++iNow;
 814     if (iNow == iVal) idVal1 = idVal[i];
 815     else if ( idVal2 == 0) idVal2 = idVal[i];
 816     else idVal3 = idVal[i];
 817   }
 818
 819   // Construct diquark if baryon.
 820   if (idVal3 != 0) idVal2 = flavSelPtr->makeDiquark( idVal2, idVal3);
 821
 822   // Done.
 823   return idVal1;
 824
 825 }
 826
 827 //*********
 828
 829 // Share lightcone momentum between two remnants in a diffractive system.
 830
 831 double BeamParticle::zShare( double mDiff, double m1, double m2) {
 832
 833   // Set up as valence in normal beam so can use xRemnant code.
 834   append(0, idVal1, 0., -3);
 835   append(0, idVal2, 0., -3);
 836   double m2Diff = mDiff*mDiff;
 837
 838   // Begin to generate z and pT until acceptable solution.
 839   double wtAcc = 0.;
 840   do {
 841     double x1 = xRemnant(0);
 842     double x2 = xRemnant(0);
 843     zRel = x1 / (x1 + x2);
 844     pxRel = diffPrimKTwidth * Rndm::gauss();
 845     pyRel = diffPrimKTwidth * Rndm::gauss();
 846
 847     // Suppress large invariant masses of remnant system.
 848     double mTS1 = m1*m1 + pxRel*pxRel + pyRel*pyRel;
 849     double mTS2 = m2*m2 + pxRel*pxRel + pyRel*pyRel;
 850     double m2Sys = mTS1 / zRel + mTS2 / (1. - zRel);
 851     wtAcc = (m2Sys < m2Diff)
 852       ? pow( 1. - m2Sys / m2Diff, diffLargeMassSuppress) : 0.;
 853   } while (wtAcc < Rndm::flat());
 854
 855   // Done.
 856   return zRel;
 857
 858 }
 859
 860 //**************************************************************************
 861
 862 } // end namespace Pythia8