PYTHIA8/pythia8130/src/FragmentationSystems.cxx

   1 // FragmentationSystems.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 // ColConfig, StringRegion and StringSystem classes.
   8
   9 #include "FragmentationSystems.h"
  10
  11 namespace Pythia8 {
  12
  13 //**************************************************************************
  14
  15 // The ColConfig class.
  16
  17 //*********
  18
  19 // Initialize and save pointers.
  20
  21 void ColConfig::init(StringFlav* flavSelPtrIn) {
  22
  23   // Save pointer.
  24   flavSelPtr    = flavSelPtrIn;
  25
  26   // Joining of nearby partons along the string.
  27   mJoin         = Settings::parm("FragmentationSystems:mJoin");
  28
  29   // For consistency ensure that mJoin is bigger than in StringRegion.
  30   mJoin         = max( mJoin, 2. * StringRegion::MJOIN);
  31
  32   // Simplification of q q q junction topology to quark - diquark one.
  33   mJoinJunction = Settings::parm("FragmentationSystems:mJoinJunction");
  34   mStringMin    = Settings::parm("HadronLevel:mStringMin");
  35
  36 }
  37
  38 //*********
  39
  40 // Insert a new colour singlet system in ascending mass order.
  41 // Calculate its properties. Join nearby partons.
  42
  43 void ColConfig::insert( vector<int>& iPartonIn, Event& event) {
  44
  45   // Find momentum and invariant mass of system, minus endpoint masses.
  46   Vec4 pSumIn;
  47   double mSumIn = 0.;
  48   bool hasJunctionIn = false;
  49   for (int i = 0; i < int(iPartonIn.size()); ++i) {
  50     if (iPartonIn[i] < 0) {
  51       hasJunctionIn = true;
  52       continue;
  53     }
  54     pSumIn += event[ iPartonIn[i] ].p();
  55     if (!event[ iPartonIn[i] ].isGluon())
  56       mSumIn += event[ iPartonIn[i] ].constituentMass();
  57   }
  58   double massIn = pSumIn.mCalc();
  59   double massExcessIn = massIn - mSumIn;
  60
  61   // Identify closed gluon loop. Assign "endpoint" masses as light quarks.
  62   bool isClosedIn = (iPartonIn[0] >= 0 && event[ iPartonIn[0] ].isGluon());
  63   if (isClosedIn) massExcessIn -= 2. * ParticleDataTable::constituentMass(1);
  64
  65   // For junction topology: join two nearby legs into a diquark.
  66   if (hasJunctionIn && joinJunction( iPartonIn, event, massExcessIn))
  67     hasJunctionIn = false;
  68
  69   // Loop while > 2 partons left and hope of finding joining pair.
  70   bool hasJoined = true;
  71   while (hasJoined && iPartonIn.size() > 2) {
  72
  73     // Look for the pair of neighbour partons (along string) with
  74     // the smallest invariant mass (subtracting quark masses).
  75     int iJoinMin = -1;
  76     double mJoinMin = 2. * mJoin;
  77     int nSize = iPartonIn.size();
  78     int nPair = (isClosedIn) ? nSize : nSize - 1;
  79     for (int i = 0; i < nPair; ++i) {
  80       // Keep three legs of junction separate.
  81       if (iPartonIn[i] < 0 || iPartonIn[(i + 1)%nSize] < 0) continue;
  82       Particle& parton1 = event[ iPartonIn[i] ];
  83       Particle& parton2 = event[ iPartonIn[(i + 1)%nSize] ];
  84       Vec4 pSumNow;
  85       pSumNow += (parton1.isGluon()) ? 0.5 * parton1.p() : parton1.p();
  86       pSumNow += (parton2.isGluon()) ? 0.5 * parton2.p() : parton2.p();
  87       double mJoinNow = pSumNow.mCalc();
  88       if (!parton1.isGluon()) mJoinNow -= parton1.m();
  89       if (!parton2.isGluon()) mJoinNow -= parton2.m();
  90       if (mJoinNow < mJoinMin) { iJoinMin = i; mJoinMin = mJoinNow; }
  91     }
  92
  93     // If sufficiently nearby then join into one new parton.
  94     // Note: error sensitivity to mJoin indicates unstable precedure??
  95     hasJoined = false;
  96     if (mJoinMin < mJoin) {
  97       int iJoin1  = iPartonIn[iJoinMin];
  98       int iJoin2  = iPartonIn[(iJoinMin + 1)%nSize];
  99       int idNew   = (event[iJoin1].isGluon()) ? event[iJoin2].id()
 100                                               : event[iJoin1].id();
 101       int colNew  = event[iJoin1].col();
 102       int acolNew = event[iJoin2].acol();
 103       if (colNew == acolNew) {
 104         colNew    = event[iJoin2].col();
 105         acolNew   = event[iJoin1].acol();
 106       }
 107       Vec4 pNew   = event[iJoin1].p() + event[iJoin2].p();
 108
 109       // Append joined parton to event record.
 110       int iNew = event.append( idNew, 73, min(iJoin1, iJoin2),
 111         max(iJoin1, iJoin2), 0, 0, colNew, acolNew, pNew, pNew.mCalc() );
 112
 113       // Mark joined partons and reduce remaining system.
 114       event[iJoin1].statusNeg();
 115       event[iJoin2].statusNeg();
 116       event[iJoin1].daughter1(iNew);
 117       event[iJoin2].daughter1(iNew);
 118       if (iJoinMin == nSize - 1) iPartonIn[0] = iNew;
 119       else {
 120         iPartonIn[iJoinMin] = iNew;
 121         for (int i = iJoinMin + 1; i < nSize - 1; ++i)
 122           iPartonIn[i] = iPartonIn[i + 1];
 123       }
 124       iPartonIn.pop_back();
 125
 126       // If joined,then loopback to look for more.
 127       hasJoined = true;
 128     }
 129   }
 130
 131   // Store new colour singlet system at the end.
 132   singlets.push_back( ColSinglet(iPartonIn, pSumIn, massIn,
 133     massExcessIn, hasJunctionIn, isClosedIn) );
 134
 135   // Now move around, so that smallest mass excesses come first.
 136   int iInsert = singlets.size() - 1;
 137   for (int iSub = singlets.size() - 2; iSub >= 0; --iSub) {
 138     if (massExcessIn > singlets[iSub].massExcess) break;
 139     singlets[iSub + 1] = singlets[iSub];
 140     iInsert = iSub;
 141   }
 142   if (iInsert < int(singlets.size()) - 1) singlets[iInsert] =
 143     ColSinglet(iPartonIn, pSumIn, massIn, massExcessIn,
 144     hasJunctionIn, isClosedIn);
 145
 146 }
 147
 148 //*********
 149
 150 // Join two legs of junction to a diquark for small invariant masses.
 151 // Note: for junction system, iPartonIn points to structure
 152 // (-code0) g...g.q0 (-code1) g...g.q1 (-code2) g...g.q2
 153
 154 bool ColConfig::joinJunction( vector<int>& iPartonIn, Event& event,
 155   double massExcessIn) {
 156
 157   // Find four-momentum and endpoint quarks and masses on the three legs.
 158   Vec4   pLeg[3];
 159   double mLeg[3];
 160   int    idAbsLeg[3];
 161   int leg = -1;
 162   for (int i = 0; i < int(iPartonIn.size()); ++ i) {
 163     if (iPartonIn[i] < 0) ++leg;
 164     else {
 165       pLeg[leg]    += event[ iPartonIn[i] ].p();
 166       mLeg[leg]     = event[ iPartonIn[i] ].m();
 167       idAbsLeg[leg] = event[ iPartonIn[i] ].idAbs();
 168     }
 169   }
 170
 171   // Calculate invariant mass of three pairs, minus endpoint masses.
 172   double m01  = (pLeg[0] + pLeg[1]).mCalc() - mLeg[0] - mLeg[1];
 173   double m02  = (pLeg[0] + pLeg[2]).mCalc() - mLeg[0] - mLeg[2];
 174   double m12  = (pLeg[1] + pLeg[2]).mCalc() - mLeg[1] - mLeg[2];
 175
 176   // Find lowest-mass pair not involving diquark.
 177   double mMin = mJoinJunction + 1.;
 178   int    legA = -1;
 179   int    legB = -1;
 180   if (m01 < mMin && idAbsLeg[0] < 9 && idAbsLeg[1] < 9) {
 181     mMin = m01;
 182     legA = 0;
 183     legB = 1;
 184   }
 185   if (m02 < mMin && idAbsLeg[0] < 9 && idAbsLeg[2] < 9) {
 186     mMin = m02;
 187     legA = 0;
 188     legB = 2;
 189   }
 190   if (m12 < mMin && idAbsLeg[1] < 9 && idAbsLeg[2] < 9) {
 191     mMin = m12;
 192     legA = 1;
 193     legB = 2;
 194   }
 195   int legC = 3 - legA - legB;
 196
 197   // Nothing to do if no two legs have small invariant mass, and
 198   // system as a whole is above MiniStringFragmentation threshold.
 199   if (mMin > mJoinJunction && massExcessIn > mStringMin) return false;
 200
 201   // Construct separate index arrays for the three legs.
 202   vector<int> iLegA, iLegB, iLegC;
 203   leg = -1;
 204   for (int i = 0; i < int(iPartonIn.size()); ++ i) {
 205     if (iPartonIn[i] < 0) ++leg;
 206     else if( leg == legA) iLegA.push_back( iPartonIn[i] );
 207     else if( leg == legB) iLegB.push_back( iPartonIn[i] );
 208     else if( leg == legC) iLegC.push_back( iPartonIn[i] );
 209   }
 210
 211   // First step: successively combine any gluons on the two legs.
 212   // (Presumably overkill; not likely to be (m)any extra gluons.)
 213   // (Do as successive binary joinings, so only need two mothers.)
 214   for (leg = 0; leg < 2; ++leg) {
 215     vector<int>& iLegNow = (leg == 0) ? iLegA : iLegB;
 216     int sizeNow = iLegNow.size();
 217     for (int i = sizeNow - 2; i >= 0; --i) {
 218       int iQ = iLegNow.back();
 219       int iG = iLegNow[i];
 220       int colNew = (event[iQ].id() > 0) ? event[iG].col() : 0;
 221       int acolNew = (event[iQ].id() < 0) ? event[iG].acol() : 0;
 222       Vec4 pNew = event[iQ].p() + event[iG].p();
 223       int iNew = event.append( event[iQ].id(), 74, min(iQ, iG),
 224         max(iQ, iG), 0, 0, colNew, acolNew, pNew, pNew.mCalc() );
 225
 226       // Mark joined partons and update iLeg end.
 227       event[iQ].statusNeg();
 228       event[iG].statusNeg();
 229       event[iQ].daughter1(iNew);
 230       event[iG].daughter1(iNew);
 231       iLegNow.back() = iNew;
 232     }
 233   }
 234
 235   // Second step: combine two quarks into a diquark.
 236   int iQA     = iLegA.back();
 237   int iQB     = iLegB.back();
 238   int idQA    = event[iQA].id();
 239   int idQB    = event[iQB].id();
 240   int idNew   = flavSelPtr->makeDiquark( idQA, idQB );
 241   // Diquark colour is opposite to parton closest to junction on third leg.
 242   int colNew  = (idNew > 0) ? 0 : event[ iLegC[0] ].acol();
 243   int acolNew = (idNew > 0) ? event[ iLegC[0] ].col() : 0;
 244   Vec4 pNew   = pLeg[legA] + pLeg[legB];
 245   int iNew    = event.append( idNew, 74, min(iQA, iQB), max( iQA, iQB),
 246      0, 0, colNew, acolNew, pNew, pNew.mCalc() );
 247
 248   // Mark joined partons and reduce remaining system.
 249   event[iQA].statusNeg();
 250   event[iQB].statusNeg();
 251   event[iQA].daughter1(iNew);
 252   event[iQB].daughter1(iNew);
 253   iPartonIn.resize(0);
 254   iPartonIn.push_back( iNew);
 255   for (int i = 0; i < int(iLegC.size()) ; ++i)
 256     iPartonIn.push_back( iLegC[i]);
 257
 258   // Remove junction from event record list, identifying by colour.
 259   int iJun = -1;
 260   for (int i = 0; i < event.sizeJunction(); ++i)
 261     for (int j = 0; j < 3; ++ j)
 262       if ( event.colJunction(i,j) == max(colNew, acolNew) ) iJun = i;
 263   if (iJun >= 0) event.eraseJunction(iJun);
 264
 265   // Done, having eliminated junction.
 266   return true;
 267
 268 }
 269
 270 //*********
 271
 272 // Collect all partons of singlet to be consecutively ordered.
 273
 274 void ColConfig::collect(int iSub, Event& event) {
 275
 276   // Partons may already have been collected, e.g. at ministring collapse.
 277   if (singlets[iSub].isCollected) return;
 278   singlets[iSub].isCollected = true;
 279
 280   // Check if partons already "by chance" happen to be ordered.
 281   bool inOrder = true;
 282   for (int i = 0; i < singlets[iSub].size() - 1; ++i) {
 283     int iFirst = singlets[iSub].iParton[i];
 284     if (iFirst < 0) continue;
 285     int iSecond = singlets[iSub].iParton[i + 1];
 286     if (iSecond < 0) iSecond = singlets[iSub].iParton[i + 2];
 287     if (iSecond != iFirst + 1) { inOrder = false; break;}
 288   }
 289   if (inOrder) return;
 290
 291   // Copy down system. Update current partons.
 292   for (int i = 0; i < singlets[iSub].size(); ++i) {
 293     int iOld = singlets[iSub].iParton[i];
 294     if (iOld < 0) continue;
 295     int iNew = event.copy(iOld, 71);
 296     singlets[iSub].iParton[i] = iNew;
 297   }
 298
 299   // Done.
 300 }
 301
 302 //*********
 303
 304 // List all currently identified singlets.
 305
 306 void ColConfig::list(ostream& os) {
 307
 308   // Header. Loop over all individual singlets.
 309   os << "\n --------  Colour Singlet Systems Listing -------------------\n";
 310   for (int iSub = 0; iSub < int(singlets.size()); ++iSub) {
 311
 312     // List all partons belonging to each singlet.
 313     os << " singlet " << iSub << " contains " ;
 314     for (int i = 0; i < singlets[iSub].size(); ++i)
 315       os << singlets[iSub].iParton[i] << " ";
 316     os << "\n";
 317
 318   // Done.
 319   }
 320 }
 321
 322 //**************************************************************************
 323
 324 // The StringRegion class.
 325
 326 // Currently a number of simplifications, in particular ??
 327 // 1) No popcorn baryon production.
 328 // 2) Simplified treatment of pT in stepping and joining.
 329
 330 //*********
 331
 332 // Constants: could be changed here if desired, but normally should not.
 333 // These are of technical nature, as described for each.
 334
 335 // If a string region is smaller thsan this it is assumed empty.
 336 const double StringRegion::MJOIN = 0.1;
 337
 338 // Avoid division by zero.
 339 const double StringRegion::TINY  = 1e-20;
 340
 341 //*********
 342
 343 // Set up four-vectors for longitudinal and transverse directions.
 344
 345 void StringRegion::setUp(Vec4 p1, Vec4 p2, bool isMassless) {
 346
 347   // Simple case: the two incoming four-vectors guaranteed massless.
 348   if (isMassless) {
 349
 350     // Calculate w2, minimum value. Lightcone directions = input.
 351     w2 = 2. * (p1 * p2);
 352     if (w2 < MJOIN*MJOIN) {isSetUp = true; isEmpty = true; return;}
 353     pPos = p1;
 354     pNeg = p2;
 355
 356   // Else allow possibility of masses for incoming partons (also gluons!).
 357   } else {
 358
 359     // Generic four-momentum combinations.
 360     double m1Sq = p1 * p1;
 361     double m2Sq = p2 * p2;
 362     double p1p2 = p1 * p2;
 363     w2 = m1Sq + 2. * p1p2 + m2Sq;
 364     double rootSq = pow2(p1p2) - m1Sq * m2Sq;
 365
 366     // If crazy kinematics (should not happen!) modify energies.
 367     if (w2 <= 0. || rootSq <= 0.) {
 368       if (m1Sq < 0.) m1Sq = 0.;
 369       p1.e( sqrt(m1Sq + p1.pAbs2()) );
 370       if (m2Sq < 0.) m2Sq = 0.;
 371       p2.e( sqrt(m2Sq + p2.pAbs2()) );
 372       p1p2 = p1 * p2;
 373       w2 = m1Sq + 2. * p1p2 + m2Sq;
 374       rootSq = pow2(p1p2) - m1Sq * m2Sq;
 375     }
 376
 377     // If still small invariant mass then empty region (e.g. in gg system).
 378     if (w2 < MJOIN*MJOIN) {isSetUp = true; isEmpty = true; return;}
 379
 380     // Find two lightconelike longitudinal four-vector directions.
 381     double root = sqrt( max(TINY, rootSq) );
 382     double k1 = 0.5 * ( (m2Sq + p1p2) / root - 1.);
 383     double k2 = 0.5 * ( (m1Sq + p1p2) / root - 1.);
 384     pPos = (1. + k1) * p1 - k2 * p2;
 385     pNeg = (1. + k2) * p2 - k1 * p1;
 386   }
 387
 388   // Find two spacelike transverse four-vector directions.
 389   // Begin by picking two sensible trial directions.
 390   Vec4 eDiff = pPos / pPos.e() - pNeg / pNeg.e();
 391   double eDx = pow2( eDiff.px() );
 392   double eDy = pow2( eDiff.py() );
 393   double eDz = pow2( eDiff.pz() );
 394   if (eDx < min(eDy, eDz)) {
 395     eX = Vec4( 1., 0., 0., 0.);
 396     eY = (eDy < eDz) ? Vec4( 0., 1., 0., 0.) : Vec4( 0., 0., 1., 0.);
 397   } else if (eDy < eDz) {
 398     eX = Vec4( 0., 1., 0., 0.);
 399     eY = (eDx < eDz) ? Vec4( 1., 0., 0., 0.) : Vec4( 0., 0., 1., 0.);
 400   } else {
 401     eX = Vec4( 0., 0., 1., 0.);
 402     eY = (eDx < eDy) ? Vec4( 1., 0., 0., 0.) : Vec4( 0., 1., 0., 0.);
 403   }
 404
 405   // Then construct orthogonal linear combinations.
 406   double pPosNeg = pPos * pNeg;
 407   double kXPos = eX * pPos / pPosNeg;
 408   double kXNeg = eX * pNeg / pPosNeg;
 409   double kXX = 1. / sqrt( 1. + 2. * kXPos * kXNeg * pPosNeg );
 410   double kYPos = eY * pPos / pPosNeg;
 411   double kYNeg = eY * pNeg / pPosNeg;
 412   double kYX = kXX * (kXPos * kYNeg + kXNeg * kYPos) * pPosNeg;
 413   double kYY = 1. / sqrt(1. + 2. * kYPos * kYNeg * pPosNeg - pow2(kYX));
 414   eX = kXX * (eX - kXNeg * pPos - kXPos * pNeg);
 415   eY = kYY * (eY - kYNeg * pPos - kYPos * pNeg - kYX * eX);
 416
 417   // Done.
 418   isSetUp = true;
 419   isEmpty = false;
 420
 421 }
 422
 423 //*********
 424
 425 // Project a four-momentum onto (x+, x-, px, py).
 426
 427 void StringRegion::project(Vec4 pIn) {
 428
 429   // Perform projections by four-vector multiplication.
 430   xPosProj = 2. * (pIn * pNeg) / w2;
 431   xNegProj = 2. * (pIn * pPos) / w2;
 432   pxProj = - (pIn * eX);
 433   pyProj = - (pIn * eY);
 434
 435 }
 436
 437 //**************************************************************************
 438
 439 // The StringSystem class.
 440
 441 //*********
 442
 443 // Set up system from parton list.
 444
 445 void StringSystem::setUp(vector<int>& iSys, Event& event) {
 446
 447   // Figure out how big the system is. (Closed gluon loops?)
 448   sizePartons = iSys.size();
 449   sizeStrings = sizePartons - 1;
 450   sizeRegions = (sizeStrings * (sizeStrings + 1)) / 2;
 451   indxReg = 2 * sizeStrings + 1;
 452   iMax = sizeStrings - 1;
 453
 454   // Reserve space for the required number of regions.
 455   system.clear();
 456   system.resize(sizeRegions);
 457
 458   // Set up the lowest-lying regions.
 459   for (int i = 0; i < sizeStrings; ++i) {
 460     Vec4 p1 = event[ iSys[i] ].p();
 461     if ( event[ iSys[i] ].isGluon() ) p1 *= 0.5;
 462     Vec4 p2 = event[ iSys[i+1] ].p();
 463     if ( event[ iSys[i+1] ].isGluon() ) p2 *= 0.5;
 464     system[ iReg(i, iMax - i) ].setUp( p1, p2, false);
 465   }
 466
 467 }
 468
 469 //**************************************************************************
 470
 471 } // end namespace Pythia8