PYTHIA8/pythia8145/src/ParticleDecays.cxx

   1 // ParticleDecays.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 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
   7 // ParticleDecays class.
   8
   9 #include "ParticleDecays.h"
  10
  11 namespace Pythia8 {
  12
  13 //==========================================================================
  14
  15 // The ParticleDecays class.
  16
  17 //--------------------------------------------------------------------------
  18
  19 // Constants: could be changed here if desired, but normally should not.
  20 // These are of technical nature, as described for each.
  21
  22 // Number of times one tries to let decay happen (for 2 nested loops).
  23 const int    ParticleDecays::NTRYDECAY   = 10;
  24
  25 // Number of times one tries to pick valid hadronic content in decay.
  26 const int    ParticleDecays::NTRYPICK    = 100;
  27
  28 // Number of times one tries to pick decay topology.
  29 const int    ParticleDecays::NTRYMEWT    = 1000;
  30
  31 // Maximal loop count in Dalitz decay treatment.
  32 const int    ParticleDecays::NTRYDALITZ  = 1000;
  33
  34 // Minimal Dalitz pair mass this factor above threshold.
  35 const double ParticleDecays::MSAFEDALITZ = 1.000001;
  36
  37 // These numbers are hardwired empirical parameters,
  38 // intended to speed up the M-generator.
  39 const double ParticleDecays::WTCORRECTION[11] = { 1., 1., 1.,
  40   2., 5., 15., 60., 250., 1250., 7000., 50000. };
  41
  42 //--------------------------------------------------------------------------
  43
  44 // Initialize and save pointers.
  45
  46 void ParticleDecays::init(Info* infoPtrIn, Settings& settings,
  47   ParticleData* particleDataPtrIn, Rndm* rndmPtrIn,
  48   Couplings* couplingsPtrIn, TimeShower* timesDecPtrIn,
  49   StringFlav* flavSelPtrIn, DecayHandler* decayHandlePtrIn,
  50   vector<int> handledParticles) {
  51
  52   // Save pointers to error messages handling and flavour generation.
  53   infoPtr         = infoPtrIn;
  54   particleDataPtr = particleDataPtrIn;
  55   rndmPtr         = rndmPtrIn;
  56   couplingsPtr    = couplingsPtrIn;
  57   flavSelPtr      = flavSelPtrIn;
  58
  59   // Save pointer to timelike shower, as needed in some few decays.
  60   timesDecPtr     = timesDecPtrIn;
  61
  62   // Save pointer for external handling of some decays.
  63   decayHandlePtr  = decayHandlePtrIn;
  64
  65   // Set which particles should be handled externally.
  66   if (decayHandlePtr != 0)
  67   for (int i = 0; i < int(handledParticles.size()); ++i)
  68     particleDataPtr->doExternalDecay(handledParticles[i], true);
  69
  70   // Safety margin in mass to avoid troubles.
  71   mSafety       = settings.parm("ParticleDecays:mSafety");
  72
  73   // Lifetime and vertex rules for determining whether decay allowed.
  74   limitTau0     = settings.flag("ParticleDecays:limitTau0");
  75   tau0Max       = settings.parm("ParticleDecays:tau0Max");
  76   limitTau      = settings.flag("ParticleDecays:limitTau");
  77   tauMax        = settings.parm("ParticleDecays:tauMax");
  78   limitRadius   = settings.flag("ParticleDecays:limitRadius");
  79   rMax          = settings.parm("ParticleDecays:rMax");
  80   limitCylinder = settings.flag("ParticleDecays:limitCylinder");
  81   xyMax         = settings.parm("ParticleDecays:xyMax");
  82   zMax          = settings.parm("ParticleDecays:zMax");
  83   limitDecay    = limitTau0 || limitTau || limitRadius || limitCylinder;
  84
  85   // B-Bbar mixing parameters.
  86   mixB          = settings.flag("ParticleDecays:mixB");
  87   xBdMix        = settings.parm("ParticleDecays:xBdMix");
  88   xBsMix        = settings.parm("ParticleDecays:xBsMix");
  89
  90   // Suppression of extra-hadron momenta in semileptonic decays.
  91   sigmaSoft     = settings.parm("ParticleDecays:sigmaSoft");
  92
  93   // Selection of multiplicity and colours in "phase space" model.
  94   multIncrease  = settings.parm("ParticleDecays:multIncrease");
  95   multRefMass   = settings.parm("ParticleDecays:multRefMass");
  96   multGoffset   = settings.parm("ParticleDecays:multGoffset");
  97   colRearrange  = settings.parm("ParticleDecays:colRearrange");
  98
  99   // Minimum energy in system (+ m_q) from StringFragmentation.
 100   stopMass      = settings.parm("StringFragmentation:stopMass");
 101
 102   // Parameters for Dalitz decay virtual gamma mass spectrum.
 103   sRhoDal       = pow2(particleDataPtr->m0(113));
 104   wRhoDal       = pow2(particleDataPtr->mWidth(113));
 105
 106   // Allow showers in decays to qqbar/gg/ggg/gammagg.
 107   doFSRinDecays = settings.flag("ParticleDecays:FSRinDecays");
 108
 109 }
 110
 111 //--------------------------------------------------------------------------
 112
 113 // Decay a particle; main method.
 114
 115 bool ParticleDecays::decay( int iDec, Event& event) {
 116
 117   // Check whether a decay is allowed, given the upcoming decay vertex.
 118   Particle& decayer = event[iDec];
 119   hasPartons  = false;
 120   keepPartons = false;
 121   if (limitDecay && !checkVertex(decayer)) return true;
 122
 123   // Fill the decaying particle in slot 0 of arrays.
 124   idDec = decayer.id();
 125   iProd.resize(0);
 126   idProd.resize(0);
 127   mProd.resize(0);
 128   iProd.push_back( iDec );
 129   idProd.push_back( idDec );
 130   mProd.push_back( decayer.m() );
 131
 132   // Check for oscillations B0 <-> B0bar or B_s0 <-> B_s0bar.
 133   bool hasOscillated = (abs(idDec) == 511 || abs(idDec) == 531)
 134     ? oscillateB(decayer) : false;
 135   if (hasOscillated) {idDec = - idDec; idProd[0] = idDec;}
 136
 137   // Particle data for decaying particle.
 138   decDataPtr = &decayer.particleDataEntry();
 139
 140   // Optionally send on to external decay program.
 141   bool doneExternally = false;
 142   if (decDataPtr->doExternalDecay()) {
 143     pProd.resize(0);
 144     pProd.push_back(decayer.p());
 145     doneExternally = decayHandlePtr->decay(idProd, mProd, pProd,
 146       iDec, event);
 147
 148     // If it worked, then store the decay products in the event record.
 149     if (doneExternally) {
 150       mult = idProd.size() - 1;
 151       int status = (hasOscillated) ? 94 : 93;
 152       for (int i = 1; i <= mult; ++i) {
 153         int iPos = event.append( idProd[i], status, iDec, 0, 0, 0,
 154         0, 0, pProd[i], mProd[i]);
 155         iProd.push_back( iPos);
 156       }
 157
 158       // Also mark mother decayed and store daughters.
 159       event[iDec].statusNeg();
 160       event[iDec].daughters( iProd[1], iProd[mult]);
 161     }
 162   }
 163
 164   // Now begin normal internal decay treatment.
 165   if (!doneExternally) {
 166
 167     // Allow up to ten tries to pick a channel.
 168     if (!decDataPtr->preparePick(idDec)) return false;
 169     bool foundChannel = false;
 170     bool hasStored    = false;
 171     for (int iTryChannel = 0; iTryChannel < NTRYDECAY; ++iTryChannel) {
 172
 173       // Remove previous failed channel.
 174       if (hasStored) event.popBack(mult);
 175       hasStored = false;
 176
 177       // Pick new channel. Read out basics.
 178       DecayChannel& channel = decDataPtr->pickChannel();
 179       meMode = channel.meMode();
 180       keepPartons = (meMode > 90 && meMode <= 100);
 181       mult = channel.multiplicity();
 182
 183       // Allow up to ten tries for each channel (e.g with different masses).
 184       bool foundMode = false;
 185       for (int iTryMode = 0; iTryMode < NTRYDECAY; ++iTryMode) {
 186         idProd.resize(1);
 187         mProd.resize(1);
 188         scale = 0.;
 189
 190         // Extract and store the decay products in local arrays.
 191         hasPartons = false;
 192         for (int i = 0; i < mult; ++i) {
 193           int idNow = channel.product(i);
 194           int idAbs = abs(idNow);
 195           if ( idAbs < 10 || idAbs == 21 || idAbs == 81 || idAbs == 82
 196             || idAbs == 83 || (idAbs > 1000 && idAbs < 10000
 197             && (idAbs/10)%10 == 0) ) hasPartons = true;
 198           if (idDec < 0 && particleDataPtr->hasAnti(idNow)) idNow = -idNow;
 199           double mNow = particleDataPtr->mass(idNow);
 200           idProd.push_back( idNow);
 201           mProd.push_back( mNow);
 202         }
 203
 204         // Decays into partons usually translate into hadrons.
 205         if (hasPartons && !keepPartons && !pickHadrons()) continue;
 206
 207         // Need to set colour flow if explicit decay to partons.
 208         cols.resize(0);
 209         acols.resize(0);
 210         for (int i = 0; i <= mult; ++i) {
 211           cols.push_back(0);
 212           acols.push_back(0);
 213         }
 214         if (hasPartons && keepPartons && !setColours(event)) continue;
 215
 216         // Check that enough phase space for decay.
 217         if (mult > 1) {
 218           double mDiff = mProd[0];
 219           for (int i = 1; i <= mult; ++i) mDiff -= mProd[i];
 220           if (mDiff < mSafety) continue;
 221         }
 222
 223         // End of inner trial loops. Check if succeeded or not.
 224         foundMode = true;
 225         break;
 226       }
 227       if (!foundMode) continue;
 228
 229       // Store decay products in the event record.
 230       int status = (hasOscillated) ? 92 : 91;
 231       for (int i = 1; i <= mult; ++i) {
 232         int iPos = event.append( idProd[i], status, iDec, 0, 0, 0,
 233           cols[i], acols[i], Vec4(0., 0., 0., 0.), mProd[i], scale);
 234         iProd.push_back( iPos);
 235       }
 236       hasStored = true;
 237
 238       // Pick mass of Dalitz decay. Temporarily change multiplicity.
 239       if ( (meMode == 11 || meMode == 12 || meMode == 13)
 240         && !dalitzMass() ) continue;
 241
 242       // Do a decay, split by multiplicity.
 243       bool decayed = false;
 244       if      (mult == 1) decayed = oneBody(event);
 245       else if (mult == 2) decayed = twoBody(event);
 246       else if (mult == 3) decayed = threeBody(event);
 247       else                decayed = mGenerator(event);
 248       if (!decayed) continue;
 249
 250       // Kinematics of gamma* -> l- l+ in Dalitz decay. Restore multiplicity.
 251       if (meMode == 11 || meMode == 12 || meMode == 13)
 252         dalitzKinematics(event);
 253
 254       // End of outer trial loops.
 255       foundChannel = true;
 256       break;
 257     }
 258
 259     // If the decay worked, then mark mother decayed and store daughters.
 260     if (foundChannel) {
 261       event[iDec].statusNeg();
 262       event[iDec].daughters( iProd[1], iProd[mult]);
 263
 264     // Else remove unused daughters and return failure.
 265     } else {
 266       if (hasStored) event.popBack(mult);
 267       infoPtr->errorMsg("Error in ParticleDecays::decay: "
 268         "failed to find workable decay channel");
 269       return false;
 270     }
 271
 272   // Now finished normal internal decay treatment.
 273   }
 274
 275   // Set decay vertex when this is displaced.
 276   if (event[iDec].hasVertex() || event[iDec].tau() > 0.) {
 277     Vec4 vDec = event[iDec].vDec();
 278     for (int i = 1; i <= mult; ++i) event[iProd[i]].vProd( vDec );
 279   }
 280
 281   // Set lifetime of daughters.
 282   for (int i = 1; i <= mult; ++i)
 283     event[iProd[i]].tau( event[iProd[i]].tau0() * rndmPtr->exp() );
 284
 285   // In a decay explicitly to partons then optionally do a shower,
 286   // and always flag that partonic system should be fragmented.
 287   if (hasPartons && keepPartons && doFSRinDecays)
 288     timesDecPtr->shower( iProd[1], iProd.back(), event, mProd[0]);
 289
 290   // Done.
 291   return true;
 292
 293 }
 294
 295 //--------------------------------------------------------------------------
 296
 297 // Check whether a decay is allowed, given the upcoming decay vertex.
 298
 299 bool ParticleDecays::checkVertex(Particle& decayer) {
 300
 301   // Check whether any of the conditions are not fulfilled.
 302   if (limitTau0 && decayer.tau0() > tau0Max) return false;
 303   if (limitTau && decayer.tau() > tauMax) return false;
 304   if (limitRadius && pow2(decayer.xDec()) + pow2(decayer.yDec())
 305     + pow2(decayer.zDec()) > pow2(rMax)) return false;
 306   if (limitCylinder && (pow2(decayer.xDec()) + pow2(decayer.yDec())
 307     > pow2(xyMax) || abs(decayer.zDec()) > zMax) ) return false;
 308
 309   // Done.
 310   return true;
 311
 312 }
 313
 314 //--------------------------------------------------------------------------
 315
 316 // Check for oscillations B0 <-> B0bar or B_s0 <-> B_s0bar.
 317
 318 bool ParticleDecays::oscillateB(Particle& decayer) {
 319
 320   // Extract relevant information and decide.
 321   double xBmix   = (abs(decayer.id()) == 511) ? xBdMix : xBsMix;
 322   double tau     = decayer.tau();
 323   double tau0    = decayer.tau0();
 324   double probosc = pow2(sin(0.5 * xBmix * tau / tau0));
 325   return (probosc > rndmPtr->flat());
 326
 327 }
 328
 329 //--------------------------------------------------------------------------
 330
 331 // Do a one-body decay. (Rare; e.g. for K0 -> K0_short.)
 332
 333 bool ParticleDecays::oneBody(Event& event) {
 334
 335   // References to the particles involved.
 336   Particle& decayer = event[iProd[0]];
 337   Particle& prod    = event[iProd[1]];
 338
 339   // Set momentum and expand mother information.
 340   prod.p( decayer.p() );
 341   prod.m( decayer.m() );
 342   prod.mother2( iProd[0] );
 343
 344   // Done.
 345   return true;
 346
 347 }
 348
 349 //--------------------------------------------------------------------------
 350
 351 // Do a two-body decay.
 352
 353 bool ParticleDecays::twoBody(Event& event) {
 354
 355   // References to the particles involved.
 356   Particle& decayer = event[iProd[0]];
 357   Particle& prod1   = event[iProd[1]];
 358   Particle& prod2   = event[iProd[2]];
 359
 360   // Masses.
 361   double m0   = mProd[0];
 362   double m1   = mProd[1];
 363   double m2   = mProd[2];
 364
 365   // Energies and absolute momentum in the rest frame.
 366   if (m1 + m2 + mSafety > m0) return false;
 367   double e1   = 0.5 * (m0*m0 + m1*m1 - m2*m2) / m0;
 368   double e2   = 0.5 * (m0*m0 + m2*m2 - m1*m1) / m0;
 369   double pAbs = 0.5 * sqrtpos( (m0 - m1 - m2) * (m0 + m1 + m2)
 370     * (m0 + m1 - m2) * (m0 - m1 + m2) ) / m0;
 371
 372   // When meMode = 2, for V -> PS2 + PS3 (V = vector, pseudoscalar),
 373   // need to check if production is PS0 -> PS1/gamma + V.
 374   int iMother = event[iProd[0]].mother1();
 375   int idSister = 0;
 376   if (meMode == 2) {
 377     if (iMother <= 0 || iMother >= iProd[0]) meMode = 0;
 378     else {
 379       int iDaughter1 = event[iMother].daughter1();
 380       int iDaughter2 = event[iMother].daughter2();
 381       if (iDaughter2 != iDaughter1 + 1) meMode = 0;
 382       else {
 383         int idMother = abs( event[iMother].id() );
 384         if (idMother <= 100 || idMother%10 !=1
 385           || (idMother/1000)%10 != 0) meMode = 0;
 386         else {
 387           int iSister = (iProd[0] == iDaughter1) ? iDaughter2 : iDaughter1;
 388           idSister = abs( event[iSister].id() );
 389           if ( (idSister <= 100 || idSister%10 !=1
 390             || (idSister/1000)%10 != 0) && idSister != 22) meMode = 0;
 391         }
 392       }
 393     }
 394   }
 395
 396   // Begin loop over matrix-element corrections.
 397   double wtME, wtMEmax;
 398   int loop = 0;
 399   do {
 400     wtME = 1.;
 401     wtMEmax = 1.;
 402     ++loop;
 403
 404     // Isotropic angles give three-momentum.
 405     double cosTheta = 2. * rndmPtr->flat() - 1.;
 406     double sinTheta = sqrt(1. - cosTheta*cosTheta);
 407     double phi      = 2. * M_PI * rndmPtr->flat();
 408     double pX       = pAbs * sinTheta * cos(phi);
 409     double pY       = pAbs * sinTheta * sin(phi);
 410     double pZ       = pAbs * cosTheta;
 411
 412     // Fill four-momenta and boost them away from mother rest frame.
 413     prod1.p(  pX,  pY,  pZ, e1);
 414     prod2.p( -pX, -pY, -pZ, e2);
 415     prod1.bst( decayer.p(), decayer.m() );
 416     prod2.bst( decayer.p(), decayer.m() );
 417
 418     // Matrix element for PS0 -> PS1 + V1 -> PS1 + PS2 + PS3 of form
 419     // cos**2(theta02) in V1 rest frame, and for PS0 -> gamma + V1
 420     // -> gamma + PS2 + PS3 of form sin**2(theta02).
 421     if (meMode == 2) {
 422       double p10 = decayer.p() * event[iMother].p();
 423       double p12 = decayer.p() * prod1.p();
 424       double p02 = event[iMother].p() * prod1.p();
 425       double s0  = pow2(event[iMother].m());
 426       double s1  = pow2(decayer.m());
 427       double s2  =  pow2(prod1.m());
 428       if (idSister != 22) wtME = pow2(p10 * p12 - s1 * p02);
 429       else wtME = s1 * (2. * p10 * p12 * p02 - s1 * p02*p02
 430         - s0 * p12*p12 - s2 * p10*p10 + s1 * s0 * s2);
 431       wtME = max( wtME, 1e-6 * s1*s1 * s0 * s2);
 432       wtMEmax = (p10*p10 - s1 * s0) * (p12*p12 - s1 * s2);
 433     }
 434
 435     // Break out of loop if no sensible ME weight.
 436     if(loop > NTRYMEWT) {
 437       infoPtr->errorMsg("ParticleDecays::twoBody: "
 438         "caught in infinite ME weight loop");
 439       wtME = abs(wtMEmax);
 440     }
 441
 442   // If rejected, try again with new invariant masses.
 443   } while ( wtME < rndmPtr->flat() * wtMEmax );
 444
 445   // Done.
 446   return true;
 447
 448 }
 449
 450 //--------------------------------------------------------------------------
 451
 452 // Do a three-body decay (except Dalitz decays).
 453
 454 bool ParticleDecays::threeBody(Event& event) {
 455
 456   // References to the particles involved.
 457   Particle& decayer = event[iProd[0]];
 458   Particle& prod1   = event[iProd[1]];
 459   Particle& prod2   = event[iProd[2]];
 460   Particle& prod3   = event[iProd[3]];
 461
 462   // Mother and sum daughter masses. Fail if too close.
 463   double m0      = mProd[0];
 464   double m1      = mProd[1];
 465   double m2      = mProd[2];
 466   double m3      = mProd[3];
 467   double mSum    = m1 + m2 + m3;
 468   double mDiff   = m0 - mSum;
 469   if (mDiff < mSafety) return false;
 470
 471   // Kinematical limits for 2+3 mass. Maximum phase-space weight.
 472   double m23Min  = m2 + m3;
 473   double m23Max  = m0 - m1;
 474   double p1Max   = 0.5 * sqrtpos( (m0 - m1 - m23Min) * (m0 + m1 + m23Min)
 475     * (m0 + m1 - m23Min) * (m0 - m1 + m23Min) ) / m0;
 476   double p23Max  = 0.5 * sqrtpos( (m23Max - m2 - m3) * (m23Max + m2 + m3)
 477     * (m23Max + m2 - m3) * (m23Max - m2 + m3) ) / m23Max;
 478   double wtPSmax = 0.5 * p1Max * p23Max;
 479
 480   // Begin loop over matrix-element corrections.
 481   double wtME, wtMEmax, wtPS, m23, p1Abs, p23Abs;
 482   do {
 483     wtME     = 1.;
 484     wtMEmax  = 1.;
 485
 486     // Pick an intermediate mass m23 flat in the allowed range.
 487     do {
 488       m23    = m23Min + rndmPtr->flat() * mDiff;
 489
 490       // Translate into relative momenta and find phase-space weight.
 491       p1Abs  = 0.5 * sqrtpos( (m0 - m1 - m23) * (m0 + m1 + m23)
 492         * (m0 + m1 - m23) * (m0 - m1 + m23) ) / m0;
 493       p23Abs = 0.5 * sqrtpos( (m23 - m2 - m3) * (m23 + m2 + m3)
 494         * (m23 + m2 - m3) * (m23 - m2 + m3) ) / m23;
 495       wtPS   = p1Abs * p23Abs;
 496
 497     // If rejected, try again with new invariant masses.
 498     } while ( wtPS < rndmPtr->flat() * wtPSmax );
 499
 500     // Set up m23 -> m2 + m3 isotropic in its rest frame.
 501     double cosTheta = 2. * rndmPtr->flat() - 1.;
 502     double sinTheta = sqrt(1. - cosTheta*cosTheta);
 503     double phi      = 2. * M_PI * rndmPtr->flat();
 504     double pX       = p23Abs * sinTheta * cos(phi);
 505     double pY       = p23Abs * sinTheta * sin(phi);
 506     double pZ       = p23Abs * cosTheta;
 507     double e2       = sqrt( m2*m2 + p23Abs*p23Abs);
 508     double e3       = sqrt( m3*m3 + p23Abs*p23Abs);
 509     prod2.p(  pX,  pY,  pZ, e2);
 510     prod3.p( -pX, -pY, -pZ, e3);
 511
 512     // Set up m0 -> m1 + m23 isotropic in its rest frame.
 513     cosTheta        = 2. * rndmPtr->flat() - 1.;
 514     sinTheta        = sqrt(1. - cosTheta*cosTheta);
 515     phi             = 2. * M_PI * rndmPtr->flat();
 516     pX              = p1Abs * sinTheta * cos(phi);
 517     pY              = p1Abs * sinTheta * sin(phi);
 518     pZ              = p1Abs * cosTheta;
 519     double e1       = sqrt( m1*m1 + p1Abs*p1Abs);
 520     double e23      = sqrt( m23*m23 + p1Abs*p1Abs);
 521     prod1.p( pX, pY, pZ, e1);
 522
 523     // Boost 2 + 3 to the 0 rest frame.
 524     Vec4 p23( -pX, -pY, -pZ, e23);
 525     prod2.bst( p23, m23 );
 526     prod3.bst( p23, m23 );
 527
 528     // Matrix-element weight for omega/phi -> pi+ pi- pi0.
 529     if (meMode == 1) {
 530       double p1p2 = prod1.p() * prod2.p();
 531       double p1p3 = prod1.p() * prod3.p();
 532       double p2p3 = prod2.p() * prod3.p();
 533       wtME = pow2(m1 * m2 * m3) - pow2(m1 * p2p3) - pow2(m2 * p1p3)
 534         - pow2(m3 * p1p2) + 2. * p1p2 * p1p3 * p2p3;
 535       wtMEmax = pow3(m0 * m0) / 150.;
 536
 537     // Effective matrix element for nu spectrum in tau -> nu + hadrons.
 538     } else if (meMode == 21) {
 539       double x1 = 2. *  prod1.e() / m0;
 540       wtME = x1 * (3. - 2. * x1);
 541       double xMax = min( 0.75, 2. * (1. - mSum / m0) );
 542       wtMEmax = xMax * (3. - 2. * xMax);
 543
 544     // Matrix element for weak decay (only semileptonic for c and b).
 545     } else if ((meMode == 22 || meMode == 23) && prod1.isLepton()) {
 546       wtME = m0 * prod1.e() * (prod2.p() * prod3.p());
 547       wtMEmax = min( pow4(m0) / 16., m0 * (m0 - m1 - m2) * (m0 - m1 - m3)
 548         * (m0 - m2 - m3) );
 549
 550     // Effective matrix element for weak decay to hadrons (B -> D, D -> K).
 551     } else if (meMode == 22 || meMode == 23) {
 552       double x1 = 2. * prod1.pAbs() / m0;
 553       wtME = x1 * (3. - 2. * x1);
 554       double xMax = min( 0.75, 2. * (1. - mSum / m0) );
 555       wtMEmax = xMax * (3. - 2. * xMax);
 556
 557     // Effective matrix element for gamma spectrum in B -> gamma + hadrons.
 558     } else if (meMode == 31) {
 559       double x1 = 2. * prod1.e() / m0;
 560       wtME = pow3(x1);
 561       double x1Max = 1. - pow2(mSum / m0);
 562       wtMEmax = pow3(x1Max);
 563
 564     // Matrix-element weight for "onium" -> g + g + g or gamma + g + g.
 565     } else if (meMode == 92) {
 566       double x1 = 2. * prod1.e() / m0;
 567       double x2 = 2. * prod2.e() / m0;
 568       double x3 = 2. * prod3.e() / m0;
 569       wtME = pow2( (1. - x1) / (x2 * x3) ) + pow2( (1. - x2) / (x1 * x3) )
 570         + pow2( (1. - x3) / (x1 * x2) );
 571       wtMEmax = 2.;
 572       // For gamma + g + g require minimum mass for g + g system.
 573       if (prod1.id() == 22 && sqrt(1. - x1) * m0 < 2. * stopMass) wtME = 0.;
 574       if (prod2.id() == 22 && sqrt(1. - x2) * m0 < 2. * stopMass) wtME = 0.;
 575       if (prod3.id() == 22 && sqrt(1. - x3) * m0 < 2. * stopMass) wtME = 0.;
 576     }
 577
 578   // If rejected, try again with new invariant masses.
 579   } while ( wtME < rndmPtr->flat() * wtMEmax );
 580
 581   // Boost 1 + 2 + 3 to the current frame.
 582   prod1.bst( decayer.p(), decayer.m() );
 583   prod2.bst( decayer.p(), decayer.m() );
 584   prod3.bst( decayer.p(), decayer.m() );
 585
 586   // Done.
 587   return true;
 588
 589 }
 590
 591 //--------------------------------------------------------------------------
 592
 593 // Do a multibody decay using the M-generator algorithm.
 594
 595 bool ParticleDecays::mGenerator(Event& event) {
 596
 597   // Mother and sum daughter masses. Fail if too close or inconsistent.
 598   double m0      = mProd[0];
 599   double mSum    = mProd[1];
 600   for (int i = 2; i <= mult; ++i) mSum += mProd[i];
 601   double mDiff   = m0 - mSum;
 602   if (mDiff < mSafety) return false;
 603
 604   // Begin setup of intermediate invariant masses.
 605   mInv.resize(0);
 606   for (int i = 0; i <= mult; ++i) mInv.push_back( mProd[i]);
 607
 608   // Calculate the maximum weight in the decay.
 609   double wtPS, wtME, wtMEmax;
 610   double wtPSmax = 1. / WTCORRECTION[mult];
 611   double mMax    = mDiff + mProd[mult];
 612   double mMin    = 0.;
 613   for (int i = mult - 1; i > 0; --i) {
 614     mMax        += mProd[i];
 615     mMin        += mProd[i+1];
 616     double mNow  = mProd[i];
 617     wtPSmax *= 0.5 * sqrtpos( (mMax - mMin - mNow) * (mMax + mMin + mNow)
 618     * (mMax + mMin - mNow) * (mMax - mMin + mNow) ) / mMax;
 619   }
 620
 621   // Begin loop over matrix-element corrections.
 622   do {
 623     wtME    = 1.;
 624     wtMEmax = 1.;
 625
 626     // Begin loop to find the set of intermediate invariant masses.
 627     do {
 628       wtPS  = 1.;
 629
 630       // Find and order random numbers in descending order.
 631       rndmOrd.resize(0);
 632       rndmOrd.push_back(1.);
 633       for (int i = 1; i < mult - 1; ++i) {
 634         double rndm = rndmPtr->flat();
 635         rndmOrd.push_back(rndm);
 636         for (int j = i - 1; j > 0; --j) {
 637           if (rndm > rndmOrd[j]) swap( rndmOrd[j], rndmOrd[j+1] );
 638           else break;
 639         }
 640       }
 641       rndmOrd.push_back(0.);
 642
 643       // Translate into intermediate masses and find weight.
 644       for (int i = mult - 1; i > 0; --i) {
 645         mInv[i] = mInv[i+1] + mProd[i] + (rndmOrd[i-1] - rndmOrd[i]) * mDiff;
 646         wtPS   *= 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
 647           * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
 648           * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
 649       }
 650
 651     // If rejected, try again with new invariant masses.
 652     } while ( wtPS < rndmPtr->flat() * wtPSmax );
 653
 654     // Perform two-particle decays in the respective rest frame.
 655     pInv.resize(mult + 1);
 656     for (int i = 1; i < mult; ++i) {
 657       double pAbs = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
 658         * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
 659         * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
 660
 661       // Isotropic angles give three-momentum.
 662       double cosTheta = 2. * rndmPtr->flat() - 1.;
 663       double sinTheta = sqrt(1. - cosTheta*cosTheta);
 664       double phi      = 2. * M_PI * rndmPtr->flat();
 665       double pX       = pAbs * sinTheta * cos(phi);
 666       double pY       = pAbs * sinTheta * sin(phi);
 667       double pZ       = pAbs * cosTheta;
 668
 669       // Calculate energies, fill four-momenta.
 670       double eHad     = sqrt( mProd[i]*mProd[i] + pAbs*pAbs);
 671       double eInv     = sqrt( mInv[i+1]*mInv[i+1] + pAbs*pAbs);
 672       event[iProd[i]].p( pX, pY, pZ, eHad);
 673       pInv[i+1].p( -pX, -pY, -pZ, eInv);
 674     }
 675
 676     // Boost decay products to the mother rest frame.
 677     event[iProd[mult]].p( pInv[mult] );
 678     for (int iFrame = mult - 1; iFrame > 1; --iFrame)
 679       for (int i = iFrame; i <= mult; ++i)
 680         event[iProd[i]].bst( pInv[iFrame], mInv[iFrame]);
 681
 682     // Effective matrix element for nu spectrum in tau -> nu + hadrons.
 683     if (meMode == 21 && event[iProd[1]].isLepton()) {
 684       double x1 = 2. * event[iProd[1]].e() / m0;
 685       wtME = x1 * (3. - 2. * x1);
 686       double xMax = min( 0.75, 2. * (1. - mSum / m0) );
 687       wtMEmax = xMax * (3. - 2. * xMax);
 688
 689     // Effective matrix element for weak decay (only semileptonic for c and b).
 690     // Particles 4 onwards should be made softer explicitly?
 691     } else if ((meMode == 22 || meMode == 23) && event[iProd[1]].isLepton()) {
 692       Vec4 pRest = event[iProd[3]].p();
 693       for (int i = 4; i <= mult; ++i) pRest += event[iProd[i]].p();
 694       wtME = m0 * event[iProd[1]].e() * (event[iProd[2]].p() * pRest);
 695       for (int i = 4; i <= mult; ++i) wtME
 696         *= exp(- event[iProd[i]].pAbs2() / pow2(sigmaSoft) );
 697       wtMEmax = pow4(m0) / 16.;
 698
 699     // Effective matrix element for weak decay to hadrons (B -> D, D -> K).
 700     } else if (meMode == 22 || meMode == 23) {
 701       double x1 = 2. * event[iProd[1]].pAbs() / m0;
 702       wtME = x1 * (3. - 2. * x1);
 703       double xMax = min( 0.75, 2. * (1. - mSum / m0) );
 704       wtMEmax = xMax * (3. - 2. * xMax);
 705
 706     // Effective matrix element for gamma spectrum in B -> gamma + hadrons.
 707     } else if (meMode == 31) {
 708       double x1 = 2. * event[iProd[1]].e() / m0;
 709       wtME = pow3(x1);
 710       double x1Max = 1. - pow2(mSum / m0);
 711       wtMEmax = pow3(x1Max);
 712     }
 713
 714   // If rejected, try again with new invariant masses.
 715   } while ( wtME < rndmPtr->flat() * wtMEmax );
 716
 717   // Boost decay products to the current frame.
 718   pInv[1].p( event[iProd[0]].p() );
 719   for (int i = 1; i <= mult; ++i) event[iProd[i]].bst( pInv[1], mInv[1] );
 720
 721   // Done.
 722   return true;
 723
 724 }
 725
 726 //--------------------------------------------------------------------------
 727
 728 // Select mass of lepton pair in a Dalitz decay.
 729
 730 bool ParticleDecays::dalitzMass() {
 731
 732   // Mother and sum daughter masses.
 733   double mSum1 = 0;
 734   for (int i = 1; i <= mult - 2; ++i) mSum1 += mProd[i];
 735   if (meMode == 13) mSum1 *= MSAFEDALITZ;
 736   double mSum2 = MSAFEDALITZ * (mProd[mult -1] + mProd[mult]);
 737   double mDiff = mProd[0] - mSum1 - mSum2;
 738
 739   // Fail if too close or inconsistent.
 740   if (mDiff < mSafety) return false;
 741   if (idProd[mult - 1] + idProd[mult] != 0
 742     || mProd[mult - 1] != mProd[mult]) {
 743     infoPtr->errorMsg("Error in ParticleDecays::dalitzMass:"
 744     " inconsistent flavour/mass assignments");
 745     return false;
 746   }
 747   if ( meMode == 13 && (idProd[1] + idProd[2] != 0
 748     || mProd[1] != mProd[2]) ) {
 749     infoPtr->errorMsg("Error in ParticleDecays::dalitzMass:"
 750     " inconsistent flavour/mass assignments");
 751     return false;
 752   }
 753
 754   // Case 1: one Dalitz pair.
 755   if (meMode == 11 || meMode == 12) {
 756
 757     // Kinematical limits for gamma* squared mass.
 758     double sGamMin = pow2(mSum2);
 759     double sGamMax = pow2(mProd[0] - mSum1);
 760     // Select virtual gamma squared mass. Guessed form for meMode == 12.
 761     double sGam, wtGam;
 762     int loop = 0;
 763     do {
 764       if (++loop > NTRYDALITZ) return false;
 765       sGam = sGamMin * pow( sGamMax / sGamMin, rndmPtr->flat() );
 766       wtGam = (1. + 0.5 * sGamMin / sGam) *  sqrt(1. - sGamMin / sGam)
 767         * pow3(1. - sGam / sGamMax) * sRhoDal * (sRhoDal + wRhoDal)
 768         / ( pow2(sGam - sRhoDal) + sRhoDal * wRhoDal );
 769     } while ( wtGam < rndmPtr->flat() );
 770
 771     // Store results in preparation for doing a one-less-body decay.
 772     --mult;
 773     mProd[mult] = sqrt(sGam);
 774
 775   // Case 2: two Dalitz pairs.
 776   } else {
 777
 778     // Kinematical limits for 1 + 2 and 3 + 4 gamma* masses.
 779     double s0 = pow2(mProd[0]);
 780     double s12Min = pow2(mSum1);
 781     double s12Max = pow2(mProd[0] - mSum2);
 782     double s34Min = pow2(mSum2);
 783     double s34Max = pow2(mProd[0] - mSum1);
 784
 785     // Select virtual gamma squared masses. Guessed form for meMode == 13.
 786     double s12, s34, wt12, wt34, wtPAbs, wtAll;
 787     int loop = 0;
 788     do {
 789       if (++loop > NTRYDALITZ) return false;
 790       s12 = s12Min * pow( s12Max / s12Min, rndmPtr->flat() );
 791       wt12 = (1. + 0.5 * s12Min / s12) *  sqrt(1. - s12Min / s12)
 792         * sRhoDal * (sRhoDal + wRhoDal)
 793         / ( pow2(s12 - sRhoDal) + sRhoDal * wRhoDal );
 794       s34 = s34Min * pow( s34Max / s34Min, rndmPtr->flat() );
 795       wt34 = (1. + 0.5 * s34Min / s34) *  sqrt(1. - s34Min / s34)
 796         * sRhoDal * (sRhoDal + wRhoDal)
 797         / ( pow2(s34 - sRhoDal) + sRhoDal * wRhoDal );
 798       wtPAbs = sqrtpos( pow2(1. - (s12 + s34)/ s0)
 799         - 4. * s12 * s34 / (s0 * s0) );
 800       wtAll = wt12 * wt34 * pow3(wtPAbs);
 801       if (wtAll > 1.) infoPtr->errorMsg(
 802         "Error in ParticleDecays::dalitzMass: weight > 1");
 803     } while (wtAll < rndmPtr->flat());
 804
 805     // Store results in preparation for doing a two-body decay.
 806     mult = 2;
 807     mProd[1] = sqrt(s12);
 808     mProd[2] = sqrt(s34);
 809   }
 810
 811   // Done.
 812   return true;
 813
 814 }
 815
 816 //--------------------------------------------------------------------------
 817
 818 // Do kinematics of gamma* -> l- l+ in Dalitz decay.
 819
 820 bool ParticleDecays::dalitzKinematics(Event& event) {
 821
 822   // Restore multiplicity.
 823   int nDal = (meMode < 13) ? 1 : 2;
 824   mult += nDal;
 825
 826   // Loop over one or two lepton pairs.
 827   for (int iDal = 0; iDal < nDal; ++iDal) {
 828
 829     // References to the particles involved.
 830     Particle& decayer = event[iProd[0]];
 831     Particle& prodA = (iDal == 0) ? event[iProd[mult - 1]]
 832       : event[iProd[1]];
 833     Particle& prodB = (iDal == 0) ? event[iProd[mult]]
 834       : event[iProd[2]];
 835
 836     // Reconstruct required rotations and boosts backwards.
 837     Vec4 pDec    = decayer.p();
 838     int  iGam    = (meMode < 13) ? mult - 1 : 2 - iDal;
 839     Vec4 pGam    = event[iProd[iGam]].p();
 840     pGam.bstback( pDec, decayer.m() );
 841     double phiGam = pGam.phi();
 842     pGam.rot( 0., -phiGam);
 843     double thetaGam = pGam.theta();
 844     pGam.rot( -thetaGam, 0.);
 845
 846     // Masses and phase space in gamma* rest frame.
 847     double mGam     = (meMode < 13) ? mProd[mult - 1] : mProd[2 - iDal];
 848     double mA       = prodA.m();
 849     double mB       = prodB.m();
 850     double mGamMin  = MSAFEDALITZ * (mA + mB);
 851     double mGamRat  = pow2(mGamMin / mGam);
 852     double pGamAbs  = 0.5 * sqrtpos( (mGam - mA - mB) * (mGam + mA + mB) );
 853
 854     // Set up decay in gamma* rest frame, reference along +z axis.
 855     double cosTheta, cos2Theta;
 856     do {
 857       cosTheta      = 2. * rndmPtr->flat() - 1.;
 858       cos2Theta     = cosTheta * cosTheta;
 859     } while ( 1. + cos2Theta + mGamRat * (1. - cos2Theta)
 860       < 2. * rndmPtr->flat() );
 861     double sinTheta = sqrt(1. - cosTheta*cosTheta);
 862     double phi      = 2. * M_PI * rndmPtr->flat();
 863     double pX       = pGamAbs * sinTheta * cos(phi);
 864     double pY       = pGamAbs * sinTheta * sin(phi);
 865     double pZ       = pGamAbs * cosTheta;
 866     double eA       = sqrt( mA*mA + pGamAbs*pGamAbs);
 867     double eB       = sqrt( mB*mB + pGamAbs*pGamAbs);
 868     prodA.p(  pX,  pY,  pZ, eA);
 869     prodB.p( -pX, -pY, -pZ, eB);
 870
 871     // Boost to lab frame.
 872     prodA.bst( pGam, mGam);
 873     prodB.bst( pGam, mGam);
 874     prodA.rot( thetaGam, phiGam);
 875     prodB.rot( thetaGam, phiGam);
 876     prodA.bst( pDec, decayer.m() );
 877     prodB.bst( pDec, decayer.m() );
 878   }
 879
 880   // Done.
 881   return true;
 882
 883 }
 884
 885 //--------------------------------------------------------------------------
 886
 887 // Translate a partonic content into a set of actual hadrons.
 888
 889 bool ParticleDecays::pickHadrons() {
 890
 891   // Find partonic decay products. Rest are known id's, mainly hadrons,
 892   // when necessary shuffled to beginning of idProd list.
 893   idPartons.resize(0);
 894   int nPartons = 0;
 895   int nKnown = 0;
 896   bool closedGLoop = false;
 897   for (int i = 1; i <= mult; ++i) {
 898     int idAbs = abs(idProd[i]);
 899     if ( idAbs < 9 || (idAbs > 1000 && idAbs < 10000 && (idAbs/10)%10 == 0)
 900       || idAbs == 81 || idAbs == 82 || idAbs == 83) {
 901       ++nPartons;
 902       idPartons.push_back(idProd[i]);
 903       if (idAbs == 83) closedGLoop = true;
 904     } else {
 905       ++nKnown;
 906       if (nPartons > 0) {
 907         idProd[nKnown] = idProd[i];
 908         mProd[nKnown] = mProd[i];
 909       }
 910     }
 911   }
 912
 913   // Replace generic spectator flavour code by the actual one.
 914   for (int i = 0; i < nPartons; ++i) {
 915     int idPart = idPartons[i];
 916     int idNew = idPart;
 917     if (idPart == 81) {
 918       int idAbs = abs(idDec);
 919       if ( (idAbs/1000)%10 == 0 ) {
 920         idNew = -(idAbs/10)%10;
 921         if ((idAbs/100)%2 == 1) idNew = -idNew;
 922       } else if ( (idAbs/100)%10 >= (idAbs/10)%10 )
 923         idNew = 100 * ((idAbs/10)%100) + 3;
 924       else idNew = 1000 * ((idAbs/10)%10) + 100 * ((idAbs/100)%10) + 1;
 925       if (idDec < 0) idNew = -idNew;
 926
 927     // Replace generic random flavour by a randomly selected one.
 928     } else if (idPart == 82 || idPart == 83) {
 929       double mFlav;
 930       do {
 931         int idDummy = -flavSelPtr->pickLightQ();
 932         FlavContainer flavDummy(idDummy, idPart - 82);
 933         do idNew = flavSelPtr->pick(flavDummy).id;
 934         while (idNew == 0);
 935         mFlav = particleDataPtr->constituentMass(idNew);
 936       } while (2. * mFlav + stopMass > mProd[0]);
 937     } else if (idPart == -82 || idPart == -83) {
 938       idNew = -idPartons[i-1];
 939     }
 940     idPartons[i] = idNew;
 941   }
 942
 943   // Determine whether fixed multiplicity or to be selected at random.
 944   int nMin = max( 2, nKnown + nPartons / 2);
 945   if (meMode == 23) nMin = 3;
 946   if (meMode > 41 && meMode <= 50) nMin = meMode - 40;
 947   if (meMode > 51 && meMode <= 60) nMin = meMode - 50;
 948   int nFix = 0;
 949   if (meMode == 0) nFix = nMin;
 950   if (meMode == 11) nFix = 3;
 951   if (meMode == 12) nFix = 4;
 952   if (meMode > 61 && meMode <= 70) nFix = meMode - 60;
 953   if (meMode > 71 && meMode <= 80) nFix = meMode - 70;
 954   if (nFix > 0 && nKnown + nPartons/2 > nFix) return false;
 955
 956   // Initial values for loop to set new hadronic content.
 957   int nFilled = nKnown + 1;
 958   int nTotal, nNew, nSpec, nLeft;
 959   double mDiff;
 960   int nTry = 0;
 961   bool diquarkClash = false;
 962   bool usedChannel  = false;
 963
 964   // Begin loop; interrupt if multiple tries fail.
 965   do {
 966     ++nTry;
 967     if (nTry > NTRYPICK) return false;
 968
 969     // Initialize variables inside new try.
 970     nFilled = nKnown + 1;
 971     idProd.resize(nFilled);
 972     mProd.resize(nFilled);
 973     nTotal = nKnown;
 974     nNew = 0;
 975     nSpec = 0;
 976     nLeft = nPartons;
 977     mDiff = mProd[0];
 978     for (int i = 1; i < nFilled; ++i) mDiff -= mProd[i];
 979     diquarkClash = false;
 980     usedChannel = false;
 981
 982     // For weak decays collapse spectators to one particle.
 983     if ( (meMode == 22 || meMode == 23) && nLeft > 1) {
 984       FlavContainer flav1( idPartons[nPartons - 2] );
 985       FlavContainer flav2( idPartons[nPartons - 1] );
 986       int idHad;
 987       do idHad = flavSelPtr->combine( flav1, flav2);
 988       while (idHad == 0);
 989       double mHad = particleDataPtr->mass(idHad);
 990       mDiff -= mHad;
 991       idProd.push_back( idHad);
 992       mProd.push_back( mHad);
 993       ++nFilled;
 994       nSpec = 1;
 995       nLeft -= 2;
 996     }
 997
 998     // If there are partons left, then determine how many hadrons to pick.
 999     if (nLeft > 0) {
1000
1001       // For B -> gamma + X use geometrical distribution.
1002       if (meMode == 31) {
1003         double geom = rndmPtr->flat();
1004         nTotal = 1;
1005         do {
1006           ++nTotal;
1007           geom *= 2.;
1008         } while (geom < 1. && nTotal < 10);
1009
1010       // Calculate mass excess and from there average multiplicity.
1011       } else if (nFix == 0) {
1012         double mDiffPS = mDiff;
1013         for (int i = 0; i < nLeft; ++i)
1014           mDiffPS -= particleDataPtr->constituentMass( idPartons[i] );
1015         double average = 0.5 * (nKnown + nSpec) + 0.25 * nPartons
1016           + multIncrease * log( max( 1.1, mDiffPS / multRefMass ) );
1017         if (closedGLoop) average += multGoffset;
1018
1019         // Pick multiplicity according to Poissonian.
1020         double value = 1.;
1021         double sum = 1.;
1022         for (int nNow = nMin + 1; nNow <= 10; ++nNow) {
1023           value *= average / nNow;
1024           sum += value;
1025         }
1026         nTotal = nMin;
1027         value = 1.;
1028         sum *= rndmPtr->flat();
1029         sum -= value;
1030         if (sum > 0.) do {
1031           ++nTotal;
1032           value *= average / nTotal;
1033           sum -= value;
1034         } while (sum > 0. && nTotal < 10);
1035
1036       // Alternatively predetermined multiplicity.
1037       } else {
1038         nTotal = nFix;
1039       }
1040       nNew = nTotal - nKnown - nSpec;
1041
1042       // Set up ends of fragmented system, as copy of idPartons.
1043       flavEnds.resize(0);
1044       for (int i = 0; i < nLeft; ++i) {
1045         flavEnds.push_back( FlavContainer(idPartons[i]) );
1046         if (abs(idPartons[i]) > 100) flavSelPtr->assignPopQ( flavEnds[i] );
1047       }
1048
1049       // Fragment off at random, but save nLeft/2 for final recombination.
1050       if (nNew > nLeft/2) {
1051         FlavContainer flavNew;
1052         int idHad;
1053         for (int i = 0; i < nNew - nLeft/2; ++i) {
1054           // When four quarks consider last one to be spectator.
1055           int iEnd = int( (nLeft - 1.) * rndmPtr->flat() );
1056           // Pick new flavour and form a new hadron.
1057           do {
1058             flavNew = flavSelPtr->pick( flavEnds[iEnd] );
1059             idHad = flavSelPtr->combine( flavEnds[iEnd], flavNew);
1060           } while (idHad == 0);
1061           // Store new hadron and endpoint flavour.
1062           idProd.push_back( idHad);
1063           flavEnds[iEnd].anti(flavNew);
1064         }
1065       }
1066
1067       // When only two quarks left, combine to form final hadron.
1068       if (nLeft == 2) {
1069         int idHad;
1070         if ( abs(flavEnds[0].id) > 8 && abs(flavEnds[1].id) > 8)
1071           diquarkClash = true;
1072         else {
1073           do idHad = flavSelPtr->combine( flavEnds[0], flavEnds[1]);
1074           while (idHad == 0);
1075           idProd.push_back( idHad);
1076         }
1077
1078       // If four quarks, decide how to pair them up.
1079       } else {
1080         int iEnd1 = 0;
1081         int iEnd2 = 1;
1082         int iEnd3 = 2;
1083         int iEnd4 = 3;
1084         if ( rndmPtr->flat() < colRearrange) iEnd2 = 3;
1085         int relColSign =
1086           ( (flavEnds[iEnd1].id > 0 && flavEnds[iEnd1].id < 9)
1087           || flavEnds[iEnd1].id < -10 ) ? 1 : -1;
1088         if ( (flavEnds[iEnd2].id < 0 && flavEnds[iEnd2].id > -9)
1089           || flavEnds[iEnd2].id > 10 ) relColSign *= -1;
1090         if (relColSign == 1) iEnd2 = 2;
1091         if (iEnd2 == 2) iEnd3 = 1;
1092         if (iEnd2 == 3) iEnd4 = 1;
1093
1094         // Then combine to get final two hadrons.
1095         int idHad;
1096         if ( abs(flavEnds[iEnd1].id) > 8 && abs(flavEnds[iEnd2].id) > 8)
1097           diquarkClash = true;
1098         else {
1099           do idHad = flavSelPtr->combine( flavEnds[iEnd1], flavEnds[iEnd2]);
1100           while (idHad == 0);
1101           idProd.push_back( idHad);
1102         }
1103         if ( abs(flavEnds[iEnd3].id) > 8 && abs(flavEnds[iEnd4].id) > 8)
1104           diquarkClash = true;
1105         else {
1106           do idHad = flavSelPtr->combine( flavEnds[iEnd3], flavEnds[iEnd4]);
1107           while (idHad == 0);
1108           idProd.push_back( idHad);
1109         }
1110       }
1111
1112       // Find masses of the new hadrons.
1113       for (int i = nFilled; i < int(idProd.size()) ; ++i) {
1114         double mHad = particleDataPtr->mass(idProd[i]);
1115         mProd.push_back( mHad);
1116         mDiff -= mHad;
1117       }
1118     }
1119
1120     // Optional: check that this decay mode is not explicitly defined.
1121     if ( (meMode > 61 && meMode <= 80) && mDiff > mSafety && !diquarkClash ) {
1122       int idMatch[10], idPNow;
1123       usedChannel = false;
1124       bool matched = false;
1125       // Loop through all channels. Done if not same multiplicity.
1126       for (int i = 0; i < decDataPtr->sizeChannels(); ++i) {
1127         DecayChannel& channel = decDataPtr->channel(i);
1128         if (channel.multiplicity() != nTotal) continue;
1129         for (int k = 0; k < nTotal; ++k) idMatch[k] = channel.product(k);
1130         // Match particles one by one until fail.
1131         // Do not distinguish K0/K0bar/K0short/K0long at this stage.
1132         for (int j = 0; j < nTotal; ++j) {
1133           matched = false;
1134           idPNow = idProd[j + 1];
1135           if (idPNow == -311 || idPNow == 130 || idPNow == 310) idPNow = 311;
1136           for (int k = 0; k < nTotal - j; ++k)
1137           if (idMatch[k] == idPNow || (idMatch[k] == -311 && idPNow == 311)) {
1138             // Compress list of unmatched when matching worked.
1139             idMatch[k] = idMatch[nTotal - 1 - j];
1140             matched = true;
1141             break;
1142           }
1143           if (!matched) break;
1144         }
1145         // If matching worked, then chosen channel to be rejected.
1146         if (matched) {usedChannel = true; break;}
1147       }
1148     }
1149
1150   // Keep on trying until enough phase space and no clash.
1151   } while (mDiff < mSafety || diquarkClash || usedChannel);
1152
1153   // Update particle multiplicity.
1154   mult = idProd.size() - 1;
1155
1156   // For Dalitz decays shuffle Dalitz pair to the end of the list.
1157   if (meMode == 11 || meMode == 12) {
1158     int iL1 = 0;
1159     int iL2 = 0;
1160     for (int i = 1; i <= mult; ++i) {
1161       if (idProd[i] ==  11 || idProd[i] ==  13 || idProd[i] ==  15) iL1 = i;
1162       if (idProd[i] == -11 || idProd[i] == -13 || idProd[i] == -15) iL2 = i;
1163     }
1164     if (iL1 > 0 && iL2 > 0) {
1165       int idL1 = idProd[iL1];
1166       int idL2 = idProd[iL2];
1167       double mL1 = mProd[iL1];
1168       double mL2 = mProd[iL2];
1169       int iMove = 0;
1170       for (int i = 1; i <= mult; ++i) if (i != iL1 && i != iL2) {
1171         ++iMove;
1172         idProd[iMove] = idProd[i];
1173         mProd[iMove] = mProd[i];
1174       }
1175       idProd[mult - 1] = idL1;
1176       idProd[mult] = idL2;
1177       mProd[mult - 1] = mL1;
1178       mProd[mult] = mL2;
1179     }
1180   }
1181
1182   // Done.
1183   return true;
1184
1185 }
1186
1187 //--------------------------------------------------------------------------
1188
1189 // Set colour flow and scale in a decay explicitly to partons.
1190
1191 bool ParticleDecays::setColours(Event& event) {
1192
1193   // Decay to q qbar (or qbar q).
1194   if (meMode == 91 && idProd[1] > 0 && idProd[1] < 9) {
1195     int newCol = event.nextColTag();
1196     cols[1] = newCol;
1197     acols[2] = newCol;
1198   } else if (meMode == 91 && idProd[1] < 0 && idProd[1] > -9) {
1199     int newCol = event.nextColTag();
1200     cols[2] = newCol;
1201     acols[1] = newCol;
1202
1203   // Decay to g g.
1204   } else if (meMode == 91 && idProd[1] == 21) {
1205     int newCol1 = event.nextColTag();
1206     int newCol2 = event.nextColTag();
1207     cols[1] = newCol1;
1208     acols[1] = newCol2;
1209     cols[2] = newCol2;
1210     acols[2] = newCol1;
1211
1212   // Decay to g g g.
1213   } else if (meMode == 92 && idProd[1] == 21 && idProd[2] == 21
1214     &&  idProd[3] == 21) {
1215     int newCol1 = event.nextColTag();
1216     int newCol2 = event.nextColTag();
1217     int newCol3 = event.nextColTag();
1218     cols[1] = newCol1;
1219     acols[1] = newCol2;
1220     cols[2] = newCol2;
1221     acols[2] = newCol3;
1222     cols[3] = newCol3;
1223     acols[3] = newCol1;
1224
1225   // Decay to g g gamma: locate which is gamma.
1226   } else if (meMode == 92) {
1227     int iGlu1 = (idProd[1] == 21) ? 1 : 3;
1228     int iGlu2 = (idProd[2] == 21) ? 2 : 3;
1229     int newCol1 = event.nextColTag();
1230     int newCol2 = event.nextColTag();
1231     cols[iGlu1] = newCol1;
1232     acols[iGlu1] = newCol2;
1233     cols[iGlu2] = newCol2;
1234     acols[iGlu2] = newCol1;
1235
1236   // Unknown decay mode means failure.
1237   } else return false;
1238
1239   // Set maximum scale to be mass of decaying particle.
1240   scale = mProd[0];
1241
1242   // Done.
1243   return true;
1244
1245 }
1246
1247 //==========================================================================
1248
1249 } // end namespace Pythia8
1250