PYTHIA8/pythia8170/src/RHadrons.cxx

   1 // RHadrons.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 RHadrons class.
   7
   8 #include "RHadrons.h"
   9
  10 namespace Pythia8 {
  11
  12 //==========================================================================
  13
  14 // The RHadrons 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 const int RHadrons::IDRHADSB[14] = {  1000512, 1000522, 1000532,
  22   1000542, 1000552, 1005113, 1005211, 1005213, 1005223, 1005311,
  23   1005313, 1005321, 1005323, 1005333 };
  24
  25 const int RHadrons::IDRHADST[14] = {  1000612, 1000622, 1000632,
  26   1000642, 1000652, 1006113, 1006211, 1006213, 1006223, 1006311,
  27   1006313, 1006321, 1006323, 1006333 };
  28
  29 // Gluino code and list of gluino R-hadron codes.
  30 const int RHadrons::IDRHADGO[38] = {  1000993, 1009113, 1009213,
  31   1009223, 1009313, 1009323, 1009333, 1009413, 1009423, 1009433,
  32   1009443, 1009513, 1009523, 1009533, 1009543, 1009553, 1091114,
  33   1092114, 1092214, 1092224, 1093114, 1093214, 1093224, 1093314,
  34   1093324, 1093334, 1094114, 1094214, 1094224, 1094314, 1094324,
  35   1094334, 1095114, 1095214, 1095224, 1095314, 1095324, 1095334 };
  36
  37 // Allow a few tries for flavour selection since failure is allowed.
  38 const int RHadrons::NTRYMAX = 10;
  39
  40 // Safety margin (in GeV) when constructing kinematics of system.
  41 const double RHadrons::MSAFETY = 0.1;
  42
  43 // Maximal energy to borrow for gluon to insert on leg in to junction.
  44 const double RHadrons::EGBORROWMAX = 4.;
  45
  46 //--------------------------------------------------------------------------
  47
  48 // Main routine to initialize R-hadron handling.
  49
  50 bool RHadrons::init( Info* infoPtrIn, Settings& settings,
  51   ParticleData* particleDataPtrIn, Rndm* rndmPtrIn) {
  52
  53   // Store input pointers for future use.
  54   infoPtr          = infoPtrIn;
  55   particleDataPtr  = particleDataPtrIn;
  56   rndmPtr          = rndmPtrIn;
  57
  58   // Flags and parameters related to R-hadron formation and decay.
  59   allowRH          = settings.flag("RHadrons:allow");
  60   maxWidthRH       = settings.parm("RHadrons:maxWidth");
  61   idRSb            = settings.mode("RHadrons:idSbottom");
  62   idRSt            = settings.mode("RHadrons:idStop");
  63   idRGo            = settings.mode("RHadrons:idGluino");
  64   setMassesRH      = settings.flag("RHadrons:setMasses");
  65   probGluinoballRH = settings.parm("RHadrons:probGluinoball");
  66   mOffsetCloudRH   = settings.parm("RHadrons:mOffsetCloud");
  67   mCollapseRH      = settings.parm("RHadrons:mCollapse");
  68   diquarkSpin1RH   = settings.parm("RHadrons:diquarkSpin1");
  69
  70   // Check whether sbottom, stop or gluino should form R-hadrons.
  71   allowRSb         = allowRH && idRSb > 0
  72     && (particleDataPtr->mWidth(idRSb) < maxWidthRH);
  73   allowRSt         = allowRH && idRSt > 0
  74     && (particleDataPtr->mWidth(idRSt) < maxWidthRH);
  75   allowRGo         = allowRH && idRGo > 0
  76     && (particleDataPtr->mWidth(idRGo) < maxWidthRH);
  77   allowSomeR       = allowRSb || allowRSt || allowRGo;
  78
  79   // Set nominal masses of sbottom R-mesons and R-baryons.
  80   if (allowRSb) {
  81     m0Sb = particleDataPtr->m0(idRSb);
  82     if (setMassesRH) {
  83       for (int i = 0; i < 14; ++i) {
  84         int idR = IDRHADSB[i];
  85         double m0RHad = m0Sb + mOffsetCloudRH;
  86         m0RHad += particleDataPtr->constituentMass( (idR%100)/10);
  87         if (i > 4)
  88         m0RHad += particleDataPtr->constituentMass( (idR%1000)/100 );
  89         particleDataPtr->m0( idR, m0RHad);
  90       }
  91     }
  92
  93     // Set widths and lifetimes of sbottom R-hadrons.
  94     double mWidthRHad = particleDataPtr->mWidth(idRSb);
  95     double tau0RHad   = particleDataPtr->tau0(  idRSb);
  96     for (int i = 0; i < 14; ++i) {
  97       particleDataPtr->mWidth( IDRHADSB[i], mWidthRHad);
  98       particleDataPtr->tau0(   IDRHADSB[i],   tau0RHad);
  99     }
 100   }
 101
 102   // Set nominal masses of stop R-mesons and R-baryons.
 103   if (allowRSt) {
 104     m0St = particleDataPtr->m0(idRSt);
 105     if (setMassesRH) {
 106       for (int i = 0; i < 14; ++i) {
 107         int idR = IDRHADST[i];
 108         double m0RHad = m0St + mOffsetCloudRH;
 109         m0RHad += particleDataPtr->constituentMass( (idR%100)/10);
 110         if (i > 4)
 111         m0RHad += particleDataPtr->constituentMass( (idR%1000)/100 );
 112         particleDataPtr->m0( idR, m0RHad);
 113       }
 114     }
 115
 116     // Set widths and lifetimes of stop R-hadrons.
 117     double mWidthRHad = particleDataPtr->mWidth(idRSt);
 118     double tau0RHad   = particleDataPtr->tau0(  idRSt);
 119     for (int i = 0; i < 14; ++i) {
 120       particleDataPtr->mWidth( IDRHADST[i], mWidthRHad);
 121       particleDataPtr->tau0(   IDRHADST[i],   tau0RHad);
 122     }
 123   }
 124
 125   // Set nominal masses of gluino R-glueballs, R-mesons and R-baryons.
 126   if (allowRGo) {
 127     m0Go = particleDataPtr->m0(idRGo);
 128     if (setMassesRH) {
 129       particleDataPtr->m0( IDRHADGO[0], m0Go + 2. * mOffsetCloudRH
 130         + particleDataPtr->constituentMass(21) );
 131       for (int i = 1; i < 38; ++i) {
 132         int idR = IDRHADGO[i];
 133         double m0RHad = m0Go + 2. * mOffsetCloudRH;
 134         m0RHad += particleDataPtr->constituentMass( (idR%1000)/100 );
 135         m0RHad += particleDataPtr->constituentMass( (idR%100)/10);
 136         if (i > 15)
 137         m0RHad += particleDataPtr->constituentMass( (idR%10000)/1000 );
 138         particleDataPtr->m0( idR, m0RHad);
 139       }
 140     }
 141
 142     // Set widths and lifetimes of gluino R-hadrons.
 143     double mWidthRHad = particleDataPtr->mWidth(idRGo);
 144     double tau0RHad   = particleDataPtr->tau0(  idRGo);
 145     for (int i = 0; i < 38; ++i) {
 146       particleDataPtr->mWidth( IDRHADGO[i], mWidthRHad);
 147       particleDataPtr->tau0(   IDRHADGO[i],   tau0RHad);
 148     }
 149   }
 150
 151   // Done.
 152   return true;
 153
 154 }
 155 //--------------------------------------------------------------------------
 156
 157 // Check if a given particle can produce R-hadrons.
 158
 159 bool RHadrons::givesRHadron( int id) {
 160   if (allowRSb && abs(id) == idRSb) return true;
 161   if (allowRSt && abs(id) == idRSt) return true;
 162   if (allowRGo && id == idRGo) return true;
 163   return false;
 164
 165 }
 166
 167 //--------------------------------------------------------------------------
 168
 169 // Produce R-hadrons by fragmenting them off from existing strings.
 170
 171 bool RHadrons::produce( ColConfig& colConfig, Event& event) {
 172
 173   // Check whether some sparticles are allowed to hadronize.
 174   if (!allowSomeR) return true;
 175
 176   // Reset arrays and values.
 177   iBefRHad.resize(0);
 178   iCreRHad.resize(0);
 179   iRHadron.resize(0);
 180   iAftRHad.resize(0);
 181   isTriplet.resize(0);
 182   nRHad = 0;
 183
 184   // Loop over event and identify hadronizing sparticles.
 185   for (int i = 0; i < event.size(); ++i)
 186    if (event[i].isFinal() && givesRHadron(event[i].id())) {
 187     iBefRHad.push_back(i);
 188     iCreRHad.push_back(i);
 189     iRHadron.push_back(0);
 190     iAftRHad.push_back(0);
 191     isTriplet.push_back(true);
 192   }
 193   nRHad = iRHadron.size();
 194
 195   // Done if no hadronizing sparticles.
 196   if (nRHad == 0) return true;
 197
 198   // Max two R-hadrons. Randomize order of processing.
 199   if (nRHad > 2) {
 200      infoPtr->errorMsg("Error in RHadrons::produce: "
 201        "cannot handle more than two R-hadrons");
 202      return false;
 203   }
 204   if (nRHad > 1 && rndmPtr->flat() > 0.5) swap( iBefRHad[0], iBefRHad[1]);
 205
 206   // Split a system with both a sparticle and a junction.
 207   iBef = iBefRHad[0];
 208   iSys = colConfig.findSinglet( iBef);
 209   systemPtr = &colConfig[iSys];
 210   if (systemPtr->hasJunction && !splitOffJunction( colConfig, event)) {
 211     infoPtr->errorMsg("Error in RHadrons::produce: "
 212       "cannot handle system with junction");
 213     return false;
 214   }
 215   if (nRHad == 2) {
 216     iBef = iBefRHad[1];
 217     iSys = colConfig.findSinglet( iBefRHad[1]);
 218     systemPtr = &colConfig[iSys];
 219     if (systemPtr->hasJunction && !splitOffJunction( colConfig, event)) {
 220       infoPtr->errorMsg("Error in RHadrons::produce: "
 221         "cannot handle system with junction");
 222       return false;
 223     }
 224   }
 225
 226   // Open up a closed gluon/gluino loop.
 227   iBef = iBefRHad[0];
 228   iSys = colConfig.findSinglet( iBef);
 229   systemPtr = &colConfig[iSys];
 230   if (systemPtr->isClosed && !openClosedLoop( colConfig, event)) {
 231     infoPtr->errorMsg("Error in RHadrons::produce: "
 232       "cannot open up closed gluon/gluino loop");
 233     return false;
 234   }
 235   if (nRHad == 2) {
 236     iBef = iBefRHad[1];
 237     iSys = colConfig.findSinglet( iBefRHad[1]);
 238     systemPtr = &colConfig[iSys];
 239     if (systemPtr->isClosed && !openClosedLoop( colConfig, event)) {
 240       infoPtr->errorMsg("Error in RHadrons::produce: "
 241         "cannot open up closed gluon/gluino loop");
 242       return false;
 243     }
 244   }
 245
 246   // Split up a colour singlet system that should give two R-hadrons.
 247   if (nRHad == 2) {
 248     int iSys1 = colConfig.findSinglet( iBefRHad[0]);
 249     int iSys2 = colConfig.findSinglet( iBefRHad[1]);
 250     if (iSys2 == iSys1) {
 251       iSys = iSys1;
 252       systemPtr = &colConfig[iSys];
 253       if ( !splitSystem( colConfig, event) ) {
 254         infoPtr->errorMsg("Error in RHadrons::produce: "
 255           "failed to handle two sparticles in same system");
 256         return false;
 257       }
 258     }
 259   }
 260
 261   // Loop over R-hadrons to be formed. Find its colour singlet system.
 262   for (iRHad = 0; iRHad < nRHad; ++iRHad) {
 263     iBef = iBefRHad[iRHad];
 264     iSys = colConfig.findSinglet( iBef);
 265     if (iSys < 0) {
 266       infoPtr->errorMsg("Error in RHadrons::produce: "
 267         "sparticle not in any colour singlet");
 268       return false;
 269     }
 270     systemPtr = &colConfig[iSys];
 271
 272     // For now don't handle systems involving junctions or loops.
 273     if (systemPtr->hasJunction) {
 274       infoPtr->errorMsg("Error in RHadrons::produce: "
 275         "cannot handle system with junction");
 276       return false;
 277     }
 278     if (systemPtr->isClosed) {
 279       infoPtr->errorMsg("Error in RHadrons::produce: "
 280         "cannot handle closed colour loop");
 281       return false;
 282     }
 283
 284     // Handle formation of R-hadron separately for gluino and squark.
 285     if (event[iBef].id() == idRGo) isTriplet[iRHad] = false;
 286     bool formed = (isTriplet[iRHad]) ? produceSquark( colConfig, event)
 287                                      : produceGluino( colConfig, event);
 288     if (!formed) return false;
 289
 290   // End of loop over R-hadrons. Done.
 291   }
 292   return true;
 293
 294 }
 295
 296 //--------------------------------------------------------------------------
 297
 298 // Decay R-hadrons by resolving them into string systems and letting the
 299 // heavy unstable particle decay as normal.
 300
 301 bool RHadrons::decay( Event& event) {
 302
 303   // Loop over R-hadrons to decay.
 304   for (iRHad = 0; iRHad < nRHad; ++iRHad) {
 305     int    iRNow  = iRHadron[iRHad];
 306     int    iRBef  = iBefRHad[iRHad];
 307     int    idRHad = event[iRNow].id();
 308     double mRHad  = event[iRNow].m();
 309     double mRBef  = event[iRBef].m();
 310     int    iR0    = 0;
 311     int    iR2    = 0;
 312
 313     // Find flavour content of squark or gluino R-hadron.
 314     pair<int,int> idPair = (isTriplet[iRHad])
 315       ? fromIdWithSquark( idRHad) : fromIdWithGluino( idRHad);
 316     int id1 = idPair.first;
 317     int id2 = idPair.second;
 318
 319     // Sharing of momentum: the squark/gluino should be restored
 320     // to original mass, but error if negative-mass spectators.
 321     double fracR = mRBef / mRHad;
 322     if (fracR >= 1.) {
 323       infoPtr->errorMsg("Error in RHadrons::decay: "
 324           "too low R-hadron mass for decay");
 325       return false;
 326     }
 327
 328     // Squark: new colour needed in the breakup.
 329     if (isTriplet[iRHad]) {
 330       int colNew = event.nextColTag();
 331       int col    = (event[iRBef].col() != 0) ? colNew : 0;
 332       int acol   = (col == 0) ? colNew : 0;
 333
 334       // Store the constituents of a squark R-hadron.
 335       iR0 = event.append( id1, 106, iRNow, 0, 0, 0, col, acol,
 336         fracR * event[iRNow].p(), fracR * mRHad, 0.);
 337       iR2 = event.append( id2, 106, iRNow, 0, 0, 0, acol, col,
 338         (1. - fracR) * event[iRNow].p(), (1. - fracR) * mRHad, 0.);
 339
 340     // Gluino: set mass sharing between two spectators.
 341     } else {
 342       double m1Eff  = particleDataPtr->constituentMass(id1) + mOffsetCloudRH;
 343       double m2Eff  = particleDataPtr->constituentMass(id2) + mOffsetCloudRH;
 344       double frac1 = (1. - fracR) * m1Eff / ( m1Eff + m2Eff);
 345       double frac2 = (1. - fracR) * m2Eff / ( m1Eff + m2Eff);
 346
 347       // Two new colours needed in the breakups.
 348       int col1 = event.nextColTag();
 349       int col2 = event.nextColTag();
 350
 351       // Store the constituents of a gluino R-hadron.
 352       iR0 = event.append( idRGo, 106, iRNow, 0, 0, 0, col2, col1,
 353         fracR * event[iRNow].p(), fracR * mRHad, 0.);
 354       event.append( id1, 106, iRNow, 0, 0, 0, col1, 0,
 355         frac1 * event[iRNow].p(), frac1 * mRHad, 0.);
 356       iR2 = event.append( id2, 106, iRNow, 0, 0, 0, 0, col2,
 357         frac2 * event[iRNow].p(), frac2 * mRHad, 0.);
 358     }
 359
 360     // Mark R-hadron as decayed and update history.
 361     event[iRNow].statusNeg();
 362     event[iRNow].daughters( iR0, iR2);
 363     iAftRHad[iRHad] = iR0;
 364
 365     // Set secondary vertex for decay products, but no lifetime.
 366     Vec4 vDec = event[iRNow].vProd() + event[iRNow].tau()
 367               * event[iR0].p() / event[iR0].m();
 368     for (int iRd = iR0; iRd <= iR2; ++iRd) event[iRd].vProd( vDec);
 369
 370   // End loop over R-hadron decays, based on velocity of squark.
 371
 372   }
 373
 374   // Done.
 375   return true;
 376
 377 }
 378
 379 //--------------------------------------------------------------------------
 380
 381 // Split a system that contains both a sparticle and a junction.
 382
 383 bool RHadrons::splitOffJunction( ColConfig& colConfig, Event& event) {
 384
 385   // Classify system into three legs, and find sparticle location.
 386   vector<int> leg1, leg2, leg3;
 387   int iLegSP = 0;
 388   int iPosSP = 0;
 389   int iLeg = 0;
 390   int iPos = 0;
 391   for (int i = 0; i < systemPtr->size(); ++i) {
 392     ++iPos;
 393     int iP = systemPtr->iParton[i];
 394     if (iP < 0) {
 395       ++iLeg;
 396       iPos = -1;
 397     } else if (iLeg == 1) leg1.push_back( iP);
 398     else if   (iLeg == 2) leg2.push_back( iP);
 399     else if   (iLeg == 3) leg3.push_back( iP);
 400     if (iP == iBef) {
 401       iLegSP = iLeg;
 402       iPosSP = iPos;
 403     }
 404   }
 405   if (iLegSP == 0) return false;
 406
 407   // Swap so leg 1 contains sparticle. If not innermost sparticle then
 408   // skip for now, and wait for this other (gluino!) to be split off.
 409   if      (iLegSP == 2) swap(leg2, leg1);
 410   else if (iLegSP == 3) swap(leg3, leg1);
 411   for (int i = 0; i < iPosSP; ++i)
 412     if (event[leg1[i]].id() != 21) return true;
 413
 414   // No gluon between sparticle and junction: find kinetic energy of system.
 415   if (iPosSP == 0) {
 416     Vec4 pSP  = event[iBef].p();
 417     Vec4 pRec = 0.;
 418     for (int i = 0; i < int(leg2.size()); ++i) pRec += event[leg2[i]].p();
 419     for (int i = 0; i < int(leg3.size()); ++i) pRec += event[leg3[i]].p();
 420     double mSys  = (pSP + pRec).mCalc();
 421     double mSP   = pSP.mCalc();
 422     double mRec  = pRec.mCalc();
 423     double eKin  = mSys - mSP - mRec;
 424
 425     // Insert new gluon that borrows part of kinetic energy.
 426     double mNewG  = EGBORROWMAX * eKin / (EGBORROWMAX + eKin) ;
 427     Vec4   pNewG  = (mNewG / mSys) * (pSP + pRec);
 428     int    newCol = event.nextColTag();
 429     bool   isCol  = (event[leg1.back()].col() > 0);
 430     int    colNG  = (isCol)? newCol :  event[iBef].acol();
 431     int    acolNG = (isCol) ? event[iBef].col() : newCol;
 432     int    iNewG  = event.append( 21, 101, iBef, 0, 0, 0, colNG, acolNG,
 433       pNewG, mNewG, 0.);
 434     leg1.insert( leg1.begin(), iNewG);
 435     ++iPosSP;
 436
 437     // Boosts for remainder systems that gave up energy.
 438     double mNewSys = mSys - mNewG;
 439     double pAbsOld = 0.5 * sqrtpos( pow2(mSys*mSys - mSP*mSP - mRec*mRec)
 440                    - pow2(2. * mSP * mRec) ) / mSys;
 441     double pAbsNew = 0.5 * sqrtpos( pow2(mNewSys*mNewSys - mSP*mSP - mRec*mRec)
 442                    - pow2(2. * mSP * mRec) ) / mNewSys;
 443     RotBstMatrix MtoCM, MfromCM, MSP, MRec;
 444     MtoCM.toCMframe( pSP, pRec);
 445     MfromCM = MtoCM;
 446     MfromCM.invert();
 447     MSP = MtoCM;
 448     MSP.bst( 0., 0., -pAbsOld / sqrt(mSP * mSP + pAbsOld * pAbsOld));
 449     MSP.bst( 0., 0.,  pAbsNew / sqrt(mSP * mSP + pAbsNew * pAbsNew));
 450     MSP.rotbst( MfromCM);
 451     MRec = MtoCM;
 452     MRec.bst( 0., 0.,  pAbsOld / sqrt(mRec * mRec + pAbsOld * pAbsOld));
 453     MRec.bst( 0., 0., -pAbsNew / sqrt(mRec * mRec + pAbsNew * pAbsNew));
 454     MRec.rotbst( MfromCM);
 455
 456     // Copy down recoling partons and boost their momenta.
 457     int iNewSP  = event.copy( iBef, 101);
 458     event[iNewSP].mother2(0);
 459     event[iBef].daughter1(iNewG);
 460     event[iNewSP].rotbst( MSP);
 461     leg1[iPosSP]   = iNewSP;
 462     if (iBefRHad[0] == iBef) iBefRHad[0] = iNewSP;
 463     else if (nRHad > 1 && iBefRHad[1] == iBef) iBefRHad[1] = iNewSP;
 464     iBef = iNewSP;
 465     for (int i = 0; i < int(leg2.size()); ++i) {
 466       int iCopy = event.copy( leg2[i], 101);
 467       event[iCopy].rotbst( MRec);
 468       if (iBefRHad[0] == leg2[i]) iBefRHad[0] = iCopy;
 469       else if (nRHad > 1 && iBefRHad[1] == leg2[i]) iBefRHad[1] = iCopy;
 470       leg2[i] = iCopy;
 471     }
 472     for (int i = 0; i < int(leg3.size()); ++i) {
 473       int iCopy = event.copy( leg3[i], 101);
 474       event[iCopy].rotbst( MRec);
 475       if (iBefRHad[0] == leg3[i]) iBefRHad[0] = iCopy;
 476       else if (nRHad > 1 && iBefRHad[1] == leg3[i]) iBefRHad[1] = iCopy;
 477       leg3[i]   = iCopy;
 478     }
 479
 480   // Now always at least one gluon between sparticle and junction.
 481   }
 482
 483   // Find gluon with largest energy in sparticle rest frame.
 484   int    iGspl = 0;
 485   double eGspl = event[leg1[0]].p() * event[iBef].p();
 486   for (int i = 1; i < iPosSP; ++i) {
 487     double eGnow = event[leg1[i]].p() * event[iBef].p();
 488     if (eGnow > eGspl) {
 489       iGspl = i;
 490       eGspl = eGnow;
 491     }
 492   }
 493   int iG = leg1[iGspl];
 494
 495   // Split this gluon into a collinear quark.antiquark pair.
 496   int idNewQ = flavSelPtr->pickLightQ();
 497   int iNewQ  = event.append(  idNewQ, 101, iG, 0, 0, 0, event[iG].col(), 0,
 498     0.5 * event[iG].p(), 0.5 * event[iG].m(), 0.);
 499   int iNewQb = event.append( -idNewQ, 101, iG, 0, 0, 0, 0, event[iG].acol(),
 500     0.5 * event[iG].p(), 0.5 * event[iG].m(), 0.);
 501   event[iG].statusNeg();
 502   event[iG].daughters( iNewQ, iNewQb);
 503   if (event[leg1.back()].col() == 0) swap( iNewQ, iNewQb);
 504
 505   // Collect two new systems after split.
 506   vector<int> iNewSys1, iNewSys2;
 507   iNewSys1.push_back( iNewQb);
 508   for (int i = iGspl + 1; i < int(leg1.size()); ++i)
 509     iNewSys1.push_back( leg1[i]);
 510   iNewSys2.push_back( -10);
 511   for (int i = 0; i < iGspl; ++i) iNewSys2.push_back( leg1[i]);
 512   iNewSys2.push_back( iNewQ);
 513   iNewSys2.push_back( -11);
 514   for (int i = 0; i < int(leg2.size()); ++i) iNewSys2.push_back( leg2[i]);
 515   iNewSys2.push_back( -12);
 516   for (int i = 0; i < int(leg3.size()); ++i) iNewSys2.push_back( leg3[i]);
 517
 518   // Remove old system and insert two new ones.
 519   colConfig.erase(iSys);
 520   colConfig.insert( iNewSys1, event);
 521   colConfig.insert( iNewSys2, event);
 522
 523   // Done.
 524   return true;
 525
 526 }
 527
 528 //--------------------------------------------------------------------------
 529
 530 // Open up a closed gluon/gluino loop.
 531
 532 bool RHadrons::openClosedLoop( ColConfig& colConfig, Event& event) {
 533
 534   // Find gluon with largest energy in gluino rest frame.
 535   int    iGspl = -1;
 536   double eGspl = 0.;
 537   for (int i = 0; i < systemPtr->size(); ++i) {
 538     int  iGNow = systemPtr->iParton[i];
 539     if (event[iGNow].id() == 21) {
 540       double eGnow = event[iGNow].p() * event[iBef].p();
 541       if (eGnow > eGspl) {
 542         iGspl = i;
 543         eGspl = eGnow;
 544       }
 545     }
 546   }
 547   if (iGspl == -1) return false;
 548
 549   // Split this gluon into a collinear quark.antiquark pair.
 550   int iG     = systemPtr->iParton[iGspl];
 551   int idNewQ = flavSelPtr->pickLightQ();
 552   int iNewQ  = event.append(  idNewQ, 101, iG, 0, 0, 0, event[iG].col(), 0,
 553     0.5 * event[iG].p(), 0.5 * event[iG].m(), 0.);
 554   int iNewQb = event.append( -idNewQ, 101, iG, 0, 0, 0, 0, event[iG].acol(),
 555     0.5 * event[iG].p(), 0.5 * event[iG].m(), 0.);
 556   event[iG].statusNeg();
 557   event[iG].daughters( iNewQ, iNewQb);
 558
 559   // Order partons in new system. Note order of colour flow.
 560   int iNext = iGspl + 1;
 561   if (iNext == systemPtr->size()) iNext = 0;
 562   if (event[ systemPtr->iParton[iNext]].acol() != event[iNewQ].col())
 563     swap( iNewQ, iNewQb);
 564   vector<int> iNewSys;
 565   iNewSys.push_back( iNewQ);
 566   for (int i = iGspl + 1; i < systemPtr->size(); ++i)
 567     iNewSys.push_back( systemPtr->iParton[i]);
 568   for (int i = 0; i < iGspl; ++i)
 569     iNewSys.push_back( systemPtr->iParton[i]);
 570   iNewSys.push_back( iNewQb);
 571
 572   // Erase the old system and insert the new one instead.
 573   colConfig.erase(iSys);
 574   colConfig.insert( iNewSys, event);
 575
 576   // Done.
 577   return true;
 578
 579 }
 580
 581 //--------------------------------------------------------------------------
 582
 583 // Split a single colour singlet that contains two sparticles.
 584 // To fix : if nLeg >= 3 && mMin large handle as in nLeg == 1??
 585
 586 bool RHadrons::splitSystem( ColConfig& colConfig, Event& event) {
 587
 588   // First and second R-hadron mother. Number of legs in between.
 589   int iFirst  = -1;
 590   int iSecond = -1;
 591   for (int i = 0; i < int(systemPtr->size()); ++i) {
 592     int  iTmp   = systemPtr->iParton[i];
 593     if ( givesRHadron( event[iTmp].id()) ) {
 594       if (iFirst == -1) iFirst  = i;
 595       else              iSecond = i;
 596     }
 597   }
 598   int nLeg = iSecond - iFirst;
 599
 600   // New flavour pair for breaking the string, and its mass.
 601   int    idNewQ = flavSelPtr->pickLightQ();
 602   double mNewQ  = particleDataPtr->constituentMass( idNewQ);
 603   vector<int> iNewSys1, iNewSys2;
 604
 605   // If sparticles are neighbours then need new q-qbar pair in between.
 606   if (nLeg == 1) {
 607
 608     // The location of the two sparticles.
 609     int i1Old = systemPtr->iParton[iFirst];
 610     int i2Old = systemPtr->iParton[iSecond];
 611
 612     // Take a fraction of momentum to give breakup quark pair.
 613     double mHat = (event[i1Old].p() + event[i2Old].p()).mCalc();
 614     double mMax = mHat - event[i1Old].m() - event[i2Old].m();
 615     if (mMax < 2. * (mNewQ + MSAFETY)) return false;
 616     double mEff = min( 2. * (mNewQ + mOffsetCloudRH), mMax - 2. * MSAFETY);
 617     double frac = mEff / mHat;
 618     Vec4   pEff = frac * (event[i1Old].p() + event[i2Old].p());
 619
 620     // New kinematics by (1) go to same mHat with bigger masses, and
 621     // (2) rescale back down to original masses and reduced mHat.
 622     Vec4 p1New, p2New;
 623     if ( !newKin( event[i1Old].p(), event[i2Old].p(),
 624       event[i1Old].m() / (1. - frac), event[i2Old].m() / (1. - frac),
 625       p1New, p2New) ) return false;
 626     p1New *= 1. - frac;
 627     p2New *= 1. - frac;
 628
 629     // Fill in new partons after branching, with correct colour/flavour sign.
 630     int i1New, i2New, i3New, i4New;
 631     int newCol = event.nextColTag();
 632     i1New = event.copy( i1Old, 101);
 633     if (event[i2Old].acol() == event[i1Old].col()) {
 634       i3New = event.append( -idNewQ, 101, i1Old, 0, 0, 0,
 635         0, event[i2Old].acol(), 0.5 * pEff, 0.5 * mEff, 0.);
 636       i2New = event.copy( i2Old, 101);
 637       event[i2New].acol( newCol);
 638       i4New = event.append(  idNewQ, 101, i2Old, 0, 0, 0,
 639         newCol, 0, 0.5 * pEff, 0.5 * mEff, 0.);
 640     } else {
 641       i3New = event.append(  idNewQ, 101, i1Old, 0, 0, 0,
 642         event[i2Old].col(), 0, 0.5 * pEff, 0.5 * mEff, 0.);
 643       i2New = event.copy( i2Old, 101);
 644       event[i2New].col( newCol);
 645       i4New = event.append( -idNewQ, 101, i2Old, 0, 0, 0,
 646         0, newCol, 0.5 * pEff, 0.5 * mEff, 0.);
 647     }
 648
 649     // Modify momenta and history. For iBotCopyId tracing asymmetric
 650     // bookkeeping where one sparticle is radiator and the other recoiler.
 651     event[i1New].p( p1New);
 652     event[i2New].p( p2New);
 653     event[i1Old].daughters( i1New, i3New);
 654     event[i1New].mother2( 0);
 655     event[i2Old].daughters( i2New, i4New);
 656     event[i2New].mother2( 0);
 657     iBefRHad[0] = i1New;
 658     iBefRHad[1] = i2New;
 659
 660     // Partons in the two new systems.
 661     for (int i = 0; i < iFirst; ++i)
 662       iNewSys1.push_back( systemPtr->iParton[i]);
 663     iNewSys1.push_back( i1New);
 664     iNewSys1.push_back( i3New);
 665     iNewSys2.push_back( i4New);
 666     iNewSys2.push_back( i2New);
 667     for (int i = iSecond + 1; i < int(systemPtr->size()); ++i)
 668       iNewSys2.push_back( systemPtr->iParton[i]);
 669
 670   // If one gluon between sparticles then split it into a
 671   // collinear q-qbar pair.
 672   } else if (nLeg == 2) {
 673
 674     // Gluon to be split and its two daughters.
 675     int iOld  = systemPtr->iParton[iFirst + 1];
 676     int i1New = event.append(  idNewQ, 101, iOld, 0, 0, 0,
 677       event[iOld].col(), 0, 0.5 * event[iOld].p(),
 678       0.5 * event[iOld].m(), 0.);
 679     int i2New = event.append( -idNewQ, 101, iOld, 0, 0, 0,
 680       0, event[iOld].acol(), 0.5 * event[iOld].p(),
 681       0.5 * event[iOld].m(), 0.);
 682     event[iOld].statusNeg();
 683     event[iOld].daughters( i1New, i2New);
 684
 685     // Partons in the two new systems.
 686     if (event[systemPtr->iParton[iFirst]].col() == event[i2New].acol())
 687       swap( i1New, i2New);
 688     for (int i = 0; i <= iFirst; ++i)
 689       iNewSys1.push_back( systemPtr->iParton[i]);
 690     iNewSys1.push_back( i1New);
 691     iNewSys2.push_back( i2New);
 692     for (int i = iSecond; i < int(systemPtr->size()); ++i)
 693       iNewSys2.push_back( systemPtr->iParton[i]);
 694
 695   // If several gluons between sparticles then find lowest-mass gluon pair
 696   // and replace it by a q-qbar pair.
 697   } else {
 698
 699     // Find lowest-mass gluon pair and adjust effective quark mass.
 700     int    iMin  = 0;
 701     int    i1Old = 0;
 702     int    i2Old = 0;
 703     double mMin  = 1e20;
 704     for (int i = iFirst + 1; i < iSecond - 1; ++i) {
 705       int    i1Tmp = systemPtr->iParton[i];
 706       int    i2Tmp = systemPtr->iParton[i + 1];
 707       double mTmp  = (event[i1Tmp].p() + event[i2Tmp].p()).mCalc();
 708       if (mTmp < mMin) {
 709         iMin  = i;
 710         i1Old = i1Tmp;
 711         i2Old = i2Tmp;
 712         mMin  = mTmp;
 713       }
 714     }
 715     double mEff = min( mNewQ + mOffsetCloudRH, 0.4 * mMin);
 716
 717     // New kinematics  by sharing mHat appropriately.
 718     Vec4 p1New, p2New;
 719     if ( !newKin( event[i1Old].p(), event[i2Old].p(),
 720       mEff, mEff, p1New, p2New, false) ) return false;
 721
 722     // Insert new partons and update history.
 723     int i1New, i2New;
 724     if (event[systemPtr->iParton[0]].acol() == 0) {
 725       i1New = event.append( -idNewQ, 101, i1Old, 0, 0, 0,
 726         0, event[i1Old].acol(), p1New, mEff, 0.);
 727       i2New = event.append(  idNewQ, 101, i2Old, 0, 0, 0,
 728         event[i2Old].col(), 0, p2New, mEff, 0.);
 729     } else {
 730       i1New = event.append(  idNewQ, 101, i1Old, 0, 0, 0,
 731         event[i1Old].col(), 0, p1New, mEff, 0.);
 732       i2New = event.append( -idNewQ, 101, i2Old, 0, 0, 0,
 733         0, event[i2Old].acol(), p2New, mEff, 0.);
 734     }
 735     event[i1Old].statusNeg();
 736     event[i2Old].statusNeg();
 737     event[i1Old].daughters( i1New, 0);
 738     event[i2Old].daughters( i2New, 0);
 739
 740     // Partons in the two new systems.
 741     for (int i = 0; i < iMin; ++i)
 742       iNewSys1.push_back( systemPtr->iParton[i]);
 743     iNewSys1.push_back( i1New);
 744     iNewSys2.push_back( i2New);
 745     for (int i = iMin + 2; i < int(systemPtr->size()); ++i)
 746       iNewSys2.push_back( systemPtr->iParton[i]);
 747   }
 748
 749   // Erase the old system and insert the two new ones instead.
 750   colConfig.erase(iSys);
 751   colConfig.insert( iNewSys1, event);
 752   colConfig.insert( iNewSys2, event);
 753
 754   // Done.
 755   return true;
 756
 757 }
 758
 759 //--------------------------------------------------------------------------
 760
 761 // Produce a R-hadron from a squark and another string end.
 762
 763 bool RHadrons::produceSquark( ColConfig& colConfig, Event& event) {
 764
 765   // Initial values.
 766   int    nBody  = 0;
 767   int    iRNow  = 0;
 768   int    iNewQ  = 0;
 769   int    iNewL  = 0;
 770
 771   // Check at which end of the string the squark is located.
 772   int    idAbsTop = event[ systemPtr->iParton[0] ].idAbs();
 773   bool   sqAtTop  = (allowRSb && idAbsTop == idRSb)
 774                  || (allowRSt && idAbsTop == idRSt);
 775
 776   // Copy down system consecutively from squark end.
 777   int    iBeg = event.size();
 778   iCreRHad[iRHad] = iBeg;
 779   if (sqAtTop) for (int i = 0; i < systemPtr->size(); ++i)
 780     event.copy( systemPtr->iParton[i], 102);
 781   else         for (int i = systemPtr->size() - 1; i >= 0; --i)
 782     event.copy( systemPtr->iParton[i], 102);
 783   int    iEnd = event.size() - 1;
 784
 785   // Input flavours. From now on H = heavy and L = light string ends.
 786   int    idOldH = event[iBeg].id();
 787   int    idOldL = event[iEnd].id();
 788
 789   // Pick new flavour to form R-hadron.
 790   FlavContainer flavOld( idOldH%10);
 791   int    idNewQ = flavSelPtr->pick(flavOld).id;
 792   int    idRHad = toIdWithSquark( idOldH, idNewQ);
 793   if (idRHad == 0) {
 794      infoPtr->errorMsg("Error in RHadrons::produceSquark: "
 795        "cannot form R-hadron code");
 796      return false;
 797   }
 798
 799   // Target mass of R-hadron and z value of fragmentation function.
 800   double mRHad  = particleDataPtr->m0(idRHad) + event[iBeg].m()
 801                 - ( (abs(idOldH) == idRSb) ? m0Sb : m0St );
 802   double z      = zSelPtr->zFrag( idOldH, idNewQ, mRHad*mRHad);
 803
 804   // Basic kinematics of string piece where break is to occur.
 805   Vec4   pOldH  = event[iBeg].p();
 806   int    iOldL  = iBeg + 1;
 807   Vec4   pOldL  = event[iOldL].p();
 808   double mOldL  = event[iOldL].m();
 809   double mNewH  = mRHad / z;
 810   double sSys   = (pOldH + pOldL).m2Calc();
 811   double sRem   = (1. - z) * (sSys - mNewH*mNewH);
 812   double sMin   = pow2(mOldL + mCollapseRH);
 813
 814   // If too low remaining mass in system then add one more parton to it.
 815   while ( ( sRem < sMin || sSys < pow2(mNewH + mOldL + MSAFETY) )
 816     && iOldL < iEnd ) {
 817     ++iOldL;
 818     pOldL      += event[iOldL].p();
 819     mOldL       = event[iOldL].m();
 820     sSys        = (pOldH + pOldL).m2Calc();
 821     sRem        = (1. - z) * (sSys - mNewH*mNewH);
 822     sMin        = pow2(mOldL + mCollapseRH);
 823   }
 824
 825   // If enough mass then split off R-hadron and reduced system.
 826   if ( sRem > sMin && sSys > pow2(mNewH + mOldL + MSAFETY) ) {
 827     Vec4 pNewH, pNewL;
 828     if ( !newKin( pOldH, pOldL, mNewH, mOldL, pNewH, pNewL) ) {
 829       infoPtr->errorMsg("Error in RHadrons::produceSquark: "
 830        "failed to construct kinematics with reduced system");
 831       return false;
 832     }
 833
 834     // Insert R-hadron with new momentum.
 835     iRNow       = event.append( idRHad, 104, iBeg, iOldL, 0, 0, 0, 0,
 836       z * pNewH, mRHad, 0.);
 837
 838     // Reduced system with new string endpoint and modified recoiler.
 839     idNewQ      = -idNewQ;
 840     bool hasCol = (idNewQ > 0 && idNewQ < 10) || idNewQ < -10;
 841     int  col    = (hasCol) ? event[iOldL].acol() : 0;
 842     int  acol   = (hasCol) ? 0 : event[iOldL].col();
 843     iNewQ       = event.append( idNewQ, 105, iBeg, iOldL, 0, 0, col, acol,
 844       (1. - z) * pNewH, (1. - z) * mNewH, 0.);
 845     iNewL       = event.copy( iOldL, 105);
 846     event[iNewL].mothers( iBeg, iOldL);
 847     event[iNewL].p( pNewL);
 848
 849     // Done with processing of split to R-hadron and reduced system.
 850     nBody = 3;
 851   }
 852
 853   // Two-body final state: form light hadron from endpoint and new flavour.
 854   if (nBody == 0) {
 855     FlavContainer flav1( idOldL);
 856     FlavContainer flav2( -idNewQ);
 857     int iTry   = 0;
 858     int idNewL = flavSelPtr->combine( flav1, flav2);
 859     while (++iTry < NTRYMAX && idNewL == 0)
 860       idNewL = flavSelPtr->combine( flav1, flav2);
 861     if (idNewL == 0) {
 862        infoPtr->errorMsg("Error in RHadrons::produceSquark: "
 863          "cannot form light hadron code");
 864        return false;
 865     }
 866     double mNewL = particleDataPtr->mass( idNewL);
 867
 868     // Check that energy enough for light hadron and R-hadron.
 869     if ( sSys > pow2(mRHad + mNewL + MSAFETY) ) {
 870       Vec4 pRHad, pNewL;
 871       if ( !newKin( pOldH, pOldL, mRHad, mNewL, pRHad, pNewL) ) {
 872         infoPtr->errorMsg("Error in RHadrons::produceSquark: "
 873          "failed to construct kinematics for two-hadron decay");
 874         return false;
 875       }
 876
 877       // Insert R-hadron and light hadron.
 878       iRNow = event.append( idRHad, 104, iBeg, iOldL, 0, 0, 0, 0,
 879         pRHad, mRHad, 0.);
 880       event.append( idNewL, 105, iBeg, iOldL, 0, 0, 0, 0,
 881         pNewL, mNewL, 0.);
 882
 883       // Done for two-body case.
 884       nBody = 2;
 885     }
 886   }
 887
 888   // Final case: for very low mass collapse to one single R-hadron.
 889   if (nBody == 0) {
 890     idRHad = toIdWithSquark( idOldH, idOldL);
 891     if (idRHad == 0) {
 892        infoPtr->errorMsg("Error in RHadrons::produceSquark: "
 893          "cannot form R-hadron code");
 894        return false;
 895     }
 896
 897     // Insert R-hadron with new momentum.
 898     iRNow = event.append( idRHad, 104, iBeg, iOldL, 0, 0, 0, 0,
 899       systemPtr->pSum, systemPtr->mass, 0.);
 900
 901     // Done with one-body case.
 902     nBody = 1;
 903   }
 904
 905   // History bookkeeping: the R-hadron and processed partons.
 906   iRHadron[iRHad] = iRNow;
 907   int iLast = event.size() - 1;
 908   for (int i = iBeg; i <= iOldL; ++i) {
 909     event[i].statusNeg();
 910     event[i].daughters( iRNow, iLast);
 911   }
 912
 913   // Remove old string system and, if needed, insert a new one.
 914   colConfig.erase(iSys);
 915   if (nBody == 3) {
 916     vector<int> iNewSys;
 917     iNewSys.push_back( iNewQ);
 918     iNewSys.push_back( iNewL);
 919     for ( int i = iOldL + 1; i <= iEnd; ++i) iNewSys.push_back( i);
 920     colConfig.insert( iNewSys, event);
 921   }
 922
 923   // Copy lifetime and vertex from sparticle to R-hadron.
 924   event[iRNow].tau( event[iBef].tau() );
 925   if (event[iBef].hasVertex()) event[iRNow].vProd( event[iBef].vProd() );
 926
 927   // Done with production of a R-hadron from a squark.
 928   return true;
 929
 930 }
 931
 932 //--------------------------------------------------------------------------
 933
 934 // Produce a R-hadron from a gluino and two string ends (or a gluon).
 935
 936 bool RHadrons::produceGluino( ColConfig& colConfig, Event& event) {
 937
 938   // Initial values.
 939   int    iGlui   = 0;
 940   int    idSave  = 0;
 941   int    idQLeap = 0;
 942   bool   isDiq1  = false;
 943   vector<int> iSide1, iSide2, iSideTmp, iNewSys1, iNewSys2;
 944   Vec4   pGlui, pSide1, pSide2;
 945
 946   // Decide whether to produce a gluinoball or not.
 947   bool isGBall = (rndmPtr->flat() < probGluinoballRH);
 948
 949   // Extract one string system on either side of the gluino.
 950   for (int i = 0; i < systemPtr->size(); ++i) {
 951     int  iTmp   = systemPtr->iParton[i];
 952     if (iGlui == 0 && event[ iTmp ].id() == idRGo) {
 953       iGlui     = iTmp;
 954       pGlui     = event[ iTmp ].p();
 955     } else if (iGlui == 0) {
 956       iSideTmp.push_back( iTmp);
 957       pSide1   += event[ iTmp ].p();
 958     } else {
 959       iSide2.push_back( iTmp);
 960       pSide2   += event[ iTmp ].p();
 961     }
 962   }
 963
 964   // Order sides from gluino outwards and with lowest relative mass first.
 965   for (int i = int(iSideTmp.size()) - 1; i >= 0; --i)
 966     iSide1.push_back( iSideTmp[i]);
 967   double m1H    = (pGlui + pSide1).mCalc() - event[ iSide1.back() ].m();
 968   double m2H    = (pGlui + pSide2).mCalc() - event[ iSide2.back() ].m();
 969   if (m2H < m1H) {
 970     swap( iSide1, iSide2);
 971     swap( pSide1, pSide2);
 972   }
 973
 974   // Begin loop over two sides of gluino, with lowest relative mass first.
 975   for (int iSide = 1; iSide < 3; ++iSide) {
 976
 977     // Begin processing of lower-mass string side.
 978     int    idRHad = 0;
 979     double mRHad  = 0.;
 980     int    nBody  = 0;
 981     int    iRNow  = 0;
 982     int    iNewQ  = 0;
 983     int    iNewL  = 0;
 984     int    statusRHad = 0;
 985
 986     // Copy down system consecutively from gluino end.
 987     int    iBeg   = event.size();
 988     if (iSide == 1) {
 989       iCreRHad[iRHad] = iBeg;
 990       event.copy( iGlui, 102);
 991       for (int i = 0; i < int(iSide1.size()); ++i)
 992         event.copy( iSide1[i], 102);
 993     } else {
 994       event.copy( iGlui, 102);
 995       for (int i = 0; i < int(iSide2.size()); ++i)
 996         event.copy( iSide2[i], 102);
 997     }
 998     int    iEnd   = event.size() - 1;
 999
1000     // Pick new flavour to help form R-hadron. Initial colour values.
1001     int    idOldL = event[iEnd].id();
1002     int    idNewQ = 21;
1003     if (!isGBall) {
1004       do {
1005         FlavContainer flavOld( idOldL);
1006         idNewQ = -flavSelPtr->pick(flavOld).id;
1007       } while (iSide == 2 && isDiq1 && abs(idNewQ) > 10);
1008       if (iSide == 1) isDiq1 = (abs(idNewQ) > 10);
1009     }
1010     bool   hasCol = (event[iEnd].col() > 0);
1011     int    colR   = 0;
1012     int    acolR  = 0;
1013
1014     // Target intermediary mass or R-hadron mass.
1015     if (iSide == 1) {
1016       idSave      = idNewQ;
1017       idRHad      = (hasCol) ? 1009002 : -1009002;
1018       if (hasCol) colR  = event[iBeg].col();
1019       else        acolR = event[iBeg].acol();
1020       statusRHad  = 103;
1021       double mNewQ = particleDataPtr->constituentMass( idNewQ);
1022       if (isGBall) mNewQ *= 0.5;
1023       mRHad       = event[iBeg].m() + mOffsetCloudRH + mNewQ;
1024     } else {
1025       idRHad      = toIdWithGluino( idSave, idNewQ);
1026       if (idRHad == 0) {
1027          infoPtr->errorMsg("Error in RHadrons::produceGluino: "
1028            "cannot form R-hadron code");
1029          return false;
1030       }
1031       statusRHad  = 104;
1032       mRHad       = particleDataPtr->m0(idRHad) + event[iBeg].m() - m0Go;
1033     }
1034
1035     // z value of fragmentation function.
1036     double z      = zSelPtr->zFrag( idRGo, idNewQ, mRHad*mRHad);
1037
1038     // Basic kinematics of string piece where break is to occur.
1039     Vec4   pOldH  = event[iBeg].p();
1040     int    iOldL  = iBeg + 1;
1041     Vec4   pOldL  = event[iOldL].p();
1042     double mOldL  = event[iOldL].m();
1043     double mNewH  = mRHad / z;
1044     double sSys   = (pOldH + pOldL).m2Calc();
1045     double sRem   = (1. - z) * (sSys - mNewH*mNewH);
1046     double sMin   = pow2(mOldL + mCollapseRH);
1047
1048     // If too low remaining mass in system then add one more parton to it.
1049     while ( ( sRem < sMin || sSys < pow2(mNewH + mOldL + MSAFETY) )
1050       && iOldL < iEnd ) {
1051       ++iOldL;
1052       pOldL      += event[iOldL].p();
1053       mOldL       = event[iOldL].m();
1054       sSys        = (pOldH + pOldL).m2Calc();
1055       sRem        = (1. - z) * (sSys - mNewH*mNewH);
1056       sMin        = pow2(mOldL + mCollapseRH);
1057     }
1058
1059     // If enough mass then split off R-hadron and reduced system.
1060     if ( sRem > sMin && sSys > pow2(mNewH + mOldL + MSAFETY) ) {
1061       Vec4 pNewH, pNewL;
1062       if ( !newKin( pOldH, pOldL, mNewH, mOldL, pNewH, pNewL) ) {
1063         infoPtr->errorMsg("Error in RHadrons::produceGluino: "
1064          "failed to construct kinematics with reduced system");
1065         return false;
1066       }
1067
1068       // Insert intermediary or R-hadron with new momentum, less colour.
1069       iRNow       = event.append( idRHad, statusRHad, iBeg, iOldL,
1070         0, 0, colR, acolR, z * pNewH, mRHad, 0.);
1071
1072       // Reduced system with new string endpoint and modified recoiler.
1073       if (!isGBall) idNewQ = -idNewQ;
1074       int  colN   = (hasCol) ? 0 : event[iOldL].acol();
1075       int  acolN  = (hasCol) ? event[iOldL].col() : 0;
1076       iNewQ       = event.append( idNewQ, 105, iBeg, iOldL, 0, 0,
1077         colN, acolN, (1. - z) * pNewH, (1. - z) * mNewH, 0.);
1078       iNewL       = event.copy( iOldL, 105);
1079       event[iNewL].mothers( iBeg, iOldL);
1080       event[iNewL].p( pNewL);
1081
1082       // Collect partons of new string system.
1083       if (iSide == 1) {
1084         iNewSys1.push_back( iNewQ);
1085         iNewSys1.push_back( iNewL);
1086         for ( int i = iOldL + 1; i <= iEnd; ++i) iNewSys1.push_back( i);
1087       } else {
1088         iNewSys2.push_back( iNewQ);
1089         iNewSys2.push_back( iNewL);
1090         for ( int i = iOldL + 1; i <= iEnd; ++i) iNewSys2.push_back( i);
1091       }
1092
1093       // Done with processing of split to R-hadron and reduced system.
1094       nBody = 3;
1095     }
1096
1097     // For side-1 low-mass glueball system reabsorb full momentum.
1098     if (nBody == 0 && isGBall && iSide == 1) {
1099       idQLeap    = event[iOldL].id();
1100       Vec4 pNewH = event[iBeg].p() + pOldL;
1101
1102       // Insert intermediary R-hadron with new momentum, less colour.
1103       iRNow      = event.append( idRHad, statusRHad, iBeg, iEnd,
1104         0, 0, colR, acolR, pNewH, pNewH.mCalc(), 0.);
1105       nBody = 1;
1106     }
1107
1108     // Two-body final state: form light hadron from endpoint and new flavour.
1109     // Also for side 2 if gluinoball formation gives quark from side 1.
1110     if (nBody == 0 && (!isGBall || (iSide == 2 && idQLeap != 0) )) {
1111       if (isGBall) idNewQ = -idQLeap;
1112       FlavContainer flav1( idOldL);
1113       FlavContainer flav2( -idNewQ);
1114       int iTry   = 0;
1115       int idNewL = flavSelPtr->combine( flav1, flav2);
1116       while (++iTry < NTRYMAX && idNewL == 0)
1117         idNewL = flavSelPtr->combine( flav1, flav2);
1118       if (idNewL == 0) {
1119          infoPtr->errorMsg("Error in RHadrons::produceGluino: "
1120            "cannot form light hadron code");
1121          return false;
1122       }
1123       double mNewL = particleDataPtr->mass( idNewL);
1124
1125       // Check that energy enough for light hadron and R-hadron.
1126       if ( sSys > pow2(mRHad + mNewL + MSAFETY) ) {
1127         Vec4 pRHad, pNewL;
1128         if ( !newKin( pOldH, pOldL, mRHad, mNewL, pRHad, pNewL) ) {
1129           infoPtr->errorMsg("Error in RHadrons::produceGluino: "
1130            "failed to construct kinematics for two-hadron decay");
1131           return false;
1132         }
1133
1134         // Insert intermediary or R-hadron and light hadron.
1135         iRNow = event.append( idRHad, statusRHad, iBeg, iOldL, 0, 0,
1136           colR, acolR, pRHad, mRHad, 0.);
1137         event.append( idNewL, 105, iBeg, iOldL, 0, 0, 0, 0,
1138           pNewL, mNewL, 0.);
1139
1140         // Done for two-body case.
1141         nBody   = 2;
1142         // For this case gluinoball should be handled as normal flavour.
1143         isGBall = false;
1144       }
1145     }
1146
1147     // Final case: for very low mass collapse to one single R-hadron.
1148     if (nBody == 0 && (!isGBall || (iSide == 2 && idQLeap != 0) )) {
1149       if (isGBall) idSave = idQLeap;
1150       if (iSide == 1) idSave = idOldL;
1151       else            idRHad = toIdWithGluino( idSave, idOldL);
1152       if (idRHad == 0) {
1153          infoPtr->errorMsg("Error in RHadrons::produceGluino: "
1154            "cannot form R-hadron code");
1155          return false;
1156       }
1157
1158       // Insert R-hadron with new momentum.
1159       iRNow = event.append( idRHad, statusRHad, iBeg, iOldL, 0, 0,
1160         colR, acolR, pOldH + pOldL, (pOldH + pOldL).mCalc(), 0.);
1161
1162       // Done with one-body case.
1163       // Even if hoped-for, it was not possible to create a gluinoball.
1164       isGBall = false;
1165     }
1166
1167     // History bookkeeping: the processed partons.
1168     int iLast = event.size() - 1;
1169     for (int i = iBeg; i <= iOldL; ++i) {
1170       event[i].statusNeg();
1171       event[i].daughters( iRNow, iLast);
1172     }
1173
1174     // End of loop over two sides of the gluino.
1175     iGlui   = iRNow;
1176   }
1177
1178   // History bookkeeping: insert R-hadron; delete old string system.
1179   if (iGlui == 0) {
1180      infoPtr->errorMsg("Error in RHadrons::produceGluino: "
1181            "could not handle gluinoball kinematics");
1182      return false;
1183   }
1184   iRHadron[iRHad] = iGlui;
1185   colConfig.erase(iSys);
1186
1187   // Non-gluinoball: insert (up to) two new string systems.
1188   if (!isGBall) {
1189     if (iNewSys1.size() > 0) colConfig.insert( iNewSys1, event);
1190     if (iNewSys2.size() > 0) colConfig.insert( iNewSys2, event);
1191
1192   // Gluinoball with enough energy in first string:
1193   // join two temporary gluons into one.
1194   } else if (idQLeap == 0) {
1195     int iG1   = iNewSys1[0];
1196     int iG2   = iNewSys2[0];
1197     int colG  = event[iG1].col()  + event[iG2].col();
1198     int acolG = event[iG1].acol() + event[iG2].acol();
1199     Vec4 pG   = event[iG1].p()    + event[iG2].p();
1200     int iG12  = event.append( 21, 107, iG1, iG2, 0, 0, colG, acolG,
1201       pG, pG.mCalc(), 0.);
1202
1203     // Temporary gluons no longer needed, but new colour to avoid warnings.
1204     event[iG1].id( 21);
1205     event[iG2].id( 21);
1206     event[iG1].statusNeg();
1207     event[iG2].statusNeg();
1208     event[iG1].daughter1( iG12);
1209     event[iG2].daughter1( iG12);
1210     int colBridge = event.nextColTag();
1211     if (event[iG1].col() == 0) {
1212       event[iG1].col(  colBridge);
1213       event[iG2].acol( colBridge);
1214     } else {
1215       event[iG1].acol( colBridge);
1216       event[iG2].col(  colBridge);
1217     }
1218
1219     // Form the remnant system from which the R-hadron has been split off.
1220     vector<int> iNewSys12;
1221     for (int i = int(iNewSys1.size()) - 1; i > 0; --i)
1222       iNewSys12.push_back( iNewSys1[i]);
1223     iNewSys12.push_back( iG12);
1224     for (int i = 1; i < int(iNewSys2.size()); ++i)
1225       iNewSys12.push_back( iNewSys2[i]);
1226     colConfig.insert( iNewSys12, event);
1227
1228   // Gluinoball where side 1 was fully eaten, and its flavour content
1229   // should leap over to the other side, to a gluon there.
1230   } else {
1231     int iG2   = iNewSys2[0];
1232     event[iG2].id( idQLeap);
1233     colConfig.insert( iNewSys2, event);
1234   }
1235
1236   // Copy lifetime and vertex from sparticle to R-hadron.
1237   event[iGlui].tau( event[iBef].tau() );
1238   if (event[iBef].hasVertex()) event[iGlui].vProd( event[iBef].vProd() );
1239
1240   // Done with production of a R-hadron from a gluino.
1241   return true;
1242
1243 }
1244
1245 //--------------------------------------------------------------------------
1246
1247 // Form a R-hadron code from a squark and a (di)quark code.
1248 // First argument is the (anti)squark, second the (anti)(di)quark.
1249
1250 int RHadrons::toIdWithSquark( int id1, int id2) {
1251
1252   // Check that physical combination; return 0 if not.
1253   int id1Abs = abs(id1);
1254   int id2Abs = abs(id2);
1255   if (id2Abs < 10 && id1 > 0 && id2 > 0) return 0;
1256   if (id2Abs < 10 && id1 < 0 && id2 < 0) return 0;
1257   if (id2Abs > 10 && id1 > 0 && id2 < 0) return 0;
1258   if (id2Abs > 10 && id1 < 0 && id2 > 0) return 0;
1259
1260   // Form R-hadron code. Flip sign for antisquark.
1261   bool isSt = (id1Abs == idRSt);
1262   int idRHad = 1000000;
1263   if (id2Abs < 10) idRHad += ((isSt) ? 600 : 500) + 10 * id2Abs + 2;
1264   else idRHad += ((isSt) ? 6000 : 5000) + 10 * (id2Abs/100) + id2Abs%10;
1265   if (id1 < 0) idRHad = -idRHad;
1266
1267   // Done.
1268   return idRHad;
1269
1270 }
1271
1272 //--------------------------------------------------------------------------
1273
1274 // Resolve a R-hadron code into a squark and a (di)quark.
1275 // Return a pair, where first is the squark and the second the (di)quark.
1276
1277 pair<int,int> RHadrons::fromIdWithSquark( int idRHad) {
1278
1279   // Find squark flavour content.
1280   int idLight = (abs(idRHad) - 1000000) / 10;
1281   int idSq    = (idLight < 100) ? idLight/10 : idLight/100;
1282   int id1     = (idSq == 6) ? idRSt : idRSb;
1283   if (idRHad < 0) id1 = -id1;
1284
1285   // Find light (di)quark flavour content.
1286   int id2     =  (idLight < 100) ? idLight%10 : idLight%100;
1287   if (id2 > 10) id2 = 100 * id2 + abs(idRHad)%10;
1288   if ((id2 < 10 && idRHad > 0) || (id2 > 10 && idRHad < 0)) id2 = -id2;
1289
1290   // Done.
1291   return make_pair( id1, id2);
1292
1293 }
1294
1295 //--------------------------------------------------------------------------
1296
1297 // Form a R-hadron code from two quark/diquark endpoints and a gluino.
1298
1299 int RHadrons::toIdWithGluino( int id1, int id2) {
1300
1301   // Check that physical combination; return 0 if not. Handle gluinoball.
1302   int id1Abs = abs(id1);
1303   int id2Abs = abs(id2);
1304   if (id1Abs == 21 && id2Abs == 21) return 1000993;
1305   int idMax  = max( id1Abs, id2Abs);
1306   int idMin  = min( id1Abs, id2Abs);
1307   if (idMin > 10) return 0;
1308   if (idMax > 10 && id1 > 0 && id2 < 0) return 0;
1309   if (idMax > 10 && id1 < 0 && id2 > 0) return 0;
1310   if (idMax < 10 && id1 > 0 && id2 > 0) return 0;
1311   if (idMax < 10 && id1 < 0 && id2 < 0) return 0;
1312
1313   // Form R-meson code. Flip sign for antiparticle.
1314   int idRHad = 0;
1315   if (idMax < 10) {
1316     idRHad = 1009003 + 100 * idMax + 10 * idMin;
1317     if (idMin != idMax && idMax%2 == 1) {
1318       if (id1Abs == idMax && id1 > 0) idRHad = -idRHad;
1319       if (id2Abs == idMax && id2 > 0) idRHad = -idRHad;
1320     }
1321     if (idMin != idMax && idMax%2 == 0) {
1322       if (id1Abs == idMax && id1 < 0) idRHad = -idRHad;
1323       if (id2Abs == idMax && id2 < 0) idRHad = -idRHad;
1324     }
1325
1326   // Form R-baryon code. Flip sign for antiparticle.
1327   } else {
1328     int idA = idMax/1000;
1329     int idB = (idMax/100)%10;
1330     int idC = idMin;
1331     if (idC > idB) swap( idB, idC);
1332     if (idB > idA) swap( idA, idB);
1333     if (idC > idB) swap( idB, idC);
1334     idRHad  = 1090004 + 1000 * idA + 100 * idB + 10 * idC;
1335     if (id1 < 0) idRHad = -idRHad;
1336   }
1337
1338   // Done.
1339   return idRHad;
1340
1341 }
1342
1343 //--------------------------------------------------------------------------
1344
1345 // Resolve a R-hadron code into two quark/diquark endpoints (and a gluino).
1346 // Return a pair, where first carries colour and second anticolour.
1347
1348 pair<int,int> RHadrons::fromIdWithGluino( int idRHad) {
1349
1350   // Find light flavour content of R-hadron.
1351   int idLight = (abs(idRHad) - 1000000) / 10;
1352   int id1, id2, idTmp, idA, idB, idC;
1353
1354   // Gluinoballs: split g into d dbar or u ubar.
1355   if (idLight < 100) {
1356     id1 = (rndmPtr->flat() < 0.5) ? 1 : 2;
1357     id2 = -id1;
1358
1359   // Gluino-meson: split into q + qbar.
1360   } else if (idLight < 1000) {
1361     id1 = (idLight / 10) % 10;
1362     id2 = -(idLight % 10);
1363     // Flip signs when first quark of down-type.
1364     if (id1%2 == 1) {
1365       idTmp = id1;
1366       id1   = -id2;
1367       id2   = -idTmp;
1368     }
1369
1370   // Gluino-baryon: split to q + qq (diquark).
1371   // Pick diquark at random, except if c or b involved.
1372   } else {
1373     idA = (idLight / 100) % 10;
1374     idB = (idLight / 10) % 10;
1375     idC = idLight % 10;
1376     double rndmQ = 3. * rndmPtr->flat();
1377     if (idA > 3) rndmQ = 0.5;
1378     if (rndmQ < 1.) {
1379       id1 = idA;
1380       id2 = 1000 * idB + 100 * idC + 3;
1381       if (idB != idC && rndmPtr->flat() > diquarkSpin1RH) id2 -= 2;
1382     } else if (rndmQ < 2.) {
1383       id1 = idB;
1384       id2 = 1000 * idA + 100 * idC + 3;
1385       if (idA != idC && rndmPtr->flat() > diquarkSpin1RH) id2 -= 2;
1386     } else {
1387       id1 = idC;
1388       id2 = 1000 * idA + 100 * idB +3;
1389       if (idA != idB && rndmPtr->flat() > diquarkSpin1RH) id2 -= 2;
1390     }
1391   }
1392
1393   // Flip signs for anti-R-hadron.
1394   if (idRHad < 0) {
1395     idTmp = id1;
1396     id1   = -id2;
1397     id2   = -idTmp;
1398   }
1399
1400   // Done.
1401   return make_pair( id1, id2);
1402
1403 }
1404
1405 //--------------------------------------------------------------------------
1406
1407 // Construct modified four-vectors to match modified masses:
1408 // minimal reshuffling of momentum along common axis.
1409 // Note that last two arguments are used to provide output!
1410
1411 bool RHadrons::newKin( Vec4 pOld1, Vec4 pOld2, double mNew1, double mNew2,
1412   Vec4& pNew1, Vec4& pNew2, bool checkMargin) {
1413
1414   // Squared masses in initial and final kinematics.
1415   double sSum   = (pOld1 + pOld2).m2Calc();
1416   double sOld1  = pOld1.m2Calc();
1417   double sOld2  = pOld2.m2Calc();
1418   double sNew1  = mNew1 * mNew1;
1419   double sNew2  = mNew2 * mNew2;
1420
1421   // Check that kinematically possible.
1422   if (checkMargin && pow2(mNew1 + mNew2 + MSAFETY) > sSum) return false;
1423
1424   // Transfer coefficients to give four-vectors with the new masses.
1425   double lamOld = sqrt( pow2(sSum - sOld1 - sOld2) - 4. * sOld1 * sOld2 );
1426   double lamNew = sqrt( pow2(sSum - sNew1 - sNew2) - 4. * sNew1 * sNew2 );
1427   double move1  = (lamNew * (sSum - sOld1 + sOld2)
1428                 -  lamOld * (sSum - sNew1 + sNew2)) / (2. * sSum * lamOld);
1429   double move2  = (lamNew * (sSum + sOld1 - sOld2)
1430                 -  lamOld * (sSum + sNew1 - sNew2)) / (2. * sSum * lamOld);
1431
1432   // Construct final vectors. Done.
1433   pNew1 = (1. + move1) * pOld1 - move2 * pOld2;
1434   pNew2 = (1. + move2) * pOld2 - move1 * pOld1;
1435   return true;
1436
1437 }
1438
1439 //==========================================================================
1440
1441 } // end namespace Pythia8