using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / src / ResonanceDecays.cxx
CommitLineData
5ad4eb21 1// ResonanceDecays.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2008 Torbjorn Sjostrand.
3// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4// Please respect the MCnet Guidelines, see GUIDELINES for details.
5
6// Function definitions (not found in the header) for
7// the ResonanceDecays class.
8
9#include "ResonanceDecays.h"
10
11namespace Pythia8 {
12
13//**************************************************************************
14
15// The ResonanceDecays class.
16// Do all resonance decays sequentially.
17
18//*********
19
20// Constants: could be changed here if desired, but normally should not.
21// These are of technical nature, as described for each.
22
23// Number of tries to pick a decay channel.
24const int ResonanceDecays::NTRYCHANNEL = 10;
25
26// Number of tries to pick a set of daughter masses.
27const int ResonanceDecays::NTRYMASSES = 10000;
28
29// Mass above threshold for allowed decays.
30const double ResonanceDecays::MSAFETY = 0.1;
31
32// When constrainted kinematics cut high-mass tail of Breit-Wigner.
33const double ResonanceDecays::WIDTHCUT = 5.;
34
35// Small number (relative to 1) to protect against roundoff errors.
36const double ResonanceDecays::TINY = 1e-10;
37
38// These numbers are hardwired empirical parameters,
39// intended to speed up the M-generator.
40const double ResonanceDecays::WTCORRECTION[11] = { 1., 1., 1.,
41 2., 5., 15., 60., 250., 1250., 7000., 50000. };
42
43//*********
44
45bool ResonanceDecays::next( Event& process) {
46
47 // Loop over all entries to find resonances that should decay.
48 int iDec = 0;
49 do {
50 Particle& decayer = process[iDec];
51 if (decayer.isFinal() && decayer.canDecay() && decayer.mayDecay()
52 && decayer.isResonance() ) {
53
54 // Fill the decaying particle in slot 0 of arrays.
55 id0 = decayer.id();
56 m0 = decayer.m();
57 idProd.resize(0);
58 mProd.resize(0);
59 idProd.push_back( id0 );
60 mProd.push_back( m0 );
61
62 // Mother flavour - relevant for gamma*/Z0 mixing. (Not always??)
63 int idIn = process[decayer.mother1()].id();
64
65 // Prepare decay selection.
66 decayer.particleData().preparePick(id0, m0, idIn);
67
68 // Pick a decay channel; allow up to ten tries.
69 bool foundChannel = false;
70 for (int iTryChannel = 0; iTryChannel < NTRYCHANNEL; ++iTryChannel) {
71
72 // Pick decay channel. Find multiplicity.
73 DecayChannel& channel = decayer.particleData().pickChannel();
74 mult = channel.multiplicity();
75
76 // Read out flavours.
77 idProd.resize(1);
78 int idNow;
79 for (int i = 1; i <= mult; ++i) {
80 idNow = channel.product(i - 1);
81 if (id0 < 0 && ParticleDataTable::hasAnti(idNow)) idNow = -idNow;
82 idProd.push_back( idNow);
83 }
84
85 // Pick masses. Pick new channel if fail.
86 if (!pickMasses()) continue;
87 foundChannel = true;
88 break;
89 }
90
91 // Failed to find acceptable decays.
92 if (!foundChannel) {
93 infoPtr->errorMsg("Error in ResonanceDecays::next:"
94 " failed to find workable decay channel");
95 return false;
96 }
97
98 // Select colours in decay.
99 if (!pickColours(iDec, process)) return false;
100
101 // Select four-momenta in decay, boosted to lab frame.
102 pProd.resize(0);
103 pProd.push_back( decayer.p() );
104 if (!pickKinematics()) return false;
105
106 // Append decay products to the process event record.
107 int iFirst = process.size();
108 for (int i = 1; i <= mult; ++i) {
109 process.append( idProd[i], 23, iDec, 0, 0, 0, cols[i], acols[i],
110 pProd[i], mProd[i], m0);
111 }
112 int iLast = process.size() - 1;
113
114 // Modify mother status and daughters.
115 decayer.status(-22);
116 decayer.daughters(iFirst, iLast);
117
118 // End of loop over all entries.
119 }
120 } while (++iDec < process.size());
121
122 // Done.
123 return true;
124
125}
126
127// Select masses of decay products.
128
129bool ResonanceDecays::pickMasses() {
130
131 // Arrays with properties of particles. Fill with dummy values for mother.
132 vector<bool> useBW;
133 vector<double> m0BW, mMinBW, mMaxBW, widthBW;
134 double mMother = mProd[0];
135 double m2Mother = mMother * mMother;
136 useBW.push_back( false );
137 m0BW.push_back( mMother );
138 mMinBW.push_back( mMother );
139 mMaxBW.push_back( mMother );
140 widthBW.push_back( 0. );
141
142 // Loop throught products for masses and widths. Set nominal mass.
143 bool useBWNow;
144 double m0Now, mMinNow, mMaxNow, widthNow;
145 for (int i = 1; i <= mult; ++i) {
146 useBWNow = ParticleDataTable::useBreitWigner( idProd[i] );
147 m0Now = ParticleDataTable::m0( idProd[i] );
148 mMinNow = ParticleDataTable::m0Min( idProd[i] );
149 mMaxNow = ParticleDataTable::m0Max( idProd[i] );
150 if (useBWNow && mMaxNow < mMinNow) mMaxNow = mMother;
151 widthNow = ParticleDataTable::mWidth( idProd[i] );
152 useBW.push_back( useBWNow );
153 m0BW.push_back( m0Now );
154 mMinBW.push_back( mMinNow );
155 mMaxBW.push_back( mMaxNow );
156 widthBW.push_back( widthNow );
157 mProd.push_back( m0Now );
158 }
159
160 // Find number of Breit-Wigners and summed (minimal) masses.
161 int nBW = 0;
162 double mSum = 0.;
163 double mSumMin = 0.;
164 for (int i = 1; i <= mult; ++i) {
165 if (useBW[i]) ++nBW;
166 mSum += max( m0BW[i], mMinBW[i]);
167 mSumMin += mMinBW[i];
168 }
169
170 // If sum of minimal masses above mother mass then give up.
171 if (mSumMin + MSAFETY > mMother) return false;
172
173 // If sum of masses below and no Breit-Wigners then done.
174 if (mSum + 0.5 * MSAFETY < mMother && nBW == 0) return true;
175
176 // Else if below then retry Breit-Wigners, with simple treshold.
177 if (mSum + MSAFETY < mMother) {
178 double wtMax = 2. * sqrtpos(1. - mSum*mSum / m2Mother);
179 double wt;
180 for (int iTryMasses = 0; iTryMasses <= NTRYMASSES; ++ iTryMasses) {
181 if (iTryMasses == NTRYMASSES) return false;
182 mSum = 0.;
183 for (int i = 1; i <= mult; ++i) {
184 if (useBW[i]) mProd[i] = ParticleDataTable::mass( idProd[i] );
185 mSum += mProd[i];
186 }
187 wt = (mSum + 0.5 * MSAFETY < mMother)
188 ? sqrtpos(1. - mSum*mSum / m2Mother) : 0.;
189 if (wt > Rndm::flat() * wtMax) break;
190 }
191 return true;
192 }
193
194 // From now on some particles will have to be forced off shell.
195
196 // Order Breit-Wigners in decreasing widths. Sum of other masses.
197 vector<int> iBW;
198 double mSum0 = 0.;
199 for (int i = 1; i <= mult; ++i) {
200 if (useBW[i]) iBW.push_back(i);
201 else mSum0 += mProd[i];
202 }
203 for (int i = 1; i < nBW; ++i) {
204 for (int j = i - 1; j >= 0; --j) {
205 if (widthBW[iBW[j+1]] > widthBW[iBW[j]]) swap (iBW[j+1], iBW[j]);
206 else break;
207 }
208 }
209
210 // Do all but broadest two in increasing-width order. Includes only one.
211 if (nBW != 2) {
212 int iMin = (nBW == 1) ? 0 : 2;
213 for (int i = nBW - 1; i >= iMin; --i) {
214 int iBWi = iBW[i];
215
216 // Find allowed mass range of current resonance.
217 double mMax = mMother - mSum0 - MSAFETY;
218 if (nBW != 1) for (int j = 0; j < i; ++j) mMax -= mMinBW[iBW[j]];
219 mMax = min( mMaxBW[iBWi], mMax );
220 double mMin = min( mMinBW[iBWi], mMax - MSAFETY);
221 if (mMin < 0.) return false;
222
223 // Parameters for Breit-Wigner choice, with constrained mass range.
224 double m2Nom = pow2( m0BW[iBWi] );
225 double m2Max = mMax * mMax;
226 double m2Min = mMin * mMin;
227 double mmWid = m0BW[iBWi] * widthBW[iBWi];
228 double atanMin = atan( (m2Min - m2Nom) / mmWid );
229 double atanMax = atan( (m2Max - m2Nom) / mmWid );
230 double atanDif = atanMax - atanMin;
231
232 // Retry mass according to Breit-Wigner, with simple threshold factor.
233 double mr1 = mSum0*mSum0 / m2Mother;
234 double mr2 = m2Min / m2Mother;
235 double wtMax = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
236 double m2Now, wt;
237 for (int iTryMasses = 0; iTryMasses <= NTRYMASSES; ++ iTryMasses) {
238 if (iTryMasses == NTRYMASSES) return false;
239 m2Now = m2Nom + mmWid * tan(atanMin + Rndm::flat() * atanDif);
240 mr2 = m2Now / m2Mother;
241 wt = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
242 if (wt > Rndm::flat() * wtMax) break;
243 }
244
245 // Prepare to iterate for more. Done for one Breit-Wigner.
246 mProd[iBWi] = sqrt(m2Now);
247 mSum0 += mProd[iBWi];
248 }
249 if (nBW == 1) return true;
250 }
251
252 // Left to do two broadest Breit-Wigners correlated, i.e. more realistic.
253 int iBW1 = iBW[0];
254 int iBW2 = iBW[1];
255 int idMother = abs(idProd[0]);
256 int idDau1 = abs(idProd[iBW1]);
257 int idDau2 = abs(idProd[iBW2]);
258
259 // In some cases known phase-space behaviour; else simple beta factor.
260 int psMode = 1 ;
261 if ( (idMother == 25 || idMother == 35) && idDau1 < 19
262 && idDau2 == idDau1 ) psMode = 3;
263 if ( (idMother == 25 || idMother == 35 || idMother == 36)
264 && (idDau1 == 23 || idDau1 == 24) && idDau2 == idDau1 ) psMode = 5;
265
266 // Find allowed mass ranges
267 double mRem = mMother - mSum0 - MSAFETY;
268 double mMax1 = min( mMaxBW[iBW1], mRem - mMinBW[iBW2] );
269 double mMin1 = min( mMinBW[iBW1], mMax1 - MSAFETY);
270 double mMax2 = min( mMaxBW[iBW2], mRem - mMinBW[iBW1] );
271 double mMin2 = min( mMinBW[iBW2], mMax2 - MSAFETY);
272
273 // At least one range must extend below half remaining mass.
274 double mMid = 0.5 * mRem;
275 bool hasMid1 = (mMin1 < 0.5 * mRem);
276 bool hasMid2 = (mMin2 < 0.5 * mRem);
277 if (!hasMid1 && !hasMid2) return false;
278
279 // Parameters for Breit-Wigner choice, with constrained mass range.
280 double m2Nom1 = pow2( m0BW[iBW1] );
281 double m2Max1 = mMax1 * mMax1;
282 double m2Min1 = mMin1 * mMin1;
283 double mmWid1 = m0BW[iBW1] * widthBW[iBW1];
284 double atanMin1 = atan( (m2Min1 - m2Nom1) / mmWid1 );
285 double atanMax1 = atan( (m2Max1 - m2Nom1) / mmWid1 );
286 double atanMid1 = (hasMid1) ? atan( (mMid*mMid - m2Nom1) / mmWid1 ) : 0.;
287 double m2Nom2 = pow2( m0BW[iBW2] );
288 double m2Max2 = mMax1 * mMax2;
289 double m2Min2 = mMin1 * mMin2;
290 double mmWid2 = m0BW[iBW2] * widthBW[iBW2];
291 double atanMin2 = atan( (m2Min2 - m2Nom2) / mmWid2 );
292 double atanMax2 = atan( (m2Max2 - m2Nom2) / mmWid2 );
293 double atanMid2 = (hasMid2) ? atan( (mMid*mMid - m2Nom2) / mmWid2 ) : 0.;
294
295 // Relative weight to pick either below half remaining mass.
296 double probLow1 = (hasMid1) ? 1. : 0.;
297 if (hasMid1 && hasMid2) probLow1 = (atanMid1 - atanMin1)
298 / ( (atanMid1 - atanMin1) + (atanMid2 - atanMin1) );
299
300 // Maximum matrix element times phase space weight.
301 double m2Rem = mRem * mRem;
302 double mr1 = m2Min1 / m2Rem;
303 double mr2 = m2Min2 / m2Rem;
304 double psMax = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
305 double wtMax = 1.;
306 if (psMode == 1) wtMax = psMax;
307 else if (psMode == 2) wtMax = psMax * psMax;
308 else if (psMode == 3) wtMax = pow3(psMax);
309 else if (psMode == 5) wtMax = psMax
310 * (pow2(1. - mr1 - mr2) + 8. * mr1 * mr2);
311
312 // Retry mass according to Breit-Wigners, with simple threshold factor.
313 double atanDif1, atanDif2, m2Now1, m2Now2, mNow1, mNow2, ps, wt;
314 for (int iTryMasses = 0; iTryMasses <= NTRYMASSES; ++ iTryMasses) {
315 if (iTryMasses == NTRYMASSES) return false;
316
317 // Pick either below half remaining mass.
318 if (Rndm::flat() < probLow1) {
319 atanDif1 = atanMid1 - atanMin1;
320 atanDif2 = atanMax2 - atanMin2;
321 } else {
322 atanDif1 = atanMax1 - atanMin1;
323 atanDif2 = atanMid2 - atanMin2;
324 }
325 m2Now1 = m2Nom1 + mmWid1 * tan(atanMin1 + Rndm::flat() * atanDif1);
326 m2Now2 = m2Nom2 + mmWid2 * tan(atanMin2 + Rndm::flat() * atanDif2);
327 mNow1 = sqrt(m2Now1);
328 mNow2 = sqrt(m2Now2);
329
330 // Threshold weight.
331 mr1 = m2Now1 / m2Rem;
332 mr2 = m2Now2 / m2Rem;
333 wt = 0.;
334 if (mNow1 + mNow2 + MSAFETY < mMother) {
335 ps = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
336 wt = 1.;
337 if (psMode == 1) wt = ps;
338 else if (psMode == 2) wt = ps * ps;
339 else if (psMode == 3) wt = pow3(ps);
340 else if (psMode == 5) wt = ps
341 * (pow2(1. - mr1 - mr2) + 8. * mr1 * mr2);
342 }
343 if (wt > Rndm::flat() * wtMax) break;
344 }
345 mProd[iBW1] = mNow1;
346 mProd[iBW2] = mNow2;
347
348 // Done.
349 return true;
350
351}
352
353//*********
354
355// Select colours of decay products.
356
357bool ResonanceDecays::pickColours(int iDec, Event& process) {
358
359 // Reset or create arrays with colour info.
360 cols.resize(0);
361 acols.resize(0);
362 vector<int> iTriplet, iAtriplet, iOctet, iDipCol, iDipAcol;
363
364 // Mother colours already known.
365 int col0 = process[iDec].col();
366 int acol0 = process[iDec].acol();
367 cols.push_back( col0);
368 acols.push_back(acol0);
369
370 // Loop through all daughters.
371 int colTypeNow;
372 for (int i = 1; i <= mult; ++i) {
373 // Daughter colours initially empty, so that all is set for singlet.
374 cols.push_back(0);
375 acols.push_back(0);
376 // Find character (singlet, triplet, antitriplet, octet) of daughters.
377 colTypeNow = ParticleDataTable::colType( idProd[i] );
378 if (colTypeNow == 0);
379 else if (colTypeNow == 1) iTriplet.push_back(i);
380 else if (colTypeNow == -1) iAtriplet.push_back(i);
381 else if (colTypeNow == 2) iOctet.push_back(i);
382 else {
383 infoPtr->errorMsg("Error in ResonanceDecays::pickColours:"
384 " unknown colour type encountered");
385 return false;
386 }
387 }
388
389 // Check excess of colours and anticolours in final over initial state.
390 int nCol = iTriplet.size();
391 if (col0 != 0) --nCol;
392 int nAcol = iAtriplet.size();
393 if (acol0 != 0) --nAcol;
394
395 // If net creation of three colours then find junction kind:
396 // mother is 1 = singlet, 3 = antitriplet, 5 = octet.
397 if (nCol - nAcol == 3) {
398 int kindJun = (col0 == 0) ? ((acol0 == 0) ? 1 : 3) : 5;
399
400 // Set colours in three junction legs and store junction.
401 int colJun[3];
402 colJun[0] = (kindJun == 1) ? process.nextColTag() : acol0;
403 colJun[1] = process.nextColTag();
404 colJun[2] = process.nextColTag();
405 process.appendJunction( kindJun, colJun[0], colJun[1], colJun[2]);
406
407 // Loop over three legs. Remove an incoming anticolour on first leg.
408 for (int leg = 0; leg < 3; ++leg) {
409 if (leg == 0 && kindJun != 1) acol0 = 0;
410
411 // Pick final-state triplets to carry these new colours.
412 else {
413 int pickT = (iTriplet.size() == 1) ? 0
414 : int( TINY + Rndm::flat() * (iTriplet.size() - TINY) );
415 int iPickT = iTriplet[pickT];
416 cols[iPickT] = colJun[leg];
417
418 // Remove matched triplet and store new colour dipole ends.
419 iTriplet[pickT] = iTriplet.back();
420 iTriplet.pop_back();
421 iDipCol.push_back(iPickT);
422 iDipAcol.push_back(0);
423 }
424 }
425
426 // Update colour counter. Done with junction.
427 nCol -= 3;
428 }
429
430 // If net creation of three anticolours then find antijunction kind:
431 // mother is 2 = singlet, 4 = triplet, 6 = octet.
432 if (nAcol - nCol == 3) {
433 int kindJun = (acol0 == 0) ? ((col0 == 0) ? 2 : 4) : 6;
434
435 // Set anticolours in three antijunction legs and store antijunction.
436 int acolJun[3];
437 acolJun[0] = (kindJun == 2) ? process.nextColTag() : col0;
438 acolJun[1] = process.nextColTag();
439 acolJun[2] = process.nextColTag();
440 process.appendJunction( kindJun, acolJun[0], acolJun[1], acolJun[2]);
441
442 // Loop over three legs. Remove an incoming colour on first leg.
443 for (int leg = 0; leg < 3; ++leg) {
444 if (leg == 0 && kindJun != 2) col0 = 0;
445
446 // Pick final-state antitriplets to carry these new anticolours.
447 else {
448 int pickA = (iAtriplet.size() == 1) ? 0
449 : int( TINY + Rndm::flat() * (iAtriplet.size() - TINY) );
450 int iPickA = iAtriplet[pickA];
451 acols[iPickA] = acolJun[leg];
452
453 // Remove matched antitriplet and store new colour dipole ends.
454 iAtriplet[pickA] = iAtriplet.back();
455 iAtriplet.pop_back();
456 iDipCol.push_back(0);
457 iDipAcol.push_back(iPickA);
458 }
459 }
460
461 // Update anticolour counter. Done with antijunction.
462 nAcol -= 3;
463 }
464
465 // If colours and anticolours do not match now then unphysical.
466 if (nCol != nAcol) {
467 infoPtr->errorMsg("Error in ResonanceDecays::pickColours:"
468 " inconsistent colour tags");
469 return false;
470 }
471
472 // Pick final-state triplet (if any) to carry initial colour.
473 if (col0 != 0 && iTriplet.size() > 0) {
474 int pickT = (iTriplet.size() == 1) ? 0
475 : int( TINY + Rndm::flat() * (iTriplet.size() - TINY) );
476 int iPickT = iTriplet[pickT];
477 cols[iPickT] = col0;
478
479 // Remove matched triplet and store new colour dipole ends.
480 col0 = 0;
481 iTriplet[pickT] = iTriplet.back();
482 iTriplet.pop_back();
483 iDipCol.push_back(iPickT);
484 iDipAcol.push_back(0);
485 }
486
487 // Pick final-state antitriplet (if any) to carry initial anticolour.
488 if (acol0 != 0 && iAtriplet.size() > 0) {
489 int pickA = (iAtriplet.size() == 1) ? 0
490 : int( TINY + Rndm::flat() * (iAtriplet.size() - TINY) );
491 int iPickA = iAtriplet[pickA];
492 acols[iPickA] = acol0;
493
494 // Remove matched antitriplet and store new colour dipole ends.
495 acol0 = 0;
496 iAtriplet[pickA] = iAtriplet.back();
497 iAtriplet.pop_back();
498 iDipCol.push_back(0);
499 iDipAcol.push_back(iPickA);
500 }
501
502 // Error checks that amount of leftover colours and anticolours match.
503 if ( (iTriplet.size() != iAtriplet.size())
504 || (col0 != 0 && acol0 == 0) || (col0 == 0 && acol0 != 0) ) {
505 infoPtr->errorMsg("Error in ResonanceDecays::pickColours:"
506 " inconsistent colour tags");
507 return false;
508 }
509
510 // Match triplets to antitriplets in the final state.
511 for (int pickT = 0; pickT < int(iTriplet.size()); ++pickT) {
512 int iPickT = iTriplet[pickT];
513 int pickA = (iAtriplet.size() == 1) ? 0
514 : int( TINY + Rndm::flat() * (iAtriplet.size() - TINY) );
515 int iPickA = iAtriplet[pickA];
516
517 // Connect pair with new colour tag.
518 cols[iPickT] = process.nextColTag();
519 acols[iPickA] = cols[iPickT];
520
521 // Remove matched antitriplet and store new colour dipole ends.
522 iAtriplet[pickT] = iAtriplet.back();
523 iAtriplet.pop_back();
524 iDipCol.push_back(iPickT);
525 iDipAcol.push_back(iPickA);
526 }
527
528 // If no octets are around then matching is done.
529 if (col0 == 0 && acol0 == 0 && iOctet.size() == 0) return true;
530
531 // If initial-state octet remains then store as (first!) new dipole.
532 if (col0 != 0) {
533 iDipCol.push_back(0);
534 iDipAcol.push_back(0);
535 }
536
537 // Now attach all final-state octets at random to existing dipoles.
538 for (int i = 0; i < int(iOctet.size()); ++i) {
539 int iOct = iOctet[i];
540
541 // If no dipole then start new one. (Happens for singlet -> octets.)
542 if (iDipCol.size() == 0) {
543 cols[iOct] = process.nextColTag();
544 acols[iOct] = cols[iOct] ;
545 iDipCol.push_back(iOct);
546 iDipAcol.push_back(iOct);
547 }
548
549 // Else attach to existing dipole picked at random.
550 else {
551 int pickDip = (iDipCol.size() == 1) ? 0
552 : int( TINY + Rndm::flat() * (iDipCol.size() - TINY) );
553
554 // Case with dipole in initial state: reattach existing colours.
555 if (iDipCol[pickDip] == 0 && iDipAcol[pickDip] == 0) {
556 cols[iOct] = col0;
557 acols[iOct] = acol0;
558 iDipAcol[pickDip] = iOct;
559 iDipCol.push_back(iOct);
560 iDipAcol.push_back(0);
561
562 // Case with dipole from colour in initial state: also new colour.
563 } else if (iDipAcol[pickDip] == 0) {
564 int iPickCol = iDipCol[pickDip];
565 cols[iOct] = cols[iPickCol];
566 acols[iOct] = process.nextColTag();
567 cols[iPickCol] = acols[iOct];
568 iDipCol[pickDip] = iOct;
569 iDipCol.push_back(iPickCol);
570 iDipAcol.push_back(iOct);
571
572 // Remaining cases with dipole from anticolour in initial state
573 // or dipole inside final state: also new colour.
574 } else {
575 int iPickAcol = iDipAcol[pickDip];
576 acols[iOct] = acols[iPickAcol];
577 cols[iOct] = process.nextColTag();
578 acols[iPickAcol] = cols[iOct];
579 iDipAcol[pickDip] = iOct;
580 iDipCol.push_back(iOct);
581 iDipAcol.push_back(iPickAcol);
582 }
583 }
584 }
585
586 // Must now have at least two dipoles (no 1 -> 8 or 8 -> 1).
587 if (iDipCol.size() < 2) {
588 infoPtr->errorMsg("Error in ResonanceDecays::pickColours:"
589 " inconsistent colour tags");
590 return false;
591 }
592
593 // Done.
594 return true;
595
596}
597
598//*********
599
600// Select decay products momenta isotropically in phase space.
601// Process-dependent angular distributions may be imposed in SigmaProcess.
602
603bool ResonanceDecays::pickKinematics() {
604
605 // Description of two-body decays as simple special case.
606 if (mult == 2) {
607
608 // Masses.
609 m0 = mProd[0];
610 double m1 = mProd[1];
611 double m2 = mProd[2];
612
613 // Energies and absolute momentum in the rest frame.
614 double e1 = 0.5 * (m0*m0 + m1*m1 - m2*m2) / m0;
615 double e2 = 0.5 * (m0*m0 + m2*m2 - m1*m1) / m0;
616 double pAbs = 0.5 * sqrtpos( (m0 - m1 - m2) * (m0 + m1 + m2)
617 * (m0 + m1 - m2) * (m0 - m1 + m2) ) / m0;
618
619 // Pick isotropic angles to give three-momentum.
620 double cosTheta = 2. * Rndm::flat() - 1.;
621 double sinTheta = sqrt(1. - cosTheta*cosTheta);
622 double phi = 2. * M_PI * Rndm::flat();
623 double pX = pAbs * sinTheta * cos(phi);
624 double pY = pAbs * sinTheta * sin(phi);
625 double pZ = pAbs * cosTheta;
626
627 // Fill four-momenta in mother rest frame and then boost to lab frame.
628 pProd.push_back( Vec4( pX, pY, pZ, e1) );
629 pProd.push_back( Vec4( -pX, -pY, -pZ, e2) );
630 pProd[1].bst( pProd[0] );
631 pProd[2].bst( pProd[0] );
632
633 // Done for two-body decay.
634 return true;
635 }
636
637 // Description of three-body decays as semi-simple special case.
638 if (mult == 3) {
639
640 // Masses.
641 m0 = mProd[0];
642 double m1 = mProd[1];
643 double m2 = mProd[2];
644 double m3 = mProd[3];
645 double mDiff = m0 - (m1 + m2 + m3);
646
647 // Kinematical limits for 2+3 mass. Maximum phase-space weight.
648 double m23Min = m2 + m3;
649 double m23Max = m0 - m1;
650 double p1Max = 0.5 * sqrtpos( (m0 - m1 - m23Min) * (m0 + m1 + m23Min)
651 * (m0 + m1 - m23Min) * (m0 - m1 + m23Min) ) / m0;
652 double p23Max = 0.5 * sqrtpos( (m23Max - m2 - m3) * (m23Max + m2 + m3)
653 * (m23Max + m2 - m3) * (m23Max - m2 + m3) ) / m23Max;
654 double wtPSmax = 0.5 * p1Max * p23Max;
655
656 // Pick an intermediate mass m23 flat in the allowed range.
657 double wtPS, m23, p1Abs, p23Abs;
658 do {
659 m23 = m23Min + Rndm::flat() * mDiff;
660
661 // Translate into relative momenta and find phase-space weight.
662 p1Abs = 0.5 * sqrtpos( (m0 - m1 - m23) * (m0 + m1 + m23)
663 * (m0 + m1 - m23) * (m0 - m1 + m23) ) / m0;
664 p23Abs = 0.5 * sqrtpos( (m23 - m2 - m3) * (m23 + m2 + m3)
665 * (m23 + m2 - m3) * (m23 - m2 + m3) ) / m23;
666 wtPS = p1Abs * p23Abs;
667
668 // If rejected, try again with new invariant masses.
669 } while ( wtPS < Rndm::flat() * wtPSmax );
670
671 // Set up m23 -> m2 + m3 isotropic in its rest frame.
672 double cosTheta = 2. * Rndm::flat() - 1.;
673 double sinTheta = sqrt(1. - cosTheta*cosTheta);
674 double phi = 2. * M_PI * Rndm::flat();
675 double pX = p23Abs * sinTheta * cos(phi);
676 double pY = p23Abs * sinTheta * sin(phi);
677 double pZ = p23Abs * cosTheta;
678 double e2 = sqrt( m2*m2 + p23Abs*p23Abs);
679 double e3 = sqrt( m3*m3 + p23Abs*p23Abs);
680 Vec4 p2( pX, pY, pZ, e2);
681 Vec4 p3( -pX, -pY, -pZ, e3);
682
683 // Set up 0 -> 1 + 23 isotropic in its rest frame.
684 cosTheta = 2. * Rndm::flat() - 1.;
685 sinTheta = sqrt(1. - cosTheta*cosTheta);
686 phi = 2. * M_PI * Rndm::flat();
687 pX = p1Abs * sinTheta * cos(phi);
688 pY = p1Abs * sinTheta * sin(phi);
689 pZ = p1Abs * cosTheta;
690 double e1 = sqrt( m1*m1 + p1Abs*p1Abs);
691 double e23 = sqrt( m23*m23 + p1Abs*p1Abs);
692 pProd.push_back( Vec4( pX, pY, pZ, e1) );
693
694 // Boost 2 + 3 to the 0 rest frame and then boost to lab frame.
695 Vec4 p23( -pX, -pY, -pZ, e23);
696 p2.bst( p23 );
697 p3.bst( p23 );
698 pProd.push_back( p2 );
699 pProd.push_back( p3 );
700 pProd[1].bst( pProd[0] );
701 pProd[2].bst( pProd[0] );
702 pProd[3].bst( pProd[0] );
703
704 // Done for three-body decay.
705 return true;
706 }
707
708 // Do a multibody decay using the M-generator algorithm.
709
710 // Mother and sum daughter masses.
711 m0 = mProd[0];
712 double mSum = mProd[1];
713 for (int i = 2; i <= mult; ++i) mSum += mProd[i];
714 double mDiff = m0 - mSum;
715
716 // Begin setup of intermediate invariant masses.
717 vector<double> mInv;
718 for (int i = 0; i <= mult; ++i) mInv.push_back( mProd[i]);
719
720 // Calculate the maximum weight in the decay.
721 double wtPSmax = 1. / WTCORRECTION[mult];
722 double mMax = mDiff + mProd[mult];
723 double mMin = 0.;
724 for (int i = mult - 1; i > 0; --i) {
725 mMax += mProd[i];
726 mMin += mProd[i+1];
727 double mNow = mProd[i];
728 wtPSmax *= 0.5 * sqrtpos( (mMax - mMin - mNow) * (mMax + mMin + mNow)
729 * (mMax + mMin - mNow) * (mMax - mMin + mNow) ) / mMax;
730 }
731
732 // Begin loop to find the set of intermediate invariant masses.
733 vector<double> rndmOrd;
734 double wtPS;
735 do {
736 wtPS = 1.;
737
738 // Find and order random numbers in descending order.
739 rndmOrd.resize(0);
740 rndmOrd.push_back(1.);
741 for (int i = 1; i < mult - 1; ++i) {
742 double rndm = Rndm::flat();
743 rndmOrd.push_back(rndm);
744 for (int j = i - 1; j > 0; --j) {
745 if (rndm > rndmOrd[j]) swap( rndmOrd[j], rndmOrd[j+1] );
746 else break;
747 }
748 }
749 rndmOrd.push_back(0.);
750
751 // Translate into intermediate masses and find weight.
752 for (int i = mult - 1; i > 0; --i) {
753 mInv[i] = mInv[i+1] + mProd[i] + (rndmOrd[i-1] - rndmOrd[i]) * mDiff;
754 wtPS *= 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
755 * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
756 * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
757 }
758
759 // If rejected, try again with new invariant masses.
760 } while ( wtPS < Rndm::flat() * wtPSmax );
761
762 // Perform two-particle decays in the respective rest frame.
763 vector<Vec4> pInv;
764 pInv.resize(mult + 1);
765 for (int i = 1; i < mult; ++i) {
766 double pAbs = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i])
767 * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i])
768 * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i];
769
770 // Isotropic angles give three-momentum.
771 double cosTheta = 2. * Rndm::flat() - 1.;
772 double sinTheta = sqrt(1. - cosTheta*cosTheta);
773 double phi = 2. * M_PI * Rndm::flat();
774 double pX = pAbs * sinTheta * cos(phi);
775 double pY = pAbs * sinTheta * sin(phi);
776 double pZ = pAbs * cosTheta;
777
778 // Calculate energies, fill four-momenta.
779 double eHad = sqrt( mProd[i]*mProd[i] + pAbs*pAbs);
780 double eInv = sqrt( mInv[i+1]*mInv[i+1] + pAbs*pAbs);
781 pProd.push_back( Vec4( pX, pY, pZ, eHad) );
782 pInv[i+1].p( -pX, -pY, -pZ, eInv);
783 }
784 pProd.push_back( pInv[mult] );
785
786 // Boost decay products to the mother rest frame and on to lab frame.
787 pInv[1] = pProd[0];
788 for (int iFrame = mult - 1; iFrame > 0; --iFrame)
789 for (int i = iFrame; i <= mult; ++i) pProd[i].bst(pInv[iFrame]);
790
791 // Done for multibody decay.
792 return true;
793
794}
795
796//**************************************************************************
797
798} // end namespace Pythia8