source: HiSusy/trunk/Pythia8/pythia8170/src/ParticleDecays.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES

File size: 43.6 KB
Line 
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
11namespace 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).
23const int    ParticleDecays::NTRYDECAY   = 10;
24
25// Number of times one tries to pick valid hadronic content in decay.
26const int    ParticleDecays::NTRYPICK    = 100;
27
28// Number of times one tries to pick decay topology.
29const int    ParticleDecays::NTRYMEWT    = 1000;
30
31// Maximal loop count in Dalitz decay treatment.
32const int    ParticleDecays::NTRYDALITZ  = 1000;
33
34// Minimal Dalitz pair mass this factor above threshold.
35const double ParticleDecays::MSAFEDALITZ = 1.000001;
36
37// These numbers are hardwired empirical parameters,
38// intended to speed up the M-generator.
39const 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
46void 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
123bool 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
328bool 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
347bool 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
363bool 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
383bool 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
484bool 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
625bool 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
760bool 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
850bool 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
919bool 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
1221bool 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
Note: See TracBrowser for help on using the repository browser.