PYTHIA8/pythia8170/src/ParticleDecays.cxx

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