1
// ResonanceDecays.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
7
// the ResonanceDecays class.
9
#include "ResonanceDecays.h"
13
//==========================================================================
15
// The ResonanceDecays class.
16
// Do all resonance decays sequentially.
18
//--------------------------------------------------------------------------
20
// Constants: could be changed here if desired, but normally should not.
21
// These are of technical nature, as described for each.
23
// Number of tries to pick a decay channel.
24
const int ResonanceDecays::NTRYCHANNEL = 10;
26
// Number of tries to pick a set of daughter masses.
27
const int ResonanceDecays::NTRYMASSES = 10000;
29
// Mass above threshold for allowed decays.
30
const double ResonanceDecays::MSAFETY = 0.1;
32
// When constrainted kinematics cut high-mass tail of Breit-Wigner.
33
const double ResonanceDecays::WIDTHCUT = 5.;
35
// Small number (relative to 1) to protect against roundoff errors.
36
const double ResonanceDecays::TINY = 1e-10;
38
// Forbid small Breit-Wigner mass range, as mapped onto atan range.
39
const double ResonanceDecays::TINYBWRANGE = 1e-8;
41
// These numbers are hardwired empirical parameters,
42
// intended to speed up the M-generator.
43
const double ResonanceDecays::WTCORRECTION[11] = { 1., 1., 1.,
44
2., 5., 15., 60., 250., 1250., 7000., 50000. };
46
//--------------------------------------------------------------------------
48
bool ResonanceDecays::next( Event& process) {
50
// Loop over all entries to find resonances that should decay.
53
Particle& decayer = process[iDec];
54
if (decayer.isFinal() && decayer.canDecay() && decayer.mayDecay()
55
&& decayer.isResonance() ) {
57
// Fill the decaying particle in slot 0 of arrays.
62
idProd.push_back( id0 );
63
mProd.push_back( m0 );
65
// Mother flavour - relevant for gamma*/Z0 mixing. (Not always??)
66
int idIn = process[decayer.mother1()].id();
68
// Prepare decay selection.
69
decayer.particleDataEntry().preparePick(id0, m0, idIn);
71
// Pick a decay channel; allow up to ten tries.
72
bool foundChannel = false;
73
for (int iTryChannel = 0; iTryChannel < NTRYCHANNEL; ++iTryChannel) {
75
// Pick decay channel. Find multiplicity.
76
DecayChannel& channel = decayer.particleDataEntry().pickChannel();
77
mult = channel.multiplicity();
82
for (int i = 1; i <= mult; ++i) {
83
idNow = channel.product(i - 1);
84
if (id0 < 0 && particleDataPtr->hasAnti(idNow)) idNow = -idNow;
85
idProd.push_back( idNow);
88
// Pick masses. Pick new channel if fail.
89
if (!pickMasses()) continue;
94
// Failed to find acceptable decays.
96
infoPtr->errorMsg("Error in ResonanceDecays::next:"
97
" failed to find workable decay channel");
101
// Select colours in decay.
102
if (!pickColours(iDec, process)) return false;
104
// Select four-momenta in decay, boosted to lab frame.
106
pProd.push_back( decayer.p() );
107
if (!pickKinematics()) return false;
109
// Append decay products to the process event record. Set lifetimes.
110
int iFirst = process.size();
111
for (int i = 1; i <= mult; ++i) {
112
process.append( idProd[i], 23, iDec, 0, 0, 0, cols[i], acols[i],
113
pProd[i], mProd[i], m0);
115
int iLast = process.size() - 1;
117
// Set decay vertex when this is displaced.
118
if (process[iDec].hasVertex() || process[iDec].tau() > 0.) {
119
Vec4 vDec = process[iDec].vDec();
120
for (int i = iFirst; i <= iLast; ++i) process[i].vProd( vDec );
123
// Set lifetime of daughters.
124
for (int i = iFirst; i <= iLast; ++i)
125
process[i].tau( process[i].tau0() * rndmPtr->exp() );
127
// Modify mother status and daughters.
129
decayer.daughters(iFirst, iLast);
131
// End of loop over all entries.
133
} while (++iDec < process.size());
140
//--------------------------------------------------------------------------
142
// Select masses of decay products.
144
bool ResonanceDecays::pickMasses() {
146
// Arrays with properties of particles. Fill with dummy values for mother.
148
vector<double> m0BW, mMinBW, mMaxBW, widthBW;
149
double mMother = mProd[0];
150
double m2Mother = mMother * mMother;
151
useBW.push_back( false );
152
m0BW.push_back( mMother );
153
mMinBW.push_back( mMother );
154
mMaxBW.push_back( mMother );
155
widthBW.push_back( 0. );
157
// Loop throught products for masses and widths. Set nominal mass.
159
double m0Now, mMinNow, mMaxNow, widthNow;
160
for (int i = 1; i <= mult; ++i) {
161
useBWNow = particleDataPtr->useBreitWigner( idProd[i] );
162
m0Now = particleDataPtr->m0( idProd[i] );
163
mMinNow = particleDataPtr->m0Min( idProd[i] );
164
mMaxNow = particleDataPtr->m0Max( idProd[i] );
165
if (useBWNow && mMaxNow < mMinNow) mMaxNow = mMother;
166
widthNow = particleDataPtr->mWidth( idProd[i] );
167
useBW.push_back( useBWNow );
168
m0BW.push_back( m0Now );
169
mMinBW.push_back( mMinNow );
170
mMaxBW.push_back( mMaxNow );
171
widthBW.push_back( widthNow );
172
mProd.push_back( m0Now );
175
// Find number of Breit-Wigners and summed (minimal) masses.
179
for (int i = 1; i <= mult; ++i) {
181
mSum += max( m0BW[i], mMinBW[i]);
182
mSumMin += mMinBW[i];
185
// If sum of minimal masses above mother mass then give up.
186
if (mSumMin + MSAFETY > mMother) return false;
188
// If sum of masses below and no Breit-Wigners then done.
189
if (mSum + 0.5 * MSAFETY < mMother && nBW == 0) return true;
191
// Else if below then retry Breit-Wigners, with simple treshold.
192
if (mSum + MSAFETY < mMother) {
193
double wtMax = 2. * sqrtpos(1. - mSum*mSum / m2Mother);
195
for (int iTryMasses = 0; iTryMasses <= NTRYMASSES; ++ iTryMasses) {
196
if (iTryMasses == NTRYMASSES) return false;
198
for (int i = 1; i <= mult; ++i) {
199
if (useBW[i]) mProd[i] = particleDataPtr->mass( idProd[i] );
202
wt = (mSum + 0.5 * MSAFETY < mMother)
203
? sqrtpos(1. - mSum*mSum / m2Mother) : 0.;
204
if (wt > rndmPtr->flat() * wtMax) break;
209
// From now on some particles will have to be forced off shell.
211
// Order Breit-Wigners in decreasing widths. Sum of other masses.
214
for (int i = 1; i <= mult; ++i) {
215
if (useBW[i]) iBW.push_back(i);
216
else mSum0 += mProd[i];
218
for (int i = 1; i < nBW; ++i) {
219
for (int j = i - 1; j >= 0; --j) {
220
if (widthBW[iBW[j+1]] > widthBW[iBW[j]]) swap (iBW[j+1], iBW[j]);
225
// Do all but broadest two in increasing-width order. Includes only one.
227
int iMin = (nBW == 1) ? 0 : 2;
228
for (int i = nBW - 1; i >= iMin; --i) {
231
// Find allowed mass range of current resonance.
232
double mMax = mMother - mSum0 - MSAFETY;
233
if (nBW != 1) for (int j = 0; j < i; ++j) mMax -= mMinBW[iBW[j]];
234
mMax = min( mMaxBW[iBWi], mMax );
235
double mMin = min( mMinBW[iBWi], mMax - MSAFETY);
236
if (mMin < 0.) return false;
238
// Parameters for Breit-Wigner choice, with constrained mass range.
239
double m2Nom = pow2( m0BW[iBWi] );
240
double m2Max = mMax * mMax;
241
double m2Min = mMin * mMin;
242
double mmWid = m0BW[iBWi] * widthBW[iBWi];
243
double atanMin = atan( (m2Min - m2Nom) / mmWid );
244
double atanMax = atan( (m2Max - m2Nom) / mmWid );
245
double atanDif = atanMax - atanMin;
247
// Fail if too narrow mass range; e.g. out in tail of Breit-Wigner.
248
if (atanDif < TINYBWRANGE) return false;
250
// Retry mass according to Breit-Wigner, with simple threshold factor.
251
double mr1 = mSum0*mSum0 / m2Mother;
252
double mr2 = m2Min / m2Mother;
253
double wtMax = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
255
for (int iTryMasses = 0; iTryMasses <= NTRYMASSES; ++ iTryMasses) {
256
if (iTryMasses == NTRYMASSES) return false;
257
m2Now = m2Nom + mmWid * tan(atanMin + rndmPtr->flat() * atanDif);
258
mr2 = m2Now / m2Mother;
259
wt = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
260
if (wt > rndmPtr->flat() * wtMax) break;
263
// Prepare to iterate for more. Done for one Breit-Wigner.
264
mProd[iBWi] = sqrt(m2Now);
265
mSum0 += mProd[iBWi];
267
if (nBW == 1) return true;
270
// Left to do two broadest Breit-Wigners correlated, i.e. more realistic.
273
int idMother = abs(idProd[0]);
274
int idDau1 = abs(idProd[iBW1]);
275
int idDau2 = abs(idProd[iBW2]);
277
// In some cases known phase-space behaviour; else simple beta factor.
279
if ( (idMother == 25 || idMother == 35) && idDau1 < 19
280
&& idDau2 == idDau1 ) psMode = 3;
281
if ( (idMother == 25 || idMother == 35 )
282
&& (idDau1 == 23 || idDau1 == 24) && idDau2 == idDau1 ) psMode = 5;
284
&& (idDau1 == 23 || idDau1 == 24) && idDau2 == idDau1 ) psMode = 6;
286
// Find allowed mass ranges. Ensure that they are not closed.
287
double mRem = mMother - mSum0 - MSAFETY;
288
double mMax1 = min( mMaxBW[iBW1], mRem - mMinBW[iBW2] );
289
double mMin1 = min( mMinBW[iBW1], mMax1 - MSAFETY);
290
double mMax2 = min( mMaxBW[iBW2], mRem - mMinBW[iBW1] );
291
double mMin2 = min( mMinBW[iBW2], mMax2 - MSAFETY);
293
// At least one range must extend below half remaining mass.
294
if (mMin1 + mMin2 > mRem) return false;
295
double mMid = 0.5 * mRem;
296
bool hasMid1 = (mMin1 < mMid);
297
bool hasMid2 = (mMin2 < mMid);
298
if (!hasMid1 && !hasMid2) return false;
300
// Parameters for Breit-Wigner choice, with constrained mass range.
301
double m2Nom1 = pow2( m0BW[iBW1] );
302
double m2Max1 = mMax1 * mMax1;
303
double m2Min1 = mMin1 * mMin1;
304
double m2Mid1 = min( mMid * mMid, m2Max1);
305
double mmWid1 = m0BW[iBW1] * widthBW[iBW1];
306
double atanMin1 = atan( (m2Min1 - m2Nom1) / mmWid1 );
307
double atanMax1 = atan( (m2Max1 - m2Nom1) / mmWid1 );
308
double atanMid1 = (hasMid1) ? atan( (m2Mid1 - m2Nom1) / mmWid1 ) : 0.;
309
double m2Nom2 = pow2( m0BW[iBW2] );
310
double m2Max2 = mMax2 * mMax2;
311
double m2Min2 = mMin2 * mMin2;
312
double m2Mid2 = min( mMid * mMid, m2Max2);
313
double mmWid2 = m0BW[iBW2] * widthBW[iBW2];
314
double atanMin2 = atan( (m2Min2 - m2Nom2) / mmWid2 );
315
double atanMax2 = atan( (m2Max2 - m2Nom2) / mmWid2 );
316
double atanMid2 = (hasMid2) ? atan( (m2Mid2 - m2Nom2) / mmWid2 ) : 0.;
318
// Relative weight to pick either below half remaining mass.
319
double probLow1 = (hasMid1) ? 1. : 0.;
320
if (hasMid1 && hasMid2) {
321
double intLow1 = (atanMid1 - atanMin1) * (atanMax2 - atanMin2);
322
double intLow2 = (atanMax1 - atanMin1) * (atanMid2 - atanMin2);
323
probLow1 = intLow1 / (intLow1 + intLow2);
326
// Maximum matrix element times phase space weight.
327
double m2Rem = mRem * mRem;
328
double mr1 = m2Min1 / m2Rem;
329
double mr2 = m2Min2 / m2Rem;
330
double psMax = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
332
if (psMode == 1) wtMax = psMax;
333
else if (psMode == 2) wtMax = psMax * psMax;
334
else if (psMode == 3) wtMax = pow3(psMax);
335
else if (psMode == 5) wtMax = psMax
336
* (pow2(1. - mr1 - mr2) + 8. * mr1 * mr2);
337
else if (psMode == 6) wtMax = pow3(psMax);
339
// Retry mass according to Breit-Wigners, with simple threshold factor.
340
double atanDif1, atanDif2, m2Now1, m2Now2, mNow1, mNow2, ps, wt;
341
for (int iTryMasses = 0; iTryMasses <= NTRYMASSES; ++ iTryMasses) {
342
if (iTryMasses == NTRYMASSES) return false;
344
// Pick either below half remaining mass.
345
bool pickLow1 = false;
346
if (rndmPtr->flat() < probLow1) {
347
atanDif1 = atanMid1 - atanMin1;
348
atanDif2 = atanMax2 - atanMin2;
351
atanDif1 = atanMax1 - atanMin1;
352
atanDif2 = atanMid2 - atanMin2;
354
m2Now1 = m2Nom1 + mmWid1 * tan(atanMin1 + rndmPtr->flat() * atanDif1);
355
m2Now2 = m2Nom2 + mmWid2 * tan(atanMin2 + rndmPtr->flat() * atanDif2);
356
mNow1 = sqrt(m2Now1);
357
mNow2 = sqrt(m2Now2);
359
// Check that intended mass ordering is fulfilled.
360
bool rejectRegion = (pickLow1) ? (mNow1 > mNow2) : (mNow2 > mNow1);
361
if (rejectRegion) continue;
364
mr1 = m2Now1 / m2Rem;
365
mr2 = m2Now2 / m2Rem;
367
if (mNow1 + mNow2 + MSAFETY < mMother) {
368
ps = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
370
if (psMode == 1) wt = ps;
371
else if (psMode == 2) wt = ps * ps;
372
else if (psMode == 3) wt = pow3(ps);
373
else if (psMode == 5) wt = ps
374
* (pow2(1. - mr1 - mr2) + 8. * mr1 * mr2);
375
else if (psMode == 6) wt = pow3(ps)*mr1*mr2;
377
if (wt > rndmPtr->flat() * wtMax) break;
387
//--------------------------------------------------------------------------
389
// Select colours of decay products.
391
bool ResonanceDecays::pickColours(int iDec, Event& process) {
393
// Reset or create arrays with colour info.
396
vector<int> iTriplet, iAtriplet, iOctet, iDipCol, iDipAcol;
398
// Mother colours already known.
399
int col0 = process[iDec].col();
400
int acol0 = process[iDec].acol();
401
int colType0 = process[iDec].colType();
402
cols.push_back( col0);
403
acols.push_back(acol0);
405
// Loop through all daughters.
407
for (int i = 1; i <= mult; ++i) {
408
// Daughter colours initially empty, so that all is set for singlet.
411
// Find character (singlet, triplet, antitriplet, octet) of daughters.
412
colTypeNow = particleDataPtr->colType( idProd[i] );
413
if (colTypeNow == 0);
414
else if (colTypeNow == 1) iTriplet.push_back(i);
415
else if (colTypeNow == -1) iAtriplet.push_back(i);
416
else if (colTypeNow == 2) iOctet.push_back(i);
417
// Add two entries for sextets;
418
else if (colTypeNow == 3) {
419
iTriplet.push_back(i);
420
iTriplet.push_back(i);
421
} else if (colTypeNow == -3) {
422
iAtriplet.push_back(i);
423
iAtriplet.push_back(i);
425
infoPtr->errorMsg("Error in ResonanceDecays::pickColours:"
426
" unknown colour type encountered");
431
// Check excess of colours and anticolours in final over initial state.
432
int nCol = iTriplet.size();
433
if (colType0 == 1 || colType0 == 2) nCol -= 1;
434
else if (colType0 == 3) nCol -= 2;
435
int nAcol = iAtriplet.size();
436
if (colType0 == -1 || colType0 == 2) nAcol -= 1;
437
else if (colType0 == -3) nAcol -= 2;
439
// If net creation of three colours then find junction kind:
440
// mother is 1 = singlet, 3 = antitriplet, 5 = octet.
441
if (nCol - nAcol == 3) {
442
int kindJun = (col0 == 0) ? ((acol0 == 0) ? 1 : 3) : 5;
444
// Set colours in three junction legs and store junction.
446
colJun[0] = (kindJun == 1) ? process.nextColTag() : acol0;
447
colJun[1] = process.nextColTag();
448
colJun[2] = process.nextColTag();
449
process.appendJunction( kindJun, colJun[0], colJun[1], colJun[2]);
451
// Loop over three legs. Remove an incoming anticolour on first leg.
452
for (int leg = 0; leg < 3; ++leg) {
453
if (leg == 0 && kindJun != 1) acol0 = 0;
455
// Pick final-state triplets to carry these new colours.
457
int pickT = (iTriplet.size() == 1) ? 0
458
: int( TINY + rndmPtr->flat() * (iTriplet.size() - TINY) );
459
int iPickT = iTriplet[pickT];
460
cols[iPickT] = colJun[leg];
462
// Remove matched triplet and store new colour dipole ends.
463
iTriplet[pickT] = iTriplet.back();
465
iDipCol.push_back(iPickT);
466
iDipAcol.push_back(0);
470
// Update colour counter. Done with junction.
474
// If net creation of three anticolours then find antijunction kind:
475
// mother is 2 = singlet, 4 = triplet, 6 = octet.
476
if (nAcol - nCol == 3) {
477
int kindJun = (acol0 == 0) ? ((col0 == 0) ? 2 : 4) : 6;
479
// Set anticolours in three antijunction legs and store antijunction.
481
acolJun[0] = (kindJun == 2) ? process.nextColTag() : col0;
482
acolJun[1] = process.nextColTag();
483
acolJun[2] = process.nextColTag();
484
process.appendJunction( kindJun, acolJun[0], acolJun[1], acolJun[2]);
486
// Loop over three legs. Remove an incoming colour on first leg.
487
for (int leg = 0; leg < 3; ++leg) {
488
if (leg == 0 && kindJun != 2) col0 = 0;
490
// Pick final-state antitriplets to carry these new anticolours.
492
int pickA = (iAtriplet.size() == 1) ? 0
493
: int( TINY + rndmPtr->flat() * (iAtriplet.size() - TINY) );
494
int iPickA = iAtriplet[pickA];
495
acols[iPickA] = acolJun[leg];
497
// Remove matched antitriplet and store new colour dipole ends.
498
iAtriplet[pickA] = iAtriplet.back();
499
iAtriplet.pop_back();
500
iDipCol.push_back(0);
501
iDipAcol.push_back(iPickA);
505
// Update anticolour counter. Done with antijunction.
509
// If colours and anticolours do not match now then unphysical.
511
infoPtr->errorMsg("Error in ResonanceDecays::pickColours:"
512
" inconsistent colour tags");
516
// Pick final-state triplet (if any) to carry initial colour.
517
if (col0 > 0 && iTriplet.size() > 0) {
518
int pickT = (iTriplet.size() == 1) ? 0
519
: int( TINY + rndmPtr->flat() * (iTriplet.size() - TINY) );
520
int iPickT = iTriplet[pickT];
523
// Remove matched triplet and store new colour dipole ends.
525
iTriplet[pickT] = iTriplet.back();
527
iDipCol.push_back(iPickT);
528
iDipAcol.push_back(0);
531
// Pick final-state antitriplet (if any) to carry initial anticolour.
532
if (acol0 > 0 && iAtriplet.size() > 0) {
533
int pickA = (iAtriplet.size() == 1) ? 0
534
: int( TINY + rndmPtr->flat() * (iAtriplet.size() - TINY) );
535
int iPickA = iAtriplet[pickA];
536
acols[iPickA] = acol0;
538
// Remove matched antitriplet and store new colour dipole ends.
540
iAtriplet[pickA] = iAtriplet.back();
541
iAtriplet.pop_back();
542
iDipCol.push_back(0);
543
iDipAcol.push_back(iPickA);
546
// Sextets: second final-state triplet (if any)
547
if (acol0 < 0 && iTriplet.size() > 0) {
548
int pickT = (iTriplet.size() == 1) ? 0
549
: int( TINY + rndmPtr->flat() * (iTriplet.size() - TINY) );
550
int iPickT = iTriplet[pickT];
551
cols[iPickT] = -acol0;
553
// Remove matched antitriplet and store new colour dipole ends.
555
iTriplet[pickT] = iTriplet.back();
557
iDipCol.push_back(iPickT);
558
iDipAcol.push_back(0);
561
// Sextets: second final-state antitriplet (if any)
562
if (col0 < 0 && iAtriplet.size() > 0) {
563
int pickA = (iAtriplet.size() == 1) ? 0
564
: int( TINY + rndmPtr->flat() * (iAtriplet.size() - TINY) );
565
int iPickA = iAtriplet[pickA];
566
acols[iPickA] = -col0;
568
// Remove matched triplet and store new colour dipole ends.
570
iAtriplet[pickA] = iAtriplet.back();
571
iAtriplet.pop_back();
572
iDipCol.push_back(0);
573
iDipAcol.push_back(iPickA);
576
// Error checks that amount of leftover colours and anticolours match.
577
if ( (iTriplet.size() != iAtriplet.size())
578
|| (col0 != 0 && acol0 == 0) || (col0 == 0 && acol0 != 0) ) {
579
infoPtr->errorMsg("Error in ResonanceDecays::pickColours:"
580
" inconsistent colour tags");
584
// Match triplets to antitriplets in the final state.
585
for (int pickT = 0; pickT < int(iTriplet.size()); ++pickT) {
586
int iPickT = iTriplet[pickT];
587
int pickA = (iAtriplet.size() == 1) ? 0
588
: int( TINY + rndmPtr->flat() * (iAtriplet.size() - TINY) );
589
int iPickA = iAtriplet[pickA];
591
// Connect pair with new colour tag.
592
cols[iPickT] = process.nextColTag();
593
acols[iPickA] = cols[iPickT];
595
// Remove matched antitriplet and store new colour dipole ends.
596
iAtriplet[pickT] = iAtriplet.back();
597
iAtriplet.pop_back();
598
iDipCol.push_back(iPickT);
599
iDipAcol.push_back(iPickA);
602
// If no octets are around then matching is done.
603
if (col0 == 0 && acol0 == 0 && iOctet.size() == 0) return true;
605
// If initial-state octet remains then store as (first!) new dipole.
607
iDipCol.push_back(0);
608
iDipAcol.push_back(0);
611
// Now attach all final-state octets at random to existing dipoles.
612
for (int i = 0; i < int(iOctet.size()); ++i) {
613
int iOct = iOctet[i];
615
// If no dipole then start new one. (Happens for singlet -> octets.)
616
if (iDipCol.size() == 0) {
617
cols[iOct] = process.nextColTag();
618
acols[iOct] = cols[iOct] ;
619
iDipCol.push_back(iOct);
620
iDipAcol.push_back(iOct);
623
// Else attach to existing dipole picked at random.
625
int pickDip = (iDipCol.size() == 1) ? 0
626
: int( TINY + rndmPtr->flat() * (iDipCol.size() - TINY) );
628
// Case with dipole in initial state: reattach existing colours.
629
if (iDipCol[pickDip] == 0 && iDipAcol[pickDip] == 0) {
632
iDipAcol[pickDip] = iOct;
633
iDipCol.push_back(iOct);
634
iDipAcol.push_back(0);
636
// Case with dipole from colour in initial state: also new colour.
637
} else if (iDipAcol[pickDip] == 0) {
638
int iPickCol = iDipCol[pickDip];
639
cols[iOct] = cols[iPickCol];
640
acols[iOct] = process.nextColTag();
641
cols[iPickCol] = acols[iOct];
642
iDipCol[pickDip] = iOct;
643
iDipCol.push_back(iPickCol);
644
iDipAcol.push_back(iOct);
646
// Remaining cases with dipole from anticolour in initial state
647
// or dipole inside final state: also new colour.
649
int iPickAcol = iDipAcol[pickDip];
650
acols[iOct] = acols[iPickAcol];
651
cols[iOct] = process.nextColTag();
652
acols[iPickAcol] = cols[iOct];
653
iDipAcol[pickDip] = iOct;
654
iDipCol.push_back(iOct);
655
iDipAcol.push_back(iPickAcol);
660
// Must now have at least two dipoles (no 1 -> 8 or 8 -> 1).
661
if (iDipCol.size() < 2) {
662
infoPtr->errorMsg("Error in ResonanceDecays::pickColours:"
663
" inconsistent colour tags");
672
//--------------------------------------------------------------------------
674
// Select decay products momenta isotropically in phase space.
675
// Process-dependent angular distributions may be imposed in SigmaProcess.
677
bool ResonanceDecays::pickKinematics() {
679
// Description of two-body decays as simple special case.
684
double m1 = mProd[1];
685
double m2 = mProd[2];
687
// Energies and absolute momentum in the rest frame.
688
double e1 = 0.5 * (m0*m0 + m1*m1 - m2*m2) / m0;
689
double e2 = 0.5 * (m0*m0 + m2*m2 - m1*m1) / m0;
690
double pAbs = 0.5 * sqrtpos( (m0 - m1 - m2) * (m0 + m1 + m2)
691
* (m0 + m1 - m2) * (m0 - m1 + m2) ) / m0;
693
// Pick isotropic angles to give three-momentum.
694
double cosTheta = 2. * rndmPtr->flat() - 1.;
695
double sinTheta = sqrt(1. - cosTheta*cosTheta);
696
double phi = 2. * M_PI * rndmPtr->flat();
697
double pX = pAbs * sinTheta * cos(phi);
698
double pY = pAbs * sinTheta * sin(phi);
699
double pZ = pAbs * cosTheta;
701
// Fill four-momenta in mother rest frame and then boost to lab frame.
702
pProd.push_back( Vec4( pX, pY, pZ, e1) );
703
pProd.push_back( Vec4( -pX, -pY, -pZ, e2) );
704
pProd[1].bst( pProd[0] );
705
pProd[2].bst( pProd[0] );
707
// Done for two-body decay.
711
// Description of three-body decays as semi-simple special case.
716
double m1 = mProd[1];
717
double m2 = mProd[2];
718
double m3 = mProd[3];
719
double mDiff = m0 - (m1 + m2 + m3);
721
// Kinematical limits for 2+3 mass. Maximum phase-space weight.
722
double m23Min = m2 + m3;
723
double m23Max = m0 - m1;
724
double p1Max = 0.5 * sqrtpos( (m0 - m1 - m23Min) * (m0 + m1 + m23Min)
725
* (m0 + m1 - m23Min) * (m0 - m1 + m23Min) ) / m0;
726
double p23Max = 0.5 * sqrtpos( (m23Max - m2 - m3) * (m23Max + m2 + m3)
727
* (m23Max + m2 - m3) * (m23Max - m2 + m3) ) / m23Max;
728
double wtPSmax = 0.5 * p1Max * p23Max;
730
// Pick an intermediate mass m23 flat in the allowed range.
731
double wtPS, m23, p1Abs, p23Abs;
733
m23 = m23Min + rndmPtr->flat() * mDiff;
735
// Translate into relative momenta and find phase-space weight.
736
p1Abs = 0.5 * sqrtpos( (m0 - m1 - m23) * (m0 + m1 + m23)
737
* (m0 + m1 - m23) * (m0 - m1 + m23) ) / m0;
738
p23Abs = 0.5 * sqrtpos( (m23 - m2 - m3) * (m23 + m2 + m3)
739
* (m23 + m2 - m3) * (m23 - m2 + m3) ) / m23;
740
wtPS = p1Abs * p23Abs;
742
// If rejected, try again with new invariant masses.
743
} while ( wtPS < rndmPtr->flat() * wtPSmax );
745
// Set up m23 -> m2 + m3 isotropic in its rest frame.
746
double cosTheta = 2. * rndmPtr->flat() - 1.;
747
double sinTheta = sqrt(1. - cosTheta*cosTheta);
748
double phi = 2. * M_PI * rndmPtr->flat();
749
double pX = p23Abs * sinTheta * cos(phi);
750
double pY = p23Abs * sinTheta * sin(phi);
751
double pZ = p23Abs * cosTheta;
752
double e2 = sqrt( m2*m2 + p23Abs*p23Abs);
753
double e3 = sqrt( m3*m3 + p23Abs*p23Abs);
754
Vec4 p2( pX, pY, pZ, e2);
755
Vec4 p3( -pX, -pY, -pZ, e3);
757
// Set up 0 -> 1 + 23 isotropic in its rest frame.
758
cosTheta = 2. * rndmPtr->flat() - 1.;
759
sinTheta = sqrt(1. - cosTheta*cosTheta);
760
phi = 2. * M_PI * rndmPtr->flat();
761
pX = p1Abs * sinTheta * cos(phi);
762
pY = p1Abs * sinTheta * sin(phi);
763
pZ = p1Abs * cosTheta;
764
double e1 = sqrt( m1*m1 + p1Abs*p1Abs);
765
double e23 = sqrt( m23*m23 + p1Abs*p1Abs);
766
pProd.push_back( Vec4( pX, pY, pZ, e1) );
768
// Boost 2 + 3 to the 0 rest frame and then boost to lab frame.
769
Vec4 p23( -pX, -pY, -pZ, e23);
772
pProd.push_back( p2 );
773
pProd.push_back( p3 );
774
pProd[1].bst( pProd[0] );
775
pProd[2].bst( pProd[0] );
776
pProd[3].bst( pProd[0] );
778
// Done for three-body decay.
782
// Do a multibody decay using the M-generator algorithm.
784
// Mother and sum daughter masses.
786
double mSum = mProd[1];
787
for (int i = 2; i <= mult; ++i) mSum += mProd[i];
788
double mDiff = m0 - mSum;
790
// Begin setup of intermediate invariant masses.
792
for (int i = 0; i <= mult; ++i) mInv.push_back( mProd[i]);
794
// Calculate the maximum weight in the decay.
795
double wtPSmax = 1. / WTCORRECTION[mult];
796
double mMax = mDiff + mProd[mult];
798
for (int i = mult - 1; i > 0; --i) {
801
double mNow = mProd[i];
802
wtPSmax *= 0.5 * sqrtpos( (mMax - mMin - mNow) * (mMax + mMin + mNow)
803
* (mMax + mMin - mNow) * (mMax - mMin + mNow) ) / mMax;
806
// Begin loop to find the set of intermediate invariant masses.
807
vector<double> rndmOrd;
812
// Find and order random numbers in descending order.
814
rndmOrd.push_back(1.);
815
for (int i = 1; i < mult - 1; ++i) {
816
double rndm = rndmPtr->flat();
817
rndmOrd.push_back(rndm);
818
for (int j = i - 1; j > 0; --j) {
819
if (rndm > rndmOrd[j]) swap( rndmOrd[j], rndmOrd[j+1] );
823
rndmOrd.push_back(0.);
825
// Translate into intermediate masses and find weight.
826
for (int i = mult - 1; i > 0; --i) {
827
mInv[i] = mInv[i+1] + mProd[i] + (rndmOrd[i-1] - rndmOrd[i]) * mDiff;
828
wtPS *= 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
829
* (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
830
* (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
833
// If rejected, try again with new invariant masses.
834
} while ( wtPS < rndmPtr->flat() * wtPSmax );
836
// Perform two-particle decays in the respective rest frame.
838
pInv.resize(mult + 1);
839
for (int i = 1; i < mult; ++i) {
840
double pAbs = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
841
* (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
842
* (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
844
// Isotropic angles give three-momentum.
845
double cosTheta = 2. * rndmPtr->flat() - 1.;
846
double sinTheta = sqrt(1. - cosTheta*cosTheta);
847
double phi = 2. * M_PI * rndmPtr->flat();
848
double pX = pAbs * sinTheta * cos(phi);
849
double pY = pAbs * sinTheta * sin(phi);
850
double pZ = pAbs * cosTheta;
852
// Calculate energies, fill four-momenta.
853
double eHad = sqrt( mProd[i]*mProd[i] + pAbs*pAbs);
854
double eInv = sqrt( mInv[i+1]*mInv[i+1] + pAbs*pAbs);
855
pProd.push_back( Vec4( pX, pY, pZ, eHad) );
856
pInv[i+1].p( -pX, -pY, -pZ, eInv);
858
pProd.push_back( pInv[mult] );
860
// Boost decay products to the mother rest frame and on to lab frame.
862
for (int iFrame = mult - 1; iFrame > 0; --iFrame)
863
for (int i = iFrame; i <= mult; ++i) pProd[i].bst(pInv[iFrame]);
865
// Done for multibody decay.
870
//==========================================================================
872
} // end namespace Pythia8
added
removed
src/ResonanceDecays.cc
