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