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