PYTHIA8/pythia8170/src/FragmentationSystems.cxx

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