PYTHIA8/pythia8130/src/StringFragmentation.cxx

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