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