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