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.
6
// Function definitions (not found in the header) for the
7
// ParticleDecays class.
9
#include "ParticleDecays.h"
13
//==========================================================================
15
// The ParticleDecays class.
17
//--------------------------------------------------------------------------
19
// Constants: could be changed here if desired, but normally should not.
20
// These are of technical nature, as described for each.
22
// Number of times one tries to let decay happen (for 2 nested loops).
23
const int ParticleDecays::NTRYDECAY = 10;
25
// Number of times one tries to pick valid hadronic content in decay.
26
const int ParticleDecays::NTRYPICK = 100;
28
// Number of times one tries to pick decay topology.
29
const int ParticleDecays::NTRYMEWT = 1000;
31
// Maximal loop count in Dalitz decay treatment.
32
const int ParticleDecays::NTRYDALITZ = 1000;
34
// Minimal Dalitz pair mass this factor above threshold.
35
const double ParticleDecays::MSAFEDALITZ = 1.000001;
37
// These numbers are hardwired empirical parameters,
38
// intended to speed up the M-generator.
39
const double ParticleDecays::WTCORRECTION[11] = { 1., 1., 1.,
40
2., 5., 15., 60., 250., 1250., 7000., 50000. };
42
//--------------------------------------------------------------------------
44
// Initialize and save pointers.
46
void ParticleDecays::init(Info* infoPtrIn, Settings& settings,
47
ParticleData* particleDataPtrIn, Rndm* rndmPtrIn,
48
Couplings* couplingsPtrIn, TimeShower* timesDecPtrIn,
49
StringFlav* flavSelPtrIn, DecayHandler* decayHandlePtrIn,
50
vector<int> handledParticles) {
52
// Save pointers to error messages handling and flavour generation.
54
particleDataPtr = particleDataPtrIn;
56
couplingsPtr = couplingsPtrIn;
57
flavSelPtr = flavSelPtrIn;
59
// Save pointer to timelike shower, as needed in some few decays.
60
timesDecPtr = timesDecPtrIn;
62
// Save pointer for external handling of some decays.
63
decayHandlePtr = decayHandlePtrIn;
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);
70
// Safety margin in mass to avoid troubles.
71
mSafety = settings.parm("ParticleDecays:mSafety");
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;
85
// B-Bbar mixing parameters.
86
mixB = settings.flag("ParticleDecays:mixB");
87
xBdMix = settings.parm("ParticleDecays:xBdMix");
88
xBsMix = settings.parm("ParticleDecays:xBsMix");
90
// Suppression of extra-hadron momenta in semileptonic decays.
91
sigmaSoft = settings.parm("ParticleDecays:sigmaSoft");
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");
100
// Minimum energy in system (+ m_q) from StringFragmentation.
101
stopMass = settings.parm("StringFragmentation:stopMass");
103
// Parameters for Dalitz decay virtual gamma mass spectrum.
104
sRhoDal = pow2(particleDataPtr->m0(113));
105
wRhoDal = pow2(particleDataPtr->mWidth(113));
107
// Allow showers in decays to qqbar/gg/ggg/gammagg.
108
doFSRinDecays = settings.flag("ParticleDecays:FSRinDecays");
110
// Use standard decays or dedicated tau decay package
111
sophisticatedTau = settings.mode("ParticleDecays:sophisticatedTau");
113
// Initialize the dedicated tau decay handler.
114
if (sophisticatedTau) tauDecayer.init(infoPtr, &settings,
115
particleDataPtr, rndmPtr, couplingsPtr);
119
//--------------------------------------------------------------------------
121
// Decay a particle; main method.
123
bool ParticleDecays::decay( int iDec, Event& event) {
125
// Check whether a decay is allowed, given the upcoming decay vertex.
126
Particle& decayer = event[iDec];
129
if (limitDecay && !checkVertex(decayer)) return true;
131
// Fill the decaying particle in slot 0 of arrays.
132
idDec = decayer.id();
136
iProd.push_back( iDec );
137
idProd.push_back( idDec );
138
mProd.push_back( decayer.m() );
140
// Check for oscillations B0 <-> B0bar or B_s0 <-> B_s0bar.
141
bool hasOscillated = (abs(idDec) == 511 || abs(idDec) == 531)
142
? oscillateB(decayer) : false;
143
if (hasOscillated) {idDec = - idDec; idProd[0] = idDec;}
145
// Particle data for decaying particle.
146
decDataPtr = &decayer.particleDataEntry();
148
// Optionally send on to external decay program.
149
bool doneExternally = false;
150
if (decDataPtr->doExternalDecay()) {
152
pProd.push_back(decayer.p());
153
doneExternally = decayHandlePtr->decay(idProd, mProd, pProd,
156
// If it worked, then store the decay products in the event record.
157
if (doneExternally) {
158
mult = idProd.size() - 1;
159
int status = (hasOscillated) ? 94 : 93;
160
for (int i = 1; i <= mult; ++i) {
161
int iPos = event.append( idProd[i], status, iDec, 0, 0, 0,
162
0, 0, pProd[i], mProd[i]);
163
iProd.push_back( iPos);
166
// Also mark mother decayed and store daughters.
167
event[iDec].statusNeg();
168
event[iDec].daughters( iProd[1], iProd[mult]);
172
// Check if the particle is tau and let the special tau decayer handle it.
173
if (decayer.idAbs() == 15 && !doneExternally && sophisticatedTau) {
174
doneExternally = tauDecayer.decay(iDec, event);
175
if (doneExternally) return true;
178
// Now begin normal internal decay treatment.
179
if (!doneExternally) {
181
// Allow up to ten tries to pick a channel.
182
if (!decDataPtr->preparePick(idDec)) return false;
183
bool foundChannel = false;
184
bool hasStored = false;
185
for (int iTryChannel = 0; iTryChannel < NTRYDECAY; ++iTryChannel) {
187
// Remove previous failed channel.
188
if (hasStored) event.popBack(mult);
191
// Pick new channel. Read out basics.
192
DecayChannel& channel = decDataPtr->pickChannel();
193
meMode = channel.meMode();
194
keepPartons = (meMode > 90 && meMode <= 100);
195
mult = channel.multiplicity();
197
// Allow up to ten tries for each channel (e.g with different masses).
198
bool foundMode = false;
199
for (int iTryMode = 0; iTryMode < NTRYDECAY; ++iTryMode) {
204
// Extract and store the decay products in local arrays.
206
for (int i = 0; i < mult; ++i) {
207
int idNow = channel.product(i);
208
int idAbs = abs(idNow);
209
if ( idAbs < 10 || idAbs == 21 || idAbs == 81 || idAbs == 82
210
|| idAbs == 83 || (idAbs > 1000 && idAbs < 10000
211
&& (idAbs/10)%10 == 0) ) hasPartons = true;
212
if (idDec < 0 && particleDataPtr->hasAnti(idNow)) idNow = -idNow;
213
double mNow = particleDataPtr->mass(idNow);
214
idProd.push_back( idNow);
215
mProd.push_back( mNow);
218
// Decays into partons usually translate into hadrons.
219
if (hasPartons && !keepPartons && !pickHadrons()) continue;
221
// Need to set colour flow if explicit decay to partons.
224
for (int i = 0; i <= mult; ++i) {
228
if (hasPartons && keepPartons && !setColours(event)) continue;
230
// Check that enough phase space for decay.
232
double mDiff = mProd[0];
233
for (int i = 1; i <= mult; ++i) mDiff -= mProd[i];
234
if (mDiff < mSafety) continue;
237
// End of inner trial loops. Check if succeeded or not.
241
if (!foundMode) continue;
243
// Store decay products in the event record.
244
int status = (hasOscillated) ? 92 : 91;
245
for (int i = 1; i <= mult; ++i) {
246
int iPos = event.append( idProd[i], status, iDec, 0, 0, 0,
247
cols[i], acols[i], Vec4(0., 0., 0., 0.), mProd[i], scale);
248
iProd.push_back( iPos);
252
// Pick mass of Dalitz decay. Temporarily change multiplicity.
253
if ( (meMode == 11 || meMode == 12 || meMode == 13)
254
&& !dalitzMass() ) continue;
256
// Do a decay, split by multiplicity.
257
bool decayed = false;
258
if (mult == 1) decayed = oneBody(event);
259
else if (mult == 2) decayed = twoBody(event);
260
else if (mult == 3) decayed = threeBody(event);
261
else decayed = mGenerator(event);
262
if (!decayed) continue;
264
// Kinematics of gamma* -> l- l+ in Dalitz decay. Restore multiplicity.
265
if (meMode == 11 || meMode == 12 || meMode == 13)
266
dalitzKinematics(event);
268
// End of outer trial loops.
273
// If the decay worked, then mark mother decayed and store daughters.
275
event[iDec].statusNeg();
276
event[iDec].daughters( iProd[1], iProd[mult]);
278
// Else remove unused daughters and return failure.
280
if (hasStored) event.popBack(mult);
281
infoPtr->errorMsg("Error in ParticleDecays::decay: "
282
"failed to find workable decay channel");
286
// Now finished normal internal decay treatment.
289
// Set decay vertex when this is displaced.
290
if (event[iDec].hasVertex() || event[iDec].tau() > 0.) {
291
Vec4 vDec = event[iDec].vDec();
292
for (int i = 1; i <= mult; ++i) event[iProd[i]].vProd( vDec );
295
// Set lifetime of daughters.
296
for (int i = 1; i <= mult; ++i)
297
event[iProd[i]].tau( event[iProd[i]].tau0() * rndmPtr->exp() );
299
// In a decay explicitly to partons then optionally do a shower,
300
// and always flag that partonic system should be fragmented.
301
if (hasPartons && keepPartons && doFSRinDecays)
302
timesDecPtr->shower( iProd[1], iProd.back(), event, mProd[0]);
304
// For Hidden Valley particles also allow leptons to shower.
305
else if (event[iDec].idAbs() > 4900000 && event[iDec].idAbs() < 5000000
306
&& doFSRinDecays && mult == 2 && event[iProd[1]].isLepton()) {
307
event[iProd[1]].scale(mProd[0]);
308
event[iProd[2]].scale(mProd[0]);
309
timesDecPtr->shower( iProd[1], iProd.back(), event, mProd[0]);
317
//--------------------------------------------------------------------------
319
// Check whether a decay is allowed, given the upcoming decay vertex.
321
bool ParticleDecays::checkVertex(Particle& decayer) {
323
// Check whether any of the conditions are not fulfilled.
324
if (limitTau0 && decayer.tau0() > tau0Max) return false;
325
if (limitTau && decayer.tau() > tauMax) return false;
326
if (limitRadius && pow2(decayer.xDec()) + pow2(decayer.yDec())
327
+ pow2(decayer.zDec()) > pow2(rMax)) return false;
328
if (limitCylinder && (pow2(decayer.xDec()) + pow2(decayer.yDec())
329
> pow2(xyMax) || abs(decayer.zDec()) > zMax) ) return false;
336
//--------------------------------------------------------------------------
338
// Check for oscillations B0 <-> B0bar or B_s0 <-> B_s0bar.
340
bool ParticleDecays::oscillateB(Particle& decayer) {
342
// Extract relevant information and decide.
343
if (!mixB) return false;
344
double xBmix = (abs(decayer.id()) == 511) ? xBdMix : xBsMix;
345
double tau = decayer.tau();
346
double tau0 = decayer.tau0();
347
double probosc = pow2(sin(0.5 * xBmix * tau / tau0));
348
return (probosc > rndmPtr->flat());
352
//--------------------------------------------------------------------------
354
// Do a one-body decay. (Rare; e.g. for K0 -> K0_short.)
356
bool ParticleDecays::oneBody(Event& event) {
358
// References to the particles involved.
359
Particle& decayer = event[iProd[0]];
360
Particle& prod = event[iProd[1]];
362
// Set momentum and expand mother information.
363
prod.p( decayer.p() );
364
prod.m( decayer.m() );
365
prod.mother2( iProd[0] );
372
//--------------------------------------------------------------------------
374
// Do a two-body decay.
376
bool ParticleDecays::twoBody(Event& event) {
378
// References to the particles involved.
379
Particle& decayer = event[iProd[0]];
380
Particle& prod1 = event[iProd[1]];
381
Particle& prod2 = event[iProd[2]];
384
double m0 = mProd[0];
385
double m1 = mProd[1];
386
double m2 = mProd[2];
388
// Energies and absolute momentum in the rest frame.
389
if (m1 + m2 + mSafety > m0) return false;
390
double e1 = 0.5 * (m0*m0 + m1*m1 - m2*m2) / m0;
391
double e2 = 0.5 * (m0*m0 + m2*m2 - m1*m1) / m0;
392
double pAbs = 0.5 * sqrtpos( (m0 - m1 - m2) * (m0 + m1 + m2)
393
* (m0 + m1 - m2) * (m0 - m1 + m2) ) / m0;
395
// When meMode = 2, for V -> PS2 + PS3 (V = vector, pseudoscalar),
396
// need to check if production is PS0 -> PS1/gamma + V.
397
int iMother = event[iProd[0]].mother1();
400
if (iMother <= 0 || iMother >= iProd[0]) meMode = 0;
402
int iDaughter1 = event[iMother].daughter1();
403
int iDaughter2 = event[iMother].daughter2();
404
if (iDaughter2 != iDaughter1 + 1) meMode = 0;
406
int idMother = abs( event[iMother].id() );
407
if (idMother <= 100 || idMother%10 !=1
408
|| (idMother/1000)%10 != 0) meMode = 0;
410
int iSister = (iProd[0] == iDaughter1) ? iDaughter2 : iDaughter1;
411
idSister = abs( event[iSister].id() );
412
if ( (idSister <= 100 || idSister%10 !=1
413
|| (idSister/1000)%10 != 0) && idSister != 22) meMode = 0;
419
// Begin loop over matrix-element corrections.
420
double wtME, wtMEmax;
427
// Isotropic angles give three-momentum.
428
double cosTheta = 2. * rndmPtr->flat() - 1.;
429
double sinTheta = sqrt(1. - cosTheta*cosTheta);
430
double phi = 2. * M_PI * rndmPtr->flat();
431
double pX = pAbs * sinTheta * cos(phi);
432
double pY = pAbs * sinTheta * sin(phi);
433
double pZ = pAbs * cosTheta;
435
// Fill four-momenta and boost them away from mother rest frame.
436
prod1.p( pX, pY, pZ, e1);
437
prod2.p( -pX, -pY, -pZ, e2);
438
prod1.bst( decayer.p(), decayer.m() );
439
prod2.bst( decayer.p(), decayer.m() );
441
// Matrix element for PS0 -> PS1 + V1 -> PS1 + PS2 + PS3 of form
442
// cos**2(theta02) in V1 rest frame, and for PS0 -> gamma + V1
443
// -> gamma + PS2 + PS3 of form sin**2(theta02).
445
double p10 = decayer.p() * event[iMother].p();
446
double p12 = decayer.p() * prod1.p();
447
double p02 = event[iMother].p() * prod1.p();
448
double s0 = pow2(event[iMother].m());
449
double s1 = pow2(decayer.m());
450
double s2 = pow2(prod1.m());
451
if (idSister != 22) wtME = pow2(p10 * p12 - s1 * p02);
452
else wtME = s1 * (2. * p10 * p12 * p02 - s1 * p02*p02
453
- s0 * p12*p12 - s2 * p10*p10 + s1 * s0 * s2);
454
wtME = max( wtME, 1e-6 * s1*s1 * s0 * s2);
455
wtMEmax = (p10*p10 - s1 * s0) * (p12*p12 - s1 * s2);
458
// Break out of loop if no sensible ME weight.
459
if(loop > NTRYMEWT) {
460
infoPtr->errorMsg("ParticleDecays::twoBody: "
461
"caught in infinite ME weight loop");
465
// If rejected, try again with new invariant masses.
466
} while ( wtME < rndmPtr->flat() * wtMEmax );
473
//--------------------------------------------------------------------------
475
// Do a three-body decay (except Dalitz decays).
477
bool ParticleDecays::threeBody(Event& event) {
479
// References to the particles involved.
480
Particle& decayer = event[iProd[0]];
481
Particle& prod1 = event[iProd[1]];
482
Particle& prod2 = event[iProd[2]];
483
Particle& prod3 = event[iProd[3]];
485
// Mother and sum daughter masses. Fail if too close.
486
double m0 = mProd[0];
487
double m1 = mProd[1];
488
double m2 = mProd[2];
489
double m3 = mProd[3];
490
double mSum = m1 + m2 + m3;
491
double mDiff = m0 - mSum;
492
if (mDiff < mSafety) return false;
494
// Kinematical limits for 2+3 mass. Maximum phase-space weight.
495
double m23Min = m2 + m3;
496
double m23Max = m0 - m1;
497
double p1Max = 0.5 * sqrtpos( (m0 - m1 - m23Min) * (m0 + m1 + m23Min)
498
* (m0 + m1 - m23Min) * (m0 - m1 + m23Min) ) / m0;
499
double p23Max = 0.5 * sqrtpos( (m23Max - m2 - m3) * (m23Max + m2 + m3)
500
* (m23Max + m2 - m3) * (m23Max - m2 + m3) ) / m23Max;
501
double wtPSmax = 0.5 * p1Max * p23Max;
503
// Begin loop over matrix-element corrections.
504
double wtME, wtMEmax, wtPS, m23, p1Abs, p23Abs;
509
// Pick an intermediate mass m23 flat in the allowed range.
511
m23 = m23Min + rndmPtr->flat() * mDiff;
513
// Translate into relative momenta and find phase-space weight.
514
p1Abs = 0.5 * sqrtpos( (m0 - m1 - m23) * (m0 + m1 + m23)
515
* (m0 + m1 - m23) * (m0 - m1 + m23) ) / m0;
516
p23Abs = 0.5 * sqrtpos( (m23 - m2 - m3) * (m23 + m2 + m3)
517
* (m23 + m2 - m3) * (m23 - m2 + m3) ) / m23;
518
wtPS = p1Abs * p23Abs;
520
// If rejected, try again with new invariant masses.
521
} while ( wtPS < rndmPtr->flat() * wtPSmax );
523
// Set up m23 -> m2 + m3 isotropic in its rest frame.
524
double cosTheta = 2. * rndmPtr->flat() - 1.;
525
double sinTheta = sqrt(1. - cosTheta*cosTheta);
526
double phi = 2. * M_PI * rndmPtr->flat();
527
double pX = p23Abs * sinTheta * cos(phi);
528
double pY = p23Abs * sinTheta * sin(phi);
529
double pZ = p23Abs * cosTheta;
530
double e2 = sqrt( m2*m2 + p23Abs*p23Abs);
531
double e3 = sqrt( m3*m3 + p23Abs*p23Abs);
532
prod2.p( pX, pY, pZ, e2);
533
prod3.p( -pX, -pY, -pZ, e3);
535
// Set up m0 -> m1 + m23 isotropic in its rest frame.
536
cosTheta = 2. * rndmPtr->flat() - 1.;
537
sinTheta = sqrt(1. - cosTheta*cosTheta);
538
phi = 2. * M_PI * rndmPtr->flat();
539
pX = p1Abs * sinTheta * cos(phi);
540
pY = p1Abs * sinTheta * sin(phi);
541
pZ = p1Abs * cosTheta;
542
double e1 = sqrt( m1*m1 + p1Abs*p1Abs);
543
double e23 = sqrt( m23*m23 + p1Abs*p1Abs);
544
prod1.p( pX, pY, pZ, e1);
546
// Boost 2 + 3 to the 0 rest frame.
547
Vec4 p23( -pX, -pY, -pZ, e23);
548
prod2.bst( p23, m23 );
549
prod3.bst( p23, m23 );
551
// Matrix-element weight for omega/phi -> pi+ pi- pi0.
553
double p1p2 = prod1.p() * prod2.p();
554
double p1p3 = prod1.p() * prod3.p();
555
double p2p3 = prod2.p() * prod3.p();
556
wtME = pow2(m1 * m2 * m3) - pow2(m1 * p2p3) - pow2(m2 * p1p3)
557
- pow2(m3 * p1p2) + 2. * p1p2 * p1p3 * p2p3;
558
wtMEmax = pow3(m0 * m0) / 150.;
560
// Effective matrix element for nu spectrum in tau -> nu + hadrons.
561
} else if (meMode == 21) {
562
double x1 = 2. * prod1.e() / m0;
563
wtME = x1 * (3. - 2. * x1);
564
double xMax = min( 0.75, 2. * (1. - mSum / m0) );
565
wtMEmax = xMax * (3. - 2. * xMax);
567
// Matrix element for weak decay (only semileptonic for c and b).
568
} else if ((meMode == 22 || meMode == 23) && prod1.isLepton()) {
569
wtME = m0 * prod1.e() * (prod2.p() * prod3.p());
570
wtMEmax = min( pow4(m0) / 16., m0 * (m0 - m1 - m2) * (m0 - m1 - m3)
573
// Effective matrix element for weak decay to hadrons (B -> D, D -> K).
574
} else if (meMode == 22 || meMode == 23) {
575
double x1 = 2. * prod1.pAbs() / m0;
576
wtME = x1 * (3. - 2. * x1);
577
double xMax = min( 0.75, 2. * (1. - mSum / m0) );
578
wtMEmax = xMax * (3. - 2. * xMax);
580
// Effective matrix element for gamma spectrum in B -> gamma + hadrons.
581
} else if (meMode == 31) {
582
double x1 = 2. * prod1.e() / m0;
584
double x1Max = 1. - pow2(mSum / m0);
585
wtMEmax = pow3(x1Max);
587
// Matrix-element weight for "onium" -> g + g + g or gamma + g + g.
588
} else if (meMode == 92) {
589
double x1 = 2. * prod1.e() / m0;
590
double x2 = 2. * prod2.e() / m0;
591
double x3 = 2. * prod3.e() / m0;
592
wtME = pow2( (1. - x1) / (x2 * x3) ) + pow2( (1. - x2) / (x1 * x3) )
593
+ pow2( (1. - x3) / (x1 * x2) );
595
// For gamma + g + g require minimum mass for g + g system.
596
if (prod1.id() == 22 && sqrt(1. - x1) * m0 < 2. * stopMass) wtME = 0.;
597
if (prod2.id() == 22 && sqrt(1. - x2) * m0 < 2. * stopMass) wtME = 0.;
598
if (prod3.id() == 22 && sqrt(1. - x3) * m0 < 2. * stopMass) wtME = 0.;
601
// If rejected, try again with new invariant masses.
602
} while ( wtME < rndmPtr->flat() * wtMEmax );
604
// Boost 1 + 2 + 3 to the current frame.
605
prod1.bst( decayer.p(), decayer.m() );
606
prod2.bst( decayer.p(), decayer.m() );
607
prod3.bst( decayer.p(), decayer.m() );
614
//--------------------------------------------------------------------------
616
// Do a multibody decay using the M-generator algorithm.
618
bool ParticleDecays::mGenerator(Event& event) {
620
// Mother and sum daughter masses. Fail if too close or inconsistent.
621
double m0 = mProd[0];
622
double mSum = mProd[1];
623
for (int i = 2; i <= mult; ++i) mSum += mProd[i];
624
double mDiff = m0 - mSum;
625
if (mDiff < mSafety) return false;
627
// Begin setup of intermediate invariant masses.
629
for (int i = 0; i <= mult; ++i) mInv.push_back( mProd[i]);
631
// Calculate the maximum weight in the decay.
632
double wtPS, wtME, wtMEmax;
633
double wtPSmax = 1. / WTCORRECTION[mult];
634
double mMax = mDiff + mProd[mult];
636
for (int i = mult - 1; i > 0; --i) {
639
double mNow = mProd[i];
640
wtPSmax *= 0.5 * sqrtpos( (mMax - mMin - mNow) * (mMax + mMin + mNow)
641
* (mMax + mMin - mNow) * (mMax - mMin + mNow) ) / mMax;
644
// Begin loop over matrix-element corrections.
649
// Begin loop to find the set of intermediate invariant masses.
653
// Find and order random numbers in descending order.
655
rndmOrd.push_back(1.);
656
for (int i = 1; i < mult - 1; ++i) {
657
double rndm = rndmPtr->flat();
658
rndmOrd.push_back(rndm);
659
for (int j = i - 1; j > 0; --j) {
660
if (rndm > rndmOrd[j]) swap( rndmOrd[j], rndmOrd[j+1] );
664
rndmOrd.push_back(0.);
666
// Translate into intermediate masses and find weight.
667
for (int i = mult - 1; i > 0; --i) {
668
mInv[i] = mInv[i+1] + mProd[i] + (rndmOrd[i-1] - rndmOrd[i]) * mDiff;
669
wtPS *= 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
670
* (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
671
* (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
674
// If rejected, try again with new invariant masses.
675
} while ( wtPS < rndmPtr->flat() * wtPSmax );
677
// Perform two-particle decays in the respective rest frame.
678
pInv.resize(mult + 1);
679
for (int i = 1; i < mult; ++i) {
680
double pAbs = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
681
* (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
682
* (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
684
// Isotropic angles give three-momentum.
685
double cosTheta = 2. * rndmPtr->flat() - 1.;
686
double sinTheta = sqrt(1. - cosTheta*cosTheta);
687
double phi = 2. * M_PI * rndmPtr->flat();
688
double pX = pAbs * sinTheta * cos(phi);
689
double pY = pAbs * sinTheta * sin(phi);
690
double pZ = pAbs * cosTheta;
692
// Calculate energies, fill four-momenta.
693
double eHad = sqrt( mProd[i]*mProd[i] + pAbs*pAbs);
694
double eInv = sqrt( mInv[i+1]*mInv[i+1] + pAbs*pAbs);
695
event[iProd[i]].p( pX, pY, pZ, eHad);
696
pInv[i+1].p( -pX, -pY, -pZ, eInv);
699
// Boost decay products to the mother rest frame.
700
event[iProd[mult]].p( pInv[mult] );
701
for (int iFrame = mult - 1; iFrame > 1; --iFrame)
702
for (int i = iFrame; i <= mult; ++i)
703
event[iProd[i]].bst( pInv[iFrame], mInv[iFrame]);
705
// Effective matrix element for nu spectrum in tau -> nu + hadrons.
706
if (meMode == 21 && event[iProd[1]].isLepton()) {
707
double x1 = 2. * event[iProd[1]].e() / m0;
708
wtME = x1 * (3. - 2. * x1);
709
double xMax = min( 0.75, 2. * (1. - mSum / m0) );
710
wtMEmax = xMax * (3. - 2. * xMax);
712
// Effective matrix element for weak decay (only semileptonic for c and b).
713
// Particles 4 onwards should be made softer explicitly?
714
} else if ((meMode == 22 || meMode == 23) && event[iProd[1]].isLepton()) {
715
Vec4 pRest = event[iProd[3]].p();
716
for (int i = 4; i <= mult; ++i) pRest += event[iProd[i]].p();
717
wtME = m0 * event[iProd[1]].e() * (event[iProd[2]].p() * pRest);
718
for (int i = 4; i <= mult; ++i) wtME
719
*= exp(- event[iProd[i]].pAbs2() / pow2(sigmaSoft) );
720
wtMEmax = pow4(m0) / 16.;
722
// Effective matrix element for weak decay to hadrons (B -> D, D -> K).
723
} else if (meMode == 22 || meMode == 23) {
724
double x1 = 2. * event[iProd[1]].pAbs() / m0;
725
wtME = x1 * (3. - 2. * x1);
726
double xMax = min( 0.75, 2. * (1. - mSum / m0) );
727
wtMEmax = xMax * (3. - 2. * xMax);
729
// Effective matrix element for gamma spectrum in B -> gamma + hadrons.
730
} else if (meMode == 31) {
731
double x1 = 2. * event[iProd[1]].e() / m0;
733
double x1Max = 1. - pow2(mSum / m0);
734
wtMEmax = pow3(x1Max);
737
// If rejected, try again with new invariant masses.
738
} while ( wtME < rndmPtr->flat() * wtMEmax );
740
// Boost decay products to the current frame.
741
pInv[1].p( event[iProd[0]].p() );
742
for (int i = 1; i <= mult; ++i) event[iProd[i]].bst( pInv[1], mInv[1] );
749
//--------------------------------------------------------------------------
751
// Select mass of lepton pair in a Dalitz decay.
753
bool ParticleDecays::dalitzMass() {
755
// Mother and sum daughter masses.
757
for (int i = 1; i <= mult - 2; ++i) mSum1 += mProd[i];
758
if (meMode == 13) mSum1 *= MSAFEDALITZ;
759
double mSum2 = MSAFEDALITZ * (mProd[mult -1] + mProd[mult]);
760
double mDiff = mProd[0] - mSum1 - mSum2;
762
// Fail if too close or inconsistent.
763
if (mDiff < mSafety) return false;
764
if (idProd[mult - 1] + idProd[mult] != 0
765
|| mProd[mult - 1] != mProd[mult]) {
766
infoPtr->errorMsg("Error in ParticleDecays::dalitzMass:"
767
" inconsistent flavour/mass assignments");
770
if ( meMode == 13 && (idProd[1] + idProd[2] != 0
771
|| mProd[1] != mProd[2]) ) {
772
infoPtr->errorMsg("Error in ParticleDecays::dalitzMass:"
773
" inconsistent flavour/mass assignments");
777
// Case 1: one Dalitz pair.
778
if (meMode == 11 || meMode == 12) {
780
// Kinematical limits for gamma* squared mass.
781
double sGamMin = pow2(mSum2);
782
double sGamMax = pow2(mProd[0] - mSum1);
783
// Select virtual gamma squared mass. Guessed form for meMode == 12.
787
if (++loop > NTRYDALITZ) return false;
788
sGam = sGamMin * pow( sGamMax / sGamMin, rndmPtr->flat() );
789
wtGam = (1. + 0.5 * sGamMin / sGam) * sqrt(1. - sGamMin / sGam)
790
* pow3(1. - sGam / sGamMax) * sRhoDal * (sRhoDal + wRhoDal)
791
/ ( pow2(sGam - sRhoDal) + sRhoDal * wRhoDal );
792
} while ( wtGam < rndmPtr->flat() );
794
// Store results in preparation for doing a one-less-body decay.
796
mProd[mult] = sqrt(sGam);
798
// Case 2: two Dalitz pairs.
801
// Kinematical limits for 1 + 2 and 3 + 4 gamma* masses.
802
double s0 = pow2(mProd[0]);
803
double s12Min = pow2(mSum1);
804
double s12Max = pow2(mProd[0] - mSum2);
805
double s34Min = pow2(mSum2);
806
double s34Max = pow2(mProd[0] - mSum1);
808
// Select virtual gamma squared masses. Guessed form for meMode == 13.
809
double s12, s34, wt12, wt34, wtPAbs, wtAll;
812
if (++loop > NTRYDALITZ) return false;
813
s12 = s12Min * pow( s12Max / s12Min, rndmPtr->flat() );
814
wt12 = (1. + 0.5 * s12Min / s12) * sqrt(1. - s12Min / s12)
815
* sRhoDal * (sRhoDal + wRhoDal)
816
/ ( pow2(s12 - sRhoDal) + sRhoDal * wRhoDal );
817
s34 = s34Min * pow( s34Max / s34Min, rndmPtr->flat() );
818
wt34 = (1. + 0.5 * s34Min / s34) * sqrt(1. - s34Min / s34)
819
* sRhoDal * (sRhoDal + wRhoDal)
820
/ ( pow2(s34 - sRhoDal) + sRhoDal * wRhoDal );
821
wtPAbs = sqrtpos( pow2(1. - (s12 + s34)/ s0)
822
- 4. * s12 * s34 / (s0 * s0) );
823
wtAll = wt12 * wt34 * pow3(wtPAbs);
824
if (wtAll > 1.) infoPtr->errorMsg(
825
"Error in ParticleDecays::dalitzMass: weight > 1");
826
} while (wtAll < rndmPtr->flat());
828
// Store results in preparation for doing a two-body decay.
830
mProd[1] = sqrt(s12);
831
mProd[2] = sqrt(s34);
839
//--------------------------------------------------------------------------
841
// Do kinematics of gamma* -> l- l+ in Dalitz decay.
843
bool ParticleDecays::dalitzKinematics(Event& event) {
845
// Restore multiplicity.
846
int nDal = (meMode < 13) ? 1 : 2;
849
// Loop over one or two lepton pairs.
850
for (int iDal = 0; iDal < nDal; ++iDal) {
852
// References to the particles involved.
853
Particle& decayer = event[iProd[0]];
854
Particle& prodA = (iDal == 0) ? event[iProd[mult - 1]]
856
Particle& prodB = (iDal == 0) ? event[iProd[mult]]
859
// Reconstruct required rotations and boosts backwards.
860
Vec4 pDec = decayer.p();
861
int iGam = (meMode < 13) ? mult - 1 : 2 - iDal;
862
Vec4 pGam = event[iProd[iGam]].p();
863
pGam.bstback( pDec, decayer.m() );
864
double phiGam = pGam.phi();
865
pGam.rot( 0., -phiGam);
866
double thetaGam = pGam.theta();
867
pGam.rot( -thetaGam, 0.);
869
// Masses and phase space in gamma* rest frame.
870
double mGam = (meMode < 13) ? mProd[mult - 1] : mProd[2 - iDal];
871
double mA = prodA.m();
872
double mB = prodB.m();
873
double mGamMin = MSAFEDALITZ * (mA + mB);
874
double mGamRat = pow2(mGamMin / mGam);
875
double pGamAbs = 0.5 * sqrtpos( (mGam - mA - mB) * (mGam + mA + mB) );
877
// Set up decay in gamma* rest frame, reference along +z axis.
878
double cosTheta, cos2Theta;
880
cosTheta = 2. * rndmPtr->flat() - 1.;
881
cos2Theta = cosTheta * cosTheta;
882
} while ( 1. + cos2Theta + mGamRat * (1. - cos2Theta)
883
< 2. * rndmPtr->flat() );
884
double sinTheta = sqrt(1. - cosTheta*cosTheta);
885
double phi = 2. * M_PI * rndmPtr->flat();
886
double pX = pGamAbs * sinTheta * cos(phi);
887
double pY = pGamAbs * sinTheta * sin(phi);
888
double pZ = pGamAbs * cosTheta;
889
double eA = sqrt( mA*mA + pGamAbs*pGamAbs);
890
double eB = sqrt( mB*mB + pGamAbs*pGamAbs);
891
prodA.p( pX, pY, pZ, eA);
892
prodB.p( -pX, -pY, -pZ, eB);
894
// Boost to lab frame.
895
prodA.bst( pGam, mGam);
896
prodB.bst( pGam, mGam);
897
prodA.rot( thetaGam, phiGam);
898
prodB.rot( thetaGam, phiGam);
899
prodA.bst( pDec, decayer.m() );
900
prodB.bst( pDec, decayer.m() );
908
//--------------------------------------------------------------------------
910
// Translate a partonic content into a set of actual hadrons.
912
bool ParticleDecays::pickHadrons() {
914
// Find partonic decay products. Rest are known id's, mainly hadrons,
915
// when necessary shuffled to beginning of idProd list.
919
bool closedGLoop = false;
920
for (int i = 1; i <= mult; ++i) {
921
int idAbs = abs(idProd[i]);
922
if ( idAbs < 9 || (idAbs > 1000 && idAbs < 10000 && (idAbs/10)%10 == 0)
923
|| idAbs == 81 || idAbs == 82 || idAbs == 83) {
925
idPartons.push_back(idProd[i]);
926
if (idAbs == 83) closedGLoop = true;
930
idProd[nKnown] = idProd[i];
931
mProd[nKnown] = mProd[i];
936
// Replace generic spectator flavour code by the actual one.
937
for (int i = 0; i < nPartons; ++i) {
938
int idPart = idPartons[i];
941
int idAbs = abs(idDec);
942
if ( (idAbs/1000)%10 == 0 ) {
943
idNew = -(idAbs/10)%10;
944
if ((idAbs/100)%2 == 1) idNew = -idNew;
945
} else if ( (idAbs/100)%10 >= (idAbs/10)%10 )
946
idNew = 100 * ((idAbs/10)%100) + 3;
947
else idNew = 1000 * ((idAbs/10)%10) + 100 * ((idAbs/100)%10) + 1;
948
if (idDec < 0) idNew = -idNew;
950
// Replace generic random flavour by a randomly selected one.
951
} else if (idPart == 82 || idPart == 83) {
954
int idDummy = -flavSelPtr->pickLightQ();
955
FlavContainer flavDummy(idDummy, idPart - 82);
956
do idNew = flavSelPtr->pick(flavDummy).id;
958
mFlav = particleDataPtr->constituentMass(idNew);
959
} while (2. * mFlav + stopMass > mProd[0]);
960
} else if (idPart == -82 || idPart == -83) {
961
idNew = -idPartons[i-1];
963
idPartons[i] = idNew;
966
// Determine whether fixed multiplicity or to be selected at random.
967
int nMin = max( 2, nKnown + nPartons / 2);
968
if (meMode == 23) nMin = 3;
969
if (meMode > 41 && meMode <= 50) nMin = meMode - 40;
970
if (meMode > 51 && meMode <= 60) nMin = meMode - 50;
972
if (meMode == 0) nFix = nMin;
973
if (meMode == 11) nFix = 3;
974
if (meMode == 12) nFix = 4;
975
if (meMode > 61 && meMode <= 70) nFix = meMode - 60;
976
if (meMode > 71 && meMode <= 80) nFix = meMode - 70;
977
if (nFix > 0 && nKnown + nPartons/2 > nFix) return false;
979
// Initial values for loop to set new hadronic content.
980
int nFilled, nTotal, nNew, nSpec, nLeft;
983
bool diquarkClash = false;
984
bool usedChannel = false;
986
// Begin loop; interrupt if multiple tries fail.
989
if (nTry > NTRYPICK) return false;
991
// Initialize variables inside new try.
992
nFilled = nKnown + 1;
993
idProd.resize(nFilled);
994
mProd.resize(nFilled);
999
for (int i = 1; i < nFilled; ++i) mDiff -= mProd[i];
1000
diquarkClash = false;
1001
usedChannel = false;
1003
// For weak decays collapse spectators to one particle.
1004
if ( (meMode == 22 || meMode == 23) && nLeft > 1) {
1005
FlavContainer flav1( idPartons[nPartons - 2] );
1006
FlavContainer flav2( idPartons[nPartons - 1] );
1008
do idHad = flavSelPtr->combine( flav1, flav2);
1010
double mHad = particleDataPtr->mass(idHad);
1012
idProd.push_back( idHad);
1013
mProd.push_back( mHad);
1019
// If there are partons left, then determine how many hadrons to pick.
1022
// For B -> gamma + X use geometrical distribution.
1024
double geom = rndmPtr->flat();
1029
} while (geom < 1. && nTotal < 10);
1031
// Calculate mass excess and from there average multiplicity.
1032
} else if (nFix == 0) {
1033
double multIncreaseNow = (meMode == 23)
1034
? multIncreaseWeak : multIncrease;
1035
double mDiffPS = mDiff;
1036
for (int i = 0; i < nLeft; ++i)
1037
mDiffPS -= particleDataPtr->constituentMass( idPartons[i] );
1038
double average = 0.5 * (nKnown + nSpec) + 0.25 * nPartons
1039
+ multIncreaseNow * log( max( 1.1, mDiffPS / multRefMass ) );
1040
if (closedGLoop) average += multGoffset;
1042
// Pick multiplicity according to Poissonian.
1045
for (int nNow = nMin + 1; nNow <= 10; ++nNow) {
1046
value *= average / nNow;
1051
sum *= rndmPtr->flat();
1055
value *= average / nTotal;
1057
} while (sum > 0. && nTotal < 10);
1059
// Alternatively predetermined multiplicity.
1063
nNew = nTotal - nKnown - nSpec;
1065
// Set up ends of fragmented system, as copy of idPartons.
1067
for (int i = 0; i < nLeft; ++i) {
1068
flavEnds.push_back( FlavContainer(idPartons[i]) );
1069
if (abs(idPartons[i]) > 100) flavSelPtr->assignPopQ( flavEnds[i] );
1072
// Fragment off at random, but save nLeft/2 for final recombination.
1073
if (nNew > nLeft/2) {
1074
FlavContainer flavNew;
1076
for (int i = 0; i < nNew - nLeft/2; ++i) {
1077
// When four quarks consider last one to be spectator.
1078
int iEnd = int( (nLeft - 1.) * rndmPtr->flat() );
1079
// Pick new flavour and form a new hadron.
1081
flavNew = flavSelPtr->pick( flavEnds[iEnd] );
1082
idHad = flavSelPtr->combine( flavEnds[iEnd], flavNew);
1083
} while (idHad == 0);
1084
// Store new hadron and endpoint flavour.
1085
idProd.push_back( idHad);
1086
flavEnds[iEnd].anti(flavNew);
1090
// When only two quarks left, combine to form final hadron.
1093
if ( abs(flavEnds[0].id) > 8 && abs(flavEnds[1].id) > 8)
1094
diquarkClash = true;
1096
do idHad = flavSelPtr->combine( flavEnds[0], flavEnds[1]);
1098
idProd.push_back( idHad);
1101
// If four quarks, decide how to pair them up.
1107
if ( rndmPtr->flat() < colRearrange) iEnd2 = 3;
1109
( (flavEnds[iEnd1].id > 0 && flavEnds[iEnd1].id < 9)
1110
|| flavEnds[iEnd1].id < -10 ) ? 1 : -1;
1111
if ( (flavEnds[iEnd2].id < 0 && flavEnds[iEnd2].id > -9)
1112
|| flavEnds[iEnd2].id > 10 ) relColSign *= -1;
1113
if (relColSign == 1) iEnd2 = 2;
1114
if (iEnd2 == 2) iEnd3 = 1;
1115
if (iEnd2 == 3) iEnd4 = 1;
1117
// Then combine to get final two hadrons.
1119
if ( abs(flavEnds[iEnd1].id) > 8 && abs(flavEnds[iEnd2].id) > 8)
1120
diquarkClash = true;
1122
do idHad = flavSelPtr->combine( flavEnds[iEnd1], flavEnds[iEnd2]);
1124
idProd.push_back( idHad);
1126
if ( abs(flavEnds[iEnd3].id) > 8 && abs(flavEnds[iEnd4].id) > 8)
1127
diquarkClash = true;
1129
do idHad = flavSelPtr->combine( flavEnds[iEnd3], flavEnds[iEnd4]);
1131
idProd.push_back( idHad);
1135
// Find masses of the new hadrons.
1136
for (int i = nFilled; i < int(idProd.size()) ; ++i) {
1137
double mHad = particleDataPtr->mass(idProd[i]);
1138
mProd.push_back( mHad);
1143
// Optional: check that this decay mode is not explicitly defined.
1144
if ( (meMode > 61 && meMode <= 80) && mDiff > mSafety && !diquarkClash ) {
1145
int idMatch[10], idPNow;
1146
usedChannel = false;
1147
bool matched = false;
1148
// Loop through all channels. Done if not same multiplicity.
1149
for (int i = 0; i < decDataPtr->sizeChannels(); ++i) {
1150
DecayChannel& channel = decDataPtr->channel(i);
1151
if (channel.multiplicity() != nTotal) continue;
1152
for (int k = 0; k < nTotal; ++k) idMatch[k] = channel.product(k);
1153
// Match particles one by one until fail.
1154
// Do not distinguish K0/K0bar/K0short/K0long at this stage.
1155
for (int j = 0; j < nTotal; ++j) {
1157
idPNow = idProd[j + 1];
1158
if (idPNow == -311 || idPNow == 130 || idPNow == 310) idPNow = 311;
1159
for (int k = 0; k < nTotal - j; ++k)
1160
if (idMatch[k] == idPNow || (idMatch[k] == -311 && idPNow == 311)) {
1161
// Compress list of unmatched when matching worked.
1162
idMatch[k] = idMatch[nTotal - 1 - j];
1166
if (!matched) break;
1168
// If matching worked, then chosen channel to be rejected.
1169
if (matched) {usedChannel = true; break;}
1173
// Keep on trying until enough phase space and no clash.
1174
} while (mDiff < mSafety || diquarkClash || usedChannel);
1176
// Update particle multiplicity.
1177
mult = idProd.size() - 1;
1179
// For Dalitz decays shuffle Dalitz pair to the end of the list.
1180
if (meMode == 11 || meMode == 12) {
1183
for (int i = 1; i <= mult; ++i) {
1184
if (idProd[i] == 11 || idProd[i] == 13 || idProd[i] == 15) iL1 = i;
1185
if (idProd[i] == -11 || idProd[i] == -13 || idProd[i] == -15) iL2 = i;
1187
if (iL1 > 0 && iL2 > 0) {
1188
int idL1 = idProd[iL1];
1189
int idL2 = idProd[iL2];
1190
double mL1 = mProd[iL1];
1191
double mL2 = mProd[iL2];
1193
for (int i = 1; i <= mult; ++i) if (i != iL1 && i != iL2) {
1195
idProd[iMove] = idProd[i];
1196
mProd[iMove] = mProd[i];
1198
idProd[mult - 1] = idL1;
1199
idProd[mult] = idL2;
1200
mProd[mult - 1] = mL1;
1210
//--------------------------------------------------------------------------
1212
// Set colour flow and scale in a decay explicitly to partons.
1214
bool ParticleDecays::setColours(Event& event) {
1216
// Decay to q qbar (or qbar q).
1217
if (meMode == 91 && idProd[1] > 0 && idProd[1] < 9) {
1218
int newCol = event.nextColTag();
1221
} else if (meMode == 91 && idProd[1] < 0 && idProd[1] > -9) {
1222
int newCol = event.nextColTag();
1227
} else if (meMode == 91 && idProd[1] == 21) {
1228
int newCol1 = event.nextColTag();
1229
int newCol2 = event.nextColTag();
1236
} else if (meMode == 92 && idProd[1] == 21 && idProd[2] == 21
1237
&& idProd[3] == 21) {
1238
int newCol1 = event.nextColTag();
1239
int newCol2 = event.nextColTag();
1240
int newCol3 = event.nextColTag();
1248
// Decay to g g gamma: locate which is gamma.
1249
} else if (meMode == 92) {
1250
int iGlu1 = (idProd[1] == 21) ? 1 : 3;
1251
int iGlu2 = (idProd[2] == 21) ? 2 : 3;
1252
int newCol1 = event.nextColTag();
1253
int newCol2 = event.nextColTag();
1254
cols[iGlu1] = newCol1;
1255
acols[iGlu1] = newCol2;
1256
cols[iGlu2] = newCol2;
1257
acols[iGlu2] = newCol1;
1259
// Unknown decay mode means failure.
1260
} else return false;
1262
// Set maximum scale to be mass of decaying particle.
1270
//==========================================================================
1272
} // end namespace Pythia8
added
removed
src/ParticleDecays.cc
