PYTHIA8/pythia8130/src/HadronLevel.cxx

   1 // HadronLevel.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 HadronLevel class.
   7
   8 #include "HadronLevel.h"
   9
  10 namespace Pythia8 {
  11
  12 //**************************************************************************
  13
  14 // The HadronLevel class.
  15
  16 //*********
  17
  18 // Constants: could be changed here if desired, but normally should not.
  19 // These are of technical nature, as described for each.
  20
  21 // For breaking J-J string, pick a Gamma by taking a step with fictitious mass.
  22 const double HadronLevel::JJSTRINGM2MAX  = 25.;
  23 const double HadronLevel::JJSTRINGM2FRAC = 0.1;
  24
  25 // Iterate junction rest frame boost until convergence or too many tries.
  26 const double HadronLevel::CONVJNREST     = 1e-5;
  27 const int HadronLevel::NTRYJNREST        = 20;
  28
  29 // Typical average transvere primary hadron mass <mThad>.
  30 const double HadronLevel::MTHAD          = 0.9;
  31
  32 //*********
  33
  34 // Find settings. Initialize HadronLevel classes as required.
  35
  36 bool HadronLevel::init(Info* infoPtrIn, TimeShower* timesDecPtr,
  37   DecayHandler* decayHandlePtr, vector<int> handledParticles) {
  38
  39   // Save pointer.
  40   infoPtr        = infoPtrIn;
  41
  42   // Main flags.
  43   doHadronize    = Settings::flag("HadronLevel:Hadronize");
  44   doDecay        = Settings::flag("HadronLevel:Decay");
  45   doBoseEinstein = Settings::flag("HadronLevel:BoseEinstein");
  46
  47   // Boundary mass between string and ministring handling.
  48   mStringMin     = Settings::parm("HadronLevel:mStringMin");
  49
  50   // For junction processing.
  51   eNormJunction  = Settings::parm("StringFragmentation:eNormJunction");
  52
  53   // Particles that should decay or not before Bose-Einstein stage.
  54   widthSepBE     = Settings::parm("BoseEinstein:widthSep");
  55
  56   // Initialize string and ministring fragmentation.
  57   stringFrag.init(infoPtr, &flavSel, &pTSel, &zSel);
  58   ministringFrag.init(infoPtr, &flavSel);
  59
  60   // Initialize particle decays.
  61   decays.init(infoPtr, timesDecPtr, &flavSel, decayHandlePtr,
  62     handledParticles);
  63
  64   // Initialize BoseEinstein.
  65   boseEinstein.init(infoPtr);
  66
  67   // Initialize auxiliary administrative classes.
  68   colConfig.init( &flavSel);
  69
  70   // Initialize auxiliary fragmentation classes.
  71   flavSel.init();
  72   pTSel.init();
  73   zSel.init();
  74
  75   // Done.
  76   return true;
  77
  78 }
  79
  80 //*********
  81
  82 // Hadronize and decay the next parton-level.
  83
  84 bool HadronLevel::next( Event& event) {
  85
  86   // Colour-octet onia states must be decayed to singlet + gluon.
  87   if (!decayOctetOnia(event)) return false;
  88
  89   // Possibility of hadronization inside decay, but then no BE second time.
  90   bool moreToDo;
  91   bool doBoseEinsteinNow = doBoseEinstein;
  92   do {
  93     moreToDo = false;
  94
  95     // First part: string fragmentation.
  96     if (doHadronize) {
  97
  98       // Find the complete colour singlet configuration of the event.
  99       if (!findSinglets( event)) return false;
 100
 101       // Process all colour singlet (sub)system
 102       for (int iSub = 0; iSub < colConfig.size(); ++iSub) {
 103
 104         // Collect sequentially all partons in a colour singlet subsystem.
 105         colConfig.collect(iSub, event);
 106
 107         // String fragmentation of each colour singlet (sub)system.
 108         if ( colConfig[iSub].massExcess > mStringMin ) {
 109           if (!stringFrag.fragment( iSub, colConfig, event)) return false;
 110
 111         // Low-mass string treated separately. Tell if diffractive system.
 112         } else {
 113           bool isDiff = infoPtr->isDiffractiveA()
 114             || infoPtr->isDiffractiveB();
 115           if (!ministringFrag.fragment( iSub, colConfig, event, isDiff))
 116             return false;
 117         }
 118       }
 119     }
 120
 121     // Second part: sequential decays of short-lived particles (incl. K0).
 122     if (doDecay) {
 123
 124       // Loop through all entries to find those that should decay.
 125       int iDec = 0;
 126       do {
 127         if ( event[iDec].isFinal() && event[iDec].canDecay()
 128           && event[iDec].mayDecay() && (event[iDec].idAbs() == 311
 129           || event[iDec].mWidth() > widthSepBE) ) {
 130           decays.decay( iDec, event);
 131           if (decays.moreToDo()) moreToDo = true;
 132         }
 133       } while (++iDec < event.size());
 134     }
 135
 136     // Third part: include Bose-Einstein effects among current particles.
 137     if (doBoseEinsteinNow) {
 138       if (!boseEinstein.shiftEvent(event)) return false;
 139       doBoseEinsteinNow = false;
 140     }
 141
 142     // Fourth part: sequential decays also of long-lived particles.
 143     if (doDecay) {
 144
 145       // Loop through all entries to find those that should decay.
 146       int iDec = 0;
 147       do {
 148         if ( event[iDec].isFinal() && event[iDec].canDecay()
 149           && event[iDec].mayDecay() ) {
 150           decays.decay( iDec, event);
 151           if (decays.moreToDo()) moreToDo = true;
 152         }
 153       } while (++iDec < event.size());
 154     }
 155
 156
 157   // Normally done first time around, but sometimes not (e.g. Upsilon).
 158   } while (moreToDo);
 159
 160   // Done.
 161   return true;
 162
 163 }
 164
 165 //*********
 166
 167 // Allow more decays if on/off switches changed.
 168 // Note: does not do sequential hadronization, e.g. for Upsilon.
 169
 170 bool HadronLevel::moreDecays( Event& event) {
 171
 172   // Colour-octet onia states must be decayed to singlet + gluon.
 173   if (!decayOctetOnia(event)) return false;
 174
 175   // Loop through all entries to find those that should decay.
 176   int iDec = 0;
 177   do {
 178     if ( event[iDec].isFinal() && event[iDec].canDecay()
 179       && event[iDec].mayDecay() ) decays.decay( iDec, event);
 180   } while (++iDec < event.size());
 181
 182   // Done.
 183   return true;
 184
 185 }
 186
 187 //*********
 188
 189 // Decay colour-octet onium states.
 190
 191 bool HadronLevel::decayOctetOnia(Event& event) {
 192
 193   // Onium states to be decayed.
 194   int idOnium[6] = { 9900443, 9900441, 9910441,
 195                      9900553, 9900551, 9910551 };
 196
 197   // Loop over particles and identify onia.
 198   for (int iDec = 0; iDec < event.size(); ++iDec)
 199   if (event[iDec].isFinal()) {
 200     int id = event[iDec].id();
 201     bool isOnium = false;
 202     for (int j = 0; j < 6; ++j) if (id == idOnium[j]) isOnium = true;
 203
 204     // Decay any onia encountered.
 205     if (isOnium) {
 206       if (!decays.decay( iDec, event)) return false;
 207
 208       // Set colour flow by hand: gluon inherits octet-onium state.
 209       int iGlu = event.size() - 1;
 210       event[iGlu].cols( event[iDec].col(), event[iDec].acol() );
 211     }
 212   }
 213
 214   // Done.
 215   return true;
 216
 217 }
 218
 219 //*********
 220
 221 // Trace colour flow in the event to form colour singlet subsystems.
 222
 223 bool HadronLevel::findSinglets(Event& event) {
 224
 225   // Find a list of final partons and of all colour ends and gluons.
 226   iColEnd.resize(0);
 227   iAcolEnd.resize(0);
 228   iColAndAcol.resize(0);
 229   for (int i = 0; i < event.size(); ++ i) if (event[i].isFinal()) {
 230     if (event[i].col() > 0 && event[i].acol() > 0) iColAndAcol.push_back(i);
 231     else if (event[i].col() > 0) iColEnd.push_back(i);
 232     else if (event[i].acol() > 0) iAcolEnd.push_back(i);
 233   }
 234
 235   // Begin arrange the partons into separate colour singlets.
 236   colConfig.clear();
 237   iPartonJun.resize(0);
 238   iPartonAntiJun.resize(0);
 239
 240   // Junctions: loop over them, and identify kind.
 241   for (int iJun = 0; iJun < event.sizeJunction(); ++iJun)
 242   if (event.remainsJunction(iJun)) {
 243     event.remainsJunction(iJun, false);
 244     int kindJun = event.kindJunction(iJun);
 245     iParton.resize(0);
 246
 247     // First kind: three (anti)colours out from junction.
 248     if (kindJun == 1 || kindJun == 2) {
 249       for (int iCol = 0; iCol < 3; ++iCol) {
 250         int indxCol = event.colJunction(iJun, iCol);
 251         iParton.push_back( -(10 + 10 * iJun + iCol) );
 252         if (kindJun == 1 && !traceFromAcol(indxCol, event, iJun, iCol))
 253           return false;
 254         if (kindJun == 2 && !traceFromCol(indxCol, event, iJun, iCol))
 255           return false;
 256       }
 257     }
 258
 259     // Keep in memory a junction hooked up with an antijunction,
 260     // else store found single-junction system.
 261     int nNeg = 0;
 262     for (int i = 0; i < int(iParton.size()); ++i) if (iParton[i] < 0)
 263       ++nNeg;
 264     if (nNeg > 3 && kindJun == 1) {
 265       for (int i = 0; i < int(iParton.size()); ++i)
 266         iPartonJun.push_back(iParton[i]);
 267     } else if (nNeg > 3 && kindJun == 2) {
 268       for (int i = 0; i < int(iParton.size()); ++i)
 269         iPartonAntiJun.push_back(iParton[i]);
 270     } else {
 271       // A junction may be eliminated by insert if two quarks are nearby.
 272       int nJunOld = event.sizeJunction();
 273       colConfig.insert(iParton, event);
 274       if (event.sizeJunction() < nJunOld) --iJun;
 275     }
 276   }
 277
 278   // Split junction-antijunction system into two, and store those.
 279   // (Only one system in extreme cases, and then second empty.)
 280   if (iPartonJun.size() > 0 && iPartonAntiJun.size() > 0) {
 281     if (!splitJunctionPair(event)) return false;
 282     colConfig.insert(iPartonJun, event);
 283     if (iPartonAntiJun.size() > 0)
 284       colConfig.insert(iPartonAntiJun, event);
 285   // Error if only one of junction and antijuction left here.
 286   } else if (iPartonJun.size() > 0 || iPartonAntiJun.size() > 0) {
 287     infoPtr->errorMsg("Error in HadronLevel::findSinglets: "
 288       "unmatched (anti)junction");
 289     return false;
 290   }
 291
 292   // Open strings: pick up each colour end and trace to its anticolor end.
 293   for (int iEnd = 0; iEnd < int(iColEnd.size()); ++iEnd) {
 294     iParton.resize(0);
 295     iParton.push_back( iColEnd[iEnd] );
 296     int indxCol = event[ iColEnd[iEnd] ].col();
 297     if (!traceFromCol(indxCol, event)) return false;
 298
 299     // Store found open string system. Analyze its properties.
 300     colConfig.insert(iParton, event);
 301   }
 302
 303   // Closed strings : begin at any gluon and trace until back at it.
 304   while (iColAndAcol.size() > 0) {
 305     iParton.resize(0);
 306     iParton.push_back( iColAndAcol[0] );
 307     int indxCol = event[ iColAndAcol[0] ].col();
 308     int indxAcol = event[ iColAndAcol[0] ].acol();
 309     iColAndAcol[0] = iColAndAcol.back();
 310     iColAndAcol.pop_back();
 311     if (!traceInLoop(indxCol, indxAcol, event)) return false;
 312
 313     // Store found closed string system. Analyze its properties.
 314     colConfig.insert(iParton, event);
 315   }
 316
 317   // Done.
 318   return true;
 319
 320 }
 321
 322 //*********
 323
 324 // Trace a colour line, from a colour to an anticolour.
 325
 326 bool HadronLevel::traceFromCol(int indxCol, Event& event, int iJun,
 327   int iCol) {
 328
 329   // Junction kind, if any.
 330   int kindJun = (iJun >= 0) ? event.kindJunction(iJun) : 0;
 331
 332   // Begin to look for a matching anticolour.
 333   int loop = 0;
 334   int loopMax = iColAndAcol.size() + 2;
 335   bool hasFound = false;
 336   do {
 337     ++loop;
 338     hasFound= false;
 339
 340     // First check list of matching anticolour ends.
 341     for (int i = 0; i < int(iAcolEnd.size()); ++i)
 342     if (event[ iAcolEnd[i] ].acol() == indxCol) {
 343       iParton.push_back( iAcolEnd[i] );
 344       indxCol = 0;
 345       iAcolEnd[i] = iAcolEnd.back();
 346       iAcolEnd.pop_back();
 347       hasFound = true;
 348       break;
 349     }
 350
 351     // Then check list of intermediate gluons.
 352     if (!hasFound)
 353     for (int i = 0; i < int(iColAndAcol.size()); ++i)
 354     if (event[ iColAndAcol[i] ].acol() == indxCol) {
 355       iParton.push_back( iColAndAcol[i] );
 356
 357       // Update to new colour. Remove gluon.
 358       indxCol = event[ iColAndAcol[i] ].col();
 359       if (kindJun > 0) event.endColJunction(iJun, iCol, indxCol);
 360       iColAndAcol[i] = iColAndAcol.back();
 361       iColAndAcol.pop_back();
 362       hasFound = true;
 363       break;
 364     }
 365
 366     // In a pinch, check list of end colours on other (anti)junction.
 367     if (!hasFound && kindJun == 2 && event.sizeJunction() > 1)
 368     for (int iAntiJun = 0; iAntiJun < event.sizeJunction(); ++iAntiJun)
 369     if (iAntiJun != iJun && event.kindJunction(iAntiJun) == 1)
 370     for (int iColAnti = 0; iColAnti < 3; ++iColAnti)
 371     if (event.endColJunction(iAntiJun, iColAnti) == indxCol) {
 372       iParton.push_back( -(10 + 10 * iAntiJun + iColAnti) );
 373       indxCol = 0;
 374       hasFound = true;
 375       break;
 376     }
 377
 378   // Keep on tracing via gluons until reached end of leg.
 379   } while (hasFound && indxCol > 0 && loop < loopMax);
 380
 381   // Something went wrong in colour tracing.
 382   if (!hasFound || loop == loopMax) {
 383     infoPtr->errorMsg("Error in HadronLevel::traceFromCol: "
 384       "colour tracing failed");
 385     return false;
 386   }
 387
 388   // Done.
 389   return true;
 390
 391 }
 392
 393 //*********
 394
 395 // Trace a colour line, from an anticolour to a colour.
 396
 397 bool HadronLevel::traceFromAcol(int indxCol, Event& event, int iJun,
 398   int iCol) {
 399
 400   // Junction kind, if any.
 401   int kindJun = (iJun >= 0) ? event.kindJunction(iJun) : 0;
 402
 403   // Begin to look for a matching colour.
 404   int loop = 0;
 405   int loopMax = iColAndAcol.size() + 2;
 406   bool hasFound = false;
 407   do {
 408     ++loop;
 409     hasFound= false;
 410
 411     // First check list of matching colour ends.
 412     for (int i = 0; i < int(iColEnd.size()); ++i)
 413     if (event[ iColEnd[i] ].col() == indxCol) {
 414       iParton.push_back( iColEnd[i] );
 415       indxCol = 0;
 416       iColEnd[i] = iColEnd.back();
 417       iColEnd.pop_back();
 418       hasFound = true;
 419       break;
 420     }
 421
 422     // Then check list of intermediate gluons.
 423     if (!hasFound)
 424     for (int i = 0; i < int(iColAndAcol.size()); ++i)
 425     if (event[ iColAndAcol[i] ].col() == indxCol) {
 426       iParton.push_back( iColAndAcol[i] );
 427
 428       // Update to new colour. Remove gluon.
 429       indxCol = event[ iColAndAcol[i] ].acol();
 430       if (kindJun > 0) event.endColJunction(iJun, iCol, indxCol);
 431       iColAndAcol[i] = iColAndAcol.back();
 432       iColAndAcol.pop_back();
 433       hasFound = true;
 434       break;
 435     }
 436
 437     // In a pinch, check list of colours on other (anti)junction.
 438     if (!hasFound && kindJun == 1 && event.sizeJunction() > 1)
 439     for (int iAntiJun = 0; iAntiJun < event.sizeJunction(); ++iAntiJun)
 440     if (iAntiJun != iJun && event.kindJunction(iAntiJun) == 2)
 441     for (int iColAnti = 0; iColAnti < 3; ++iColAnti)
 442     if (event.endColJunction(iAntiJun, iColAnti) == indxCol) {
 443       iParton.push_back( -(10 + 10 * iAntiJun + iColAnti) );
 444       indxCol = 0;
 445       hasFound = true;
 446       break;
 447     }
 448
 449   // Keep on tracing via gluons until reached end of leg.
 450   } while (hasFound && indxCol > 0 && loop < loopMax);
 451
 452   // Something went wrong in colour tracing.
 453   if (!hasFound || loop == loopMax) {
 454     infoPtr->errorMsg("Error in HadronLevel::traceFromAcol: "
 455       "colour tracing failed");
 456     return false;
 457   }
 458
 459   // Done.
 460   return true;
 461
 462 }
 463
 464 //*********
 465
 466 // Trace a colour loop, from a colour back to the anticolour of the same.
 467
 468 bool HadronLevel::traceInLoop(int indxCol, int indxAcol, Event& event) {
 469
 470   // Move around until back where begun.
 471   int loop = 0;
 472   int loopMax = iColAndAcol.size() + 2;
 473   bool hasFound = false;
 474   do {
 475     ++loop;
 476     hasFound= false;
 477
 478     // Check list of gluons.
 479     for (int i = 0; i < int(iColAndAcol.size()); ++i)
 480       if (event[ iColAndAcol[i] ].acol() == indxCol) {
 481         iParton.push_back( iColAndAcol[i] );
 482         indxCol = event[ iColAndAcol[i] ].col();
 483         iColAndAcol[i] = iColAndAcol.back();
 484         iColAndAcol.pop_back();
 485         hasFound = true;
 486         break;
 487     }
 488   } while (hasFound && indxCol != indxAcol && loop < loopMax);
 489
 490   // Something went wrong in colour tracing.
 491   if (!hasFound || loop == loopMax) {
 492     infoPtr->errorMsg("Error in HadronLevel::traceInLoop: "
 493       "colour tracing failed");
 494     return false;
 495   }
 496
 497   // Done.
 498   return true;
 499
 500 }
 501
 502 //*********
 503
 504 // Split junction-antijunction system into two, or simplify other way.
 505
 506 bool HadronLevel::splitJunctionPair(Event& event) {
 507
 508   // Construct separate index arrays for the three junction legs.
 509   int identJun = (-iPartonJun[0])/10;
 510   iJunLegA.resize(0);
 511   iJunLegB.resize(0);
 512   iJunLegC.resize(0);
 513   int leg = -1;
 514   for (int i = 0; i < int(iPartonJun.size()); ++ i) {
 515     if ( (-iPartonJun[i])/10 == identJun) ++leg;
 516     if (leg == 0) iJunLegA.push_back( iPartonJun[i] );
 517     else if (leg == 1) iJunLegB.push_back( iPartonJun[i] );
 518     else iJunLegC.push_back( iPartonJun[i] );
 519   }
 520
 521   // Construct separate index arrays for the three antijunction legs.
 522   int identAnti = (-iPartonAntiJun[0])/10;
 523   iAntiLegA.resize(0);
 524   iAntiLegB.resize(0);
 525   iAntiLegC.resize(0);
 526   leg = -1;
 527   for (int i = 0; i < int(iPartonAntiJun.size()); ++ i) {
 528     if ( (-iPartonAntiJun[i])/10 == identAnti) ++leg;
 529     if (leg == 0) iAntiLegA.push_back( iPartonAntiJun[i] );
 530     else if (leg == 1) iAntiLegB.push_back( iPartonAntiJun[i] );
 531     else iAntiLegC.push_back( iPartonAntiJun[i] );
 532   }
 533
 534   // Find interjunction legs, i.e. between junction and antijunction.
 535   int nMatch = 0;
 536   int legJun[3], legAnti[3], nGluLeg[3];
 537   if (iJunLegA.back() < 0) { legJun[nMatch] = 0;
 538     legAnti[nMatch] = (-iJunLegA.back())%10; ++nMatch;}
 539   if (iJunLegB.back() < 0) { legJun[nMatch] = 1;
 540     legAnti[nMatch] = (-iJunLegB.back())%10; ++nMatch;}
 541   if (iJunLegC.back() < 0) { legJun[nMatch] = 2;
 542     legAnti[nMatch] = (-iJunLegC.back())%10; ++nMatch;}
 543
 544   // Loop over interjunction legs.
 545   for (int iMatch = 0; iMatch < nMatch; ++iMatch) {
 546     vector<int>& iJunLeg = (legJun[iMatch] == 0) ? iJunLegA
 547       : ( (legJun[iMatch] == 1) ? iJunLegB : iJunLegC );
 548     vector<int>& iAntiLeg = (legAnti[iMatch] == 0) ? iAntiLegA
 549       : ( (legAnti[iMatch] == 1) ? iAntiLegB : iAntiLegC );
 550
 551     // Find number of gluons on each. Do nothing for now if none.
 552     nGluLeg[iMatch] = iJunLeg.size() + iAntiLeg.size() - 4;
 553     if (nGluLeg[iMatch] == 0) continue;
 554
 555     // Else pick up the gluons on the interjunction leg in order.
 556     iGluLeg.resize(0);
 557     for (int i = 1; i < int(iJunLeg.size()) - 1; ++i)
 558       iGluLeg.push_back( iJunLeg[i] );
 559     for (int i = int(iAntiLeg.size()) - 2; i > 0; --i)
 560       iGluLeg.push_back( iAntiLeg[i] );
 561
 562     // Remove those gluons from the junction/antijunction leg lists.
 563     iJunLeg.resize(1);
 564     iAntiLeg.resize(1);
 565
 566    // Pick a new quark at random; for simplicity no diquarks.
 567     int idQ = flavSel.pickLightQ();
 568     int colQ, acolQ;
 569
 570     // If one gluon on leg, split it into a collinear q-qbar pair.
 571     if (iGluLeg.size() == 1) {
 572
 573       // Store the new q qbar pair, sharing gluon colour and momentum.
 574       colQ = event[ iGluLeg[0] ].col();
 575       acolQ = event[ iGluLeg[0] ].acol();
 576       Vec4 pQ = 0.5 * event[ iGluLeg[0] ].p();
 577       double mQ = 0.5 * event[ iGluLeg[0] ].m();
 578       int iQ = event.append( idQ, 75, iGluLeg[0], 0, 0, 0, colQ, 0, pQ, mQ );
 579       int iQbar = event.append( -idQ, 75, iGluLeg[0], 0, 0, 0, 0, acolQ,
 580         pQ, mQ );
 581
 582       // Mark split gluon and update junction and antijunction legs.
 583       event[ iGluLeg[0] ].statusNeg();
 584       event[ iGluLeg[0] ].daughters( iQ, iQbar);
 585       iJunLeg.push_back(iQ);
 586       iAntiLeg.push_back(iQbar);
 587
 588     // If several gluons on the string, decide which g-g region to split up.
 589     } else {
 590
 591       // Evaluate mass-squared for all adjacent gluon pairs.
 592       m2Pair.resize(0);
 593       double m2Sum = 0.;
 594       for (int i = 0; i < int(iGluLeg.size()) - 1; ++i) {
 595         double m2Now = 0.5 * event[ iGluLeg[i] ].p()
 596           * event[ iGluLeg[i + 1] ].p();
 597         m2Pair.push_back(m2Now);
 598         m2Sum += m2Now;
 599       }
 600
 601       // Pick breakup region with probability proportional to mass-squared.
 602       double m2Reg = m2Sum * Rndm::flat();
 603       int iReg = -1;
 604       do m2Reg -= m2Pair[++iReg];
 605       while (m2Reg > 0. && iReg < int(iGluLeg.size()) - 1);
 606       m2Reg = m2Pair[iReg];
 607
 608       // Pick breaking point of string in chosen region (symmetrically).
 609       double m2Temp = min( JJSTRINGM2MAX, JJSTRINGM2FRAC * m2Reg);
 610       double xPos = 0.5;
 611       double xNeg = 0.5;
 612       do {
 613         double zTemp = zSel.zFrag( idQ, 0, m2Temp);
 614         xPos = 1. - zTemp;
 615         xNeg = m2Temp / (zTemp * m2Reg);
 616       } while (xNeg > 1.);
 617       if (Rndm::flat() > 0.5) swap(xPos, xNeg);
 618
 619       // Pick up two "mother" gluons of breakup. Mark them decayed.
 620       Particle& gJun = event[ iGluLeg[iReg] ];
 621       Particle& gAnti = event[ iGluLeg[iReg + 1] ];
 622       gJun.statusNeg();
 623       gAnti.statusNeg();
 624       int dau1 = event.size();
 625       gJun.daughters(dau1, dau1 + 3);
 626       gAnti.daughters(dau1, dau1 + 3);
 627       int mother1 = min( iGluLeg[iReg], iGluLeg[iReg + 1]);
 628       int mother2 = max( iGluLeg[iReg], iGluLeg[iReg + 1]);
 629
 630       // Can keep one of old colours but need one new so unambiguous.
 631       colQ = gJun.acol();
 632       acolQ = event.nextColTag();
 633
 634       // Store copied gluons with reduced momenta.
 635       int iGjun = event.append( 21, 75, mother1, mother2, 0, 0,
 636         gJun.col(), gJun.acol(), (1. - 0.5 * xPos) * gJun.p(),
 637         (1. - 0.5 * xPos) * gJun.m());
 638       int iGanti = event.append( 21, 75, mother1, mother2, 0, 0,
 639         acolQ, gAnti.acol(), (1. - 0.5 * xNeg) * gAnti.p(),
 640         (1. - 0.5 * xNeg) * gAnti.m());
 641
 642       // Store the new q qbar pair with remaining momenta.
 643       int iQ = event.append( idQ, 75, mother1, mother2, 0, 0,
 644         colQ, 0, 0.5 * xNeg * gAnti.p(), 0.5 * xNeg * gAnti.m() );
 645       int iQbar = event.append( -idQ, 75, mother1, mother2, 0, 0,
 646         0, acolQ, 0.5 * xPos * gJun.p(), 0.5 * xPos * gJun.m() );
 647
 648       // Update junction and antijunction legs with gluons and quarks.
 649       for (int i = 0; i < iReg; ++i)
 650         iJunLeg.push_back( iGluLeg[i] );
 651       iJunLeg.push_back(iGjun);
 652       iJunLeg.push_back(iQ);
 653       for (int i = int(iGluLeg.size()) - 1; i > iReg + 1; --i)
 654         iAntiLeg.push_back( iGluLeg[i] );
 655       iAntiLeg.push_back(iGanti);
 656       iAntiLeg.push_back(iQbar);
 657     }
 658
 659     // Update end colours for both g -> q qbar and g g -> g g q qbar.
 660     event.endColJunction(identJun - 1, legJun[iMatch], colQ);
 661     event.endColJunction(identAnti - 1, legAnti[iMatch], acolQ);
 662   }
 663
 664   // Update list of interjunction legs after splittings above.
 665   int iMatchUp = 0;
 666   while (iMatchUp < nMatch) {
 667     if (nGluLeg[iMatchUp] > 0) {
 668       for (int i = iMatchUp; i < nMatch - 1; ++i) {
 669         legJun[i] = legJun[i + 1];
 670         legAnti[i] = legAnti[i + 1];
 671         nGluLeg[i] = nGluLeg[i + 1];
 672       } --nMatch;
 673     } else ++iMatchUp;
 674   }
 675
 676   // Should not ever have three empty interjunction legs.
 677   if (nMatch == 3) {
 678     infoPtr->errorMsg("Error in HadronLevel::splitJunctionPair: "
 679       "three empty junction-junction legs");
 680     return false;
 681   }
 682
 683   // If two legs are empty, then collapse system to a single string.
 684   if (nMatch == 2) {
 685     int legJunLeft = 3 - legJun[0] - legJun[1];
 686     int legAntiLeft = 3 - legAnti[0] - legAnti[1];
 687     vector<int>& iJunLeg = (legJunLeft == 0) ? iJunLegA
 688       : ( (legJunLeft == 1) ? iJunLegB : iJunLegC );
 689     vector<int>& iAntiLeg = (legAntiLeft == 0) ? iAntiLegA
 690       : ( (legAntiLeft == 1) ? iAntiLegB : iAntiLegC );
 691     iPartonJun.resize(0);
 692     for (int i = int(iJunLeg.size()) - 1; i > 0; --i)
 693       iPartonJun.push_back( iJunLeg[i] );
 694     for (int i = 1; i < int(iAntiLeg.size()); ++i)
 695       iPartonJun.push_back( iAntiLeg[i] );
 696
 697     // Other string system empty. Remove junctions from their list. Done.
 698     iPartonAntiJun.resize(0);
 699     event.eraseJunction( max(identJun, identAnti) - 1);
 700     event.eraseJunction( min(identJun, identAnti) - 1);
 701     return true;
 702   }
 703
 704   // If one leg is empty then, depending on string length, either
 705   // (a) annihilate junction and antijunction into two simple strings, or
 706   // (b) split the empty leg by borrowing energy from nearby legs.
 707   if (nMatch == 1) {
 708
 709     // Identify the two external legs of either junction.
 710     vector<int>& iJunLeg0 = (legJun[0] == 0) ? iJunLegB : iJunLegA;
 711     vector<int>& iJunLeg1 = (legJun[0] == 2) ? iJunLegB : iJunLegC;
 712     vector<int>& iAntiLeg0 = (legAnti[0] == 0) ? iAntiLegB : iAntiLegA;
 713     vector<int>& iAntiLeg1 = (legAnti[0] == 2) ? iAntiLegB : iAntiLegC;
 714
 715     // Simplified procedure: mainly study first parton on each leg.
 716     Vec4 pJunLeg0 = event[ iJunLeg0[1] ].p();
 717     Vec4 pJunLeg1 = event[ iJunLeg1[1] ].p();
 718     Vec4 pAntiLeg0 = event[ iAntiLeg0[1] ].p();
 719     Vec4 pAntiLeg1 = event[ iAntiLeg1[1] ].p();
 720
 721     // Starting frame hopefully intermediate to two junction directions.
 722     Vec4 pStart = pJunLeg0 / pJunLeg0.e() + pJunLeg1 / pJunLeg1.e()
 723       + pAntiLeg0 / pAntiLeg0.e() + pAntiLeg1 / pAntiLeg1.e();
 724
 725     // Loop over iteration to junction/antijunction rest frames (JRF/ARF).
 726     RotBstMatrix MtoJRF, MtoARF;
 727     Vec4 pInJRF[3], pInARF[3];
 728     for (int iJun = 0; iJun < 2; ++iJun) {
 729       int offset = (iJun == 0) ? 0 : 2;
 730
 731       // Iterate from system rest frame towards the junction rest frame.
 732       RotBstMatrix MtoRF, Mstep;
 733       MtoRF.bstback(pStart);
 734       Vec4 pInRF[4];
 735       int iter = 0;
 736       do {
 737         ++iter;
 738
 739         // Find rest-frame momenta on the three sides of the junction.
 740         // Only consider first parton on each leg, for simplicity.
 741         pInRF[0 + offset] = pJunLeg0;
 742         pInRF[1 + offset] = pJunLeg1;
 743         pInRF[2 - offset] = pAntiLeg0;
 744         pInRF[3 - offset] = pAntiLeg1;
 745         for (int i = 0; i < 4; ++i) pInRF[i].rotbst(MtoRF);
 746
 747         // For third side add both legs beyond other junction, weighted.
 748         double wt2 = 1. - exp( -pInRF[2].e() / eNormJunction);
 749         double wt3 = 1. - exp( -pInRF[3].e() / eNormJunction);
 750         pInRF[2] = wt2 * pInRF[2] + wt3 * pInRF[3];
 751
 752         // Find new junction rest frame from the set of momenta.
 753         Mstep = stringFrag.junctionRestFrame( pInRF[0], pInRF[1], pInRF[2]);
 754         MtoRF.rotbst( Mstep );
 755       } while (iter < 3 || (Mstep.deviation() > CONVJNREST
 756         && iter < NTRYJNREST) );
 757
 758       // Store final boost and rest-frame (weighted) momenta.
 759       if (iJun == 0) {
 760         MtoJRF = MtoRF;
 761         for (int i = 0; i < 3; ++i) pInJRF[i] = pInRF[i];
 762       } else {
 763         MtoARF = MtoRF;
 764         for (int i = 0; i < 3; ++i) pInARF[i] = pInRF[i];
 765       }
 766     }
 767
 768     // Opposite operations: boost from JRF/ARF to original system.
 769     RotBstMatrix MfromJRF = MtoJRF;
 770     MfromJRF.invert();
 771     RotBstMatrix MfromARF = MtoARF;
 772     MfromARF.invert();
 773
 774     // Velocity vectors of junctions and momentum of legs in lab frame.
 775     Vec4 vJun(0., 0., 0., 1.);
 776     vJun.rotbst(MfromJRF);
 777     Vec4 vAnti(0., 0., 0., 1.);
 778     vAnti.rotbst(MfromARF);
 779     Vec4 pLabJ[3], pLabA[3];
 780     for (int i = 0; i < 3; ++i) {
 781       pLabJ[i] = pInJRF[i];
 782       pLabJ[i].rotbst(MfromJRF);
 783       pLabA[i] = pInARF[i];
 784       pLabA[i].rotbst(MfromARF);
 785     }
 786
 787     // Calculate Lambda-measure length of three possible topologies.
 788     double vJvA = vJun * vAnti;
 789     double vJvAe2y = vJvA + sqrt(vJvA*vJvA - 1.);
 790     double LambdaJA = (2. * pInJRF[0].e()) * (2. * pInJRF[1].e())
 791       * (2. * pInARF[0].e()) * (2. * pInARF[1].e()) * vJvAe2y;
 792     double Lambda00 = (2. * pLabJ[0] * pLabA[0])
 793       * (2. * pLabJ[1] * pLabA[1]);
 794     double Lambda01 = (2. * pLabJ[0] * pLabA[1])
 795       * (2. * pLabJ[1] * pLabA[0]);
 796
 797     // Case when either topology without junctions is the shorter one.
 798     if (LambdaJA > min( Lambda00, Lambda01)) {
 799       vector<int>& iAntiMatch0 = (Lambda00 < Lambda01)
 800         ? iAntiLeg0 : iAntiLeg1;
 801       vector<int>& iAntiMatch1 = (Lambda00 < Lambda01)
 802         ? iAntiLeg1 : iAntiLeg0;
 803
 804       // Define two quark-antiquark strings.
 805       iPartonJun.resize(0);
 806       for (int i = int(iJunLeg0.size()) - 1; i > 0; --i)
 807         iPartonJun.push_back( iJunLeg0[i] );
 808       for (int i = 1; i < int(iAntiMatch0.size()); ++i)
 809         iPartonJun.push_back( iAntiMatch0[i] );
 810       iPartonAntiJun.resize(0);
 811       for (int i = int(iJunLeg1.size()) - 1; i > 0; --i)
 812         iPartonAntiJun.push_back( iJunLeg1[i] );
 813       for (int i = 1; i < int(iAntiMatch1.size()); ++i)
 814         iPartonAntiJun.push_back( iAntiMatch1[i] );
 815
 816       // Remove junctions from their list. Done.
 817       event.eraseJunction( max(identJun, identAnti) - 1);
 818       event.eraseJunction( min(identJun, identAnti) - 1);
 819       return true;
 820     }
 821
 822     // Case where junction and antijunction to be separated.
 823     // Shuffle (p+/p-)  momentum of order <mThad> between systems,
 824     // times 2/3 for 120 degree in JRF, times 1/2 for two legs,
 825     // but not more than half of what nearest parton carries.
 826     double eShift = MTHAD / (3. * sqrt(vJvAe2y));
 827     double fracJ0 = min(0.5, eShift / pInJRF[0].e());
 828     double fracJ1 = min(0.5, eShift / pInJRF[0].e());
 829     Vec4 pFromJun = fracJ0 * pJunLeg0 + fracJ1 * pJunLeg1;
 830     double fracA0 = min(0.5, eShift / pInARF[0].e());
 831     double fracA1 = min(0.5, eShift / pInARF[0].e());
 832     Vec4 pFromAnti = fracA0 * pAntiLeg0 + fracA1 * pAntiLeg1;
 833
 834     // Copy partons with scaled-down momenta and update legs.
 835     int iNew = event.copy(iJunLeg0[1], 76);
 836     event[iNew].rescale5(1. - fracJ0);
 837     iJunLeg0[1] = iNew;
 838     iNew = event.copy(iJunLeg1[1], 76);
 839     event[iNew].rescale5(1. - fracJ1);
 840     iJunLeg1[1] = iNew;
 841     iNew = event.copy(iAntiLeg0[1], 76);
 842     event[iNew].rescale5(1. - fracA0);
 843     iAntiLeg0[1] = iNew;
 844     iNew = event.copy(iAntiLeg1[1], 76);
 845     event[iNew].rescale5(1. - fracA1);
 846     iAntiLeg1[1] = iNew;
 847
 848    // Pick a new quark at random; for simplicity no diquarks.
 849     int idQ = flavSel.pickLightQ();
 850
 851     // Update junction colours for new quark and antiquark.
 852     int colQ = event.nextColTag();
 853     int acolQ = event.nextColTag();
 854     event.endColJunction(identJun - 1, legJun[0], colQ);
 855     event.endColJunction(identAnti - 1, legAnti[0], acolQ);
 856
 857     // Store new quark and antiquark with momentum from other junction.
 858     int mother1 = min(iJunLeg0[1], iJunLeg1[1]);
 859     int mother2 = max(iJunLeg0[1], iJunLeg1[1]);
 860     int iNewJ = event.append( idQ, 76, mother1, mother2, 0, 0,
 861       colQ, 0, pFromAnti, pFromAnti.mCalc() );
 862     mother1 = min(iAntiLeg0[1], iAntiLeg1[1]);
 863     mother2 = max(iAntiLeg0[1], iAntiLeg1[1]);
 864     int iNewA = event.append( -idQ, 76, mother1, mother2, 0, 0,
 865       0, acolQ, pFromJun, pFromJun.mCalc() );
 866
 867     // Bookkeep new quark and antiquark on third legs.
 868     if (legJun[0] == 0) iJunLegA[1] = iNewJ;
 869     else if (legJun[0] == 1) iJunLegB[1] = iNewJ;
 870     else iJunLegC[1] = iNewJ;
 871     if (legAnti[0] == 0) iAntiLegA[1] = iNewA;
 872     else if (legAnti[0] == 1) iAntiLegB[1] = iNewA;
 873     else iAntiLegC[1] = iNewA;
 874
 875   // Done with splitting junction from antijunction.
 876   }
 877
 878   // Put together new junction parton list.
 879   iPartonJun.resize(0);
 880   for (int i = 0; i < int(iJunLegA.size()); ++i)
 881     iPartonJun.push_back( iJunLegA[i] );
 882   for (int i = 0; i < int(iJunLegB.size()); ++i)
 883     iPartonJun.push_back( iJunLegB[i] );
 884   for (int i = 0; i < int(iJunLegC.size()); ++i)
 885     iPartonJun.push_back( iJunLegC[i] );
 886
 887   // Put together new antijunction parton list.
 888   iPartonAntiJun.resize(0);
 889   for (int i = 0; i < int(iAntiLegA.size()); ++i)
 890     iPartonAntiJun.push_back( iAntiLegA[i] );
 891   for (int i = 0; i < int(iAntiLegB.size()); ++i)
 892     iPartonAntiJun.push_back( iAntiLegB[i] );
 893   for (int i = 0; i < int(iAntiLegC.size()); ++i)
 894     iPartonAntiJun.push_back( iAntiLegC[i] );
 895
 896   // Now the two junction systems are separated and can be stored.
 897   return true;
 898
 899 }
 900
 901 //**************************************************************************
 902
 903 } // end namespace Pythia8
 904