]>
Commit | Line | Data |
---|---|---|
1 | // TauDecays.cc is a part of the PYTHIA event generator. | |
2 | // Copyright (C) 2012 Philip Ilten, 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 TauDecays class. | |
7 | ||
8 | #include "TauDecays.h" | |
9 | ||
10 | namespace Pythia8 { | |
11 | ||
12 | //========================================================================== | |
13 | ||
14 | // The TauDecays class. | |
15 | ||
16 | //-------------------------------------------------------------------------- | |
17 | ||
18 | // Constants: could be changed here if desired, but normally should not. | |
19 | // These are of technical nature, as described for each. | |
20 | ||
21 | // Number of times to try a decay channel. | |
22 | const int TauDecays::NTRYCHANNEL = 10; | |
23 | ||
24 | // Number of times to try a decay sampling. | |
25 | const int TauDecays::NTRYDECAY = 10000; | |
26 | //const int TauDecays::NTRYDECAY = 100000; | |
27 | ||
28 | // These numbers are hardwired empirical parameters, | |
29 | // intended to speed up the M-generator. | |
30 | const double TauDecays::WTCORRECTION[11] = { 1., 1., 1., | |
31 | 2., 5., 15., 60., 250., 1250., 7000., 50000. }; | |
32 | ||
33 | //-------------------------------------------------------------------------- | |
34 | ||
35 | // Initialize the TauDecays class with the necessary pointers to info, | |
36 | // particle data, random numbers, and Standard Model couplings. | |
37 | // Additionally, the necessary matrix elements are initialized with the | |
38 | // Standard Model couplings, and particle data pointers. | |
39 | ||
40 | void TauDecays::init(Info* infoPtrIn, Settings* settingsPtrIn, | |
41 | ParticleData* particleDataPtrIn, Rndm* rndmPtrIn, | |
42 | Couplings* couplingsPtrIn) { | |
43 | ||
44 | // Set the pointers. | |
45 | infoPtr = infoPtrIn; | |
46 | settingsPtr = settingsPtrIn; | |
47 | particleDataPtr = particleDataPtrIn; | |
48 | rndmPtr = rndmPtrIn; | |
49 | couplingsPtr = couplingsPtrIn; | |
50 | ||
51 | // Initialize the hard matrix elements. | |
52 | hmeTwoFermions2W2TwoFermions .initPointers(particleDataPtr, couplingsPtr); | |
53 | hmeTwoFermions2Z2TwoFermions .initPointers(particleDataPtr, couplingsPtr); | |
54 | hmeTwoFermions2Gamma2TwoFermions .initPointers(particleDataPtr, | |
55 | couplingsPtr); | |
56 | hmeTwoFermions2GammaZ2TwoFermions.initPointers(particleDataPtr, | |
57 | couplingsPtr); | |
58 | hmeZ2TwoFermions .initPointers(particleDataPtr, couplingsPtr); | |
59 | hmeHiggsEven2TwoFermions .initPointers(particleDataPtr, couplingsPtr); | |
60 | hmeHiggsOdd2TwoFermions .initPointers(particleDataPtr, couplingsPtr); | |
61 | hmeHiggsCharged2TwoFermions .initPointers(particleDataPtr, couplingsPtr); | |
62 | hmeUnpolarized .initPointers(particleDataPtr, couplingsPtr); | |
63 | ||
64 | // Initialize the tau decay matrix elements. | |
65 | hmeTau2Meson .initPointers(particleDataPtr, couplingsPtr); | |
66 | hmeTau2TwoLeptons .initPointers(particleDataPtr, couplingsPtr); | |
67 | hmeTau2TwoMesonsViaVector .initPointers(particleDataPtr, couplingsPtr); | |
68 | hmeTau2TwoMesonsViaVectorScalar.initPointers(particleDataPtr, couplingsPtr); | |
69 | hmeTau2ThreePions .initPointers(particleDataPtr, couplingsPtr); | |
70 | hmeTau2ThreeMesonsWithKaons .initPointers(particleDataPtr, couplingsPtr); | |
71 | hmeTau2ThreeMesonsGeneric .initPointers(particleDataPtr, couplingsPtr); | |
72 | hmeTau2TwoPionsGamma .initPointers(particleDataPtr, couplingsPtr); | |
73 | hmeTau2FourPions .initPointers(particleDataPtr, couplingsPtr); | |
74 | hmeTau2FivePions .initPointers(particleDataPtr, couplingsPtr); | |
75 | hmeTau2PhaseSpace .initPointers(particleDataPtr, couplingsPtr); | |
76 | ||
77 | // User selected tau decay mode. | |
78 | tauModeSave = settingsPtr->mode("ParticleDecays:sophisticatedTau"); | |
79 | ||
80 | // User selected tau decay mother. | |
81 | tauMotherSave = settingsPtr->mode("ParticleDecays:tauMother"); | |
82 | ||
83 | // User selected tau polarization. | |
84 | polSave = settingsPtr->parm("ParticleDecays:tauPolarization"); | |
85 | ||
86 | } | |
87 | ||
88 | //-------------------------------------------------------------------------- | |
89 | ||
90 | // Main method of the TauDecays class. Pass the index of the tau requested | |
91 | // to be decayed along with the event record in which the tau exists. The | |
92 | // tau is then decayed with proper spin correlations as well any partner. | |
93 | // The children of the decays are written to the event record, and if the | |
94 | // decays were succesful, a return value of true is supplied. | |
95 | ||
96 | bool TauDecays::decay(int idxOut1, Event& event) { | |
97 | ||
98 | // User selected tau decay mode, mother, and polarization. | |
99 | tauMode = tauModeSave; | |
100 | tauMother = tauMotherSave; | |
101 | polarization = polSave; | |
102 | ||
103 | // Set the first outgoing particle of the hard process. | |
104 | out1 = HelicityParticle(event[idxOut1]); | |
105 | out1.idx = idxOut1; | |
106 | ||
107 | // Begin PS April 2012. | |
108 | // Check if this tau already has helicity information (eg from LHEF). | |
109 | bool hasHelicity = false; | |
110 | double helicityNow = 0.; | |
111 | if (tauMode >= 1 && abs(out1.pol()) <= 1.001) { | |
112 | hasHelicity = true; | |
113 | helicityNow = out1.pol(); | |
114 | } | |
115 | // End PS April 2012. | |
116 | ||
117 | // Find the mediator of the hard process. Create temporary copy. | |
118 | int idxMediator = out1.mother1(); | |
119 | int idxFirstOut1 = idxOut1; | |
120 | while(idxMediator > 0 && event[idxMediator].id() == out1.id()) { | |
121 | idxFirstOut1 = idxMediator; | |
122 | idxMediator = event[idxMediator].mother1(); | |
123 | } | |
124 | Particle medTmp = event[idxMediator]; | |
125 | ||
126 | // Find and set up the incoming particles of the hard process. | |
127 | int idxIn1 = medTmp.mother1(); | |
128 | int idxIn2 = medTmp.mother2(); | |
129 | while(idxIn1 > 0 && event[idxIn1].id() == medTmp.id()) { | |
130 | idxIn1 = event[idxIn1].mother1(); | |
131 | idxIn2 = event[idxIn2].mother2(); | |
132 | } | |
133 | in1 = HelicityParticle(event[idxIn1]); | |
134 | in1.idx = idxIn1; | |
135 | in1.direction = -1; | |
136 | in2 = HelicityParticle(event[idxIn2]); | |
137 | in2.idx = idxIn2; | |
138 | in2.direction = -1; | |
139 | ||
140 | // Find and set up the second outgoing particle of the hard process. | |
141 | int idxOut2 = (medTmp.daughter1() == idxFirstOut1) | |
142 | ? medTmp.daughter2() : medTmp.daughter1(); | |
143 | while (idxOut2 > 0 && event[idxOut2].daughter1() != 0) { | |
144 | idxOut2 = event[idxOut2].daughter1(); | |
145 | } | |
146 | out2 = HelicityParticle(event[idxOut2]); | |
147 | out2.idx = idxOut2; | |
148 | ||
149 | // Set up the mediator. Special case for dipole shower, | |
150 | // where a massless photon can branch to a tau pair. | |
151 | if (medTmp.id() == 22 && out2.idAbs() == 15 | |
152 | && medTmp.m() < out1.m() + out2.m()) { | |
153 | Vec4 pTmp = out1.p() + out2.p(); | |
154 | medTmp.p( pTmp); | |
155 | medTmp.m( pTmp.mCalc() ); | |
156 | } | |
157 | mediator = HelicityParticle(medTmp); | |
158 | mediator.idx = idxMediator; | |
159 | mediator.direction = -1; | |
160 | ||
161 | // Set the particles vector. | |
162 | particles.clear(); | |
163 | particles.push_back(in1); | |
164 | particles.push_back(in2); | |
165 | particles.push_back(out1); | |
166 | particles.push_back(out2); | |
167 | ||
168 | // Set the hard matrix element. | |
169 | // Polarized tau (decayed one by one). | |
170 | if (hasHelicity) { | |
171 | correlated = false; | |
172 | ||
173 | // Produced from a W. | |
174 | } else if (abs(mediator.id()) == 24) { | |
175 | // Produced from quarks: s-channel. | |
176 | if (abs(in1.id()) <= 18 && abs(in2.id()) <= 18) | |
177 | hardME = hmeTwoFermions2W2TwoFermions.initChannel(particles); | |
178 | // Produced from quarks: t-channel. | |
179 | else if (abs(in1.id()) <= 18 || abs(in2.id()) <= 18) { | |
180 | bool fermion = (abs(in1.id()) <= 18) ? 0 : 1; | |
181 | particles[!fermion] | |
182 | = (event[particles[fermion].daughter1()].id() == mediator.id()) | |
183 | ? HelicityParticle(event[particles[fermion].daughter2()]) | |
184 | : HelicityParticle(event[particles[fermion].daughter1()]); | |
185 | particles[!fermion].direction = 1; | |
186 | if (abs(particles[!fermion].id()) <= 18) | |
187 | hardME = hmeTwoFermions2W2TwoFermions.initChannel(particles); | |
188 | else { | |
189 | infoPtr->errorMsg("Warning in TauDecays::decay: unknown " | |
190 | "tau production, assuming unpolarized and uncorrelated"); | |
191 | hardME = hmeUnpolarized.initChannel(particles); | |
192 | } | |
193 | // Unknown W production: assume negative helicity. | |
194 | } else if (tauMode == 1) { | |
195 | tauMode = 3; | |
196 | polarization = -1; | |
197 | } | |
198 | correlated = false; | |
199 | ||
200 | // Produced from a photon. | |
201 | } else if (abs(mediator.id()) == 22 && abs(in1.id()) <= 18) { | |
202 | particles.push_back(mediator); | |
203 | hardME = hmeTwoFermions2Gamma2TwoFermions.initChannel(particles); | |
204 | correlated = true; | |
205 | ||
206 | // Produced from a photon/Z. | |
207 | } else if (abs(mediator.id()) == 23 && abs(in1.id()) <= 18) { | |
208 | particles.push_back(mediator); | |
209 | if (settingsPtr->mode("WeakZ0:gmZmode") == 0) | |
210 | hardME = hmeTwoFermions2GammaZ2TwoFermions.initChannel(particles); | |
211 | else if (settingsPtr->mode("WeakZ0:gmZmode") == 1) | |
212 | hardME = hmeTwoFermions2Gamma2TwoFermions.initChannel(particles); | |
213 | else if (settingsPtr->mode("WeakZ0:gmZmode") == 2) | |
214 | hardME = hmeTwoFermions2Z2TwoFermions.initChannel(particles); | |
215 | correlated = true; | |
216 | ||
217 | // Unkown Z production: assume unpolarized Z. | |
218 | } else if (abs(mediator.id()) == 23) { | |
219 | particles[1] = mediator; | |
220 | hardME = hmeZ2TwoFermions.initChannel(particles); | |
221 | correlated = true; | |
222 | ||
223 | // Produced from a CP even Higgs. | |
224 | } else if (abs(mediator.id()) == 25 || abs(mediator.id()) == 35) { | |
225 | hardME = hmeHiggsEven2TwoFermions.initChannel(particles); | |
226 | correlated = true; | |
227 | ||
228 | // Produced from a CP odd Higgs. | |
229 | } else if (abs(mediator.id()) == 36) { | |
230 | hardME = hmeHiggsOdd2TwoFermions.initChannel(particles); | |
231 | correlated = true; | |
232 | ||
233 | // Produced from a charged Higgs. | |
234 | } else if (abs(mediator.id()) == 37) { | |
235 | hardME = hmeHiggsCharged2TwoFermions.initChannel(particles); | |
236 | correlated = false; | |
237 | ||
238 | // Produced from a D or B hadron decay with a single tau. | |
239 | } else if ((abs(mediator.id()) == 411 || abs(mediator.id()) == 431 | |
240 | || abs(mediator.id()) == 511 || abs(mediator.id()) == 521 | |
241 | || abs(mediator.id()) == 531 || abs(mediator.id()) == 541 | |
242 | || (abs(mediator.id()) > 5100 && abs(mediator.id()) < 5600) ) | |
243 | && abs(out2.id()) == 16) { | |
244 | int idBmother = (mediator.id() > 0) ? -5 : 5; | |
245 | if (abs(mediator.id()) > 5100) idBmother = -idBmother; | |
246 | particles[0] = HelicityParticle( idBmother, 0, 0, 0, 0, 0, 0, 0, | |
247 | 0., 0., 0., 0., 0., 0., particleDataPtr); | |
248 | particles[1] = HelicityParticle( -idBmother, 0, 0, 0, 0, 0, 0, 0, | |
249 | 0., 0., 0., 0., 0., 0., particleDataPtr); | |
250 | particles[0].idx = 0; | |
251 | particles[1].idx = 1; | |
252 | ||
253 | // D or B meson decays into neutrino + tau + meson. | |
254 | if (mediator.daughter1() + 2 == mediator.daughter2()) { | |
255 | particles[0].p(mediator.p()); | |
256 | particles[1].direction = 1; | |
257 | particles[1].id(-particles[1].id()); | |
258 | particles[1].p(particles[0].p() - particles[2].p() - particles[3].p()); | |
259 | } | |
260 | ||
261 | // D or B meson decays into neutrino + tau. | |
262 | else { | |
263 | particles[0].p(mediator.p()/2); | |
264 | particles[1].p(mediator.p()/2); | |
265 | } | |
266 | hardME = hmeTwoFermions2W2TwoFermions.initChannel(particles); | |
267 | correlated = false; | |
268 | ||
269 | // Produced from a virtual photon with correlated taus. | |
270 | } else if (abs(out1.id()) == 15 && abs(out2.id()) == 15) { | |
271 | particles.push_back(mediator); | |
272 | particles[0] = HelicityParticle(-1, 0, 0, 0, 0, 0, 0, 0, | |
273 | mediator.p()/2, 0., 0., particleDataPtr); | |
274 | particles[1] = HelicityParticle(1, 0, 0, 0, 0, 0, 0, 0, | |
275 | mediator.p()/2, 0., 0., particleDataPtr); | |
276 | particles[0].direction = -1; | |
277 | particles[1].direction = -1; | |
278 | particles[0].idx = 0; | |
279 | particles[1].idx = 0; | |
280 | hardME = hmeTwoFermions2Gamma2TwoFermions.initChannel(particles); | |
281 | correlated = true; | |
282 | ||
283 | // Produced from an unknown process, assume unpolarized and uncorrelated. | |
284 | } else { | |
285 | if (tauMode <= 1) | |
286 | infoPtr->errorMsg("Warning in TauDecays::decay: unknown " | |
287 | "tau production, assuming unpolarized and uncorrelated"); | |
288 | hardME = hmeUnpolarized.initChannel(particles); | |
289 | correlated = false; | |
290 | } | |
291 | ||
292 | // Pick the first tau to decay. | |
293 | HelicityParticle* tau; | |
294 | int idx; | |
295 | if (correlated) idx = (rndmPtr->flat() < 0.5) ? 2 : 3; | |
296 | else idx = (abs(particles[2].id()) == 15) ? 2 : 3; | |
297 | tau = &particles[idx]; | |
298 | ||
299 | // Calculate the density matrix and decay the tau. | |
300 | if ( (tauMode == 2 && abs(mediator.id()) == tauMother) || tauMode == 3 ) { | |
301 | tau->rho[0][0] = (tau->id() > 0) | |
302 | ? (1 - polarization) / 2 : (1 + polarization) / 2; | |
303 | tau->rho[1][1] = (tau->id() > 0) | |
304 | ? (1 + polarization) / 2 : (1 - polarization) / 2; | |
305 | correlated = false; | |
306 | } | |
307 | ||
308 | // Begin PS April 2012. | |
309 | // Else use tau helicity provided by event record (LHEF). | |
310 | else if (hasHelicity) { | |
311 | tau->rho[0][0] = (1. - helicityNow) / 2.; | |
312 | tau->rho[1][1] = (1. + helicityNow) / 2.; | |
313 | } | |
314 | // End PS April 2012. | |
315 | ||
316 | // Else compute density matrix according to matrix element. | |
317 | else | |
318 | hardME->calculateRho(idx, particles); | |
319 | vector<HelicityParticle> children = createChildren(*tau); | |
320 | if (children.size() == 0) return false; | |
321 | ||
322 | // Decay the first tau. | |
323 | bool accepted = false; | |
324 | int tries = 0; | |
325 | while (!accepted) { | |
326 | isotropicDecay(children); | |
327 | double decayWeight = decayME->decayWeight(children); | |
328 | double decayWeightMax = decayME->decayWeightMax(children); | |
329 | accepted = (rndmPtr->flat() < decayWeight / decayWeightMax); | |
330 | if (decayWeight > decayWeightMax) | |
331 | infoPtr->errorMsg("Warning in TauDecays::decay: maximum " | |
332 | "decay weight exceeded in tau decay"); | |
333 | if (tries > NTRYDECAY) { | |
334 | infoPtr->errorMsg("Warning in TauDecays::decay: maximum " | |
335 | "number of decay attempts exceeded"); | |
336 | break; | |
337 | } | |
338 | ++tries; | |
339 | } | |
340 | writeDecay(event,children); | |
341 | ||
342 | // If a correlated second tau exists, decay that tau as well. | |
343 | if (correlated) { | |
344 | idx = (idx == 2) ? 3 : 2; | |
345 | // Calculate the first tau decay matrix. | |
346 | decayME->calculateD(children); | |
347 | // Update the decay matrix for the tau. | |
348 | (*tau).D = children[0].D; | |
349 | // Switch the taus. | |
350 | tau = &particles[idx]; | |
351 | // Calculate second tau's density matrix. | |
352 | hardME->calculateRho(idx, particles); | |
353 | ||
354 | // Decay the second tau. | |
355 | children.clear(); | |
356 | children = createChildren(*tau); | |
357 | if (children.size() == 0) return false; | |
358 | accepted = false; | |
359 | tries = 0; | |
360 | while (!accepted) { | |
361 | isotropicDecay(children); | |
362 | double decayWeight = decayME->decayWeight(children); | |
363 | double decayWeightMax = decayME->decayWeightMax(children); | |
364 | accepted = (rndmPtr->flat() < decayWeight / decayWeightMax); | |
365 | if (decayWeight > decayWeightMax) | |
366 | infoPtr->errorMsg("Warning in TauDecays::decay: maximum " | |
367 | "decay weight exceeded in correlated tau decay"); | |
368 | if (tries > NTRYDECAY) { | |
369 | infoPtr->errorMsg("Warning in TauDecays::decay: maximum " | |
370 | "number of decay attempts exceeded"); | |
371 | break; | |
372 | } | |
373 | ++tries; | |
374 | } | |
375 | writeDecay(event,children); | |
376 | } | |
377 | ||
378 | // Done. | |
379 | return true; | |
380 | ||
381 | } | |
382 | ||
383 | //-------------------------------------------------------------------------- | |
384 | ||
385 | // Given a HelicityParticle parent, select the decay channel and return | |
386 | // a vector of HelicityParticles containing the children, with the parent | |
387 | // particle duplicated in the first entry of the vector. | |
388 | ||
389 | vector<HelicityParticle> TauDecays::createChildren(HelicityParticle parent) { | |
390 | ||
391 | // Initial values. | |
392 | int meMode = 0; | |
393 | vector<HelicityParticle> children; | |
394 | ||
395 | // Set the parent as incoming. | |
396 | parent.direction = -1; | |
397 | ||
398 | // Setup decay data for the decaying particle. | |
399 | ParticleDataEntry decayData = parent.particleDataEntry(); | |
400 | ||
401 | // Initialize the decay data. | |
402 | if (!decayData.preparePick(parent.id())) return children; | |
403 | ||
404 | // Try to pick a decay channel. | |
405 | bool decayed = false; | |
406 | int decayTries = 0; | |
407 | while (!decayed && decayTries < NTRYCHANNEL) { | |
408 | ||
409 | // Pick a decay channel. | |
410 | DecayChannel decayChannel = decayData.pickChannel(); | |
411 | meMode = decayChannel.meMode(); | |
412 | int decayMult = decayChannel.multiplicity(); | |
413 | ||
414 | // Select children masses. | |
415 | bool allowed = false; | |
416 | int channelTries = 0; | |
417 | while (!allowed && channelTries < NTRYCHANNEL) { | |
418 | children.resize(0); | |
419 | children.push_back(parent); | |
420 | for (int i = 0; i < decayMult; ++i) { | |
421 | // Grab child ID. | |
422 | int childId = decayChannel.product(i); | |
423 | // Flip sign for anti-particle decay. | |
424 | if (parent.id() < 0 && particleDataPtr->hasAnti(childId)) | |
425 | childId = -childId; | |
426 | // Grab child mass. | |
427 | double childMass = particleDataPtr->mass(childId); | |
428 | // Push back the child into the children vector. | |
429 | children.push_back( HelicityParticle(childId, 91, parent.idx, 0, 0, | |
430 | 0, 0, 0, 0., 0., 0., 0.,childMass, 0., particleDataPtr) ); | |
431 | } | |
432 | ||
433 | // Check there is enough phase space for decay. | |
434 | if (decayMult > 1) { | |
435 | double massDiff = parent.m(); | |
436 | for (int i = 0; i < decayMult; ++i) massDiff -= children[i].m(); | |
437 | // For now we just check kinematically available. | |
438 | if (massDiff > 0) { | |
439 | allowed = true; | |
440 | decayed = true; | |
441 | } | |
442 | } | |
443 | ||
444 | // End pick a channel. | |
445 | ++channelTries; | |
446 | } | |
447 | ++decayTries; | |
448 | } | |
449 | ||
450 | // Set the decay matrix element. | |
451 | // Two body decays. | |
452 | if (children.size() == 3) { | |
453 | if (meMode == 1521) | |
454 | decayME = hmeTau2Meson.initChannel(children); | |
455 | else decayME = hmeTau2PhaseSpace.initChannel(children); | |
456 | } | |
457 | ||
458 | // Three body decays. | |
459 | else if (children.size() == 4) { | |
460 | // Leptonic decay. | |
461 | if (meMode == 1531) | |
462 | decayME = hmeTau2TwoLeptons.initChannel(children); | |
463 | // Two meson decay via vector meson. | |
464 | else if (meMode == 1532) | |
465 | decayME = hmeTau2TwoMesonsViaVector.initChannel(children); | |
466 | // Two meson decay via vector or scalar meson. | |
467 | else if (meMode == 1533) | |
468 | decayME = hmeTau2TwoMesonsViaVectorScalar.initChannel(children); | |
469 | // Flat phase space. | |
470 | else decayME = hmeTau2PhaseSpace.initChannel(children); | |
471 | } | |
472 | ||
473 | // Four body decays. | |
474 | else if (children.size() == 5) { | |
475 | // Three pion CLEO decay. | |
476 | if (meMode == 1541) | |
477 | decayME = hmeTau2ThreePions.initChannel(children); | |
478 | // Three meson decay with one or more kaons decay. | |
479 | else if (meMode == 1542) | |
480 | decayME = hmeTau2ThreeMesonsWithKaons.initChannel(children); | |
481 | // Generic three meson decay. | |
482 | else if (meMode == 1543) | |
483 | decayME = hmeTau2ThreeMesonsGeneric.initChannel(children); | |
484 | // Two pions and photon decay. | |
485 | else if (meMode == 1544) | |
486 | decayME = hmeTau2TwoPionsGamma.initChannel(children); | |
487 | // Flat phase space. | |
488 | else decayME = hmeTau2PhaseSpace.initChannel(children); | |
489 | } | |
490 | ||
491 | // Five body decays. | |
492 | else if (children.size() == 6) { | |
493 | // Four pion Novosibirsk current. | |
494 | if (meMode == 1551) | |
495 | decayME = hmeTau2FourPions.initChannel(children); | |
496 | // Flat phase space. | |
497 | else decayME = hmeTau2PhaseSpace.initChannel(children); | |
498 | } | |
499 | ||
500 | // Six body decays. | |
501 | else if (children.size() == 7) { | |
502 | // Four pion Novosibirsk current. | |
503 | if (meMode == 1561) | |
504 | decayME = hmeTau2FivePions.initChannel(children); | |
505 | // Flat phase space. | |
506 | else decayME = hmeTau2PhaseSpace.initChannel(children); | |
507 | } | |
508 | ||
509 | // Flat phase space. | |
510 | else decayME = hmeTau2PhaseSpace.initChannel(children); | |
511 | ||
512 | // Done. | |
513 | return children; | |
514 | } | |
515 | ||
516 | //-------------------------------------------------------------------------- | |
517 | ||
518 | // N-body decay using the M-generator algorithm described in | |
519 | // "Monte Carlo Phase Space" by F. James in CERN 68-15, May 1968. Taken | |
520 | // from ParticleDecays::mGenerator but modified to handle spin particles. | |
521 | // Given a vector of HelicityParticles where the first particle is | |
522 | // the mother, the remaining particles are decayed isotropically. | |
523 | ||
524 | void TauDecays::isotropicDecay(vector<HelicityParticle>& children) { | |
525 | ||
526 | // Mother and sum daughter masses. | |
527 | int decayMult = children.size() - 1; | |
528 | double m0 = children[0].m(); | |
529 | double mSum = children[1].m(); | |
530 | for (int i = 2; i <= decayMult; ++i) mSum += children[i].m(); | |
531 | double mDiff = m0 - mSum; | |
532 | ||
533 | // Begin setup of intermediate invariant masses. | |
534 | vector<double> mInv; | |
535 | for (int i = 0; i <= decayMult; ++i) mInv.push_back( children[i].m()); | |
536 | ||
537 | // Calculate the maximum weight in the decay. | |
538 | double wtPS; | |
539 | double wtPSmax = 1. / WTCORRECTION[decayMult]; | |
540 | double mMax = mDiff + children[decayMult].m(); | |
541 | double mMin = 0.; | |
542 | for (int i = decayMult - 1; i > 0; --i) { | |
543 | mMax += children[i].m(); | |
544 | mMin += children[i+1].m(); | |
545 | double mNow = children[i].m(); | |
546 | wtPSmax *= 0.5 * sqrtpos( (mMax - mMin - mNow) * (mMax + mMin + mNow) | |
547 | * (mMax + mMin - mNow) * (mMax - mMin + mNow) ) / mMax; | |
548 | } | |
549 | ||
550 | // Begin loop to find the set of intermediate invariant masses. | |
551 | vector<double> rndmOrd; | |
552 | do { | |
553 | wtPS = 1.; | |
554 | ||
555 | // Find and order random numbers in descending order. | |
556 | rndmOrd.clear(); | |
557 | rndmOrd.push_back(1.); | |
558 | for (int i = 1; i < decayMult - 1; ++i) { | |
559 | double random = rndmPtr->flat(); | |
560 | rndmOrd.push_back(random); | |
561 | for (int j = i - 1; j > 0; --j) { | |
562 | if (random > rndmOrd[j]) swap( rndmOrd[j], rndmOrd[j+1] ); | |
563 | else break; | |
564 | } | |
565 | } | |
566 | rndmOrd.push_back(0.); | |
567 | ||
568 | // Translate into intermediate masses and find weight. | |
569 | for (int i = decayMult - 1; i > 0; --i) { | |
570 | mInv[i] = mInv[i+1] + children[i].m() | |
571 | + (rndmOrd[i-1] - rndmOrd[i]) * mDiff; | |
572 | wtPS *= 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - children[i].m()) | |
573 | * (mInv[i] + mInv[i+1] + children[i].m()) | |
574 | * (mInv[i] + mInv[i+1] - children[i].m()) | |
575 | * (mInv[i] - mInv[i+1] + children[i].m()) ) / mInv[i]; | |
576 | } | |
577 | ||
578 | // If rejected, try again with new invariant masses. | |
579 | } while ( wtPS < rndmPtr->flat() * wtPSmax ); | |
580 | ||
581 | // Perform two-particle decays in the respective rest frame. | |
582 | vector<Vec4> pInv(decayMult + 1); | |
583 | for (int i = 1; i < decayMult; ++i) { | |
584 | double pAbs = 0.5 * sqrtpos( (mInv[i] - mInv[i+1] - children[i].m()) | |
585 | * (mInv[i] + mInv[i+1] + children[i].m()) | |
586 | * (mInv[i] + mInv[i+1] - children[i].m()) | |
587 | * (mInv[i] - mInv[i+1] + children[i].m()) ) / mInv[i]; | |
588 | ||
589 | // Isotropic angles give three-momentum. | |
590 | double cosTheta = 2. * rndmPtr->flat() - 1.; | |
591 | double sinTheta = sqrt(1. - cosTheta*cosTheta); | |
592 | double phi = 2. * M_PI * rndmPtr->flat(); | |
593 | double pX = pAbs * sinTheta * cos(phi); | |
594 | double pY = pAbs * sinTheta * sin(phi); | |
595 | double pZ = pAbs * cosTheta; | |
596 | ||
597 | // Calculate energies, fill four-momenta. | |
598 | double eHad = sqrt( children[i].m()*children[i].m() + pAbs*pAbs); | |
599 | double eInv = sqrt( mInv[i+1]*mInv[i+1] + pAbs*pAbs); | |
600 | children[i].p( pX, pY, pZ, eHad); | |
601 | pInv[i+1].p( -pX, -pY, -pZ, eInv); | |
602 | } | |
603 | ||
604 | // Boost decay products to the mother rest frame. | |
605 | children[decayMult].p( pInv[decayMult] ); | |
606 | for (int iFrame = decayMult - 1; iFrame > 1; --iFrame) | |
607 | for (int i = iFrame; i <= decayMult; ++i) | |
608 | children[i].bst( pInv[iFrame], mInv[iFrame]); | |
609 | ||
610 | // Boost decay products to the current frame. | |
611 | pInv[1].p( children[0].p() ); | |
612 | for (int i = 1; i <= decayMult; ++i) children[i].bst( pInv[1], mInv[1] ); | |
613 | ||
614 | // Done. | |
615 | return; | |
616 | } | |
617 | ||
618 | //-------------------------------------------------------------------------- | |
619 | ||
620 | // Write the vector of HelicityParticles to the event record, excluding the | |
621 | // first particle. Set the lifetime and production vertex of the particles | |
622 | // and mark the first particle of the vector as decayed. | |
623 | ||
624 | void TauDecays::writeDecay(Event& event, vector<HelicityParticle>& children) { | |
625 | ||
626 | // Set additional information and append children to event. | |
627 | int decayMult = children.size() - 1; | |
628 | Vec4 decayVertex = children[0].vDec(); | |
629 | for (int i = 1; i <= decayMult; i++) { | |
630 | // Set child lifetime. | |
631 | children[i].tau(children[i].tau0() * rndmPtr->exp()); | |
632 | // Set child production vertex. | |
633 | children[i].vProd(decayVertex); | |
634 | // Append child to record. | |
635 | children[i].idx = event.append(children[i]); | |
636 | } | |
637 | ||
638 | // Mark the parent as decayed and set children. | |
639 | event[children[0].idx].statusNeg(); | |
640 | event[children[0].idx].daughters(children[1].idx, children[decayMult].idx); | |
641 | ||
642 | } | |
643 | ||
644 | //========================================================================== | |
645 | ||
646 | } // end namespace Pythia8 |