PYTHIA8/pythia8170/src/ResonanceDecays.cxx

   1 // ResonanceDecays.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
   7 // the ResonanceDecays class.
   8
   9 #include "ResonanceDecays.h"
  10
  11 namespace Pythia8 {
  12
  13 //==========================================================================
  14
  15 // The ResonanceDecays class.
  16 // Do all resonance decays sequentially.
  17
  18 //--------------------------------------------------------------------------
  19
  20 // Constants: could be changed here if desired, but normally should not.
  21 // These are of technical nature, as described for each.
  22
  23 // Number of tries to pick a decay channel.
  24 const int    ResonanceDecays::NTRYCHANNEL = 10;
  25
  26 // Number of tries to pick a set of daughter masses.
  27 const int    ResonanceDecays::NTRYMASSES  = 10000;
  28
  29 // Mass above threshold for allowed decays.
  30 const double ResonanceDecays::MSAFETY     = 0.1;
  31
  32 // When constrainted kinematics cut high-mass tail of Breit-Wigner.
  33 const double ResonanceDecays::WIDTHCUT    = 5.;
  34
  35 // Small number (relative to 1) to protect against roundoff errors.
  36 const double ResonanceDecays::TINY        = 1e-10;
  37
  38 // Forbid small Breit-Wigner mass range, as mapped onto atan range.
  39 const double ResonanceDecays::TINYBWRANGE = 1e-8;
  40
  41 // These numbers are hardwired empirical parameters,
  42 // intended to speed up the M-generator.
  43 const double ResonanceDecays::WTCORRECTION[11] = { 1., 1., 1.,
  44   2., 5., 15., 60., 250., 1250., 7000., 50000. };
  45
  46 //--------------------------------------------------------------------------
  47
  48 bool ResonanceDecays::next( Event& process) {
  49
  50   // Loop over all entries to find resonances that should decay.
  51   int iDec = 0;
  52   do {
  53     Particle& decayer = process[iDec];
  54     if (decayer.isFinal() && decayer.canDecay() && decayer.mayDecay()
  55     && decayer.isResonance() ) {
  56
  57       // Fill the decaying particle in slot 0 of arrays.
  58       id0    = decayer.id();
  59       m0     = decayer.m();
  60       idProd.resize(0);
  61       mProd.resize(0);
  62       idProd.push_back( id0 );
  63       mProd.push_back( m0 );
  64
  65       // Mother flavour - relevant for gamma*/Z0 mixing. (Not always??)
  66       int idIn = process[decayer.mother1()].id();
  67
  68       // Prepare decay selection.
  69       if (!decayer.particleDataEntry().preparePick(id0, m0, idIn)) {
  70         ostringstream osWarn;
  71         osWarn << "for id = " << id0;
  72         infoPtr->errorMsg("Error in ResonanceDecays::next:"
  73           " no open decay channel", osWarn.str());
  74         return false;
  75       }
  76
  77       // Pick a decay channel; allow up to ten tries.
  78       bool foundChannel = false;
  79       for (int iTryChannel = 0; iTryChannel < NTRYCHANNEL; ++iTryChannel) {
  80
  81         // Pick decay channel. Find multiplicity.
  82         DecayChannel& channel = decayer.particleDataEntry().pickChannel();
  83         mult = channel.multiplicity();
  84
  85         // Read out flavours.
  86         idProd.resize(1);
  87         int idNow;
  88         for (int i = 1; i <= mult; ++i) {
  89           idNow = channel.product(i - 1);
  90           if (id0 < 0 && particleDataPtr->hasAnti(idNow)) idNow = -idNow;
  91           idProd.push_back( idNow);
  92         }
  93
  94         // Pick masses. Pick new channel if fail.
  95         if (!pickMasses()) continue;
  96         foundChannel = true;
  97         break;
  98       }
  99
 100       // Failed to find acceptable decays.
 101       if (!foundChannel) {
 102         ostringstream osWarn;
 103         osWarn << "for id = " << id0;
 104         infoPtr->errorMsg("Error in ResonanceDecays::next:"
 105           " failed to find workable decay channel", osWarn.str());
 106         return false;
 107       }
 108
 109       // Select colours in decay.
 110       if (!pickColours(iDec, process)) return false;
 111
 112       // Select four-momenta in decay, boosted to lab frame.
 113       pProd.resize(0);
 114       pProd.push_back( decayer.p() );
 115       if (!pickKinematics()) return false;
 116
 117       // Append decay products to the process event record. Set lifetimes.
 118       int iFirst = process.size();
 119         for (int i = 1; i <= mult; ++i) {
 120           process.append( idProd[i], 23, iDec, 0, 0, 0, cols[i], acols[i],
 121             pProd[i], mProd[i], m0);
 122         }
 123       int iLast = process.size() - 1;
 124
 125       // Set decay vertex when this is displaced.
 126       if (process[iDec].hasVertex() || process[iDec].tau() > 0.) {
 127         Vec4 vDec = process[iDec].vDec();
 128         for (int i = iFirst; i <= iLast; ++i) process[i].vProd( vDec );
 129       }
 130
 131       // Set lifetime of daughters.
 132       for (int i = iFirst; i <= iLast; ++i)
 133         process[i].tau( process[i].tau0() * rndmPtr->exp() );
 134
 135       // Modify mother status and daughters.
 136       decayer.status(-22);
 137       decayer.daughters(iFirst, iLast);
 138
 139     // End of loop over all entries.
 140     }
 141   } while (++iDec < process.size());
 142
 143   // Done.
 144   return true;
 145
 146 }
 147
 148 //--------------------------------------------------------------------------
 149
 150 // Select masses of decay products.
 151
 152 bool ResonanceDecays::pickMasses() {
 153
 154   // Arrays with properties of particles. Fill with dummy values for mother.
 155   vector<bool>   useBW;
 156   vector<double> m0BW, mMinBW, mMaxBW, widthBW;
 157   double mMother  = mProd[0];
 158   double m2Mother = mMother * mMother;
 159   useBW.push_back( false );
 160   m0BW.push_back( mMother );
 161   mMinBW.push_back( mMother );
 162   mMaxBW.push_back( mMother );
 163   widthBW.push_back( 0. );
 164
 165   // Loop throught products for masses and widths. Set nominal mass.
 166   bool   useBWNow;
 167   double m0Now, mMinNow, mMaxNow, widthNow;
 168   for (int i = 1; i <= mult; ++i) {
 169     useBWNow  = particleDataPtr->useBreitWigner( idProd[i] );
 170     m0Now     = particleDataPtr->m0( idProd[i] );
 171     mMinNow   = particleDataPtr->m0Min( idProd[i] );
 172     mMaxNow   = particleDataPtr->m0Max( idProd[i] );
 173     if (useBWNow && mMaxNow < mMinNow) mMaxNow = mMother;
 174     widthNow  = particleDataPtr->mWidth( idProd[i] );
 175     useBW.push_back( useBWNow );
 176     m0BW.push_back( m0Now );
 177     mMinBW.push_back( mMinNow );
 178     mMaxBW.push_back( mMaxNow );
 179     widthBW.push_back( widthNow );
 180     mProd.push_back( m0Now );
 181   }
 182
 183   // Find number of Breit-Wigners and summed (minimal) masses.
 184   int    nBW     = 0;
 185   double mSum    = 0.;
 186   double mSumMin = 0.;
 187   for (int i = 1; i <= mult; ++i) {
 188     if (useBW[i]) ++nBW;
 189     mSum        += max( m0BW[i], mMinBW[i]);
 190     mSumMin     += mMinBW[i];
 191   }
 192
 193   // If sum of minimal masses above mother mass then give up.
 194   if (mSumMin + MSAFETY > mMother) return false;
 195
 196   // If sum of masses below and no Breit-Wigners then done.
 197   if (mSum + 0.5 * MSAFETY < mMother && nBW == 0) return true;
 198
 199   // Else if below then retry Breit-Wigners, with simple treshold.
 200   if (mSum + MSAFETY < mMother) {
 201     double wtMax = 2. * sqrtpos(1. - mSum*mSum / m2Mother);
 202     double wt;
 203     for (int iTryMasses = 0; iTryMasses <= NTRYMASSES; ++ iTryMasses) {
 204       if (iTryMasses == NTRYMASSES) return false;
 205       mSum = 0.;
 206       for (int i = 1; i <= mult; ++i) {
 207         if (useBW[i])  mProd[i] = particleDataPtr->mass( idProd[i] );
 208         mSum += mProd[i];
 209       }
 210       wt = (mSum + 0.5 * MSAFETY < mMother)
 211          ? sqrtpos(1. - mSum*mSum / m2Mother) : 0.;
 212       if (wt > rndmPtr->flat() * wtMax) break;
 213     }
 214     return true;
 215   }
 216
 217   // From now on some particles will have to be forced off shell.
 218
 219   // Order Breit-Wigners in decreasing widths. Sum of other masses.
 220   vector<int> iBW;
 221   double mSum0 = 0.;
 222   for (int i = 1; i <= mult; ++i) {
 223     if (useBW[i]) iBW.push_back(i);
 224     else          mSum0 += mProd[i];
 225   }
 226   for (int i = 1; i < nBW; ++i) {
 227     for (int j = i - 1; j >= 0; --j) {
 228       if (widthBW[iBW[j+1]] > widthBW[iBW[j]]) swap (iBW[j+1], iBW[j]);
 229       else break;
 230     }
 231   }
 232
 233   // Do all but broadest two in increasing-width order. Includes only one.
 234   if (nBW != 2) {
 235     int iMin = (nBW == 1) ? 0 : 2;
 236     for (int i = nBW - 1; i >= iMin; --i) {
 237       int iBWi = iBW[i];
 238
 239       // Find allowed mass range of current resonance.
 240       double mMax    = mMother - mSum0 - MSAFETY;
 241       if (nBW  != 1) for (int j = 0; j < i; ++j) mMax -= mMinBW[iBW[j]];
 242       mMax           = min( mMaxBW[iBWi], mMax );
 243       double mMin    = min( mMinBW[iBWi], mMax - MSAFETY);
 244       if (mMin < 0.) return false;
 245
 246       // Parameters for Breit-Wigner choice, with constrained mass range.
 247       double m2Nom   = pow2( m0BW[iBWi] );
 248       double m2Max   = mMax * mMax;
 249       double m2Min   = mMin * mMin;
 250       double mmWid   = m0BW[iBWi] * widthBW[iBWi];
 251       double atanMin = atan( (m2Min - m2Nom) / mmWid );
 252       double atanMax = atan( (m2Max - m2Nom) / mmWid );
 253       double atanDif = atanMax - atanMin;
 254
 255       // Fail if too narrow mass range; e.g. out in tail of Breit-Wigner.
 256       if (atanDif < TINYBWRANGE) return false;
 257
 258       // Retry mass according to Breit-Wigner, with simple threshold factor.
 259       double mr1     = mSum0*mSum0 / m2Mother;
 260       double mr2     = m2Min / m2Mother;
 261       double wtMax   = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
 262       double m2Now, wt;
 263       for (int iTryMasses = 0; iTryMasses <= NTRYMASSES; ++ iTryMasses) {
 264         if (iTryMasses == NTRYMASSES) return false;
 265         m2Now = m2Nom + mmWid * tan(atanMin + rndmPtr->flat() * atanDif);
 266         mr2   = m2Now / m2Mother;
 267         wt    = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
 268         if (wt > rndmPtr->flat() * wtMax) break;
 269       }
 270
 271       // Prepare to iterate for more. Done for one Breit-Wigner.
 272       mProd[iBWi] = sqrt(m2Now);
 273       mSum0        += mProd[iBWi];
 274     }
 275     if (nBW == 1) return true;
 276   }
 277
 278   // Left to do two broadest Breit-Wigners correlated, i.e. more realistic.
 279   int iBW1        = iBW[0];
 280   int iBW2        = iBW[1];
 281   int idMother    = abs(idProd[0]);
 282   int idDau1      = abs(idProd[iBW1]);
 283   int idDau2      = abs(idProd[iBW2]);
 284
 285   // In some cases known phase-space behaviour; else simple beta factor.
 286   int psMode      = 1 ;
 287   if ( (idMother == 25 || idMother == 35) && idDau1 < 19
 288     && idDau2 == idDau1 ) psMode = 3;
 289   if ( (idMother == 25 || idMother == 35 )
 290     && (idDau1 == 23 || idDau1 == 24) && idDau2 == idDau1 ) psMode = 5;
 291   if ( idMother == 36
 292     && (idDau1 == 23 || idDau1 == 24) && idDau2 == idDau1 ) psMode = 6;
 293
 294   // Find allowed mass ranges. Ensure that they are not closed.
 295   double mRem     = mMother - mSum0 - MSAFETY;
 296   double mMax1    = min( mMaxBW[iBW1], mRem - mMinBW[iBW2] );
 297   double mMin1    = min( mMinBW[iBW1], mMax1 - MSAFETY);
 298   double mMax2    = min( mMaxBW[iBW2], mRem - mMinBW[iBW1] );
 299   double mMin2    = min( mMinBW[iBW2], mMax2 - MSAFETY);
 300
 301   // At least one range must extend below half remaining mass.
 302   if (mMin1 + mMin2 > mRem) return false;
 303   double mMid     = 0.5 * mRem;
 304   bool   hasMid1  = (mMin1 < mMid);
 305   bool   hasMid2  = (mMin2 < mMid);
 306   if (!hasMid1 && !hasMid2) return false;
 307
 308   // Parameters for Breit-Wigner choice, with constrained mass range.
 309   double m2Nom1   = pow2( m0BW[iBW1] );
 310   double m2Max1   = mMax1 * mMax1;
 311   double m2Min1   = mMin1 * mMin1;
 312   double m2Mid1   = min( mMid * mMid, m2Max1);
 313   double mmWid1   = m0BW[iBW1] * widthBW[iBW1];
 314   double atanMin1 = atan( (m2Min1 - m2Nom1) / mmWid1 );
 315   double atanMax1 = atan( (m2Max1 - m2Nom1) / mmWid1 );
 316   double atanMid1 = (hasMid1) ? atan( (m2Mid1 - m2Nom1) / mmWid1 ) : 0.;
 317   double m2Nom2   = pow2( m0BW[iBW2] );
 318   double m2Max2   = mMax2 * mMax2;
 319   double m2Min2   = mMin2 * mMin2;
 320   double m2Mid2   = min( mMid * mMid, m2Max2);
 321   double mmWid2   = m0BW[iBW2] * widthBW[iBW2];
 322   double atanMin2 = atan( (m2Min2 - m2Nom2) / mmWid2 );
 323   double atanMax2 = atan( (m2Max2 - m2Nom2) / mmWid2 );
 324   double atanMid2 = (hasMid2) ? atan( (m2Mid2 - m2Nom2) / mmWid2 ) : 0.;
 325
 326   // Relative weight to pick either below half remaining mass.
 327   double probLow1 = (hasMid1) ? 1. : 0.;
 328   if (hasMid1 && hasMid2) {
 329     double intLow1 = (atanMid1 - atanMin1) * (atanMax2 - atanMin2);
 330     double intLow2 = (atanMax1 - atanMin1) * (atanMid2 - atanMin2);
 331     probLow1 = intLow1 / (intLow1 + intLow2);
 332   }
 333
 334   // Maximum matrix element times phase space weight.
 335   double m2Rem    = mRem * mRem;
 336   double mr1      = m2Min1 / m2Rem;
 337   double mr2      = m2Min2 / m2Rem;
 338   double psMax    = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
 339   double wtMax   = 1.;
 340   if      (psMode == 1) wtMax = psMax;
 341   else if (psMode == 2) wtMax = psMax * psMax;
 342   else if (psMode == 3) wtMax = pow3(psMax);
 343   else if (psMode == 5) wtMax = psMax
 344     * (pow2(1. - mr1 - mr2) + 8. * mr1 * mr2);
 345   else if (psMode == 6) wtMax = pow3(psMax);
 346
 347   // Retry mass according to Breit-Wigners, with simple threshold factor.
 348   double atanDif1, atanDif2, m2Now1, m2Now2, mNow1, mNow2, ps, wt;
 349   for (int iTryMasses = 0; iTryMasses <= NTRYMASSES; ++ iTryMasses) {
 350     if (iTryMasses == NTRYMASSES) return false;
 351
 352     // Pick either below half remaining mass.
 353     bool pickLow1 = false;
 354     if (rndmPtr->flat() < probLow1) {
 355       atanDif1 = atanMid1 - atanMin1;
 356       atanDif2 = atanMax2 - atanMin2;
 357       pickLow1 = true;
 358     } else {
 359       atanDif1 = atanMax1 - atanMin1;
 360       atanDif2 = atanMid2 - atanMin2;
 361     }
 362     m2Now1 = m2Nom1 + mmWid1 * tan(atanMin1 + rndmPtr->flat() * atanDif1);
 363     m2Now2 = m2Nom2 + mmWid2 * tan(atanMin2 + rndmPtr->flat() * atanDif2);
 364     mNow1  = sqrt(m2Now1);
 365     mNow2  = sqrt(m2Now2);
 366
 367     // Check that intended mass ordering is fulfilled.
 368     bool rejectRegion = (pickLow1) ? (mNow1 > mNow2) : (mNow2 > mNow1);
 369     if (rejectRegion) continue;
 370
 371     // Threshold weight.
 372     mr1    = m2Now1 / m2Rem;
 373     mr2    = m2Now2 / m2Rem;
 374     wt     = 0.;
 375     if (mNow1 + mNow2 + MSAFETY < mMother) {
 376       ps   = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
 377       wt   = 1.;
 378       if      (psMode == 1) wt = ps;
 379       else if (psMode == 2) wt = ps * ps;
 380       else if (psMode == 3) wt = pow3(ps);
 381       else if (psMode == 5) wt = ps
 382         * (pow2(1. - mr1 - mr2) + 8. * mr1 * mr2);
 383       else if (psMode == 6) wt = pow3(ps)*mr1*mr2;
 384     }
 385     if (wt > rndmPtr->flat() * wtMax) break;
 386   }
 387   mProd[iBW1] = mNow1;
 388   mProd[iBW2] = mNow2;
 389
 390   // Done.
 391   return true;
 392
 393 }
 394
 395 //--------------------------------------------------------------------------
 396
 397 // Select colours of decay products.
 398
 399 bool ResonanceDecays::pickColours(int iDec, Event& process) {
 400
 401   // Reset or create arrays with colour info.
 402   cols.resize(0);
 403   acols.resize(0);
 404   vector<int> iTriplet, iAtriplet, iOctet, iDipCol, iDipAcol;
 405
 406   // Mother colours already known.
 407   int col0     = process[iDec].col();
 408   int acol0    = process[iDec].acol();
 409   int colType0 = process[iDec].colType();
 410   cols.push_back(  col0);
 411   acols.push_back(acol0);
 412
 413   // Loop through all daughters.
 414   int colTypeNow;
 415   for (int i = 1; i <= mult; ++i) {
 416     // Daughter colours initially empty, so that all is set for singlet.
 417     cols.push_back(0);
 418     acols.push_back(0);
 419     // Find character (singlet, triplet, antitriplet, octet) of daughters.
 420     colTypeNow = particleDataPtr->colType( idProd[i] );
 421     if      (colTypeNow ==  0);
 422     else if (colTypeNow ==  1) iTriplet.push_back(i);
 423     else if (colTypeNow == -1) iAtriplet.push_back(i);
 424     else if (colTypeNow ==  2) iOctet.push_back(i);
 425     // Add two entries for sextets;
 426     else if (colTypeNow ==  3) {
 427       iTriplet.push_back(i);
 428       iTriplet.push_back(i);
 429     } else if (colTypeNow == -3) {
 430       iAtriplet.push_back(i);
 431       iAtriplet.push_back(i);
 432     } else {
 433       infoPtr->errorMsg("Error in ResonanceDecays::pickColours:"
 434         " unknown colour type encountered");
 435       return false;
 436     }
 437   }
 438
 439   // Check excess of colours and anticolours in final over initial state.
 440   int nCol = iTriplet.size();
 441   if (colType0 == 1 || colType0 == 2) nCol -= 1;
 442   else if (colType0 == 3) nCol -= 2;
 443   int nAcol = iAtriplet.size();
 444   if (colType0 == -1 || colType0 == 2) nAcol -= 1;
 445   else if (colType0 == -3) nAcol -= 2;
 446
 447   // If net creation of three colours then find junction kind:
 448   // mother is 1 = singlet, triplet, or sextet (no incoming RPV tags)
 449   //           3 = antitriplet, octet, or anti-sextet (acol0 = incoming RPV tag)
 450   //           5 = not applicable to decays (needs two incoming RPV tags)
 451   if (nCol - nAcol == 3) {
 452     int kindJun = (colType0 == 0 || colType0 == 1 || colType0 == 3) ? 1 : 3;
 453
 454     // Set colours in three junction legs and store junction.
 455     int colJun[3];
 456     colJun[0] = (kindJun == 1) ? process.nextColTag() : acol0;
 457     colJun[1] = process.nextColTag();
 458     colJun[2] = process.nextColTag();
 459     process.appendJunction( kindJun, colJun[0], colJun[1], colJun[2]);
 460
 461     // Loop over three legs. Remove an incoming anticolour on first leg.
 462     for (int leg = 0; leg < 3; ++leg) {
 463       if (leg == 0 && kindJun != 1) acol0 = 0;
 464
 465       // Pick final-state triplets to carry these new colours.
 466       else {
 467         int pickT    = (iTriplet.size() == 1) ? 0
 468           : int( TINY + rndmPtr->flat() * (iTriplet.size() - TINY) );
 469         int iPickT   = iTriplet[pickT];
 470         cols[iPickT] = colJun[leg];
 471
 472         // Remove matched triplet and store new colour dipole ends.
 473         iTriplet[pickT] = iTriplet.back();
 474         iTriplet.pop_back();
 475         iDipCol.push_back(iPickT);
 476         iDipAcol.push_back(0);
 477       }
 478     }
 479
 480     // Update colour counter. Done with junction.
 481     nCol -= 3;
 482   }
 483
 484   // If net creation of three anticolours then find antijunction kind:
 485   // mother is 2 = singlet, antitriplet, or antisextet (no incoming RPV tags)
 486   //           4 = triplet, octet, or sextet (col0 = incoming RPV tag)
 487   //           6 = not applicable to decays (needs two incoming RPV tags)
 488   if (nAcol - nCol == 3) {
 489     int kindJun = (colType0 == 0 || colType0 == -1 || colType0 == -3) ? 2 : 4;
 490
 491     // Set anticolours in three antijunction legs and store antijunction.
 492     int acolJun[3];
 493     acolJun[0] = (kindJun == 2) ? process.nextColTag() : col0;
 494     acolJun[1] = process.nextColTag();
 495     acolJun[2] = process.nextColTag();
 496     process.appendJunction( kindJun, acolJun[0], acolJun[1], acolJun[2]);
 497
 498     // Loop over three legs. Remove an incoming colour on first leg.
 499     for (int leg = 0; leg < 3; ++leg) {
 500       if (leg == 0 && kindJun != 2) col0 = 0;
 501
 502       // Pick final-state antitriplets to carry these new anticolours.
 503       else {
 504         int pickA     = (iAtriplet.size() == 1) ? 0
 505           : int( TINY + rndmPtr->flat() * (iAtriplet.size() - TINY) );
 506         int iPickA    = iAtriplet[pickA];
 507         acols[iPickA] = acolJun[leg];
 508
 509         // Remove matched antitriplet and store new colour dipole ends.
 510         iAtriplet[pickA] = iAtriplet.back();
 511         iAtriplet.pop_back();
 512         iDipCol.push_back(0);
 513         iDipAcol.push_back(iPickA);
 514       }
 515     }
 516
 517     // Update anticolour counter. Done with antijunction.
 518     nAcol -= 3;
 519   }
 520
 521   // If colours and anticolours do not match now then unphysical.
 522   if (nCol != nAcol) {
 523     infoPtr->errorMsg("Error in ResonanceDecays::pickColours:"
 524       " inconsistent colour tags");
 525     return false;
 526   }
 527
 528   // Pick final-state triplet (if any) to carry initial colour.
 529   if (col0 > 0 && iTriplet.size() > 0) {
 530     int pickT    = (iTriplet.size() == 1) ? 0
 531       : int( TINY + rndmPtr->flat() * (iTriplet.size() - TINY) );
 532     int iPickT = iTriplet[pickT];
 533     cols[iPickT] = col0;
 534
 535     // Remove matched triplet and store new colour dipole ends.
 536     col0 = 0;
 537     iTriplet[pickT] = iTriplet.back();
 538     iTriplet.pop_back();
 539     iDipCol.push_back(iPickT);
 540     iDipAcol.push_back(0);
 541   }
 542
 543   // Pick final-state antitriplet (if any) to carry initial anticolour.
 544   if (acol0 > 0 && iAtriplet.size() > 0) {
 545     int pickA = (iAtriplet.size() == 1) ? 0
 546       : int( TINY + rndmPtr->flat() * (iAtriplet.size() - TINY) );
 547     int iPickA = iAtriplet[pickA];
 548     acols[iPickA] = acol0;
 549
 550     // Remove matched antitriplet and store new colour dipole ends.
 551     acol0 = 0;
 552     iAtriplet[pickA] = iAtriplet.back();
 553     iAtriplet.pop_back();
 554     iDipCol.push_back(0);
 555     iDipAcol.push_back(iPickA);
 556   }
 557
 558   // Sextets: second final-state triplet (if any)
 559   if (acol0 < 0 && iTriplet.size() > 0) {
 560     int pickT = (iTriplet.size() == 1) ? 0
 561       : int( TINY + rndmPtr->flat() * (iTriplet.size() - TINY) );
 562     int iPickT = iTriplet[pickT];
 563     cols[iPickT] = -acol0;
 564
 565     // Remove matched antitriplet and store new colour dipole ends.
 566     acol0 = 0;
 567     iTriplet[pickT] = iTriplet.back();
 568     iTriplet.pop_back();
 569     iDipCol.push_back(iPickT);
 570     iDipAcol.push_back(0);
 571   }
 572
 573   // Sextets: second final-state antitriplet (if any)
 574   if (col0 < 0 && iAtriplet.size() > 0) {
 575     int pickA    = (iAtriplet.size() == 1) ? 0
 576       : int( TINY + rndmPtr->flat() * (iAtriplet.size() - TINY) );
 577     int iPickA = iAtriplet[pickA];
 578     acols[iPickA] = -col0;
 579
 580     // Remove matched triplet and store new colour dipole ends.
 581     col0 = 0;
 582     iAtriplet[pickA] = iAtriplet.back();
 583     iAtriplet.pop_back();
 584     iDipCol.push_back(0);
 585     iDipAcol.push_back(iPickA);
 586   }
 587
 588   // Error checks that amount of leftover colours and anticolours match.
 589   if ( (iTriplet.size() != iAtriplet.size())
 590     || (col0 != 0 && acol0 == 0) || (col0 == 0 && acol0 != 0) ) {
 591     infoPtr->errorMsg("Error in ResonanceDecays::pickColours:"
 592       " inconsistent colour tags");
 593     return false;
 594   }
 595
 596   // Match triplets to antitriplets in the final state.
 597   for (int pickT = 0; pickT < int(iTriplet.size()); ++pickT) {
 598     int iPickT = iTriplet[pickT];
 599     int pickA  = (iAtriplet.size() == 1) ? 0
 600       : int( TINY + rndmPtr->flat() * (iAtriplet.size() - TINY) );
 601     int iPickA = iAtriplet[pickA];
 602
 603     // Connect pair with new colour tag.
 604     cols[iPickT]  = process.nextColTag();
 605     acols[iPickA] = cols[iPickT];
 606
 607     // Remove matched antitriplet and store new colour dipole ends.
 608     iAtriplet[pickT] = iAtriplet.back();
 609     iAtriplet.pop_back();
 610     iDipCol.push_back(iPickT);
 611     iDipAcol.push_back(iPickA);
 612   }
 613
 614   // If no octets are around then matching is done.
 615   if (col0 == 0 && acol0 == 0 && iOctet.size() == 0) return true;
 616
 617   // If initial-state octet remains then store as (first!) new dipole.
 618   if (col0 != 0) {
 619     iDipCol.push_back(0);
 620     iDipAcol.push_back(0);
 621   }
 622
 623   // Now attach all final-state octets at random to existing dipoles.
 624   for (int i = 0; i < int(iOctet.size()); ++i) {
 625     int iOct = iOctet[i];
 626
 627     // If no dipole then start new one. (Happens for singlet -> octets.)
 628     if (iDipCol.size() == 0) {
 629       cols[iOct]  = process.nextColTag();
 630       acols[iOct] = cols[iOct] ;
 631       iDipCol.push_back(iOct);
 632       iDipAcol.push_back(iOct);
 633     }
 634
 635     // Else attach to existing dipole picked at random.
 636     else {
 637       int pickDip = (iDipCol.size() == 1) ? 0
 638         : int( TINY + rndmPtr->flat() * (iDipCol.size() - TINY) );
 639
 640       // Case with dipole in initial state: reattach existing colours.
 641       if (iDipCol[pickDip] == 0 && iDipAcol[pickDip] == 0) {
 642         cols[iOct]        = col0;
 643         acols[iOct]       = acol0;
 644         iDipAcol[pickDip] = iOct;
 645         iDipCol.push_back(iOct);
 646         iDipAcol.push_back(0);
 647
 648       // Case with dipole from colour in initial state: also new colour.
 649       } else if (iDipAcol[pickDip] == 0) {
 650         int iPickCol      = iDipCol[pickDip];
 651         cols[iOct]        = cols[iPickCol];
 652         acols[iOct]       = process.nextColTag();
 653         cols[iPickCol]    = acols[iOct];
 654         iDipCol[pickDip]  = iOct;
 655         iDipCol.push_back(iPickCol);
 656         iDipAcol.push_back(iOct);
 657
 658       // Remaining cases with dipole from anticolour in initial state
 659       // or dipole inside final state: also new colour.
 660       } else {
 661         int iPickAcol     = iDipAcol[pickDip];
 662         acols[iOct]       = acols[iPickAcol];
 663         cols[iOct]        = process.nextColTag();
 664         acols[iPickAcol]  = cols[iOct];
 665         iDipAcol[pickDip] = iOct;
 666         iDipCol.push_back(iOct);
 667         iDipAcol.push_back(iPickAcol);
 668       }
 669     }
 670   }
 671
 672   // Must now have at least two dipoles (no 1 -> 8 or 8 -> 1).
 673   if (iDipCol.size() < 2) {
 674     infoPtr->errorMsg("Error in ResonanceDecays::pickColours:"
 675       " inconsistent colour tags");
 676     return false;
 677   }
 678
 679   // Done.
 680   return true;
 681
 682 }
 683
 684 //--------------------------------------------------------------------------
 685
 686 // Select decay products momenta isotropically in phase space.
 687 // Process-dependent angular distributions may be imposed in SigmaProcess.
 688
 689 bool ResonanceDecays::pickKinematics() {
 690
 691   // Description of two-body decays as simple special case.
 692   if (mult == 2) {
 693
 694     // Masses.
 695     m0          = mProd[0];
 696     double m1   = mProd[1];
 697     double m2   = mProd[2];
 698
 699     // Energies and absolute momentum in the rest frame.
 700     double e1   = 0.5 * (m0*m0 + m1*m1 - m2*m2) / m0;
 701     double e2   = 0.5 * (m0*m0 + m2*m2 - m1*m1) / m0;
 702     double pAbs = 0.5 * sqrtpos( (m0 - m1 - m2) * (m0 + m1 + m2)
 703       * (m0 + m1 - m2) * (m0 - m1 + m2) ) / m0;
 704
 705     // Pick isotropic angles to give three-momentum.
 706     double cosTheta = 2. * rndmPtr->flat() - 1.;
 707     double sinTheta = sqrt(1. - cosTheta*cosTheta);
 708     double phi      = 2. * M_PI * rndmPtr->flat();
 709     double pX       = pAbs * sinTheta * cos(phi);
 710     double pY       = pAbs * sinTheta * sin(phi);
 711     double pZ       = pAbs * cosTheta;
 712
 713     // Fill four-momenta in mother rest frame and then boost to lab frame.
 714     pProd.push_back( Vec4(  pX,  pY,  pZ, e1) );
 715     pProd.push_back( Vec4( -pX, -pY, -pZ, e2) );
 716     pProd[1].bst( pProd[0] );
 717     pProd[2].bst( pProd[0] );
 718
 719     // Done for two-body decay.
 720     return true;
 721   }
 722
 723   // Description of three-body decays as semi-simple special case.
 724   if (mult == 3) {
 725
 726     // Masses.
 727     m0             = mProd[0];
 728     double m1      = mProd[1];
 729     double m2      = mProd[2];
 730     double m3      = mProd[3];
 731     double mDiff   = m0 - (m1 + m2 + m3);
 732
 733     // Kinematical limits for 2+3 mass. Maximum phase-space weight.
 734     double m23Min  = m2 + m3;
 735     double m23Max  = m0 - m1;
 736     double p1Max   = 0.5 * sqrtpos( (m0 - m1 - m23Min) * (m0 + m1 + m23Min)
 737       * (m0 + m1 - m23Min) * (m0 - m1 + m23Min) ) / m0;
 738     double p23Max  = 0.5 * sqrtpos( (m23Max - m2 - m3) * (m23Max + m2 + m3)
 739       * (m23Max + m2 - m3) * (m23Max - m2 + m3) ) / m23Max;
 740     double wtPSmax = 0.5 * p1Max * p23Max;
 741
 742     // Pick an intermediate mass m23 flat in the allowed range.
 743     double wtPS, m23, p1Abs, p23Abs;
 744     do {
 745       m23 = m23Min + rndmPtr->flat() * mDiff;
 746
 747       // Translate into relative momenta and find phase-space weight.
 748       p1Abs  = 0.5 * sqrtpos( (m0 - m1 - m23) * (m0 + m1 + m23)
 749         * (m0 + m1 - m23) * (m0 - m1 + m23) ) / m0;
 750       p23Abs = 0.5 * sqrtpos( (m23 - m2 - m3) * (m23 + m2 + m3)
 751         * (m23 + m2 - m3) * (m23 - m2 + m3) ) / m23;
 752       wtPS   = p1Abs * p23Abs;
 753
 754     // If rejected, try again with new invariant masses.
 755     } while ( wtPS < rndmPtr->flat() * wtPSmax );
 756
 757     // Set up m23 -> m2 + m3 isotropic in its rest frame.
 758     double cosTheta = 2. * rndmPtr->flat() - 1.;
 759     double sinTheta = sqrt(1. - cosTheta*cosTheta);
 760     double phi      = 2. * M_PI * rndmPtr->flat();
 761     double pX       = p23Abs * sinTheta * cos(phi);
 762     double pY       = p23Abs * sinTheta * sin(phi);
 763     double pZ       = p23Abs * cosTheta;
 764     double e2       = sqrt( m2*m2 + p23Abs*p23Abs);
 765     double e3       = sqrt( m3*m3 + p23Abs*p23Abs);
 766     Vec4 p2(  pX,  pY,  pZ, e2);
 767     Vec4 p3( -pX, -pY, -pZ, e3);
 768
 769     // Set up 0 -> 1 + 23 isotropic in its rest frame.
 770     cosTheta        = 2. * rndmPtr->flat() - 1.;
 771     sinTheta        = sqrt(1. - cosTheta*cosTheta);
 772     phi             = 2. * M_PI * rndmPtr->flat();
 773     pX              = p1Abs * sinTheta * cos(phi);
 774     pY              = p1Abs * sinTheta * sin(phi);
 775     pZ              = p1Abs * cosTheta;
 776     double e1       = sqrt( m1*m1 + p1Abs*p1Abs);
 777     double e23      = sqrt( m23*m23 + p1Abs*p1Abs);
 778     pProd.push_back( Vec4( pX, pY, pZ, e1) );
 779
 780     // Boost 2 + 3 to the 0 rest frame and then boost to lab frame.
 781     Vec4 p23( -pX, -pY, -pZ, e23);
 782     p2.bst( p23 );
 783     p3.bst( p23 );
 784     pProd.push_back( p2 );
 785     pProd.push_back( p3 );
 786     pProd[1].bst( pProd[0] );
 787     pProd[2].bst( pProd[0] );
 788     pProd[3].bst( pProd[0] );
 789
 790     // Done for three-body decay.
 791     return true;
 792   }
 793
 794   // Do a multibody decay using the M-generator algorithm.
 795
 796   // Mother and sum daughter masses.
 797   m0             = mProd[0];
 798   double mSum    = mProd[1];
 799   for (int i = 2; i <= mult; ++i) mSum += mProd[i];
 800   double mDiff   = m0 - mSum;
 801
 802   // Begin setup of intermediate invariant masses.
 803   vector<double> mInv;
 804   for (int i = 0; i <= mult; ++i) mInv.push_back( mProd[i]);
 805
 806   // Calculate the maximum weight in the decay.
 807   double wtPSmax = 1. / WTCORRECTION[mult];
 808   double mMax    = mDiff + mProd[mult];
 809   double mMin    = 0.;
 810   for (int i = mult - 1; i > 0; --i) {
 811     mMax        += mProd[i];
 812     mMin        += mProd[i+1];
 813     double mNow  = mProd[i];
 814     wtPSmax *= 0.5 * sqrtpos( (mMax - mMin - mNow) * (mMax + mMin + mNow)
 815     * (mMax + mMin - mNow) * (mMax - mMin + mNow) ) / mMax;
 816   }
 817
 818   // Begin loop to find the set of intermediate invariant masses.
 819   vector<double> rndmOrd;
 820   double wtPS;
 821   do {
 822     wtPS  = 1.;
 823
 824     // Find and order random numbers in descending order.
 825     rndmOrd.resize(0);
 826     rndmOrd.push_back(1.);
 827     for (int i = 1; i < mult - 1; ++i) {
 828       double rndm = rndmPtr->flat();
 829       rndmOrd.push_back(rndm);
 830       for (int j = i - 1; j > 0; --j) {
 831         if (rndm > rndmOrd[j]) swap( rndmOrd[j], rndmOrd[j+1] );
 832         else break;
 833       }
 834     }
 835     rndmOrd.push_back(0.);
 836
 837     // Translate into intermediate masses and find weight.
 838     for (int i = mult - 1; i > 0; --i) {
 839       mInv[i] = mInv[i+1] + mProd[i] + (rndmOrd[i-1] - rndmOrd[i]) * mDiff;
 840       wtPS   *= 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
 841         * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
 842         * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
 843     }
 844
 845   // If rejected, try again with new invariant masses.
 846   } while ( wtPS < rndmPtr->flat() * wtPSmax );
 847
 848   // Perform two-particle decays in the respective rest frame.
 849   vector<Vec4> pInv;
 850   pInv.resize(mult + 1);
 851   for (int i = 1; i < mult; ++i) {
 852     double pAbs = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
 853       * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
 854       * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
 855
 856     // Isotropic angles give three-momentum.
 857     double cosTheta = 2. * rndmPtr->flat() - 1.;
 858     double sinTheta = sqrt(1. - cosTheta*cosTheta);
 859     double phi      = 2. * M_PI * rndmPtr->flat();
 860     double pX       = pAbs * sinTheta * cos(phi);
 861     double pY       = pAbs * sinTheta * sin(phi);
 862     double pZ       = pAbs * cosTheta;
 863
 864     // Calculate energies, fill four-momenta.
 865     double eHad     = sqrt( mProd[i]*mProd[i] + pAbs*pAbs);
 866     double eInv     = sqrt( mInv[i+1]*mInv[i+1] + pAbs*pAbs);
 867     pProd.push_back( Vec4( pX, pY, pZ, eHad) );
 868     pInv[i+1].p( -pX, -pY, -pZ, eInv);
 869   }
 870   pProd.push_back( pInv[mult] );
 871
 872   // Boost decay products to the mother rest frame and on to lab frame.
 873   pInv[1] = pProd[0];
 874   for (int iFrame = mult - 1; iFrame > 0; --iFrame)
 875     for (int i = iFrame; i <= mult; ++i) pProd[i].bst(pInv[iFrame]);
 876
 877   // Done for multibody decay.
 878   return true;
 879
 880 }
 881
 882 //==========================================================================
 883
 884 } // end namespace Pythia8