1 | // ParticleDecays.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 the |
---|
7 | // ParticleDecays class. |
---|
8 | |
---|
9 | #include "ParticleDecays.h" |
---|
10 | |
---|
11 | namespace Pythia8 { |
---|
12 | |
---|
13 | //========================================================================== |
---|
14 | |
---|
15 | // The ParticleDecays class. |
---|
16 | |
---|
17 | //-------------------------------------------------------------------------- |
---|
18 | |
---|
19 | // Constants: could be changed here if desired, but normally should not. |
---|
20 | // These are of technical nature, as described for each. |
---|
21 | |
---|
22 | // Number of times one tries to let decay happen (for 2 nested loops). |
---|
23 | const int ParticleDecays::NTRYDECAY = 10; |
---|
24 | |
---|
25 | // Number of times one tries to pick valid hadronic content in decay. |
---|
26 | const int ParticleDecays::NTRYPICK = 100; |
---|
27 | |
---|
28 | // Number of times one tries to pick decay topology. |
---|
29 | const int ParticleDecays::NTRYMEWT = 1000; |
---|
30 | |
---|
31 | // Maximal loop count in Dalitz decay treatment. |
---|
32 | const int ParticleDecays::NTRYDALITZ = 1000; |
---|
33 | |
---|
34 | // Minimal Dalitz pair mass this factor above threshold. |
---|
35 | const double ParticleDecays::MSAFEDALITZ = 1.000001; |
---|
36 | |
---|
37 | // These numbers are hardwired empirical parameters, |
---|
38 | // intended to speed up the M-generator. |
---|
39 | const double ParticleDecays::WTCORRECTION[11] = { 1., 1., 1., |
---|
40 | 2., 5., 15., 60., 250., 1250., 7000., 50000. }; |
---|
41 | |
---|
42 | //-------------------------------------------------------------------------- |
---|
43 | |
---|
44 | // Initialize and save pointers. |
---|
45 | |
---|
46 | void ParticleDecays::init(Info* infoPtrIn, Settings& settings, |
---|
47 | ParticleData* particleDataPtrIn, Rndm* rndmPtrIn, |
---|
48 | Couplings* couplingsPtrIn, TimeShower* timesDecPtrIn, |
---|
49 | StringFlav* flavSelPtrIn, DecayHandler* decayHandlePtrIn, |
---|
50 | vector<int> handledParticles) { |
---|
51 | |
---|
52 | // Save pointers to error messages handling and flavour generation. |
---|
53 | infoPtr = infoPtrIn; |
---|
54 | particleDataPtr = particleDataPtrIn; |
---|
55 | rndmPtr = rndmPtrIn; |
---|
56 | couplingsPtr = couplingsPtrIn; |
---|
57 | flavSelPtr = flavSelPtrIn; |
---|
58 | |
---|
59 | // Save pointer to timelike shower, as needed in some few decays. |
---|
60 | timesDecPtr = timesDecPtrIn; |
---|
61 | |
---|
62 | // Save pointer for external handling of some decays. |
---|
63 | decayHandlePtr = decayHandlePtrIn; |
---|
64 | |
---|
65 | // Set which particles should be handled externally. |
---|
66 | if (decayHandlePtr != 0) |
---|
67 | for (int i = 0; i < int(handledParticles.size()); ++i) |
---|
68 | particleDataPtr->doExternalDecay(handledParticles[i], true); |
---|
69 | |
---|
70 | // Safety margin in mass to avoid troubles. |
---|
71 | mSafety = settings.parm("ParticleDecays:mSafety"); |
---|
72 | |
---|
73 | // Lifetime and vertex rules for determining whether decay allowed. |
---|
74 | limitTau0 = settings.flag("ParticleDecays:limitTau0"); |
---|
75 | tau0Max = settings.parm("ParticleDecays:tau0Max"); |
---|
76 | limitTau = settings.flag("ParticleDecays:limitTau"); |
---|
77 | tauMax = settings.parm("ParticleDecays:tauMax"); |
---|
78 | limitRadius = settings.flag("ParticleDecays:limitRadius"); |
---|
79 | rMax = settings.parm("ParticleDecays:rMax"); |
---|
80 | limitCylinder = settings.flag("ParticleDecays:limitCylinder"); |
---|
81 | xyMax = settings.parm("ParticleDecays:xyMax"); |
---|
82 | zMax = settings.parm("ParticleDecays:zMax"); |
---|
83 | limitDecay = limitTau0 || limitTau || limitRadius || limitCylinder; |
---|
84 | |
---|
85 | // B-Bbar mixing parameters. |
---|
86 | mixB = settings.flag("ParticleDecays:mixB"); |
---|
87 | xBdMix = settings.parm("ParticleDecays:xBdMix"); |
---|
88 | xBsMix = settings.parm("ParticleDecays:xBsMix"); |
---|
89 | |
---|
90 | // Suppression of extra-hadron momenta in semileptonic decays. |
---|
91 | sigmaSoft = settings.parm("ParticleDecays:sigmaSoft"); |
---|
92 | |
---|
93 | // Selection of multiplicity and colours in "phase space" model. |
---|
94 | multIncrease = settings.parm("ParticleDecays:multIncrease"); |
---|
95 | multIncreaseWeak = settings.parm("ParticleDecays:multIncreaseWeak"); |
---|
96 | multRefMass = settings.parm("ParticleDecays:multRefMass"); |
---|
97 | multGoffset = settings.parm("ParticleDecays:multGoffset"); |
---|
98 | colRearrange = settings.parm("ParticleDecays:colRearrange"); |
---|
99 | |
---|
100 | // Minimum energy in system (+ m_q) from StringFragmentation. |
---|
101 | stopMass = settings.parm("StringFragmentation:stopMass"); |
---|
102 | |
---|
103 | // Parameters for Dalitz decay virtual gamma mass spectrum. |
---|
104 | sRhoDal = pow2(particleDataPtr->m0(113)); |
---|
105 | wRhoDal = pow2(particleDataPtr->mWidth(113)); |
---|
106 | |
---|
107 | // Allow showers in decays to qqbar/gg/ggg/gammagg. |
---|
108 | doFSRinDecays = settings.flag("ParticleDecays:FSRinDecays"); |
---|
109 | |
---|
110 | // Use standard decays or dedicated tau decay package |
---|
111 | sophisticatedTau = settings.mode("ParticleDecays:sophisticatedTau"); |
---|
112 | |
---|
113 | // Initialize the dedicated tau decay handler. |
---|
114 | if (sophisticatedTau) tauDecayer.init(infoPtr, &settings, |
---|
115 | particleDataPtr, rndmPtr, couplingsPtr); |
---|
116 | |
---|
117 | } |
---|
118 | |
---|
119 | //-------------------------------------------------------------------------- |
---|
120 | |
---|
121 | // Decay a particle; main method. |
---|
122 | |
---|
123 | bool ParticleDecays::decay( int iDec, Event& event) { |
---|
124 | |
---|
125 | // Check whether a decay is allowed, given the upcoming decay vertex. |
---|
126 | Particle& decayer = event[iDec]; |
---|
127 | hasPartons = false; |
---|
128 | keepPartons = false; |
---|
129 | if (limitDecay && !checkVertex(decayer)) return true; |
---|
130 | |
---|
131 | // Do not allow resonance decays (beyond handling capability). |
---|
132 | if (decayer.isResonance()) { |
---|
133 | infoPtr->errorMsg("Warning in ParticleDecays::decay: " |
---|
134 | "resonance left undecayed"); |
---|
135 | return true; |
---|
136 | } |
---|
137 | |
---|
138 | // Fill the decaying particle in slot 0 of arrays. |
---|
139 | idDec = decayer.id(); |
---|
140 | iProd.resize(0); |
---|
141 | idProd.resize(0); |
---|
142 | mProd.resize(0); |
---|
143 | iProd.push_back( iDec ); |
---|
144 | idProd.push_back( idDec ); |
---|
145 | mProd.push_back( decayer.m() ); |
---|
146 | |
---|
147 | // Check for oscillations B0 <-> B0bar or B_s0 <-> B_s0bar. |
---|
148 | bool hasOscillated = (abs(idDec) == 511 || abs(idDec) == 531) |
---|
149 | ? oscillateB(decayer) : false; |
---|
150 | if (hasOscillated) {idDec = - idDec; idProd[0] = idDec;} |
---|
151 | |
---|
152 | // Particle data for decaying particle. |
---|
153 | decDataPtr = &decayer.particleDataEntry(); |
---|
154 | |
---|
155 | // Optionally send on to external decay program. |
---|
156 | bool doneExternally = false; |
---|
157 | if (decDataPtr->doExternalDecay()) { |
---|
158 | pProd.resize(0); |
---|
159 | pProd.push_back(decayer.p()); |
---|
160 | doneExternally = decayHandlePtr->decay(idProd, mProd, pProd, |
---|
161 | iDec, event); |
---|
162 | |
---|
163 | // If it worked, then store the decay products in the event record. |
---|
164 | if (doneExternally) { |
---|
165 | mult = idProd.size() - 1; |
---|
166 | int status = (hasOscillated) ? 94 : 93; |
---|
167 | for (int i = 1; i <= mult; ++i) { |
---|
168 | int iPos = event.append( idProd[i], status, iDec, 0, 0, 0, |
---|
169 | 0, 0, pProd[i], mProd[i]); |
---|
170 | iProd.push_back( iPos); |
---|
171 | } |
---|
172 | |
---|
173 | // Also mark mother decayed and store daughters. |
---|
174 | event[iDec].statusNeg(); |
---|
175 | event[iDec].daughters( iProd[1], iProd[mult]); |
---|
176 | } |
---|
177 | } |
---|
178 | |
---|
179 | // Check if the particle is tau and let the special tau decayer handle it. |
---|
180 | if (decayer.idAbs() == 15 && !doneExternally && sophisticatedTau) { |
---|
181 | doneExternally = tauDecayer.decay(iDec, event); |
---|
182 | if (doneExternally) return true; |
---|
183 | } |
---|
184 | |
---|
185 | // Now begin normal internal decay treatment. |
---|
186 | if (!doneExternally) { |
---|
187 | |
---|
188 | // Allow up to ten tries to pick a channel. |
---|
189 | if (!decDataPtr->preparePick(idDec, decayer.m())) return false; |
---|
190 | bool foundChannel = false; |
---|
191 | bool hasStored = false; |
---|
192 | for (int iTryChannel = 0; iTryChannel < NTRYDECAY; ++iTryChannel) { |
---|
193 | |
---|
194 | // Remove previous failed channel. |
---|
195 | if (hasStored) event.popBack(mult); |
---|
196 | hasStored = false; |
---|
197 | |
---|
198 | // Pick new channel. Read out basics. |
---|
199 | DecayChannel& channel = decDataPtr->pickChannel(); |
---|
200 | meMode = channel.meMode(); |
---|
201 | keepPartons = (meMode > 90 && meMode <= 100); |
---|
202 | mult = channel.multiplicity(); |
---|
203 | |
---|
204 | // Allow up to ten tries for each channel (e.g with different masses). |
---|
205 | bool foundMode = false; |
---|
206 | for (int iTryMode = 0; iTryMode < NTRYDECAY; ++iTryMode) { |
---|
207 | idProd.resize(1); |
---|
208 | mProd.resize(1); |
---|
209 | scale = 0.; |
---|
210 | |
---|
211 | // Extract and store the decay products in local arrays. |
---|
212 | hasPartons = false; |
---|
213 | for (int i = 0; i < mult; ++i) { |
---|
214 | int idNow = channel.product(i); |
---|
215 | int idAbs = abs(idNow); |
---|
216 | if ( idAbs < 10 || idAbs == 21 || idAbs == 81 || idAbs == 82 |
---|
217 | || idAbs == 83 || (idAbs > 1000 && idAbs < 10000 |
---|
218 | && (idAbs/10)%10 == 0) ) hasPartons = true; |
---|
219 | if (idDec < 0 && particleDataPtr->hasAnti(idNow)) idNow = -idNow; |
---|
220 | double mNow = particleDataPtr->mass(idNow); |
---|
221 | idProd.push_back( idNow); |
---|
222 | mProd.push_back( mNow); |
---|
223 | } |
---|
224 | |
---|
225 | // Decays into partons usually translate into hadrons. |
---|
226 | if (hasPartons && !keepPartons && !pickHadrons()) continue; |
---|
227 | |
---|
228 | // Need to set colour flow if explicit decay to partons. |
---|
229 | cols.resize(0); |
---|
230 | acols.resize(0); |
---|
231 | for (int i = 0; i <= mult; ++i) { |
---|
232 | cols.push_back(0); |
---|
233 | acols.push_back(0); |
---|
234 | } |
---|
235 | if (hasPartons && keepPartons && !setColours(event)) continue; |
---|
236 | |
---|
237 | // Check that enough phase space for decay. |
---|
238 | if (mult > 1) { |
---|
239 | double mDiff = mProd[0]; |
---|
240 | for (int i = 1; i <= mult; ++i) mDiff -= mProd[i]; |
---|
241 | if (mDiff < mSafety) continue; |
---|
242 | } |
---|
243 | |
---|
244 | // End of inner trial loops. Check if succeeded or not. |
---|
245 | foundMode = true; |
---|
246 | break; |
---|
247 | } |
---|
248 | if (!foundMode) continue; |
---|
249 | |
---|
250 | // Store decay products in the event record. |
---|
251 | int status = (hasOscillated) ? 92 : 91; |
---|
252 | for (int i = 1; i <= mult; ++i) { |
---|
253 | int iPos = event.append( idProd[i], status, iDec, 0, 0, 0, |
---|
254 | cols[i], acols[i], Vec4(0., 0., 0., 0.), mProd[i], scale); |
---|
255 | iProd.push_back( iPos); |
---|
256 | } |
---|
257 | hasStored = true; |
---|
258 | |
---|
259 | // Pick mass of Dalitz decay. Temporarily change multiplicity. |
---|
260 | if ( (meMode == 11 || meMode == 12 || meMode == 13) |
---|
261 | && !dalitzMass() ) continue; |
---|
262 | |
---|
263 | // Do a decay, split by multiplicity. |
---|
264 | bool decayed = false; |
---|
265 | if (mult == 1) decayed = oneBody(event); |
---|
266 | else if (mult == 2) decayed = twoBody(event); |
---|
267 | else if (mult == 3) decayed = threeBody(event); |
---|
268 | else decayed = mGenerator(event); |
---|
269 | if (!decayed) continue; |
---|
270 | |
---|
271 | // Kinematics of gamma* -> l- l+ in Dalitz decay. Restore multiplicity. |
---|
272 | if (meMode == 11 || meMode == 12 || meMode == 13) |
---|
273 | dalitzKinematics(event); |
---|
274 | |
---|
275 | // End of outer trial loops. |
---|
276 | foundChannel = true; |
---|
277 | break; |
---|
278 | } |
---|
279 | |
---|
280 | // If the decay worked, then mark mother decayed and store daughters. |
---|
281 | if (foundChannel) { |
---|
282 | event[iDec].statusNeg(); |
---|
283 | event[iDec].daughters( iProd[1], iProd[mult]); |
---|
284 | |
---|
285 | // Else remove unused daughters and return failure. |
---|
286 | } else { |
---|
287 | if (hasStored) event.popBack(mult); |
---|
288 | infoPtr->errorMsg("Error in ParticleDecays::decay: " |
---|
289 | "failed to find workable decay channel"); |
---|
290 | return false; |
---|
291 | } |
---|
292 | |
---|
293 | // Now finished normal internal decay treatment. |
---|
294 | } |
---|
295 | |
---|
296 | // Set decay vertex when this is displaced. |
---|
297 | if (event[iDec].hasVertex() || event[iDec].tau() > 0.) { |
---|
298 | Vec4 vDec = event[iDec].vDec(); |
---|
299 | for (int i = 1; i <= mult; ++i) event[iProd[i]].vProd( vDec ); |
---|
300 | } |
---|
301 | |
---|
302 | // Set lifetime of daughters. |
---|
303 | for (int i = 1; i <= mult; ++i) |
---|
304 | event[iProd[i]].tau( event[iProd[i]].tau0() * rndmPtr->exp() ); |
---|
305 | |
---|
306 | // In a decay explicitly to partons then optionally do a shower, |
---|
307 | // and always flag that partonic system should be fragmented. |
---|
308 | if (hasPartons && keepPartons && doFSRinDecays) |
---|
309 | timesDecPtr->shower( iProd[1], iProd.back(), event, mProd[0]); |
---|
310 | |
---|
311 | // For Hidden Valley particles also allow leptons to shower. |
---|
312 | else if (event[iDec].idAbs() > 4900000 && event[iDec].idAbs() < 5000000 |
---|
313 | && doFSRinDecays && mult == 2 && event[iProd[1]].isLepton()) { |
---|
314 | event[iProd[1]].scale(mProd[0]); |
---|
315 | event[iProd[2]].scale(mProd[0]); |
---|
316 | timesDecPtr->shower( iProd[1], iProd.back(), event, mProd[0]); |
---|
317 | } |
---|
318 | |
---|
319 | // Done. |
---|
320 | return true; |
---|
321 | |
---|
322 | } |
---|
323 | |
---|
324 | //-------------------------------------------------------------------------- |
---|
325 | |
---|
326 | // Check whether a decay is allowed, given the upcoming decay vertex. |
---|
327 | |
---|
328 | bool ParticleDecays::checkVertex(Particle& decayer) { |
---|
329 | |
---|
330 | // Check whether any of the conditions are not fulfilled. |
---|
331 | if (limitTau0 && decayer.tau0() > tau0Max) return false; |
---|
332 | if (limitTau && decayer.tau() > tauMax) return false; |
---|
333 | if (limitRadius && pow2(decayer.xDec()) + pow2(decayer.yDec()) |
---|
334 | + pow2(decayer.zDec()) > pow2(rMax)) return false; |
---|
335 | if (limitCylinder && (pow2(decayer.xDec()) + pow2(decayer.yDec()) |
---|
336 | > pow2(xyMax) || abs(decayer.zDec()) > zMax) ) return false; |
---|
337 | |
---|
338 | // Done. |
---|
339 | return true; |
---|
340 | |
---|
341 | } |
---|
342 | |
---|
343 | //-------------------------------------------------------------------------- |
---|
344 | |
---|
345 | // Check for oscillations B0 <-> B0bar or B_s0 <-> B_s0bar. |
---|
346 | |
---|
347 | bool ParticleDecays::oscillateB(Particle& decayer) { |
---|
348 | |
---|
349 | // Extract relevant information and decide. |
---|
350 | if (!mixB) return false; |
---|
351 | double xBmix = (abs(decayer.id()) == 511) ? xBdMix : xBsMix; |
---|
352 | double tau = decayer.tau(); |
---|
353 | double tau0 = decayer.tau0(); |
---|
354 | double probosc = pow2(sin(0.5 * xBmix * tau / tau0)); |
---|
355 | return (probosc > rndmPtr->flat()); |
---|
356 | |
---|
357 | } |
---|
358 | |
---|
359 | //-------------------------------------------------------------------------- |
---|
360 | |
---|
361 | // Do a one-body decay. (Rare; e.g. for K0 -> K0_short.) |
---|
362 | |
---|
363 | bool ParticleDecays::oneBody(Event& event) { |
---|
364 | |
---|
365 | // References to the particles involved. |
---|
366 | Particle& decayer = event[iProd[0]]; |
---|
367 | Particle& prod = event[iProd[1]]; |
---|
368 | |
---|
369 | // Set momentum and expand mother information. |
---|
370 | prod.p( decayer.p() ); |
---|
371 | prod.m( decayer.m() ); |
---|
372 | prod.mother2( iProd[0] ); |
---|
373 | |
---|
374 | // Done. |
---|
375 | return true; |
---|
376 | |
---|
377 | } |
---|
378 | |
---|
379 | //-------------------------------------------------------------------------- |
---|
380 | |
---|
381 | // Do a two-body decay. |
---|
382 | |
---|
383 | bool ParticleDecays::twoBody(Event& event) { |
---|
384 | |
---|
385 | // References to the particles involved. |
---|
386 | Particle& decayer = event[iProd[0]]; |
---|
387 | Particle& prod1 = event[iProd[1]]; |
---|
388 | Particle& prod2 = event[iProd[2]]; |
---|
389 | |
---|
390 | // Masses. |
---|
391 | double m0 = mProd[0]; |
---|
392 | double m1 = mProd[1]; |
---|
393 | double m2 = mProd[2]; |
---|
394 | |
---|
395 | // Energies and absolute momentum in the rest frame. |
---|
396 | if (m1 + m2 + mSafety > m0) return false; |
---|
397 | double e1 = 0.5 * (m0*m0 + m1*m1 - m2*m2) / m0; |
---|
398 | double e2 = 0.5 * (m0*m0 + m2*m2 - m1*m1) / m0; |
---|
399 | double pAbs = 0.5 * sqrtpos( (m0 - m1 - m2) * (m0 + m1 + m2) |
---|
400 | * (m0 + m1 - m2) * (m0 - m1 + m2) ) / m0; |
---|
401 | |
---|
402 | // When meMode = 2, for V -> PS2 + PS3 (V = vector, pseudoscalar), |
---|
403 | // need to check if production is PS0 -> PS1/gamma + V. |
---|
404 | int iMother = event[iProd[0]].mother1(); |
---|
405 | int idSister = 0; |
---|
406 | if (meMode == 2) { |
---|
407 | if (iMother <= 0 || iMother >= iProd[0]) meMode = 0; |
---|
408 | else { |
---|
409 | int iDaughter1 = event[iMother].daughter1(); |
---|
410 | int iDaughter2 = event[iMother].daughter2(); |
---|
411 | if (iDaughter2 != iDaughter1 + 1) meMode = 0; |
---|
412 | else { |
---|
413 | int idMother = abs( event[iMother].id() ); |
---|
414 | if (idMother <= 100 || idMother%10 !=1 |
---|
415 | || (idMother/1000)%10 != 0) meMode = 0; |
---|
416 | else { |
---|
417 | int iSister = (iProd[0] == iDaughter1) ? iDaughter2 : iDaughter1; |
---|
418 | idSister = abs( event[iSister].id() ); |
---|
419 | if ( (idSister <= 100 || idSister%10 !=1 |
---|
420 | || (idSister/1000)%10 != 0) && idSister != 22) meMode = 0; |
---|
421 | } |
---|
422 | } |
---|
423 | } |
---|
424 | } |
---|
425 | |
---|
426 | // Begin loop over matrix-element corrections. |
---|
427 | double wtME, wtMEmax; |
---|
428 | int loop = 0; |
---|
429 | do { |
---|
430 | wtME = 1.; |
---|
431 | wtMEmax = 1.; |
---|
432 | ++loop; |
---|
433 | |
---|
434 | // Isotropic angles give three-momentum. |
---|
435 | double cosTheta = 2. * rndmPtr->flat() - 1.; |
---|
436 | double sinTheta = sqrt(1. - cosTheta*cosTheta); |
---|
437 | double phi = 2. * M_PI * rndmPtr->flat(); |
---|
438 | double pX = pAbs * sinTheta * cos(phi); |
---|
439 | double pY = pAbs * sinTheta * sin(phi); |
---|
440 | double pZ = pAbs * cosTheta; |
---|
441 | |
---|
442 | // Fill four-momenta and boost them away from mother rest frame. |
---|
443 | prod1.p( pX, pY, pZ, e1); |
---|
444 | prod2.p( -pX, -pY, -pZ, e2); |
---|
445 | prod1.bst( decayer.p(), decayer.m() ); |
---|
446 | prod2.bst( decayer.p(), decayer.m() ); |
---|
447 | |
---|
448 | // Matrix element for PS0 -> PS1 + V1 -> PS1 + PS2 + PS3 of form |
---|
449 | // cos**2(theta02) in V1 rest frame, and for PS0 -> gamma + V1 |
---|
450 | // -> gamma + PS2 + PS3 of form sin**2(theta02). |
---|
451 | if (meMode == 2) { |
---|
452 | double p10 = decayer.p() * event[iMother].p(); |
---|
453 | double p12 = decayer.p() * prod1.p(); |
---|
454 | double p02 = event[iMother].p() * prod1.p(); |
---|
455 | double s0 = pow2(event[iMother].m()); |
---|
456 | double s1 = pow2(decayer.m()); |
---|
457 | double s2 = pow2(prod1.m()); |
---|
458 | if (idSister != 22) wtME = pow2(p10 * p12 - s1 * p02); |
---|
459 | else wtME = s1 * (2. * p10 * p12 * p02 - s1 * p02*p02 |
---|
460 | - s0 * p12*p12 - s2 * p10*p10 + s1 * s0 * s2); |
---|
461 | wtME = max( wtME, 1e-6 * s1*s1 * s0 * s2); |
---|
462 | wtMEmax = (p10*p10 - s1 * s0) * (p12*p12 - s1 * s2); |
---|
463 | } |
---|
464 | |
---|
465 | // Break out of loop if no sensible ME weight. |
---|
466 | if(loop > NTRYMEWT) { |
---|
467 | infoPtr->errorMsg("ParticleDecays::twoBody: " |
---|
468 | "caught in infinite ME weight loop"); |
---|
469 | wtME = abs(wtMEmax); |
---|
470 | } |
---|
471 | |
---|
472 | // If rejected, try again with new invariant masses. |
---|
473 | } while ( wtME < rndmPtr->flat() * wtMEmax ); |
---|
474 | |
---|
475 | // Done. |
---|
476 | return true; |
---|
477 | |
---|
478 | } |
---|
479 | |
---|
480 | //-------------------------------------------------------------------------- |
---|
481 | |
---|
482 | // Do a three-body decay (except Dalitz decays). |
---|
483 | |
---|
484 | bool ParticleDecays::threeBody(Event& event) { |
---|
485 | |
---|
486 | // References to the particles involved. |
---|
487 | Particle& decayer = event[iProd[0]]; |
---|
488 | Particle& prod1 = event[iProd[1]]; |
---|
489 | Particle& prod2 = event[iProd[2]]; |
---|
490 | Particle& prod3 = event[iProd[3]]; |
---|
491 | |
---|
492 | // Mother and sum daughter masses. Fail if too close. |
---|
493 | double m0 = mProd[0]; |
---|
494 | double m1 = mProd[1]; |
---|
495 | double m2 = mProd[2]; |
---|
496 | double m3 = mProd[3]; |
---|
497 | double mSum = m1 + m2 + m3; |
---|
498 | double mDiff = m0 - mSum; |
---|
499 | if (mDiff < mSafety) return false; |
---|
500 | |
---|
501 | // Kinematical limits for 2+3 mass. Maximum phase-space weight. |
---|
502 | double m23Min = m2 + m3; |
---|
503 | double m23Max = m0 - m1; |
---|
504 | double p1Max = 0.5 * sqrtpos( (m0 - m1 - m23Min) * (m0 + m1 + m23Min) |
---|
505 | * (m0 + m1 - m23Min) * (m0 - m1 + m23Min) ) / m0; |
---|
506 | double p23Max = 0.5 * sqrtpos( (m23Max - m2 - m3) * (m23Max + m2 + m3) |
---|
507 | * (m23Max + m2 - m3) * (m23Max - m2 + m3) ) / m23Max; |
---|
508 | double wtPSmax = 0.5 * p1Max * p23Max; |
---|
509 | |
---|
510 | // Begin loop over matrix-element corrections. |
---|
511 | double wtME, wtMEmax, wtPS, m23, p1Abs, p23Abs; |
---|
512 | do { |
---|
513 | wtME = 1.; |
---|
514 | wtMEmax = 1.; |
---|
515 | |
---|
516 | // Pick an intermediate mass m23 flat in the allowed range. |
---|
517 | do { |
---|
518 | m23 = m23Min + rndmPtr->flat() * mDiff; |
---|
519 | |
---|
520 | // Translate into relative momenta and find phase-space weight. |
---|
521 | p1Abs = 0.5 * sqrtpos( (m0 - m1 - m23) * (m0 + m1 + m23) |
---|
522 | * (m0 + m1 - m23) * (m0 - m1 + m23) ) / m0; |
---|
523 | p23Abs = 0.5 * sqrtpos( (m23 - m2 - m3) * (m23 + m2 + m3) |
---|
524 | * (m23 + m2 - m3) * (m23 - m2 + m3) ) / m23; |
---|
525 | wtPS = p1Abs * p23Abs; |
---|
526 | |
---|
527 | // If rejected, try again with new invariant masses. |
---|
528 | } while ( wtPS < rndmPtr->flat() * wtPSmax ); |
---|
529 | |
---|
530 | // Set up m23 -> m2 + m3 isotropic in its rest frame. |
---|
531 | double cosTheta = 2. * rndmPtr->flat() - 1.; |
---|
532 | double sinTheta = sqrt(1. - cosTheta*cosTheta); |
---|
533 | double phi = 2. * M_PI * rndmPtr->flat(); |
---|
534 | double pX = p23Abs * sinTheta * cos(phi); |
---|
535 | double pY = p23Abs * sinTheta * sin(phi); |
---|
536 | double pZ = p23Abs * cosTheta; |
---|
537 | double e2 = sqrt( m2*m2 + p23Abs*p23Abs); |
---|
538 | double e3 = sqrt( m3*m3 + p23Abs*p23Abs); |
---|
539 | prod2.p( pX, pY, pZ, e2); |
---|
540 | prod3.p( -pX, -pY, -pZ, e3); |
---|
541 | |
---|
542 | // Set up m0 -> m1 + m23 isotropic in its rest frame. |
---|
543 | cosTheta = 2. * rndmPtr->flat() - 1.; |
---|
544 | sinTheta = sqrt(1. - cosTheta*cosTheta); |
---|
545 | phi = 2. * M_PI * rndmPtr->flat(); |
---|
546 | pX = p1Abs * sinTheta * cos(phi); |
---|
547 | pY = p1Abs * sinTheta * sin(phi); |
---|
548 | pZ = p1Abs * cosTheta; |
---|
549 | double e1 = sqrt( m1*m1 + p1Abs*p1Abs); |
---|
550 | double e23 = sqrt( m23*m23 + p1Abs*p1Abs); |
---|
551 | prod1.p( pX, pY, pZ, e1); |
---|
552 | |
---|
553 | // Boost 2 + 3 to the 0 rest frame. |
---|
554 | Vec4 p23( -pX, -pY, -pZ, e23); |
---|
555 | prod2.bst( p23, m23 ); |
---|
556 | prod3.bst( p23, m23 ); |
---|
557 | |
---|
558 | // Matrix-element weight for omega/phi -> pi+ pi- pi0. |
---|
559 | if (meMode == 1) { |
---|
560 | double p1p2 = prod1.p() * prod2.p(); |
---|
561 | double p1p3 = prod1.p() * prod3.p(); |
---|
562 | double p2p3 = prod2.p() * prod3.p(); |
---|
563 | wtME = pow2(m1 * m2 * m3) - pow2(m1 * p2p3) - pow2(m2 * p1p3) |
---|
564 | - pow2(m3 * p1p2) + 2. * p1p2 * p1p3 * p2p3; |
---|
565 | wtMEmax = pow3(m0 * m0) / 150.; |
---|
566 | |
---|
567 | // Effective matrix element for nu spectrum in tau -> nu + hadrons. |
---|
568 | } else if (meMode == 21) { |
---|
569 | double x1 = 2. * prod1.e() / m0; |
---|
570 | wtME = x1 * (3. - 2. * x1); |
---|
571 | double xMax = min( 0.75, 2. * (1. - mSum / m0) ); |
---|
572 | wtMEmax = xMax * (3. - 2. * xMax); |
---|
573 | |
---|
574 | // Matrix element for weak decay (only semileptonic for c and b). |
---|
575 | } else if ((meMode == 22 || meMode == 23) && prod1.isLepton()) { |
---|
576 | wtME = m0 * prod1.e() * (prod2.p() * prod3.p()); |
---|
577 | wtMEmax = min( pow4(m0) / 16., m0 * (m0 - m1 - m2) * (m0 - m1 - m3) |
---|
578 | * (m0 - m2 - m3) ); |
---|
579 | |
---|
580 | // Effective matrix element for weak decay to hadrons (B -> D, D -> K). |
---|
581 | } else if (meMode == 22 || meMode == 23) { |
---|
582 | double x1 = 2. * prod1.pAbs() / m0; |
---|
583 | wtME = x1 * (3. - 2. * x1); |
---|
584 | double xMax = min( 0.75, 2. * (1. - mSum / m0) ); |
---|
585 | wtMEmax = xMax * (3. - 2. * xMax); |
---|
586 | |
---|
587 | // Effective matrix element for gamma spectrum in B -> gamma + hadrons. |
---|
588 | } else if (meMode == 31) { |
---|
589 | double x1 = 2. * prod1.e() / m0; |
---|
590 | wtME = pow3(x1); |
---|
591 | double x1Max = 1. - pow2(mSum / m0); |
---|
592 | wtMEmax = pow3(x1Max); |
---|
593 | |
---|
594 | // Matrix-element weight for "onium" -> g + g + g or gamma + g + g. |
---|
595 | } else if (meMode == 92) { |
---|
596 | double x1 = 2. * prod1.e() / m0; |
---|
597 | double x2 = 2. * prod2.e() / m0; |
---|
598 | double x3 = 2. * prod3.e() / m0; |
---|
599 | wtME = pow2( (1. - x1) / (x2 * x3) ) + pow2( (1. - x2) / (x1 * x3) ) |
---|
600 | + pow2( (1. - x3) / (x1 * x2) ); |
---|
601 | wtMEmax = 2.; |
---|
602 | // For gamma + g + g require minimum mass for g + g system. |
---|
603 | if (prod1.id() == 22 && sqrt(1. - x1) * m0 < 2. * stopMass) wtME = 0.; |
---|
604 | if (prod2.id() == 22 && sqrt(1. - x2) * m0 < 2. * stopMass) wtME = 0.; |
---|
605 | if (prod3.id() == 22 && sqrt(1. - x3) * m0 < 2. * stopMass) wtME = 0.; |
---|
606 | } |
---|
607 | |
---|
608 | // If rejected, try again with new invariant masses. |
---|
609 | } while ( wtME < rndmPtr->flat() * wtMEmax ); |
---|
610 | |
---|
611 | // Boost 1 + 2 + 3 to the current frame. |
---|
612 | prod1.bst( decayer.p(), decayer.m() ); |
---|
613 | prod2.bst( decayer.p(), decayer.m() ); |
---|
614 | prod3.bst( decayer.p(), decayer.m() ); |
---|
615 | |
---|
616 | // Done. |
---|
617 | return true; |
---|
618 | |
---|
619 | } |
---|
620 | |
---|
621 | //-------------------------------------------------------------------------- |
---|
622 | |
---|
623 | // Do a multibody decay using the M-generator algorithm. |
---|
624 | |
---|
625 | bool ParticleDecays::mGenerator(Event& event) { |
---|
626 | |
---|
627 | // Mother and sum daughter masses. Fail if too close or inconsistent. |
---|
628 | double m0 = mProd[0]; |
---|
629 | double mSum = mProd[1]; |
---|
630 | for (int i = 2; i <= mult; ++i) mSum += mProd[i]; |
---|
631 | double mDiff = m0 - mSum; |
---|
632 | if (mDiff < mSafety) return false; |
---|
633 | |
---|
634 | // Begin setup of intermediate invariant masses. |
---|
635 | mInv.resize(0); |
---|
636 | for (int i = 0; i <= mult; ++i) mInv.push_back( mProd[i]); |
---|
637 | |
---|
638 | // Calculate the maximum weight in the decay. |
---|
639 | double wtPS, wtME, wtMEmax; |
---|
640 | double wtPSmax = 1. / WTCORRECTION[mult]; |
---|
641 | double mMax = mDiff + mProd[mult]; |
---|
642 | double mMin = 0.; |
---|
643 | for (int i = mult - 1; i > 0; --i) { |
---|
644 | mMax += mProd[i]; |
---|
645 | mMin += mProd[i+1]; |
---|
646 | double mNow = mProd[i]; |
---|
647 | wtPSmax *= 0.5 * sqrtpos( (mMax - mMin - mNow) * (mMax + mMin + mNow) |
---|
648 | * (mMax + mMin - mNow) * (mMax - mMin + mNow) ) / mMax; |
---|
649 | } |
---|
650 | |
---|
651 | // Begin loop over matrix-element corrections. |
---|
652 | do { |
---|
653 | wtME = 1.; |
---|
654 | wtMEmax = 1.; |
---|
655 | |
---|
656 | // Begin loop to find the set of intermediate invariant masses. |
---|
657 | do { |
---|
658 | wtPS = 1.; |
---|
659 | |
---|
660 | // Find and order random numbers in descending order. |
---|
661 | rndmOrd.resize(0); |
---|
662 | rndmOrd.push_back(1.); |
---|
663 | for (int i = 1; i < mult - 1; ++i) { |
---|
664 | double rndm = rndmPtr->flat(); |
---|
665 | rndmOrd.push_back(rndm); |
---|
666 | for (int j = i - 1; j > 0; --j) { |
---|
667 | if (rndm > rndmOrd[j]) swap( rndmOrd[j], rndmOrd[j+1] ); |
---|
668 | else break; |
---|
669 | } |
---|
670 | } |
---|
671 | rndmOrd.push_back(0.); |
---|
672 | |
---|
673 | // Translate into intermediate masses and find weight. |
---|
674 | for (int i = mult - 1; i > 0; --i) { |
---|
675 | mInv[i] = mInv[i+1] + mProd[i] + (rndmOrd[i-1] - rndmOrd[i]) * mDiff; |
---|
676 | wtPS *= 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i]) |
---|
677 | * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i]) |
---|
678 | * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i]; |
---|
679 | } |
---|
680 | |
---|
681 | // If rejected, try again with new invariant masses. |
---|
682 | } while ( wtPS < rndmPtr->flat() * wtPSmax ); |
---|
683 | |
---|
684 | // Perform two-particle decays in the respective rest frame. |
---|
685 | pInv.resize(mult + 1); |
---|
686 | for (int i = 1; i < mult; ++i) { |
---|
687 | double pAbs = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - mProd[i]) |
---|
688 | * (mInv[i] + mInv[i+1] + mProd[i]) * (mInv[i] + mInv[i+1] - mProd[i]) |
---|
689 | * (mInv[i] - mInv[i+1] + mProd[i]) ) / mInv[i]; |
---|
690 | |
---|
691 | // Isotropic angles give three-momentum. |
---|
692 | double cosTheta = 2. * rndmPtr->flat() - 1.; |
---|
693 | double sinTheta = sqrt(1. - cosTheta*cosTheta); |
---|
694 | double phi = 2. * M_PI * rndmPtr->flat(); |
---|
695 | double pX = pAbs * sinTheta * cos(phi); |
---|
696 | double pY = pAbs * sinTheta * sin(phi); |
---|
697 | double pZ = pAbs * cosTheta; |
---|
698 | |
---|
699 | // Calculate energies, fill four-momenta. |
---|
700 | double eHad = sqrt( mProd[i]*mProd[i] + pAbs*pAbs); |
---|
701 | double eInv = sqrt( mInv[i+1]*mInv[i+1] + pAbs*pAbs); |
---|
702 | event[iProd[i]].p( pX, pY, pZ, eHad); |
---|
703 | pInv[i+1].p( -pX, -pY, -pZ, eInv); |
---|
704 | } |
---|
705 | |
---|
706 | // Boost decay products to the mother rest frame. |
---|
707 | event[iProd[mult]].p( pInv[mult] ); |
---|
708 | for (int iFrame = mult - 1; iFrame > 1; --iFrame) |
---|
709 | for (int i = iFrame; i <= mult; ++i) |
---|
710 | event[iProd[i]].bst( pInv[iFrame], mInv[iFrame]); |
---|
711 | |
---|
712 | // Effective matrix element for nu spectrum in tau -> nu + hadrons. |
---|
713 | if (meMode == 21 && event[iProd[1]].isLepton()) { |
---|
714 | double x1 = 2. * event[iProd[1]].e() / m0; |
---|
715 | wtME = x1 * (3. - 2. * x1); |
---|
716 | double xMax = min( 0.75, 2. * (1. - mSum / m0) ); |
---|
717 | wtMEmax = xMax * (3. - 2. * xMax); |
---|
718 | |
---|
719 | // Effective matrix element for weak decay (only semileptonic for c and b). |
---|
720 | // Particles 4 onwards should be made softer explicitly? |
---|
721 | } else if ((meMode == 22 || meMode == 23) && event[iProd[1]].isLepton()) { |
---|
722 | Vec4 pRest = event[iProd[3]].p(); |
---|
723 | for (int i = 4; i <= mult; ++i) pRest += event[iProd[i]].p(); |
---|
724 | wtME = m0 * event[iProd[1]].e() * (event[iProd[2]].p() * pRest); |
---|
725 | for (int i = 4; i <= mult; ++i) wtME |
---|
726 | *= exp(- event[iProd[i]].pAbs2() / pow2(sigmaSoft) ); |
---|
727 | wtMEmax = pow4(m0) / 16.; |
---|
728 | |
---|
729 | // Effective matrix element for weak decay to hadrons (B -> D, D -> K). |
---|
730 | } else if (meMode == 22 || meMode == 23) { |
---|
731 | double x1 = 2. * event[iProd[1]].pAbs() / m0; |
---|
732 | wtME = x1 * (3. - 2. * x1); |
---|
733 | double xMax = min( 0.75, 2. * (1. - mSum / m0) ); |
---|
734 | wtMEmax = xMax * (3. - 2. * xMax); |
---|
735 | |
---|
736 | // Effective matrix element for gamma spectrum in B -> gamma + hadrons. |
---|
737 | } else if (meMode == 31) { |
---|
738 | double x1 = 2. * event[iProd[1]].e() / m0; |
---|
739 | wtME = pow3(x1); |
---|
740 | double x1Max = 1. - pow2(mSum / m0); |
---|
741 | wtMEmax = pow3(x1Max); |
---|
742 | } |
---|
743 | |
---|
744 | // If rejected, try again with new invariant masses. |
---|
745 | } while ( wtME < rndmPtr->flat() * wtMEmax ); |
---|
746 | |
---|
747 | // Boost decay products to the current frame. |
---|
748 | pInv[1].p( event[iProd[0]].p() ); |
---|
749 | for (int i = 1; i <= mult; ++i) event[iProd[i]].bst( pInv[1], mInv[1] ); |
---|
750 | |
---|
751 | // Done. |
---|
752 | return true; |
---|
753 | |
---|
754 | } |
---|
755 | |
---|
756 | //-------------------------------------------------------------------------- |
---|
757 | |
---|
758 | // Select mass of lepton pair in a Dalitz decay. |
---|
759 | |
---|
760 | bool ParticleDecays::dalitzMass() { |
---|
761 | |
---|
762 | // Mother and sum daughter masses. |
---|
763 | double mSum1 = 0; |
---|
764 | for (int i = 1; i <= mult - 2; ++i) mSum1 += mProd[i]; |
---|
765 | if (meMode == 13) mSum1 *= MSAFEDALITZ; |
---|
766 | double mSum2 = MSAFEDALITZ * (mProd[mult -1] + mProd[mult]); |
---|
767 | double mDiff = mProd[0] - mSum1 - mSum2; |
---|
768 | |
---|
769 | // Fail if too close or inconsistent. |
---|
770 | if (mDiff < mSafety) return false; |
---|
771 | if (idProd[mult - 1] + idProd[mult] != 0 |
---|
772 | || mProd[mult - 1] != mProd[mult]) { |
---|
773 | infoPtr->errorMsg("Error in ParticleDecays::dalitzMass:" |
---|
774 | " inconsistent flavour/mass assignments"); |
---|
775 | return false; |
---|
776 | } |
---|
777 | if ( meMode == 13 && (idProd[1] + idProd[2] != 0 |
---|
778 | || mProd[1] != mProd[2]) ) { |
---|
779 | infoPtr->errorMsg("Error in ParticleDecays::dalitzMass:" |
---|
780 | " inconsistent flavour/mass assignments"); |
---|
781 | return false; |
---|
782 | } |
---|
783 | |
---|
784 | // Case 1: one Dalitz pair. |
---|
785 | if (meMode == 11 || meMode == 12) { |
---|
786 | |
---|
787 | // Kinematical limits for gamma* squared mass. |
---|
788 | double sGamMin = pow2(mSum2); |
---|
789 | double sGamMax = pow2(mProd[0] - mSum1); |
---|
790 | // Select virtual gamma squared mass. Guessed form for meMode == 12. |
---|
791 | double sGam, wtGam; |
---|
792 | int loop = 0; |
---|
793 | do { |
---|
794 | if (++loop > NTRYDALITZ) return false; |
---|
795 | sGam = sGamMin * pow( sGamMax / sGamMin, rndmPtr->flat() ); |
---|
796 | wtGam = (1. + 0.5 * sGamMin / sGam) * sqrt(1. - sGamMin / sGam) |
---|
797 | * pow3(1. - sGam / sGamMax) * sRhoDal * (sRhoDal + wRhoDal) |
---|
798 | / ( pow2(sGam - sRhoDal) + sRhoDal * wRhoDal ); |
---|
799 | } while ( wtGam < rndmPtr->flat() ); |
---|
800 | |
---|
801 | // Store results in preparation for doing a one-less-body decay. |
---|
802 | --mult; |
---|
803 | mProd[mult] = sqrt(sGam); |
---|
804 | |
---|
805 | // Case 2: two Dalitz pairs. |
---|
806 | } else { |
---|
807 | |
---|
808 | // Kinematical limits for 1 + 2 and 3 + 4 gamma* masses. |
---|
809 | double s0 = pow2(mProd[0]); |
---|
810 | double s12Min = pow2(mSum1); |
---|
811 | double s12Max = pow2(mProd[0] - mSum2); |
---|
812 | double s34Min = pow2(mSum2); |
---|
813 | double s34Max = pow2(mProd[0] - mSum1); |
---|
814 | |
---|
815 | // Select virtual gamma squared masses. Guessed form for meMode == 13. |
---|
816 | double s12, s34, wt12, wt34, wtPAbs, wtAll; |
---|
817 | int loop = 0; |
---|
818 | do { |
---|
819 | if (++loop > NTRYDALITZ) return false; |
---|
820 | s12 = s12Min * pow( s12Max / s12Min, rndmPtr->flat() ); |
---|
821 | wt12 = (1. + 0.5 * s12Min / s12) * sqrt(1. - s12Min / s12) |
---|
822 | * sRhoDal * (sRhoDal + wRhoDal) |
---|
823 | / ( pow2(s12 - sRhoDal) + sRhoDal * wRhoDal ); |
---|
824 | s34 = s34Min * pow( s34Max / s34Min, rndmPtr->flat() ); |
---|
825 | wt34 = (1. + 0.5 * s34Min / s34) * sqrt(1. - s34Min / s34) |
---|
826 | * sRhoDal * (sRhoDal + wRhoDal) |
---|
827 | / ( pow2(s34 - sRhoDal) + sRhoDal * wRhoDal ); |
---|
828 | wtPAbs = sqrtpos( pow2(1. - (s12 + s34)/ s0) |
---|
829 | - 4. * s12 * s34 / (s0 * s0) ); |
---|
830 | wtAll = wt12 * wt34 * pow3(wtPAbs); |
---|
831 | if (wtAll > 1.) infoPtr->errorMsg( |
---|
832 | "Error in ParticleDecays::dalitzMass: weight > 1"); |
---|
833 | } while (wtAll < rndmPtr->flat()); |
---|
834 | |
---|
835 | // Store results in preparation for doing a two-body decay. |
---|
836 | mult = 2; |
---|
837 | mProd[1] = sqrt(s12); |
---|
838 | mProd[2] = sqrt(s34); |
---|
839 | } |
---|
840 | |
---|
841 | // Done. |
---|
842 | return true; |
---|
843 | |
---|
844 | } |
---|
845 | |
---|
846 | //-------------------------------------------------------------------------- |
---|
847 | |
---|
848 | // Do kinematics of gamma* -> l- l+ in Dalitz decay. |
---|
849 | |
---|
850 | bool ParticleDecays::dalitzKinematics(Event& event) { |
---|
851 | |
---|
852 | // Restore multiplicity. |
---|
853 | int nDal = (meMode < 13) ? 1 : 2; |
---|
854 | mult += nDal; |
---|
855 | |
---|
856 | // Loop over one or two lepton pairs. |
---|
857 | for (int iDal = 0; iDal < nDal; ++iDal) { |
---|
858 | |
---|
859 | // References to the particles involved. |
---|
860 | Particle& decayer = event[iProd[0]]; |
---|
861 | Particle& prodA = (iDal == 0) ? event[iProd[mult - 1]] |
---|
862 | : event[iProd[1]]; |
---|
863 | Particle& prodB = (iDal == 0) ? event[iProd[mult]] |
---|
864 | : event[iProd[2]]; |
---|
865 | |
---|
866 | // Reconstruct required rotations and boosts backwards. |
---|
867 | Vec4 pDec = decayer.p(); |
---|
868 | int iGam = (meMode < 13) ? mult - 1 : 2 - iDal; |
---|
869 | Vec4 pGam = event[iProd[iGam]].p(); |
---|
870 | pGam.bstback( pDec, decayer.m() ); |
---|
871 | double phiGam = pGam.phi(); |
---|
872 | pGam.rot( 0., -phiGam); |
---|
873 | double thetaGam = pGam.theta(); |
---|
874 | pGam.rot( -thetaGam, 0.); |
---|
875 | |
---|
876 | // Masses and phase space in gamma* rest frame. |
---|
877 | double mGam = (meMode < 13) ? mProd[mult - 1] : mProd[2 - iDal]; |
---|
878 | double mA = prodA.m(); |
---|
879 | double mB = prodB.m(); |
---|
880 | double mGamMin = MSAFEDALITZ * (mA + mB); |
---|
881 | double mGamRat = pow2(mGamMin / mGam); |
---|
882 | double pGamAbs = 0.5 * sqrtpos( (mGam - mA - mB) * (mGam + mA + mB) ); |
---|
883 | |
---|
884 | // Set up decay in gamma* rest frame, reference along +z axis. |
---|
885 | double cosTheta, cos2Theta; |
---|
886 | do { |
---|
887 | cosTheta = 2. * rndmPtr->flat() - 1.; |
---|
888 | cos2Theta = cosTheta * cosTheta; |
---|
889 | } while ( 1. + cos2Theta + mGamRat * (1. - cos2Theta) |
---|
890 | < 2. * rndmPtr->flat() ); |
---|
891 | double sinTheta = sqrt(1. - cosTheta*cosTheta); |
---|
892 | double phi = 2. * M_PI * rndmPtr->flat(); |
---|
893 | double pX = pGamAbs * sinTheta * cos(phi); |
---|
894 | double pY = pGamAbs * sinTheta * sin(phi); |
---|
895 | double pZ = pGamAbs * cosTheta; |
---|
896 | double eA = sqrt( mA*mA + pGamAbs*pGamAbs); |
---|
897 | double eB = sqrt( mB*mB + pGamAbs*pGamAbs); |
---|
898 | prodA.p( pX, pY, pZ, eA); |
---|
899 | prodB.p( -pX, -pY, -pZ, eB); |
---|
900 | |
---|
901 | // Boost to lab frame. |
---|
902 | prodA.bst( pGam, mGam); |
---|
903 | prodB.bst( pGam, mGam); |
---|
904 | prodA.rot( thetaGam, phiGam); |
---|
905 | prodB.rot( thetaGam, phiGam); |
---|
906 | prodA.bst( pDec, decayer.m() ); |
---|
907 | prodB.bst( pDec, decayer.m() ); |
---|
908 | } |
---|
909 | |
---|
910 | // Done. |
---|
911 | return true; |
---|
912 | |
---|
913 | } |
---|
914 | |
---|
915 | //-------------------------------------------------------------------------- |
---|
916 | |
---|
917 | // Translate a partonic content into a set of actual hadrons. |
---|
918 | |
---|
919 | bool ParticleDecays::pickHadrons() { |
---|
920 | |
---|
921 | // Find partonic decay products. Rest are known id's, mainly hadrons, |
---|
922 | // when necessary shuffled to beginning of idProd list. |
---|
923 | idPartons.resize(0); |
---|
924 | int nPartons = 0; |
---|
925 | int nKnown = 0; |
---|
926 | bool closedGLoop = false; |
---|
927 | for (int i = 1; i <= mult; ++i) { |
---|
928 | int idAbs = abs(idProd[i]); |
---|
929 | if ( idAbs < 9 || (idAbs > 1000 && idAbs < 10000 && (idAbs/10)%10 == 0) |
---|
930 | || idAbs == 81 || idAbs == 82 || idAbs == 83) { |
---|
931 | ++nPartons; |
---|
932 | idPartons.push_back(idProd[i]); |
---|
933 | if (idAbs == 83) closedGLoop = true; |
---|
934 | } else { |
---|
935 | ++nKnown; |
---|
936 | if (nPartons > 0) { |
---|
937 | idProd[nKnown] = idProd[i]; |
---|
938 | mProd[nKnown] = mProd[i]; |
---|
939 | } |
---|
940 | } |
---|
941 | } |
---|
942 | |
---|
943 | // Replace generic spectator flavour code by the actual one. |
---|
944 | for (int i = 0; i < nPartons; ++i) { |
---|
945 | int idPart = idPartons[i]; |
---|
946 | int idNew = idPart; |
---|
947 | if (idPart == 81) { |
---|
948 | int idAbs = abs(idDec); |
---|
949 | if ( (idAbs/1000)%10 == 0 ) { |
---|
950 | idNew = -(idAbs/10)%10; |
---|
951 | if ((idAbs/100)%2 == 1) idNew = -idNew; |
---|
952 | } else if ( (idAbs/100)%10 >= (idAbs/10)%10 ) |
---|
953 | idNew = 100 * ((idAbs/10)%100) + 3; |
---|
954 | else idNew = 1000 * ((idAbs/10)%10) + 100 * ((idAbs/100)%10) + 1; |
---|
955 | if (idDec < 0) idNew = -idNew; |
---|
956 | |
---|
957 | // Replace generic random flavour by a randomly selected one. |
---|
958 | } else if (idPart == 82 || idPart == 83) { |
---|
959 | double mFlav; |
---|
960 | do { |
---|
961 | int idDummy = -flavSelPtr->pickLightQ(); |
---|
962 | FlavContainer flavDummy(idDummy, idPart - 82); |
---|
963 | do idNew = flavSelPtr->pick(flavDummy).id; |
---|
964 | while (idNew == 0); |
---|
965 | mFlav = particleDataPtr->constituentMass(idNew); |
---|
966 | } while (2. * mFlav + stopMass > mProd[0]); |
---|
967 | } else if (idPart == -82 || idPart == -83) { |
---|
968 | idNew = -idPartons[i-1]; |
---|
969 | } |
---|
970 | idPartons[i] = idNew; |
---|
971 | } |
---|
972 | |
---|
973 | // Determine whether fixed multiplicity or to be selected at random. |
---|
974 | int nMin = max( 2, nKnown + nPartons / 2); |
---|
975 | if (meMode == 23) nMin = 3; |
---|
976 | if (meMode > 41 && meMode <= 50) nMin = meMode - 40; |
---|
977 | if (meMode > 51 && meMode <= 60) nMin = meMode - 50; |
---|
978 | int nFix = 0; |
---|
979 | if (meMode == 0) nFix = nMin; |
---|
980 | if (meMode == 11) nFix = 3; |
---|
981 | if (meMode == 12) nFix = 4; |
---|
982 | if (meMode > 61 && meMode <= 70) nFix = meMode - 60; |
---|
983 | if (meMode > 71 && meMode <= 80) nFix = meMode - 70; |
---|
984 | if (nFix > 0 && nKnown + nPartons/2 > nFix) return false; |
---|
985 | |
---|
986 | // Initial values for loop to set new hadronic content. |
---|
987 | int nFilled, nTotal, nNew, nSpec, nLeft; |
---|
988 | double mDiff; |
---|
989 | int nTry = 0; |
---|
990 | bool diquarkClash = false; |
---|
991 | bool usedChannel = false; |
---|
992 | |
---|
993 | // Begin loop; interrupt if multiple tries fail. |
---|
994 | do { |
---|
995 | ++nTry; |
---|
996 | if (nTry > NTRYPICK) return false; |
---|
997 | |
---|
998 | // Initialize variables inside new try. |
---|
999 | nFilled = nKnown + 1; |
---|
1000 | idProd.resize(nFilled); |
---|
1001 | mProd.resize(nFilled); |
---|
1002 | nTotal = nKnown; |
---|
1003 | nSpec = 0; |
---|
1004 | nLeft = nPartons; |
---|
1005 | mDiff = mProd[0]; |
---|
1006 | for (int i = 1; i < nFilled; ++i) mDiff -= mProd[i]; |
---|
1007 | diquarkClash = false; |
---|
1008 | usedChannel = false; |
---|
1009 | |
---|
1010 | // For weak decays collapse spectators to one particle. |
---|
1011 | if ( (meMode == 22 || meMode == 23) && nLeft > 1) { |
---|
1012 | FlavContainer flav1( idPartons[nPartons - 2] ); |
---|
1013 | FlavContainer flav2( idPartons[nPartons - 1] ); |
---|
1014 | int idHad; |
---|
1015 | do idHad = flavSelPtr->combine( flav1, flav2); |
---|
1016 | while (idHad == 0); |
---|
1017 | double mHad = particleDataPtr->mass(idHad); |
---|
1018 | mDiff -= mHad; |
---|
1019 | idProd.push_back( idHad); |
---|
1020 | mProd.push_back( mHad); |
---|
1021 | ++nFilled; |
---|
1022 | nSpec = 1; |
---|
1023 | nLeft -= 2; |
---|
1024 | } |
---|
1025 | |
---|
1026 | // If there are partons left, then determine how many hadrons to pick. |
---|
1027 | if (nLeft > 0) { |
---|
1028 | |
---|
1029 | // For B -> gamma + X use geometrical distribution. |
---|
1030 | if (meMode == 31) { |
---|
1031 | double geom = rndmPtr->flat(); |
---|
1032 | nTotal = 1; |
---|
1033 | do { |
---|
1034 | ++nTotal; |
---|
1035 | geom *= 2.; |
---|
1036 | } while (geom < 1. && nTotal < 10); |
---|
1037 | |
---|
1038 | // Calculate mass excess and from there average multiplicity. |
---|
1039 | } else if (nFix == 0) { |
---|
1040 | double multIncreaseNow = (meMode == 23) |
---|
1041 | ? multIncreaseWeak : multIncrease; |
---|
1042 | double mDiffPS = mDiff; |
---|
1043 | for (int i = 0; i < nLeft; ++i) |
---|
1044 | mDiffPS -= particleDataPtr->constituentMass( idPartons[i] ); |
---|
1045 | double average = 0.5 * (nKnown + nSpec) + 0.25 * nPartons |
---|
1046 | + multIncreaseNow * log( max( 1.1, mDiffPS / multRefMass ) ); |
---|
1047 | if (closedGLoop) average += multGoffset; |
---|
1048 | |
---|
1049 | // Pick multiplicity according to Poissonian. |
---|
1050 | double value = 1.; |
---|
1051 | double sum = 1.; |
---|
1052 | for (int nNow = nMin + 1; nNow <= 10; ++nNow) { |
---|
1053 | value *= average / nNow; |
---|
1054 | sum += value; |
---|
1055 | } |
---|
1056 | nTotal = nMin; |
---|
1057 | value = 1.; |
---|
1058 | sum *= rndmPtr->flat(); |
---|
1059 | sum -= value; |
---|
1060 | if (sum > 0.) do { |
---|
1061 | ++nTotal; |
---|
1062 | value *= average / nTotal; |
---|
1063 | sum -= value; |
---|
1064 | } while (sum > 0. && nTotal < 10); |
---|
1065 | |
---|
1066 | // Alternatively predetermined multiplicity. |
---|
1067 | } else { |
---|
1068 | nTotal = nFix; |
---|
1069 | } |
---|
1070 | nNew = nTotal - nKnown - nSpec; |
---|
1071 | |
---|
1072 | // Set up ends of fragmented system, as copy of idPartons. |
---|
1073 | flavEnds.resize(0); |
---|
1074 | for (int i = 0; i < nLeft; ++i) { |
---|
1075 | flavEnds.push_back( FlavContainer(idPartons[i]) ); |
---|
1076 | if (abs(idPartons[i]) > 100) flavSelPtr->assignPopQ( flavEnds[i] ); |
---|
1077 | } |
---|
1078 | |
---|
1079 | // Fragment off at random, but save nLeft/2 for final recombination. |
---|
1080 | if (nNew > nLeft/2) { |
---|
1081 | FlavContainer flavNew; |
---|
1082 | int idHad; |
---|
1083 | for (int i = 0; i < nNew - nLeft/2; ++i) { |
---|
1084 | // When four quarks consider last one to be spectator. |
---|
1085 | int iEnd = int( (nLeft - 1.) * rndmPtr->flat() ); |
---|
1086 | // Pick new flavour and form a new hadron. |
---|
1087 | do { |
---|
1088 | flavNew = flavSelPtr->pick( flavEnds[iEnd] ); |
---|
1089 | idHad = flavSelPtr->combine( flavEnds[iEnd], flavNew); |
---|
1090 | } while (idHad == 0); |
---|
1091 | // Store new hadron and endpoint flavour. |
---|
1092 | idProd.push_back( idHad); |
---|
1093 | flavEnds[iEnd].anti(flavNew); |
---|
1094 | } |
---|
1095 | } |
---|
1096 | |
---|
1097 | // When only two quarks left, combine to form final hadron. |
---|
1098 | if (nLeft == 2) { |
---|
1099 | int idHad; |
---|
1100 | if ( abs(flavEnds[0].id) > 8 && abs(flavEnds[1].id) > 8) |
---|
1101 | diquarkClash = true; |
---|
1102 | else { |
---|
1103 | do idHad = flavSelPtr->combine( flavEnds[0], flavEnds[1]); |
---|
1104 | while (idHad == 0); |
---|
1105 | idProd.push_back( idHad); |
---|
1106 | } |
---|
1107 | |
---|
1108 | // If four quarks, decide how to pair them up. |
---|
1109 | } else { |
---|
1110 | int iEnd1 = 0; |
---|
1111 | int iEnd2 = 1; |
---|
1112 | int iEnd3 = 2; |
---|
1113 | int iEnd4 = 3; |
---|
1114 | if ( rndmPtr->flat() < colRearrange) iEnd2 = 3; |
---|
1115 | int relColSign = |
---|
1116 | ( (flavEnds[iEnd1].id > 0 && flavEnds[iEnd1].id < 9) |
---|
1117 | || flavEnds[iEnd1].id < -10 ) ? 1 : -1; |
---|
1118 | if ( (flavEnds[iEnd2].id < 0 && flavEnds[iEnd2].id > -9) |
---|
1119 | || flavEnds[iEnd2].id > 10 ) relColSign *= -1; |
---|
1120 | if (relColSign == 1) iEnd2 = 2; |
---|
1121 | if (iEnd2 == 2) iEnd3 = 1; |
---|
1122 | if (iEnd2 == 3) iEnd4 = 1; |
---|
1123 | |
---|
1124 | // Then combine to get final two hadrons. |
---|
1125 | int idHad; |
---|
1126 | if ( abs(flavEnds[iEnd1].id) > 8 && abs(flavEnds[iEnd2].id) > 8) |
---|
1127 | diquarkClash = true; |
---|
1128 | else { |
---|
1129 | do idHad = flavSelPtr->combine( flavEnds[iEnd1], flavEnds[iEnd2]); |
---|
1130 | while (idHad == 0); |
---|
1131 | idProd.push_back( idHad); |
---|
1132 | } |
---|
1133 | if ( abs(flavEnds[iEnd3].id) > 8 && abs(flavEnds[iEnd4].id) > 8) |
---|
1134 | diquarkClash = true; |
---|
1135 | else { |
---|
1136 | do idHad = flavSelPtr->combine( flavEnds[iEnd3], flavEnds[iEnd4]); |
---|
1137 | while (idHad == 0); |
---|
1138 | idProd.push_back( idHad); |
---|
1139 | } |
---|
1140 | } |
---|
1141 | |
---|
1142 | // Find masses of the new hadrons. |
---|
1143 | for (int i = nFilled; i < int(idProd.size()) ; ++i) { |
---|
1144 | double mHad = particleDataPtr->mass(idProd[i]); |
---|
1145 | mProd.push_back( mHad); |
---|
1146 | mDiff -= mHad; |
---|
1147 | } |
---|
1148 | } |
---|
1149 | |
---|
1150 | // Optional: check that this decay mode is not explicitly defined. |
---|
1151 | if ( (meMode > 61 && meMode <= 80) && mDiff > mSafety && !diquarkClash ) { |
---|
1152 | int idMatch[10], idPNow; |
---|
1153 | usedChannel = false; |
---|
1154 | bool matched = false; |
---|
1155 | // Loop through all channels. Done if not same multiplicity. |
---|
1156 | for (int i = 0; i < decDataPtr->sizeChannels(); ++i) { |
---|
1157 | DecayChannel& channel = decDataPtr->channel(i); |
---|
1158 | if (channel.multiplicity() != nTotal) continue; |
---|
1159 | for (int k = 0; k < nTotal; ++k) idMatch[k] = channel.product(k); |
---|
1160 | // Match particles one by one until fail. |
---|
1161 | // Do not distinguish K0/K0bar/K0short/K0long at this stage. |
---|
1162 | for (int j = 0; j < nTotal; ++j) { |
---|
1163 | matched = false; |
---|
1164 | idPNow = idProd[j + 1]; |
---|
1165 | if (idPNow == -311 || idPNow == 130 || idPNow == 310) idPNow = 311; |
---|
1166 | for (int k = 0; k < nTotal - j; ++k) |
---|
1167 | if (idMatch[k] == idPNow || (idMatch[k] == -311 && idPNow == 311)) { |
---|
1168 | // Compress list of unmatched when matching worked. |
---|
1169 | idMatch[k] = idMatch[nTotal - 1 - j]; |
---|
1170 | matched = true; |
---|
1171 | break; |
---|
1172 | } |
---|
1173 | if (!matched) break; |
---|
1174 | } |
---|
1175 | // If matching worked, then chosen channel to be rejected. |
---|
1176 | if (matched) {usedChannel = true; break;} |
---|
1177 | } |
---|
1178 | } |
---|
1179 | |
---|
1180 | // Keep on trying until enough phase space and no clash. |
---|
1181 | } while (mDiff < mSafety || diquarkClash || usedChannel); |
---|
1182 | |
---|
1183 | // Update particle multiplicity. |
---|
1184 | mult = idProd.size() - 1; |
---|
1185 | |
---|
1186 | // For Dalitz decays shuffle Dalitz pair to the end of the list. |
---|
1187 | if (meMode == 11 || meMode == 12) { |
---|
1188 | int iL1 = 0; |
---|
1189 | int iL2 = 0; |
---|
1190 | for (int i = 1; i <= mult; ++i) { |
---|
1191 | if (idProd[i] == 11 || idProd[i] == 13 || idProd[i] == 15) iL1 = i; |
---|
1192 | if (idProd[i] == -11 || idProd[i] == -13 || idProd[i] == -15) iL2 = i; |
---|
1193 | } |
---|
1194 | if (iL1 > 0 && iL2 > 0) { |
---|
1195 | int idL1 = idProd[iL1]; |
---|
1196 | int idL2 = idProd[iL2]; |
---|
1197 | double mL1 = mProd[iL1]; |
---|
1198 | double mL2 = mProd[iL2]; |
---|
1199 | int iMove = 0; |
---|
1200 | for (int i = 1; i <= mult; ++i) if (i != iL1 && i != iL2) { |
---|
1201 | ++iMove; |
---|
1202 | idProd[iMove] = idProd[i]; |
---|
1203 | mProd[iMove] = mProd[i]; |
---|
1204 | } |
---|
1205 | idProd[mult - 1] = idL1; |
---|
1206 | idProd[mult] = idL2; |
---|
1207 | mProd[mult - 1] = mL1; |
---|
1208 | mProd[mult] = mL2; |
---|
1209 | } |
---|
1210 | } |
---|
1211 | |
---|
1212 | // Done. |
---|
1213 | return true; |
---|
1214 | |
---|
1215 | } |
---|
1216 | |
---|
1217 | //-------------------------------------------------------------------------- |
---|
1218 | |
---|
1219 | // Set colour flow and scale in a decay explicitly to partons. |
---|
1220 | |
---|
1221 | bool ParticleDecays::setColours(Event& event) { |
---|
1222 | |
---|
1223 | // Decay to q qbar (or qbar q). |
---|
1224 | if (meMode == 91 && idProd[1] > 0 && idProd[1] < 9) { |
---|
1225 | int newCol = event.nextColTag(); |
---|
1226 | cols[1] = newCol; |
---|
1227 | acols[2] = newCol; |
---|
1228 | } else if (meMode == 91 && idProd[1] < 0 && idProd[1] > -9) { |
---|
1229 | int newCol = event.nextColTag(); |
---|
1230 | cols[2] = newCol; |
---|
1231 | acols[1] = newCol; |
---|
1232 | |
---|
1233 | // Decay to g g. |
---|
1234 | } else if (meMode == 91 && idProd[1] == 21) { |
---|
1235 | int newCol1 = event.nextColTag(); |
---|
1236 | int newCol2 = event.nextColTag(); |
---|
1237 | cols[1] = newCol1; |
---|
1238 | acols[1] = newCol2; |
---|
1239 | cols[2] = newCol2; |
---|
1240 | acols[2] = newCol1; |
---|
1241 | |
---|
1242 | // Decay to g g g. |
---|
1243 | } else if (meMode == 92 && idProd[1] == 21 && idProd[2] == 21 |
---|
1244 | && idProd[3] == 21) { |
---|
1245 | int newCol1 = event.nextColTag(); |
---|
1246 | int newCol2 = event.nextColTag(); |
---|
1247 | int newCol3 = event.nextColTag(); |
---|
1248 | cols[1] = newCol1; |
---|
1249 | acols[1] = newCol2; |
---|
1250 | cols[2] = newCol2; |
---|
1251 | acols[2] = newCol3; |
---|
1252 | cols[3] = newCol3; |
---|
1253 | acols[3] = newCol1; |
---|
1254 | |
---|
1255 | // Decay to g g gamma: locate which is gamma. |
---|
1256 | } else if (meMode == 92) { |
---|
1257 | int iGlu1 = (idProd[1] == 21) ? 1 : 3; |
---|
1258 | int iGlu2 = (idProd[2] == 21) ? 2 : 3; |
---|
1259 | int newCol1 = event.nextColTag(); |
---|
1260 | int newCol2 = event.nextColTag(); |
---|
1261 | cols[iGlu1] = newCol1; |
---|
1262 | acols[iGlu1] = newCol2; |
---|
1263 | cols[iGlu2] = newCol2; |
---|
1264 | acols[iGlu2] = newCol1; |
---|
1265 | |
---|
1266 | // Unknown decay mode means failure. |
---|
1267 | } else return false; |
---|
1268 | |
---|
1269 | // Set maximum scale to be mass of decaying particle. |
---|
1270 | scale = mProd[0]; |
---|
1271 | |
---|
1272 | // Done. |
---|
1273 | return true; |
---|
1274 | |
---|
1275 | } |
---|
1276 | |
---|
1277 | //========================================================================== |
---|
1278 | |
---|
1279 | } // end namespace Pythia8 |
---|
1280 | |
---|