PYTHIA8/pythia8170/src/StringFragmentation.cxx

   1 // StringFragmentation.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 StringEnd and
   7 // StringFragmentation classes.
   8
   9 #include "StringFragmentation.h"
  10
  11 namespace Pythia8 {
  12
  13 //==========================================================================
  14
  15 // The StringEnd class.
  16
  17 //--------------------------------------------------------------------------
  18
  19 // Constants: could be changed here if desired, but normally should not.
  20
  21 // Avoid unphysical solutions to equation system.
  22 const double StringEnd::TINY = 1e-6;
  23
  24 // Assume two (eX, eY) regions are related if pT2 differs by less.
  25 const double StringEnd::PT2SAME = 0.01;
  26
  27 //--------------------------------------------------------------------------
  28
  29 // Set up initial endpoint values from input.
  30
  31 void StringEnd::setUp(bool fromPosIn, int iEndIn, int idOldIn, int iMaxIn,
  32   double pxIn, double pyIn, double GammaIn, double xPosIn, double xNegIn) {
  33
  34   // Simple transcription from input.
  35   fromPos  = fromPosIn;
  36   iEnd     = iEndIn;
  37   iMax     = iMaxIn;
  38   flavOld  = FlavContainer(idOldIn);
  39   pxOld    = pxIn;
  40   pyOld    = pyIn;
  41   GammaOld = GammaIn;
  42   iPosOld  = (fromPos) ? 0 : iMax;
  43   iNegOld  = (fromPos) ? iMax : 0;
  44   xPosOld  = xPosIn;
  45   xNegOld  = xNegIn;
  46
  47 }
  48
  49 //--------------------------------------------------------------------------
  50
  51 // Fragment off one hadron from the string system, in flavour and pT.
  52
  53 void StringEnd::newHadron() {
  54
  55   // Pick new flavour and form a new hadron.
  56   do {
  57     flavNew = flavSelPtr->pick( flavOld);
  58     idHad   = flavSelPtr->combine( flavOld, flavNew);
  59   } while (idHad == 0);
  60
  61   // Pick its transverse momentum.
  62   pair<double, double> pxy = pTSelPtr->pxy();
  63   pxNew = pxy.first;
  64   pyNew = pxy.second;
  65   pxHad = pxOld + pxNew;
  66   pyHad = pyOld + pyNew;
  67
  68   // Pick its mass and thereby define its transverse mass.
  69   mHad   = particleDataPtr->mass(idHad);
  70   mT2Had = pow2(mHad) + pow2(pxHad) + pow2(pyHad);
  71
  72 }
  73
  74 //--------------------------------------------------------------------------
  75
  76 // Fragment off one hadron from the string system, in momentum space,
  77 // by taking steps from positive end.
  78
  79 Vec4 StringEnd::kinematicsHadron( StringSystem& system) {
  80
  81   // Pick fragmentation step z and calculate new Gamma.
  82   zHad = zSelPtr->zFrag( flavOld.id, flavNew.id, mT2Had);
  83   GammaNew = (1. - zHad) * (GammaOld + mT2Had / zHad);
  84
  85   // Set up references that are direction-neutral;
  86   // ...Dir for direction of iteration and ...Inv for its inverse.
  87   int&    iDirOld = (fromPos) ? iPosOld : iNegOld;
  88   int&    iInvOld = (fromPos) ? iNegOld : iPosOld;
  89   int&    iDirNew = (fromPos) ? iPosNew : iNegNew;
  90   int&    iInvNew = (fromPos) ? iNegNew : iPosNew;
  91   double& xDirOld = (fromPos) ? xPosOld : xNegOld;
  92   double& xInvOld = (fromPos) ? xNegOld : xPosOld;
  93   double& xDirNew = (fromPos) ? xPosNew : xNegNew;
  94   double& xInvNew = (fromPos) ? xNegNew : xPosNew;
  95   double& xDirHad = (fromPos) ? xPosHad : xNegHad;
  96   double& xInvHad = (fromPos) ? xNegHad : xPosHad;
  97
  98   // Start search for new breakup in the old region.
  99   iDirNew = iDirOld;
 100   iInvNew = iInvOld;
 101   Vec4 pTNew;
 102
 103   // Each step corresponds to trying a new string region.
 104   for (int iStep = 0; ; ++iStep) {
 105
 106     // Referance to current string region.
 107     StringRegion& region = system.region( iPosNew, iNegNew);
 108
 109     // Now begin special section for rapid processing of low region.
 110     if (iStep == 0 && iPosOld + iNegOld == iMax) {
 111
 112       // A first step within a low region is easy.
 113       if (mT2Had < zHad * xDirOld * (1. - xInvOld) * region.w2) {
 114
 115         // Translate into x coordinates.
 116         xDirHad = zHad * xDirOld;
 117         xInvHad = mT2Had / (xDirHad * region.w2);
 118         xDirNew = xDirOld - xDirHad;
 119         xInvNew = xInvOld + xInvHad;
 120
 121         // Find and return four-momentum of the produced particle.
 122         return region.pHad( xPosHad, xNegHad, pxHad, pyHad);
 123
 124       // A first step out of a low region also OK, if there are more regions.
 125       // Negative energy signals failure, i.e. in last region.
 126       } else {
 127         --iInvNew;
 128         if (iInvNew < 0) return Vec4(0., 0., 0., -1.);
 129
 130         // Momentum taken by stepping out of region. Continue to next region.
 131         xInvHad = 1. - xInvOld;
 132         xDirHad = 0.;
 133         pSoFar  = region.pHad( xPosHad, xNegHad, pxOld, pyOld);
 134         continue;
 135       }
 136
 137     // Else, for first step, take into account starting pT.
 138     } else if (iStep == 0) {
 139       pSoFar = region.pHad( 0., 0., pxOld, pyOld);
 140       pTNew  = region.pHad( 0., 0., pxNew, pyNew);
 141     }
 142
 143     // Now begin normal treatment of nontrivial regions.
 144     // Set up four-vectors in a region not visited before.
 145     if (!region.isSetUp) region.setUp(
 146       system.regionLowPos(iPosNew).pPos,
 147       system.regionLowNeg(iNegNew).pNeg, true);
 148
 149     // If new region is vanishingly small, continue immediately to next.
 150     // Negative energy signals failure to do this, i.e. moved too low.
 151     if (region.isEmpty) {
 152       xDirHad = (iDirNew == iDirOld) ? xDirOld : 1.;
 153       xInvHad = 0.;
 154       pSoFar += region.pHad( xPosHad, xNegHad, 0., 0.);
 155       ++iDirNew;
 156       if (iDirNew + iInvNew > iMax) return Vec4(0., 0., 0., -1.);
 157       continue;
 158     }
 159
 160     // Reexpress pTNew w.r.t. base vectors in new region, if possible.
 161     // Recall minus sign from normalization e_x^2 = e_y^2 = -1.
 162     double pxNewTemp = -pTNew * region.eX;
 163     double pyNewTemp = -pTNew * region.eY;
 164     if (abs( pxNewTemp * pxNewTemp + pyNewTemp * pyNewTemp
 165       - pxNew * pxNew - pyNew * pyNew) < PT2SAME) {
 166       pxNew = pxNewTemp;
 167       pyNew = pyNewTemp;
 168     }
 169
 170     // Four-momentum taken so far, including new pT.
 171     Vec4 pTemp = pSoFar + region.pHad( 0., 0., pxNew, pyNew);
 172
 173     // Derive coefficients for m2 expression.
 174     // cM2 * x+ + cM3 * x- + cM4 * x+ * x- = m^2 - cM1;
 175     double cM1 = pTemp.m2Calc();
 176     double cM2 = 2. * (pTemp * region.pPos);
 177     double cM3 = 2. * (pTemp * region.pNeg);
 178     double cM4 = region.w2;
 179     if (!fromPos) swap( cM2, cM3);
 180
 181     // Derive coefficients for Gamma expression.
 182     // cGam2 * x+ + cGam3 * x- + cGam4 * x+ * x- = Gamma_new - cGam1;
 183     double cGam1 = 0.;
 184     double cGam2 = 0.;
 185     double cGam3 = 0.;
 186     double cGam4 = 0.;
 187     for (int iInv = iInvNew; iInv <= iMax - iDirNew; ++iInv) {
 188       double xInv = 1.;
 189       if (iInv == iInvNew) xInv = (iInvNew == iInvOld) ? xInvOld : 0.;
 190       for (int iDir = iDirNew; iDir <= iMax - iInv; ++iDir) {
 191         double xDir = (iDir == iDirOld) ? xDirOld : 1.;
 192         int iPos = (fromPos) ? iDir : iInv;
 193         int iNeg = (fromPos) ? iInv : iDir;
 194         StringRegion& regionGam =  system.region( iPos, iNeg);
 195         if (!regionGam.isSetUp) regionGam.setUp(
 196           system.regionLowPos(iPos).pPos,
 197           system.regionLowNeg(iNeg).pNeg, true);
 198         double w2 = regionGam.w2;
 199         cGam1 += xDir * xInv * w2;
 200         if (iDir == iDirNew) cGam2 -= xInv * w2;
 201         if (iInv == iInvNew) cGam3 += xDir * w2;
 202         if (iDir == iDirNew && iInv == iInvNew) cGam4 -= w2;
 203       }
 204     }
 205
 206     // Solve (m2, Gamma) equation system => r2 * x-^2 + r1 * x- + r0 = 0.
 207     double cM0   = pow2(mHad) - cM1;
 208     double cGam0 = GammaNew - cGam1;
 209     double r2    = cM3 * cGam4 - cM4 * cGam3;
 210     double r1    = cM4 * cGam0 - cM0 * cGam4 + cM3 * cGam2 - cM2 * cGam3;
 211     double r0    = cM2 * cGam0 - cM0 * cGam2;
 212     double root  = sqrtpos( r1*r1 - 4. * r2 * r0 );
 213     if (abs(r2) < TINY || root < TINY) return Vec4(0., 0., 0., -1.);
 214     xInvHad      = 0.5 * (root / abs(r2) - r1 / r2);
 215     xDirHad      = (cM0 - cM3 * xInvHad) / (cM2 + cM4 * xInvHad);
 216
 217     // Define position of new trial vertex.
 218     xDirNew = (iDirNew == iDirOld) ? xDirOld - xDirHad : 1. - xDirHad;
 219     xInvNew = (iInvNew == iInvOld) ? xInvOld + xInvHad : xInvHad;
 220
 221     // Step up to new region if new x- > 1.
 222     if (xInvNew > 1.) {
 223       xInvHad = (iInvNew == iInvOld) ? 1. - xInvOld : 1.;
 224       xDirHad = 0.;
 225       pSoFar += region.pHad( xPosHad, xNegHad, 0., 0.);
 226       --iInvNew;
 227       if (iInvNew < 0) return Vec4(0., 0., 0., -1.);
 228       continue;
 229
 230     // Step down to new region if new x+ < 0.
 231     } else if (xDirNew < 0.) {
 232       xDirHad = (iDirNew == iDirOld) ? xDirOld : 1.;
 233       xInvHad = 0.;
 234       pSoFar += region.pHad( xPosHad, xNegHad, 0., 0.);
 235       ++iDirNew;
 236       if (iDirNew + iInvNew > iMax) return Vec4(0., 0., 0., -1.);
 237       continue;
 238     }
 239
 240     // Else we have found the correct region, and can return four-momentum.
 241     return pSoFar + region.pHad( xPosHad, xNegHad, pxNew, pyNew);
 242
 243   // End of "infinite" loop of stepping to new region.
 244   }
 245
 246 }
 247
 248 //--------------------------------------------------------------------------
 249
 250 // Update string end information after a hadron has been removed.
 251
 252 void StringEnd::update() {
 253
 254   flavOld.anti(flavNew);
 255   iPosOld  = iPosNew;
 256   iNegOld  = iNegNew;
 257   pxOld    = -pxNew;
 258   pyOld    = -pyNew;
 259   GammaOld = GammaNew;
 260   xPosOld  = xPosNew;
 261   xNegOld  = xNegNew;
 262
 263 }
 264
 265 //==========================================================================
 266
 267 // The StringFragmentation class.
 268
 269 //--------------------------------------------------------------------------
 270
 271 // Constants: could be changed here if desired, but normally should not.
 272 // These are of technical nature, as described for each.
 273
 274 // Maximum number of tries to (flavour-, energy) join the two string ends.
 275 const int    StringFragmentation::NTRYFLAV      = 10;
 276 const int    StringFragmentation::NTRYJOIN      = 30;
 277
 278 // The last few times gradually increase the stop mass to make it easier.
 279 const int    StringFragmentation::NSTOPMASS     = 15;
 280 const double StringFragmentation::FACSTOPMASS   = 1.05;
 281
 282 // For closed string, pick a Gamma by taking a step with fictitious mass.
 283 const double StringFragmentation::CLOSEDM2MAX   = 25.;
 284 const double StringFragmentation::CLOSEDM2FRAC  = 0.1;
 285
 286 // Do not allow too large argument to exp function.
 287 const double StringFragmentation::EXPMAX        = 50.;
 288
 289 // Matching criterion that p+ and p- not the same (can happen in gg loop).
 290 const double StringFragmentation::MATCHPOSNEG   = 1e-6;
 291
 292 // For pull on junction, do not trace too far down each leg.
 293 const double StringFragmentation::EJNWEIGHTMAX  = 10.;
 294
 295 // Iterate junction rest frame boost until convergence or too many tries.
 296 const double StringFragmentation::CONVJNREST   = 1e-5;
 297 const int    StringFragmentation::NTRYJNREST   = 20;
 298
 299 // Fail and try again when two legs combined to diquark (3 loops).
 300 const int    StringFragmentation::NTRYJNMATCH   = 20;
 301 const double StringFragmentation::EEXTRAJNMATCH = 0.5;
 302 const double StringFragmentation::MDIQUARKMIN   = -2.0;
 303
 304 // Consider junction-leg parton as massless if m2 tiny.
 305 const double StringFragmentation::M2MAXJRF      = 1e-4;
 306
 307 // Iterate junction rest frame equation until convergence or too many tries.
 308 const double StringFragmentation::CONVJRFEQ     = 1e-12;
 309 const int    StringFragmentation::NTRYJRFEQ     = 40;
 310
 311 //--------------------------------------------------------------------------
 312
 313 // Initialize and save pointers.
 314
 315 void StringFragmentation::init(Info* infoPtrIn, Settings& settings,
 316   ParticleData* particleDataPtrIn, Rndm* rndmPtrIn, StringFlav* flavSelPtrIn,
 317   StringPT* pTSelPtrIn, StringZ* zSelPtrIn) {
 318
 319   // Save pointers.
 320   infoPtr         = infoPtrIn;
 321   particleDataPtr = particleDataPtrIn;
 322   rndmPtr         = rndmPtrIn;
 323   flavSelPtr      = flavSelPtrIn;
 324   pTSelPtr        = pTSelPtrIn;
 325   zSelPtr         = zSelPtrIn;
 326
 327   // Initialize the StringFragmentation class.
 328   stopMass        = zSelPtr->stopMass();
 329   stopNewFlav     = zSelPtr->stopNewFlav();
 330   stopSmear       = zSelPtr->stopSmear();
 331   eNormJunction   = settings.parm("StringFragmentation:eNormJunction");
 332   eBothLeftJunction
 333      = settings.parm("StringFragmentation:eBothLeftJunction");
 334   eMaxLeftJunction
 335     = settings.parm("StringFragmentation:eMaxLeftJunction");
 336   eMinLeftJunction
 337     = settings.parm("StringFragmentation:eMinLeftJunction");
 338
 339   // Joining of nearby partons along the string.
 340   mJoin           = settings.parm("FragmentationSystems:mJoin");
 341
 342   // Initialize the b parameter of the z spectrum, used when joining jets.
 343   bLund           = zSelPtr->bAreaLund();
 344
 345   // Initialize the hadrons instance of an event record.
 346   hadrons.init( "(string fragmentation)", particleDataPtr);
 347
 348   // Send on pointers to the two StringEnd instances.
 349   posEnd.init( particleDataPtr, flavSelPtr, pTSelPtr, zSelPtr);
 350   negEnd.init( particleDataPtr, flavSelPtr, pTSelPtr, zSelPtr);
 351
 352 }
 353
 354 //--------------------------------------------------------------------------
 355
 356 // Perform the fragmentation.
 357
 358 bool StringFragmentation::fragment( int iSub, ColConfig& colConfig,
 359   Event& event) {
 360
 361   // Find partons and their total four-momentum.
 362   iParton            = colConfig[iSub].iParton;
 363   iPos               = iParton[0];
 364   if (iPos < 0) iPos = iParton[1];
 365   int idPos          = event[iPos].id();
 366   iNeg               = iParton.back();
 367   int idNeg          = event[iNeg].id();
 368   pSum               = colConfig[iSub].pSum;
 369
 370   // Reset the local event record.
 371   hadrons.clear();
 372
 373   // For closed gluon string: pick first breakup region.
 374   isClosed = colConfig[iSub].isClosed;
 375   if (isClosed) iParton = findFirstRegion(iParton, event);
 376
 377   // For junction topology: fragment off two of the string legs.
 378   // Then iParton overwritten to remaining leg + leftover diquark.
 379   pJunctionHadrons = 0.;
 380   hasJunction = colConfig[iSub].hasJunction;
 381   if (hasJunction && !fragmentToJunction(event)) return false;
 382   int junctionHadrons = hadrons.size();
 383   if (hasJunction) {
 384     idPos  = event[ iParton[0] ].id();
 385     idNeg  = event.back().id();
 386     pSum  -= pJunctionHadrons;
 387   }
 388
 389   // Set up kinematics of string evolution ( = motion).
 390   system.setUp(iParton, event);
 391   stopMassNow = stopMass;
 392   int nExtraJoin = 0;
 393
 394   // Fallback loop, when joining in the middle fails.  Bailout if stuck.
 395   for ( int iTry = 0; ; ++iTry) {
 396     if (iTry > NTRYJOIN) {
 397       infoPtr->errorMsg("Error in StringFragmentation::fragment: "
 398         "stuck in joining");
 399       if (hasJunction) ++nExtraJoin;
 400       if (nExtraJoin > 0) event.popBack(nExtraJoin);
 401       return false;
 402     }
 403
 404     // After several failed tries join some (extra) nearby partons.
 405     if (iTry == NTRYJOIN / 3) nExtraJoin = extraJoin( 2., event);
 406     if (iTry == 2 * NTRYJOIN / 3) nExtraJoin += extraJoin( 4., event);
 407
 408     // After several failed tries gradually allow larger stop mass.
 409     if (iTry > NTRYJOIN - NSTOPMASS) stopMassNow *= FACSTOPMASS;
 410
 411     // Set up flavours of two string ends, and reset other info.
 412     setStartEnds(idPos, idNeg, system);
 413     pRem = pSum;
 414
 415     // Begin fragmentation loop, interleaved from the two ends.
 416     bool fromPos;
 417     for ( ; ; ) {
 418
 419       // Take a step either from the positive or the negative end.
 420       fromPos           = (rndmPtr->flat() < 0.5);
 421       StringEnd& nowEnd = (fromPos) ? posEnd : negEnd;
 422
 423       // Construct trial hadron and check that energy remains.
 424       nowEnd.newHadron();
 425       if ( energyUsedUp(fromPos) ) break;
 426
 427       // Construct kinematics of the new hadron and store it.
 428       Vec4 pHad = nowEnd.kinematicsHadron(system);
 429       int statusHad = (fromPos) ? 83 : 84;
 430       hadrons.append( nowEnd.idHad, statusHad, iPos, iNeg,
 431         0, 0, 0, 0, pHad, nowEnd.mHad);
 432       if (pHad.e() < 0.) break;
 433
 434       // Update string end and remaining momentum.
 435       nowEnd.update();
 436       pRem -= pHad;
 437
 438     // End of fragmentation loop.
 439     }
 440
 441     // When done, join in the middle. If this works, then really done.
 442     if ( finalTwo(fromPos) ) break;
 443
 444     // Else remove produced particles (except from first two junction legs)
 445     // and start all over.
 446     int newHadrons = hadrons.size() - junctionHadrons;
 447     hadrons.popBack(newHadrons);
 448   }
 449
 450   // Junctions & extra joins: remove fictitious end, restore original partons.
 451   if (hasJunction) ++nExtraJoin;
 452   if (nExtraJoin > 0) {
 453     event.popBack(nExtraJoin);
 454     iParton = colConfig[iSub].iParton;
 455   }
 456
 457   // Store the hadrons in the normal event record, ordered from one end.
 458   store(event);
 459
 460   // Done.
 461   return true;
 462
 463 }
 464
 465 //--------------------------------------------------------------------------
 466
 467 // Find region where to put first string break for closed gluon loop.
 468
 469 vector<int> StringFragmentation::findFirstRegion(vector<int>& iPartonIn,
 470   Event& event) {
 471
 472   // Evaluate mass-squared for all adjacent gluon pairs.
 473   vector<double> m2Pair;
 474   double m2Sum = 0.;
 475   int size = iPartonIn.size();
 476   for (int i = 0; i < size; ++i) {
 477     double m2Now = 0.5 * event[ iPartonIn[i] ].p()
 478       * event[ iPartonIn[(i + 1)%size] ].p();
 479     m2Pair.push_back(m2Now);
 480     m2Sum += m2Now;
 481   }
 482
 483   // Pick breakup region with probability proportional to mass-squared.
 484   double m2Reg = m2Sum * rndmPtr->flat();
 485   int iReg = -1;
 486   do m2Reg -= m2Pair[++iReg];
 487   while (m2Reg > 0. && iReg < size - 1);
 488
 489   // Create reordered parton list, with breakup string region duplicated.
 490   vector<int> iPartonOut;
 491   for (int i = 0; i < size + 2; ++i)
 492     iPartonOut.push_back( iPartonIn[(i + iReg)%size] );
 493
 494   // Done.
 495   return iPartonOut;
 496
 497 }
 498
 499 //--------------------------------------------------------------------------
 500
 501 // Set flavours and momentum position for initial string endpoints.
 502
 503 void StringFragmentation::setStartEnds( int idPos, int idNeg,
 504   StringSystem systemNow) {
 505
 506   // Variables characterizing string endpoints: defaults for open string.
 507   double px          = 0.;
 508   double py          = 0.;
 509   double Gamma       = 0.;
 510   double xPosFromPos = 1.;
 511   double xNegFromPos = 0.;
 512   double xPosFromNeg = 0.;
 513   double xNegFromNeg = 1.;
 514
 515   // For closed gluon loop need to pick an initial flavour.
 516   if (isClosed) {
 517     do {
 518       int idTry = flavSelPtr->pickLightQ();
 519       FlavContainer flavTry(idTry, 1);
 520       flavTry = flavSelPtr->pick( flavTry);
 521       flavTry = flavSelPtr->pick( flavTry);
 522       idPos   = flavTry.id;
 523       idNeg   = -idPos;
 524     } while (idPos == 0);
 525
 526     // Also need pT and breakup vertex position in region.
 527    pair<double, double> pxy = pTSelPtr->pxy();
 528    px = pxy.first;
 529    py = pxy.second;
 530     double m2Region = systemNow.regionLowPos(0).w2;
 531     double m2Temp   = min( CLOSEDM2MAX, CLOSEDM2FRAC * m2Region);
 532     do {
 533       double zTemp = zSelPtr->zFrag( idPos, idNeg, m2Temp);
 534       xPosFromPos  = 1. - zTemp;
 535       xNegFromPos  = m2Temp / (zTemp * m2Region);
 536     } while (xNegFromPos > 1.);
 537     Gamma = xPosFromPos * xNegFromPos * m2Region;
 538     xPosFromNeg = xPosFromPos;
 539     xNegFromNeg = xNegFromPos;
 540   }
 541
 542   // Initialize two string endpoints.
 543   posEnd.setUp(  true, iPos, idPos, systemNow.iMax,  px,  py,
 544     Gamma, xPosFromPos, xNegFromPos);
 545   negEnd.setUp( false, iNeg, idNeg, systemNow.iMax, -px, -py,
 546     Gamma, xPosFromNeg, xNegFromNeg);
 547
 548   // For closed gluon loop can allow popcorn on one side but not both.
 549   if (isClosed) {
 550     flavSelPtr->assignPopQ(posEnd.flavOld);
 551     flavSelPtr->assignPopQ(negEnd.flavOld);
 552     if (rndmPtr->flat() < 0.5) posEnd.flavOld.nPop = 0;
 553     else                    negEnd.flavOld.nPop = 0;
 554     posEnd.flavOld.rank = 1;
 555     negEnd.flavOld.rank = 1;
 556   }
 557
 558   // Done.
 559
 560 }
 561
 562 //--------------------------------------------------------------------------
 563
 564 // Check remaining energy-momentum whether it is OK to continue.
 565
 566 bool StringFragmentation::energyUsedUp(bool fromPos) {
 567
 568   // If remaining negative energy then abort right away.
 569   if (pRem.e() < 0.) return true;
 570
 571   // Calculate W2_minimum and done if remaining W2 is below it.
 572   double wMin = stopMassNow
 573     + particleDataPtr->constituentMass(posEnd.flavOld.id)
 574     + particleDataPtr->constituentMass(negEnd.flavOld.id);
 575   if (fromPos) wMin += stopNewFlav
 576     * particleDataPtr->constituentMass(posEnd.flavNew.id);
 577   else         wMin += stopNewFlav
 578     * particleDataPtr->constituentMass(negEnd.flavNew.id);
 579   wMin *= 1. + (2. * rndmPtr->flat() - 1.) * stopSmear;
 580   w2Rem = pRem.m2Calc();
 581   if (w2Rem < pow2(wMin)) return true;
 582
 583   // Else still enough energy left to continue iteration.
 584   return false;
 585
 586 }
 587
 588 //--------------------------------------------------------------------------
 589
 590 // Produce the final two partons to complete the system.
 591
 592 bool StringFragmentation::finalTwo(bool fromPos) {
 593
 594   // Check whether we went too far in p+-.
 595   if (pRem.e() < 0.  || w2Rem < 0. || (hadrons.size() > 0
 596     && hadrons.back().e() < 0.) ) return false;
 597   if ( posEnd.iPosOld > negEnd.iPosOld || negEnd.iNegOld > posEnd.iNegOld)
 598     return false;
 599   if ( posEnd.iPosOld == negEnd.iPosOld && posEnd.xPosOld < negEnd.xPosOld)
 600     return false;
 601   if ( posEnd.iNegOld == negEnd.iNegOld && posEnd.xNegOld > negEnd.xNegOld)
 602     return false;
 603
 604   // Construct the final hadron from the leftover flavours.
 605   // Impossible to join two diquarks. Also break if stuck for other reason.
 606   FlavContainer flav1 = (fromPos) ? posEnd.flavNew.anti() : posEnd.flavOld;
 607   FlavContainer flav2 = (fromPos) ? negEnd.flavOld : negEnd.flavNew.anti();
 608   if (flav1.isDiquark() && flav2.isDiquark()) return false;
 609   int idHad;
 610   for (int iTry = 0; iTry < NTRYFLAV; ++iTry) {
 611     idHad = flavSelPtr->combine( flav1, flav2);
 612     if (idHad != 0) break;
 613   }
 614   if (idHad == 0) return false;
 615
 616   // Store the final particle and its new pT, and construct its mass.
 617   if (fromPos) {
 618     negEnd.idHad = idHad;
 619     negEnd.pxNew = -posEnd.pxNew;
 620     negEnd.pyNew = -posEnd.pyNew;
 621     negEnd.mHad  = particleDataPtr->mass(idHad);
 622   } else {
 623     posEnd.idHad = idHad;
 624     posEnd.pxNew = -negEnd.pxNew;
 625     posEnd.pyNew = -negEnd.pyNew;
 626     posEnd.mHad  = particleDataPtr->mass(idHad);
 627   }
 628
 629   // String region in which to do the joining.
 630   StringRegion region = finalRegion();
 631   if (region.isEmpty) return false;
 632
 633   // Project remaining momentum along longitudinal and transverse directions.
 634   region.project( pRem);
 635   double pxRem   = region.px() - posEnd.pxOld - negEnd.pxOld;
 636   double pyRem   = region.py() - posEnd.pyOld - negEnd.pyOld;
 637   double xPosRem = region.xPos();
 638   double xNegRem = region.xNeg();
 639
 640   // Share extra pT kick evenly between final two hadrons.
 641   posEnd.pxOld += 0.5 * pxRem;
 642   posEnd.pyOld += 0.5 * pyRem;
 643   negEnd.pxOld += 0.5 * pxRem;
 644   negEnd.pyOld += 0.5 * pyRem;
 645
 646   // Construct new pT and mT of the final two particles.
 647   posEnd.pxHad  = posEnd.pxOld + posEnd.pxNew;
 648   posEnd.pyHad  = posEnd.pyOld + posEnd.pyNew;
 649   posEnd.mT2Had = pow2(posEnd.mHad) + pow2(posEnd.pxHad)
 650     + pow2(posEnd.pyHad);
 651   negEnd.pxHad  = negEnd.pxOld + negEnd.pxNew;
 652   negEnd.pyHad  = negEnd.pyOld + negEnd.pyNew;
 653   negEnd.mT2Had = pow2(negEnd.mHad) + pow2(negEnd.pxHad)
 654     + pow2(negEnd.pyHad);
 655
 656   // Construct remaining system transverse mass.
 657   double wT2Rem = w2Rem + pow2( posEnd.pxHad + negEnd.pxHad)
 658     + pow2( posEnd.pyHad + negEnd.pyHad);
 659
 660   // Check that kinematics possible.
 661   if ( sqrt(wT2Rem) < sqrt(posEnd.mT2Had) + sqrt(negEnd.mT2Had) )
 662     return false;
 663   double lambda2 = pow2( wT2Rem - posEnd.mT2Had - negEnd.mT2Had)
 664     - 4. * posEnd.mT2Had * negEnd.mT2Had;
 665   if (lambda2 <= 0.) return false;
 666
 667   // Construct kinematics, as viewed in the transverse rest frame.
 668   double lambda = sqrt(lambda2);
 669   double probReverse = 1. / (1. + exp( min( EXPMAX, bLund * lambda) ) );
 670   double xpzPos = 0.5 * lambda/ wT2Rem;
 671   if (probReverse > rndmPtr->flat()) xpzPos = -xpzPos;
 672   double xmDiff = (posEnd.mT2Had - negEnd.mT2Had) / wT2Rem;
 673   double xePos  = 0.5 * (1. + xmDiff);
 674   double xeNeg  = 0.5 * (1. - xmDiff );
 675
 676   // Translate this into kinematics in the string frame.
 677   Vec4 pHadPos = region.pHad( (xePos + xpzPos) *  xPosRem,
 678     (xePos - xpzPos) *  xNegRem, posEnd.pxHad, posEnd.pyHad);
 679   Vec4 pHadNeg = region.pHad( (xeNeg - xpzPos) *  xPosRem,
 680     (xeNeg + xpzPos) *  xNegRem, negEnd.pxHad, negEnd.pyHad);
 681
 682   // Add produced particles to the event record.
 683   hadrons.append( posEnd.idHad, 83, posEnd.iEnd, negEnd.iEnd,
 684     0, 0, 0, 0, pHadPos, posEnd.mHad);
 685   hadrons.append( negEnd.idHad, 84, posEnd.iEnd, negEnd.iEnd,
 686     0, 0, 0, 0, pHadNeg, negEnd.mHad);
 687
 688   // It worked.
 689   return true;
 690
 691 }
 692
 693 //--------------------------------------------------------------------------
 694
 695 //  Construct a special joining region for the final two hadrons.
 696
 697 StringRegion StringFragmentation::finalRegion() {
 698
 699   // Simple case when both string ends are in the same region.
 700   if (posEnd.iPosOld == negEnd.iPosOld && posEnd.iNegOld == negEnd.iNegOld)
 701     return system.region( posEnd.iPosOld, posEnd.iNegOld);
 702
 703   // Start out with empty region. (Empty used for error returns.)
 704   StringRegion region;
 705
 706   // Add up all remaining p+.
 707   Vec4 pPosJoin;
 708   if ( posEnd.iPosOld == negEnd.iPosOld) {
 709     double xPosJoin = posEnd.xPosOld - negEnd.xPosOld;
 710     if (xPosJoin < 0.) return region;
 711     pPosJoin = system.regionLowPos(posEnd.iPosOld).pHad( xPosJoin, 0., 0., 0.);
 712   } else {
 713     for (int iPosNow = posEnd.iPosOld; iPosNow <= negEnd.iPosOld; ++iPosNow) {
 714       if (iPosNow == posEnd.iPosOld) pPosJoin
 715         += system.regionLowPos(iPosNow).pHad( posEnd.xPosOld, 0., 0., 0.);
 716       else if (iPosNow == negEnd.iPosOld) pPosJoin
 717         += system.regionLowPos(iPosNow).pHad( 1. - negEnd.xPosOld, 0., 0., 0.);
 718       else pPosJoin += system.regionLowPos(iPosNow).pHad( 1., 0., 0., 0.);
 719     }
 720   }
 721
 722   // Add up all remaining p-.
 723   Vec4 pNegJoin;
 724   if ( negEnd.iNegOld == posEnd.iNegOld) {
 725     double xNegJoin = negEnd.xNegOld - posEnd.xNegOld;
 726     if (xNegJoin < 0.) return region;
 727     pNegJoin = system.regionLowNeg(negEnd.iNegOld).pHad( 0., xNegJoin, 0., 0.);
 728   } else {
 729     for (int iNegNow = negEnd.iNegOld; iNegNow <= posEnd.iNegOld; ++iNegNow) {
 730       if (iNegNow == negEnd.iNegOld) pNegJoin
 731         += system.regionLowNeg(iNegNow).pHad( 0., negEnd.xNegOld, 0., 0.);
 732       else if (iNegNow == posEnd.iNegOld) pNegJoin
 733         += system.regionLowNeg(iNegNow).pHad( 0., 1. - posEnd.xNegOld, 0., 0.);
 734       else pNegJoin += system.regionLowNeg(iNegNow).pHad( 0., 1., 0., 0.);
 735     }
 736   }
 737
 738   // For a closed gluon loop pPosJoin == pNegJoin and the above does not work.
 739   // So reshuffle; "perfect" for g g systems, OK in general.
 740   Vec4 pTest = pPosJoin - pNegJoin;
 741   if ( abs(pTest.px()) + abs(pTest.py()) + abs(pTest.pz()) + abs(pTest.e())
 742     < MATCHPOSNEG * (pPosJoin.e() + pNegJoin.e()) ) {
 743     Vec4 delta
 744       = system.regionLowPos(posEnd.iPosOld + 1).pHad( 1., 0., 0., 0.)
 745       - system.regionLowNeg(negEnd.iNegOld + 1).pHad( 0., 1., 0., 0.);
 746     pPosJoin -= delta;
 747     pNegJoin += delta;
 748   }
 749
 750   // Construct a new region from remaining p+ and p-.
 751   region.setUp( pPosJoin, pNegJoin);
 752   if (region.isEmpty) return region;
 753
 754   // Project the existing pTold vectors onto the new directions.
 755   Vec4 pTposOld = system.region( posEnd.iPosOld, posEnd.iNegOld).pHad(
 756     0., 0.,  posEnd.pxOld, posEnd.pyOld);
 757   region.project( pTposOld);
 758   posEnd.pxOld = region.px();
 759   posEnd.pyOld = region.py();
 760   Vec4 pTnegOld = system.region( negEnd.iPosOld, negEnd.iNegOld).pHad(
 761       0., 0.,  negEnd.pxOld, negEnd.pyOld);
 762   region.project( pTnegOld);
 763   negEnd.pxOld = region.px();
 764   negEnd.pyOld = region.py();
 765
 766   // Done.
 767   return region;
 768
 769 }
 770
 771 //--------------------------------------------------------------------------
 772
 773 // Store the hadrons in the normal event record, ordered from one end.
 774
 775 void StringFragmentation::store(Event& event) {
 776
 777   // Starting position.
 778   int iFirst = event.size();
 779
 780   // Copy straight over from first two junction legs.
 781   if (hasJunction) {
 782     for (int i = 0; i < hadrons.size(); ++i)
 783     if (hadrons[i].status() == 85 || hadrons[i].status() == 86)
 784       event.append( hadrons[i] );
 785   }
 786
 787   // Loop downwards, copying all from the positive end.
 788   for (int i = 0; i < hadrons.size(); ++i)
 789     if (hadrons[i].status() == 83) event.append( hadrons[i] );
 790
 791   // Loop upwards, copying all from the negative end.
 792   for (int i = hadrons.size() - 1; i >= 0 ; --i)
 793     if (hadrons[i].status() == 84) event.append( hadrons[i] );
 794   int iLast = event.size() - 1;
 795
 796   // Set decay vertex when this is displaced.
 797   if (event[posEnd.iEnd].hasVertex()) {
 798     Vec4 vDec = event[posEnd.iEnd].vDec();
 799     for (int i = iFirst; i <= iLast; ++i) event[i].vProd( vDec );
 800   }
 801
 802   // Set lifetime of hadrons.
 803   for (int i = iFirst; i <= iLast; ++i)
 804     event[i].tau( event[i].tau0() * rndmPtr->exp() );
 805
 806   // Mark original partons as hadronized and set their daughter range.
 807   for (int i = 0; i < int(iParton.size()); ++i)
 808   if (iParton[i] >= 0) {
 809     event[ iParton[i] ].statusNeg();
 810     event[ iParton[i] ].daughters(iFirst, iLast);
 811   }
 812
 813 }
 814
 815 //--------------------------------------------------------------------------
 816
 817 // Fragment off two of the string legs in to a junction.
 818
 819 bool StringFragmentation::fragmentToJunction(Event& event) {
 820
 821   // Identify range of partons on the three legs.
 822   // (Each leg begins with an iParton[i] = -(10 + 10*junctionNumber + leg),
 823   // and partons then appear ordered from the junction outwards.)
 824   int legBeg[3] = { 0, 0, 0};
 825   int legEnd[3] = { 0, 0, 0};
 826   int leg = -1;
 827   // PS (4/10/2011) Protect against invalid systems
 828   if (iParton[0] > 0) {
 829     infoPtr->errorMsg("Error in StringFragmentation::fragment"
 830                       "ToJunction: iParton[0] not a valid junctionNumber");
 831     return false;
 832   }
 833   for (int i = 0; i < int(iParton.size()); ++i) {
 834     if (iParton[i] < 0) {
 835       if (leg == 2) {
 836         infoPtr->errorMsg("Error in StringFragmentation::fragment"
 837                           "ToJunction: unprocessed multi-junction system");
 838         return false;
 839       }
 840       legBeg[++leg] = i + 1;
 841     }
 842     else legEnd[leg] = i;
 843   }
 844
 845   // Iterate from system rest frame towards the junction rest frame (JRF).
 846   RotBstMatrix MtoJRF, Mstep;
 847   MtoJRF.bstback(pSum);
 848   Vec4 pWTinJRF[3];
 849   int iter = 0;
 850   double errInCM = 0.;
 851   do {
 852     ++iter;
 853
 854     // Find weighted sum of momenta on the three sides of the junction.
 855     for (leg = 0; leg < 3; ++ leg) {
 856       pWTinJRF[leg] = 0.;
 857       double eWeight = 0.;
 858       for (int i = legBeg[leg]; i <= legEnd[leg]; ++i) {
 859         Vec4 pTemp = event[ iParton[i] ].p();
 860         pTemp.rotbst(MtoJRF);
 861         pWTinJRF[leg] += pTemp * exp(-eWeight);
 862         eWeight += pTemp.e() / eNormJunction;
 863         if (eWeight > EJNWEIGHTMAX) break;
 864       }
 865     }
 866
 867     // Store original deviation from 120 degree topology.
 868     if (iter == 1) errInCM = pow2(costheta(pWTinJRF[0], pWTinJRF[1]) + 0.5)
 869       + pow2(costheta(pWTinJRF[0], pWTinJRF[2]) + 0.5)
 870       + pow2(costheta(pWTinJRF[1], pWTinJRF[2]) + 0.5);
 871
 872     // Find new JRF from the set of weighted momenta.
 873     Mstep = junctionRestFrame( pWTinJRF[0], pWTinJRF[1], pWTinJRF[2]);
 874     // Fortran code will not take full step after the first few
 875     // iterations. How implement this in terms of an M matrix??
 876     MtoJRF.rotbst( Mstep );
 877   } while (iter < 3 || (Mstep.deviation() > CONVJNREST && iter < NTRYJNREST) );
 878
 879   // If final deviation from 120 degrees is bigger than in CM then revert.
 880   double errInJRF = pow2(costheta(pWTinJRF[0], pWTinJRF[1]) + 0.5)
 881     + pow2(costheta(pWTinJRF[0], pWTinJRF[2]) + 0.5)
 882     + pow2(costheta(pWTinJRF[1], pWTinJRF[2]) + 0.5);
 883   if (errInJRF > errInCM) {
 884     infoPtr->errorMsg("Warning in StringFragmentation::fragmentTo"
 885       "Junction: bad convergence junction rest frame");
 886     MtoJRF.reset();
 887     MtoJRF.bstback(pSum);
 888   }
 889
 890   // Opposite operation: boost from JRF to original system.
 891   RotBstMatrix MfromJRF = MtoJRF;
 892   MfromJRF.invert();
 893
 894   // Sum leg four-momenta in original frame and in JRF.
 895   Vec4 pInLeg[3], pInJRF[3];
 896   for (leg = 0; leg < 3; ++leg) {
 897     pInLeg[leg] = 0.;
 898     for (int i = legBeg[leg]; i <= legEnd[leg]; ++i)
 899       pInLeg[leg] += event[ iParton[i] ].p();
 900     pInJRF[leg] = pInLeg[leg];
 901     pInJRF[leg].rotbst(MtoJRF);
 902   }
 903
 904   // Pick the two legs with lowest energy in JRF.
 905   int legMin = (pInJRF[0].e() < pInJRF[1].e()) ? 0 : 1;
 906   int legMax = 1 - legMin;
 907   if (pInJRF[2].e() < min(pInJRF[0].e(), pInJRF[1].e()) ) legMin = 2;
 908   else if (pInJRF[2].e() > max(pInJRF[0].e(), pInJRF[1].e()) ) legMax = 2;
 909   int legMid = 3 - legMin - legMax;
 910
 911   // Save info on which status codes belong with the three legs.
 912   int iJunction = (-iParton[0]) / 10 - 1;
 913   event.statusJunction( iJunction, legMin, 85);
 914   event.statusJunction( iJunction, legMid, 86);
 915   event.statusJunction( iJunction, legMax, 83);
 916
 917   // Temporarily copy the partons on the low-energy legs, into the JRF,
 918   // in reverse order, so (anti)quark leg end first.
 919   vector<int> iPartonMin;
 920   for (int i = legEnd[legMin]; i >= legBeg[legMin]; --i) {
 921     int iNew = event.append( event[ iParton[i] ] );
 922     event[iNew].rotbst(MtoJRF);
 923     iPartonMin.push_back( iNew );
 924   }
 925   vector<int> iPartonMid;
 926   for (int i = legEnd[legMid]; i >= legBeg[legMid]; --i) {
 927     int iNew = event.append( event[ iParton[i] ] );
 928     event[iNew].rotbst(MtoJRF);
 929     iPartonMid.push_back( iNew );
 930   }
 931
 932   // Find final weighted sum of momenta on each of the two legs.
 933   double eWeight = 0.;
 934   pWTinJRF[legMin] = 0.;
 935   for (int i = iPartonMin.size() - 1; i >= 0; --i) {
 936     pWTinJRF[legMin] += event[ iPartonMin[i] ].p() * exp(-eWeight);
 937     eWeight += event[ iPartonMin[i] ].e() / eNormJunction;
 938     if (eWeight > EJNWEIGHTMAX) break;
 939   }
 940   eWeight = 0.;
 941   pWTinJRF[legMid] = 0.;
 942   for (int i = iPartonMid.size() - 1; i >= 0; --i) {
 943     pWTinJRF[legMid] += event[ iPartonMid[i] ].p() * exp(-eWeight);
 944     eWeight += event[ iPartonMid[i] ].e() / eNormJunction;
 945     if (eWeight > EJNWEIGHTMAX) break;
 946   }
 947
 948   // Define fictitious opposing partons in JRF and store as string ends.
 949   Vec4 pOppose = pWTinJRF[legMin];
 950   pOppose.flip3();
 951   int idOppose = (rndmPtr->flat() > 0.5) ? 2 : 1;
 952   if (event[ iPartonMin[0] ].col() > 0) idOppose = -idOppose;
 953   int iOppose = event.append( idOppose, 77, 0, 0, 0, 0, 0, 0,
 954     pOppose, 0.);
 955   iPartonMin.push_back( iOppose);
 956   pOppose = pWTinJRF[legMid];
 957   pOppose.flip3();
 958   idOppose = (rndmPtr->flat() > 0.5) ? 2 : 1;
 959   if (event[ iPartonMid[0] ].col() > 0) idOppose = -idOppose;
 960   iOppose = event.append( idOppose, 77, 0, 0, 0, 0, 0, 0,
 961     pOppose, 0.);
 962   iPartonMid.push_back( iOppose);
 963
 964   // Set up kinematics of string evolution in low-energy temporary systems.
 965   systemMin.setUp(iPartonMin, event);
 966   systemMid.setUp(iPartonMid, event);
 967
 968   // Outer fallback loop, when too little energy left for third leg.
 969   int idMin = 0;
 970   int idMid = 0;
 971   Vec4 pDiquark;
 972   for ( int iTryOuter = 0; ; ++iTryOuter) {
 973
 974     // Middle fallback loop, when much unused energy in leg remnants.
 975     double eLeftMin = 0.;
 976     double eLeftMid = 0.;
 977     for ( int iTryMiddle = 0; ; ++iTryMiddle) {
 978
 979       // Loop over the two lowest-energy legs.
 980       for (int legLoop = 0; legLoop < 2; ++ legLoop) {
 981         int legNow = (legLoop == 0) ? legMin : legMid;
 982
 983         // Read in properties specific to this leg.
 984         StringSystem& systemNow = (legLoop == 0) ? systemMin : systemMid;
 985         int idPos = (legLoop == 0) ? event[ iPartonMin[0] ].id()
 986           : event[ iPartonMid[0] ].id();
 987         idOppose = (legLoop == 0) ? event[ iPartonMin.back() ].id()
 988           : event[ iPartonMid.back() ].id();
 989         double eInJRF = pInJRF[legNow].e();
 990         int statusHad = (legLoop == 0) ? 85 : 86;
 991
 992         // Inner fallback loop, when a diquark comes in to junction.
 993         double eUsed = 0.;
 994         for ( int iTryInner = 0; ; ++iTryInner) {
 995           if (iTryInner > 2 * NTRYJNMATCH) {
 996             infoPtr->errorMsg("Error in StringFragmentation::fragment"
 997               "ToJunction: caught in junction flavour loop");
 998             event.popBack( iPartonMin.size() + iPartonMid.size() );
 999             return false;
1000           }
1001           bool needBaryon = (abs(idPos) > 10 && iTryInner > NTRYJNMATCH);
1002           double eExtra   = (iTryInner > NTRYJNMATCH) ? EEXTRAJNMATCH : 0.;
1003
1004           // Set up two string ends, and begin fragmentation loop.
1005           setStartEnds(idPos, idOppose, systemNow);
1006           eUsed = 0.;
1007           int nHadrons = 0;
1008           bool noNegE = true;
1009           for ( ; ; ++nHadrons) {
1010
1011             // Construct trial hadron from positive end.
1012             posEnd.newHadron();
1013             Vec4 pHad = posEnd.kinematicsHadron(systemNow);
1014
1015             // Negative energy signals failure in construction.
1016             if (pHad.e() < 0. ) { noNegE = false; break; }
1017
1018             // Break if passed system midpoint ( = junction) in energy.
1019             // Exceptions: small systems, and/or with diquark end.
1020             bool delayedBreak = false;
1021             if (eUsed + pHad.e() + eExtra > eInJRF) {
1022               if (nHadrons > 0 || !needBaryon) break;
1023               delayedBreak = true;
1024             }
1025
1026             // Else construct kinematics of the new hadron and store it.
1027             hadrons.append( posEnd.idHad, statusHad, iPos, iNeg,
1028               0, 0, 0, 0, pHad, posEnd.mHad);
1029
1030             // Update string end and remaining momentum.
1031             posEnd.update();
1032             eUsed += pHad.e();
1033
1034             // Delayed break in small systems, and/or with diquark end.
1035             if (delayedBreak) {
1036               ++nHadrons;
1037               break;
1038             }
1039           }
1040
1041           // End of fragmentation loop. Inner loopback if ends on a diquark.
1042           if ( noNegE && abs(posEnd.flavOld.id) < 10 ) break;
1043           hadrons.popBack(nHadrons);
1044         }
1045
1046         // End of one-leg fragmentation. Store end quark and remnant energy.
1047         if (legNow == legMin) {
1048           idMin = posEnd.flavOld.id;
1049           eLeftMin = eInJRF - eUsed;
1050         } else {
1051           idMid = posEnd.flavOld.id;
1052           eLeftMid = eInJRF - eUsed;
1053         }
1054       }
1055
1056       // End of both-leg fragmentation.
1057       // Middle loopback if too much energy left.
1058       double eTrial = eBothLeftJunction + rndmPtr->flat() * eMaxLeftJunction;
1059       if (iTryMiddle > NTRYJNMATCH
1060         || ( min( eLeftMin, eLeftMid) < eBothLeftJunction
1061         && max( eLeftMin, eLeftMid) < eTrial ) ) break;
1062       hadrons.clear();
1063     }
1064
1065     // Boost hadrons away from the JRF to the original frame.
1066     for (int i = 0; i < hadrons.size(); ++i) {
1067       hadrons[i].rotbst(MfromJRF);
1068       // Recalculate energy to compensate for numerical precision loss
1069       // in iterative calculation of MfromJRF.
1070       hadrons[i].e( hadrons[i].eCalc() );
1071       pJunctionHadrons += hadrons[i].p();
1072     }
1073
1074     // Outer loopback if combined diquark mass too negative
1075     // or too little energy left in third leg.
1076     pDiquark = pInLeg[legMin] + pInLeg[legMid] - pJunctionHadrons;
1077     double m2Left = m2( pInLeg[legMax], pDiquark);
1078     if (iTryOuter >  NTRYJNMATCH || (pDiquark.mCalc() > MDIQUARKMIN
1079       && m2Left > eMinLeftJunction * pInLeg[legMax].e()) ) break;
1080     hadrons.clear();
1081     pJunctionHadrons = 0.;
1082   }
1083
1084   // Now found solution; no more loopback. Remove temporary parton copies.
1085   event.popBack( iPartonMin.size() + iPartonMid.size() );
1086
1087   // Construct and store an effective diquark string end from the
1088   // two remnant quark ends, for temporary usage.
1089   int    idDiquark = flavSelPtr->makeDiquark( idMin, idMid);
1090   double mDiquark  = pDiquark.mCalc();
1091   int    iDiquark  = event.append( idDiquark, 78, 0, 0, 0, 0, 0, 0,
1092     pDiquark, mDiquark);
1093
1094   // Find the partons on the last leg, again in reverse order.
1095   vector<int> iPartonMax;
1096   for (int i = legEnd[legMax]; i >= legBeg[legMax]; --i)
1097     iPartonMax.push_back( iParton[i] );
1098   iPartonMax.push_back( iDiquark );
1099
1100   // Recluster gluons nearby to diquark end when taken too much energy.
1101   int iPsize       = iPartonMax.size();
1102   double m0Diquark = event[iDiquark].m0();
1103   while (iPsize > 2) {
1104     Vec4 pGluNear = event[ iPartonMax[iPsize - 2] ].p();
1105     if ( pDiquark.mCalc() > 0.
1106       && (pDiquark + 0.5 * pGluNear).mCalc() > m0Diquark + mJoin ) break;
1107     pDiquark += pGluNear;
1108     event[iDiquark].p( pDiquark );
1109     event[iDiquark].m( pDiquark.mCalc() );
1110     iPartonMax.pop_back();
1111     --iPsize;
1112     iPartonMax[iPsize - 1] = iDiquark;
1113   }
1114
1115   // Modify parton list to remaining leg + remnant of the first two.
1116   iParton = iPartonMax;
1117
1118   // Done.
1119   return true;
1120 }
1121
1122 //--------------------------------------------------------------------------
1123
1124 // Find the boost matrix to the rest frame of a junction,
1125 // given the three respective endpoint four-momenta.
1126
1127 RotBstMatrix StringFragmentation::junctionRestFrame(Vec4& p0, Vec4& p1,
1128   Vec4& p2) {
1129
1130   // Calculate masses and other invariants.
1131   Vec4 pSumJun  = p0 + p1 + p2;
1132   double sHat   = pSumJun.m2Calc();
1133   double pp[3][3];
1134   pp[0][0]      = p0.m2Calc();
1135   pp[1][1]      = p1.m2Calc();
1136   pp[2][2]      = p2.m2Calc();
1137   pp[0][1] = pp[1][0] = p0 * p1;
1138   pp[0][2] = pp[2][0] = p0 * p2;
1139   pp[1][2] = pp[2][1] = p1 * p2;
1140
1141   // Requirement (eiMax)_j = pi*pj/mj < (eiMax)_k = pi*pk/mk, used below,
1142   // here rewritten as pi*pj * mk < pi*pk * mj and squared.
1143   double eMax01 = pow2(pp[0][1]) * pp[2][2];
1144   double eMax02 = pow2(pp[0][2]) * pp[1][1];
1145   double eMax12 = pow2(pp[1][2]) * pp[0][0];
1146
1147   // Initially pick i to be the most massive parton. but allow other tries.
1148   int i = (pp[1][1] > pp[0][0]) ? 1 : 0;
1149   if (pp[2][2] > max(pp[0][0], pp[1][1])) i = 2;
1150   int j, k;
1151   double ei     = 0.;
1152   double ej     = 0.;
1153   double ek     = 0.;
1154   for (int iTry = 0; iTry < 3; ++iTry) {
1155
1156     // Pick j to give minimal eiMax, and k the third vector.
1157     if (i == 0) j = (eMax02 < eMax01) ? 2 : 1;
1158     else if (i == 1) j = (eMax12 < eMax01) ? 2 : 0;
1159     else j = (eMax12 < eMax02) ? 1 : 0;
1160     k = 3 - i - j;
1161
1162     // Alternative names according to i, j, k conventions.
1163     double m2i  = pp[i][i];
1164     double m2j  = pp[j][j];
1165     double m2k  = pp[k][k];
1166     double pipj = pp[i][j];
1167     double pipk = pp[i][k];
1168     double pjpk = pp[j][k];
1169
1170     // Trivial to find new parton energies if all three partons are massless.
1171     if (m2i < M2MAXJRF) {
1172       ei        = sqrt( 2. * pipk * pipj / (3. * pjpk) );
1173       ej        = sqrt( 2. * pjpk * pipj / (3. * pipk) );
1174       ek        = sqrt( 2. * pipk * pjpk / (3. * pipj) );
1175
1176     // Else find three-momentum range for parton i and values at extremes.
1177     } else {
1178       // Minimum when i is at rest.
1179       double piMin = 0.;
1180       double eiMin = sqrt(m2i);
1181       double ejMin = pipj / eiMin;
1182       double ekMin = pipk / eiMin;
1183       double pjMin = sqrtpos( ejMin*ejMin - m2j );
1184       double pkMin = sqrtpos( ekMin*ekMin - m2k );
1185       double fMin  = ejMin * ekMin + 0.5 * pjMin * pkMin - pjpk;
1186       // Maximum estimated when j + k is at rest, alternatively j at rest.
1187       double eiMax = (pipj + pipk) / sqrt(m2j + m2k + 2. * pjpk);
1188       if (m2j > M2MAXJRF) eiMax = min( eiMax, pipj / sqrt(m2j) );
1189       double piMax = sqrtpos( eiMax*eiMax - m2i );
1190       double temp  = eiMax*eiMax - 0.25 *piMax*piMax;
1191       double pjMax = (eiMax * sqrtpos( pipj*pipj - m2j * temp )
1192         - 0.5 * piMax * pipj) / temp;
1193       double pkMax = (eiMax * sqrtpos( pipk*pipk - m2k * temp )
1194         - 0.5 * piMax * pipk) / temp;
1195       double ejMax = sqrt(pjMax*pjMax + m2j);
1196       double ekMax = sqrt(pkMax*pkMax + m2k);
1197       double fMax  = ejMax * ekMax + 0.5 * pjMax * pkMax - pjpk;
1198
1199       // If unexpected values at upper endpoint then pick another parton.
1200       if (fMax > 0.) {
1201         int iPrel = (i + 1)%3;
1202         if (pp[iPrel][iPrel] > M2MAXJRF) {i = iPrel; continue;}
1203         ++iTry;
1204         iPrel = (i + 2)%3;
1205         if (iTry < 3 && pp[iPrel][iPrel] > M2MAXJRF) {i = iPrel; continue;}
1206       }
1207
1208       // Start binary + linear search to find solution inside range.
1209       int iterMin = 0;
1210       int iterMax = 0;
1211       double pi   = 0.5 * (piMin + piMax);
1212       for (int iter = 0; iter < NTRYJRFEQ; ++iter) {
1213
1214         // Derive momentum of other two partons and distance to root.
1215         ei = sqrt(pi*pi + m2i);
1216         temp = ei*ei - 0.25 * pi*pi;
1217         double pj = (ei * sqrtpos( pipj*pipj - m2j * temp )
1218           - 0.5 * pi * pipj) / temp;
1219         double pk = (ei * sqrtpos( pipk*pipk - m2k * temp )
1220           - 0.5 * pi * pipk) / temp;
1221         ej = sqrt(pj*pj + m2j);
1222         ek = sqrt(pk*pk + m2k);
1223         double fNow = ej * ek + 0.5 * pj * pk - pjpk;
1224
1225         // Replace lower or upper bound by new value.
1226         if (fNow > 0.) { ++iterMin; piMin = pi; fMin = fNow;}
1227         else {++iterMax; piMax = pi; fMax = fNow;}
1228
1229         // Pick next i momentum to explore, hopefully closer to root.
1230         if (2 * iter < NTRYJRFEQ
1231           && (iterMin < 2 || iterMax < 2 || 4 * iter < NTRYJRFEQ))
1232           { pi = 0.5 * (piMin + piMax); continue;}
1233         if (fMin < 0. || fMax > 0. || abs(fNow) < CONVJRFEQ * sHat) break;
1234         pi = piMin + (piMax - piMin) * fMin / (fMin - fMax);
1235       }
1236
1237     // If arrived here then either succeeded or exhausted possibilities.
1238     } break;
1239   }
1240
1241   // Now we know the energies in the junction rest frame.
1242   double eNew[3] = { 0., 0., 0.};
1243   eNew[i] = ei;
1244   eNew[j] = ej;
1245   eNew[k] = ek;
1246
1247   // Boost (copy of) partons to their rest frame.
1248   RotBstMatrix Mmove;
1249   Vec4 p0cm = p0;
1250   Vec4 p1cm = p1;
1251   Vec4 p2cm = p2;
1252   Mmove.bstback(pSumJun);
1253   p0cm.rotbst(Mmove);
1254   p1cm.rotbst(Mmove);
1255   p2cm.rotbst(Mmove);
1256
1257   // Construct difference vectors and the boost to junction rest frame.
1258   Vec4 pDir01      = p0cm / p0cm.e() - p1cm / p1cm.e();
1259   Vec4 pDir02      = p0cm / p0cm.e() - p2cm / p2cm.e();
1260   double pDiff01   = pDir01.pAbs2();
1261   double pDiff02   = pDir02.pAbs2();
1262   double pDiff0102 = dot3(pDir01, pDir02);
1263   double eDiff01   = eNew[0] / p0cm.e() - eNew[1] / p1cm.e();
1264   double eDiff02   = eNew[0] / p0cm.e() - eNew[2] / p2cm.e();
1265   double denom     = pDiff01 * pDiff02 - pDiff0102*pDiff0102;
1266   double coef01    = (eDiff01 * pDiff02 - eDiff02 * pDiff0102) / denom;
1267   double coef02    = (eDiff02 * pDiff01 - eDiff01 * pDiff0102) / denom;
1268   Vec4 vJunction   = coef01 * pDir01 + coef02 * pDir02;
1269   vJunction.e( sqrt(1. + vJunction.pAbs2()) );
1270
1271   // Add two boosts, giving final result.
1272   Mmove.bst(vJunction);
1273   return Mmove;
1274
1275 }
1276
1277 //--------------------------------------------------------------------------
1278
1279 // When string fragmentation has failed several times,
1280 // try to join some more nearby partons.
1281
1282 int StringFragmentation::extraJoin(double facExtra, Event& event) {
1283
1284   // Keep on looping while pairs found below joining threshold.
1285   int nJoin  = 0;
1286   int iPsize = iParton.size();
1287   while (iPsize > 2) {
1288
1289     // Look for the pair of neighbour partons (along string) with
1290     // the smallest invariant mass (subtracting quark masses).
1291     int iJoinMin    = -1;
1292     double mJoinMin = 2. * facExtra * mJoin;
1293     for (int i = 0; i < iPsize - 1; ++i) {
1294       Particle& parton1 = event[ iParton[i] ];
1295       Particle& parton2 = event[ iParton[i + 1] ];
1296       Vec4 pSumNow;
1297       pSumNow += (parton2.isGluon()) ? 0.5 * parton1.p() : parton1.p();
1298       pSumNow += (parton2.isGluon()) ? 0.5 * parton2.p() : parton2.p();
1299       double mJoinNow = pSumNow.mCalc();
1300       if (!parton1.isGluon()) mJoinNow -= parton1.m0();
1301       if (!parton2.isGluon()) mJoinNow -= parton2.m0();
1302       if (mJoinNow < mJoinMin) { iJoinMin = i; mJoinMin = mJoinNow; }
1303     }
1304
1305     // Decide whether to join, if not finished.
1306     if (iJoinMin < 0 || mJoinMin > facExtra * mJoin) return nJoin;
1307     ++nJoin;
1308
1309     // Create new joined parton.
1310     int iJoin1  = iParton[iJoinMin];
1311     int iJoin2  = iParton[iJoinMin + 1];
1312     int idNew   = (event[iJoin1].isGluon()) ? event[iJoin2].id()
1313                                             : event[iJoin1].id();
1314     int colNew  = event[iJoin1].col();
1315     int acolNew = event[iJoin2].acol();
1316     if (colNew == acolNew) {
1317       colNew    = event[iJoin2].col();
1318       acolNew   = event[iJoin1].acol();
1319     }
1320     Vec4 pNew   = event[iJoin1].p() + event[iJoin2].p();
1321
1322     // Append joined parton to event record and reduce parton list.
1323     int iNew = event.append( idNew, 73, min(iJoin1, iJoin2),
1324       max(iJoin1, iJoin2), 0, 0, colNew, acolNew, pNew, pNew.mCalc() );
1325     iParton[iJoinMin] = iNew;
1326     for (int i = iJoinMin + 1; i < iPsize - 1; ++i)
1327       iParton[i] = iParton[i + 1];
1328     iParton.pop_back();
1329     --iPsize;
1330
1331   // Done.
1332   }
1333   return nJoin;
1334 }
1335
1336 //==========================================================================
1337
1338 } // end namespace Pythia8