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 |
---|