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