Line data Source code
1 : // HelicityMatrixElements.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 physics classes
7 : // used in tau decays.
8 :
9 : #include "Pythia8/HelicityMatrixElements.h"
10 :
11 : namespace Pythia8 {
12 :
13 : //==========================================================================
14 :
15 : // The HelicityMatrixElements class.
16 :
17 : //--------------------------------------------------------------------------
18 :
19 : // Initialize the helicity matrix element.
20 :
21 : void HelicityMatrixElement::initPointers(ParticleData* particleDataPtrIn,
22 : Couplings* couplingsPtrIn, Settings* settingsPtrIn) {
23 :
24 0 : particleDataPtr = particleDataPtrIn;
25 0 : couplingsPtr = couplingsPtrIn;
26 0 : settingsPtr = settingsPtrIn;
27 0 : for(int i = 0; i <= 5; i++)
28 0 : gamma.push_back(GammaMatrix(i));
29 :
30 0 : }
31 :
32 : //--------------------------------------------------------------------------
33 :
34 : // Initialize the channel for the helicity matrix element.
35 :
36 : HelicityMatrixElement* HelicityMatrixElement::initChannel(
37 : vector<HelicityParticle>& p) {
38 :
39 0 : pID.clear();
40 0 : pM.clear();
41 0 : for(int i = 0; i < static_cast<int>(p.size()); i++) {
42 0 : pID.push_back(p[i].id());
43 0 : pM.push_back(p[i].m());
44 : }
45 0 : initConstants();
46 0 : return this;
47 :
48 : }
49 :
50 : //--------------------------------------------------------------------------
51 :
52 : // Calculate a particle's decay matrix.
53 :
54 : void HelicityMatrixElement::calculateD(vector<HelicityParticle>& p) {
55 :
56 : // Reset the D matrix to zero.
57 0 : for (int i = 0; i < p[0].spinStates(); i++) {
58 0 : for (int j = 0; j < p[0].spinStates(); j++) {
59 0 : p[0].D[i][j] = 0;
60 : }
61 : }
62 :
63 : // Initialize the wave functions.
64 0 : initWaves(p);
65 :
66 : // Create the helicity vectors.
67 0 : vector<int> h1(p.size(),0);
68 0 : vector<int> h2(p.size(),0);
69 :
70 : // Call the recursive sub-method.
71 0 : calculateD(p, h1, h2, 0);
72 :
73 : // Normalize the decay matrix.
74 0 : p[0].normalize(p[0].D);
75 :
76 0 : }
77 :
78 : //--------------------------------------------------------------------------
79 :
80 : // Recursive sub-method for calculating a particle's decay matrix.
81 :
82 : void HelicityMatrixElement::calculateD(vector<HelicityParticle>& p,
83 : vector<int>& h1, vector<int>& h2, unsigned int i) {
84 :
85 0 : if (i < p.size()) {
86 0 : for (h1[i] = 0; h1[i] < p[i].spinStates(); h1[i]++) {
87 0 : for (h2[i] = 0; h2[i] < p[i].spinStates(); h2[i]++) {
88 0 : calculateD(p, h1, h2, i+1);
89 : }
90 : }
91 : }
92 : else {
93 0 : p[0].D[h1[0]][h2[0]] += calculateME(h1) * conj(calculateME(h2)) *
94 0 : calculateProductD(p, h1, h2);
95 : }
96 :
97 0 : }
98 :
99 : //--------------------------------------------------------------------------
100 :
101 : // Calculate a particle's helicity density matrix.
102 :
103 : void HelicityMatrixElement::calculateRho(unsigned int idx,
104 : vector<HelicityParticle>& p) {
105 :
106 : // Reset the rho matrix to zero.
107 0 : for (int i = 0; i < p[idx].spinStates(); i++) {
108 0 : for (int j = 0; j < p[idx].spinStates(); j++) {
109 0 : p[idx].rho[i][j] = 0;
110 : }
111 : }
112 :
113 : // Initialize the wave functions.
114 0 : initWaves(p);
115 :
116 : // Create the helicity vectors.
117 0 : vector<int> h1(p.size(),0);
118 0 : vector<int> h2(p.size(),0);
119 :
120 : // Call the recursive sub-method.
121 0 : calculateRho(idx, p, h1, h2, 0);
122 :
123 : // Normalize the density matrix.
124 0 : p[idx].normalize(p[idx].rho);
125 :
126 0 : }
127 :
128 : //--------------------------------------------------------------------------
129 :
130 : // Recursive sub-method for calculating a particle's helicity density matrix.
131 :
132 : void HelicityMatrixElement::calculateRho(unsigned int idx,
133 : vector<HelicityParticle>& p, vector<int>& h1, vector<int>& h2,
134 : unsigned int i) {
135 :
136 0 : if (i < p.size()) {
137 0 : for (h1[i] = 0; h1[i] < p[i].spinStates(); h1[i]++) {
138 0 : for (h2[i] = 0; h2[i] < p[i].spinStates(); h2[i]++) {
139 0 : calculateRho(idx, p, h1, h2, i+1);
140 : }
141 : }
142 : }
143 : else {
144 : // Calculate rho from a hard process.
145 0 : if (p[1].direction < 0)
146 0 : p[idx].rho[h1[idx]][h2[idx]] += p[0].rho[h1[0]][h2[0]] *
147 0 : p[1].rho[h1[1]][h2[1]] * calculateME(h1)*conj(calculateME(h2)) *
148 0 : calculateProductD(idx, 2, p, h1, h2);
149 : // Calculate rho from a decay.
150 : else
151 0 : p[idx].rho[h1[idx]][h2[idx]] += p[0].rho[h1[0]][h2[0]] *
152 0 : calculateME(h1)*conj(calculateME(h2)) *
153 0 : calculateProductD(idx, 1, p, h1, h2);
154 : return;
155 : }
156 :
157 0 : }
158 :
159 : //--------------------------------------------------------------------------
160 :
161 : // Calculate a decay's weight.
162 :
163 : double HelicityMatrixElement::decayWeight(vector<HelicityParticle>& p) {
164 :
165 0 : complex weight = complex(0,0);
166 :
167 : // Initialize the wave functions.
168 0 : initWaves(p);
169 :
170 : // Create the helicity vectors.
171 0 : vector<int> h1(p.size(),0);
172 0 : vector<int> h2(p.size(),0);
173 :
174 : // Call the recursive sub-method.
175 0 : decayWeight(p, h1, h2, weight, 0);
176 :
177 0 : return real(weight);
178 :
179 0 : }
180 :
181 : //--------------------------------------------------------------------------
182 :
183 : // Recursive sub-method for calculating a decay's weight.
184 :
185 : void HelicityMatrixElement::decayWeight(vector<HelicityParticle>& p,
186 : vector<int>& h1, vector<int>& h2, complex& weight, unsigned int i) {
187 :
188 0 : if (i < p.size()) {
189 0 : for (h1[i] = 0; h1[i] < p[i].spinStates(); h1[i]++) {
190 0 : for (h2[i] = 0; h2[i] < p[i].spinStates(); h2[i]++) {
191 0 : decayWeight(p, h1, h2, weight, i+1);
192 : }
193 : }
194 : }
195 : else {
196 0 : weight += p[0].rho[h1[0]][h2[0]] * calculateME(h1) *
197 0 : conj(calculateME(h2)) * calculateProductD(p, h1, h2);
198 : }
199 :
200 0 : }
201 :
202 : //--------------------------------------------------------------------------
203 :
204 : // Calculate the product of the decay matrices (hard process).
205 :
206 : complex HelicityMatrixElement::calculateProductD(unsigned int idx,
207 : unsigned int start, vector<HelicityParticle>& p,
208 : vector<int>& h1, vector<int>& h2) {
209 :
210 0 : complex answer(1,0);
211 0 : for (unsigned int i = start; i < p.size(); i++) {
212 0 : if (i != idx) {
213 0 : answer *= p[i].D[h1[i]][h2[i]];
214 0 : }
215 : }
216 0 : return answer;
217 :
218 : }
219 :
220 : //--------------------------------------------------------------------------
221 :
222 : // Calculate the product of the decay matrices (decay process).
223 :
224 : complex HelicityMatrixElement::calculateProductD(
225 : vector<HelicityParticle>& p, vector<int>& h1, vector<int>& h2) {
226 :
227 0 : complex answer(1,0);
228 0 : for (unsigned int i = 1; i < p.size(); i++) {
229 0 : answer *= p[i].D[h1[i]][h2[i]];
230 : }
231 0 : return answer;
232 :
233 : }
234 :
235 : //--------------------------------------------------------------------------
236 :
237 : // Initialize a fermion line.
238 :
239 : void HelicityMatrixElement::setFermionLine(int position,
240 : HelicityParticle& p0, HelicityParticle& p1) {
241 :
242 0 : vector< Wave4 > u0, u1;
243 :
244 : // First particle is incoming and particle, or outgoing and anti-particle.
245 0 : if (p0.id()*p0.direction < 0) {
246 0 : pMap[position] = position; pMap[position+1] = position+1;
247 0 : for (int h = 0; h < p0.spinStates(); h++) u0.push_back(p0.wave(h));
248 0 : for (int h = 0; h < p1.spinStates(); h++) u1.push_back(p1.waveBar(h));
249 0 : }
250 : // First particle is outgoing and particle, or incoming and anti-particle.
251 : else {
252 0 : pMap[position] = position+1; pMap[position+1] = position;
253 0 : for (int h = 0; h < p0.spinStates(); h++) u1.push_back(p0.waveBar(h));
254 0 : for (int h = 0; h < p1.spinStates(); h++) u0.push_back(p1.wave(h));
255 : }
256 0 : u.push_back(u0); u.push_back(u1);
257 :
258 0 : }
259 :
260 : //--------------------------------------------------------------------------
261 :
262 : // Return a fixed width Breit-Wigner.
263 :
264 : complex HelicityMatrixElement::breitWigner(double s, double M, double G) {
265 :
266 0 : return (-M*M + complex(0, 1) * M * G) / (s - M*M + complex(0, 1) * M * G);
267 :
268 : }
269 :
270 : //--------------------------------------------------------------------------
271 :
272 : // Return an s-wave BreitWigner.
273 :
274 : complex HelicityMatrixElement::sBreitWigner(double m0, double m1, double s,
275 : double M, double G) {
276 :
277 0 : double gs = sqrtpos((s - pow2(m0+m1)) * (s - pow2(m0-m1))) / (2*sqrtpos(s));
278 0 : double gM = sqrtpos((M*M - pow2(m0+m1)) * (M*M - pow2(m0-m1))) / (2*M);
279 0 : return M*M / (M*M - s - complex(0,1)*G*M*M/sqrtpos(s)*(gs/gM));
280 :
281 0 : }
282 :
283 : //--------------------------------------------------------------------------
284 :
285 : // Return a p-wave BreitWigner.
286 :
287 : complex HelicityMatrixElement::pBreitWigner(double m0, double m1, double s,
288 : double M, double G) {
289 :
290 0 : double gs = sqrtpos((s - pow2(m0+m1)) * (s - pow2(m0-m1))) / (2*sqrtpos(s));
291 0 : double gM = sqrtpos((M*M - pow2(m0+m1)) * (M*M - pow2(m0-m1))) / (2*M);
292 0 : return M*M / (M*M - s - complex(0,1)*G*M*M/sqrtpos(s)*pow3(gs/gM));
293 :
294 0 : }
295 :
296 : //--------------------------------------------------------------------------
297 :
298 : // Return a d-wave BreitWigner.
299 :
300 : complex HelicityMatrixElement::dBreitWigner(double m0, double m1, double s,
301 : double M, double G) {
302 :
303 0 : double gs = sqrtpos((s - pow2(m0+m1)) * (s - pow2(m0-m1))) / (2*sqrtpos(s));
304 0 : double gM = sqrtpos((M*M - pow2(m0+m1)) * (M*M - pow2(m0-m1))) / (2*M);
305 0 : return M*M / (M*M - s - complex(0,1)*G*M*M/sqrtpos(s)*pow5(gs/gM));
306 :
307 0 : }
308 :
309 : //==========================================================================
310 :
311 : // Helicity matrix element for two fermions -> W/W' -> two
312 : // fermions. This matrix element handles s-channel hard processes in
313 : // addition to t-channel, assuming the first two particles are a
314 : // fermion line and the second two particles are a fermion line. This
315 : // matrix element is not scaled with respect to W/W' propagator energy as
316 : // currently this matrix element is used only for calculating helicity
317 : // density matrices.
318 :
319 : //--------------------------------------------------------------------------
320 :
321 : // Initialize the constants for the helicity matrix element.
322 :
323 : void HMETwoFermions2W2TwoFermions::initConstants() {
324 :
325 : // Set the constants for the W'.
326 0 : if (abs(pID[4]) == 34 && settingsPtr) {
327 0 : if (abs(pID[0]) < 11) {
328 0 : p0CA = settingsPtr->parm("Wprime:aq");
329 0 : p0CV = settingsPtr->parm("Wprime:vq");
330 0 : } else {
331 0 : p0CA = settingsPtr->parm("Wprime:al");
332 0 : p0CV = settingsPtr->parm("Wprime:vl");
333 : }
334 0 : if (abs(pID[2]) < 11) {
335 0 : p2CA = settingsPtr->parm("Wprime:aq");
336 0 : p2CV = settingsPtr->parm("Wprime:vq");
337 0 : } else {
338 0 : p2CA = settingsPtr->parm("Wprime:al");
339 0 : p2CV = settingsPtr->parm("Wprime:vl");
340 : }
341 :
342 : // The default constants (SM W).
343 : } else {
344 0 : p0CA = -1; p0CV = 1;
345 0 : p2CA = -1; p2CV = 1;
346 : }
347 :
348 0 : }
349 :
350 : //--------------------------------------------------------------------------
351 :
352 : // Initialize spinors for the helicity matrix element.
353 :
354 : void HMETwoFermions2W2TwoFermions::initWaves(vector<HelicityParticle>& p) {
355 :
356 0 : u.clear();
357 0 : pMap.resize(4);
358 0 : setFermionLine(0,p[0],p[1]);
359 0 : setFermionLine(2,p[2],p[3]);
360 :
361 0 : }
362 :
363 : //--------------------------------------------------------------------------
364 :
365 : // Return element for the helicity matrix element.
366 :
367 : complex HMETwoFermions2W2TwoFermions::calculateME(vector<int> h) {
368 :
369 0 : complex answer(0,0);
370 0 : for (int mu = 0; mu <= 3; mu++) {
371 0 : answer += (u[1][h[pMap[1]]] * gamma[mu] * (p0CV + p0CA * gamma[5])
372 0 : * u[0][h[pMap[0]]]) * gamma[4](mu,mu) * (u[3][h[pMap[3]]]
373 0 : * gamma[mu] * (p2CV + p2CA * gamma[5]) * u[2][h[pMap[2]]]);
374 : }
375 0 : return answer;
376 :
377 0 : }
378 :
379 : //==========================================================================
380 :
381 : // Helicity matrix element for two fermions -> photon/Z/Z' -> two fermions.
382 :
383 : // Note that there is a double contraction in the Z matrix element, which can
384 : // be very time consuming. If the two incoming fermions are oriented along
385 : // the z-axis, their helicities must be opposite for a non-zero matrix element
386 : // term. Consequently, this check is made to help speed up the matrix element.
387 :
388 : // p0Q: charge of the incoming fermion line
389 : // p2Q: charge of the outgoing fermion line
390 : // s: center of mass energy
391 : // sin2W: sine of the Weinberg angle
392 : // cos2W: cosine of the Weinberg angle
393 : // zM: on-shell mass of the Z
394 : // zG: on-shell width of the Z
395 : // p0CA: axial coupling of particle 0 to the Z
396 : // p2CA: axial coupling of particle 2 to the Z
397 : // p0CV: vector coupling of particle 0 to the Z
398 : // p2CV: vector coupling of particle 2 to the Z
399 : // zaxis: true if the incoming fermions are oriented along the z-axis
400 :
401 : //--------------------------------------------------------------------------
402 :
403 : // Initialize the constants for the helicity matrix element.
404 :
405 : void HMETwoFermions2GammaZ2TwoFermions::initConstants() {
406 :
407 : // Set the Weinberg angle.
408 0 : sin2W = couplingsPtr->sin2thetaW();
409 0 : cos2W = couplingsPtr->cos2thetaW();
410 : // Set the on-shell Z/Z' mass and width.
411 0 : zG = particleDataPtr->mWidth(23);
412 0 : zM = particleDataPtr->m0(23);
413 0 : zpG = particleDataPtr->mWidth(32);
414 0 : zpM = particleDataPtr->m0(32);
415 : // Set the Z vector and axial couplings to the fermions.
416 0 : p0CAZ = couplingsPtr->af(abs(pID[0]));
417 0 : p2CAZ = couplingsPtr->af(abs(pID[2]));
418 0 : p0CVZ = couplingsPtr->vf(abs(pID[0]));
419 0 : p2CVZ = couplingsPtr->vf(abs(pID[2]));
420 : // Turn off the gamma/Z/Z' channels.
421 0 : includeGamma = false;
422 0 : includeZ = false;
423 0 : includeZp = false;
424 0 : if (settingsPtr) {
425 : // Set the Z' vector and axial couplings to the fermions.
426 0 : p0CAZp = zpCoupling(pID[0], "a");
427 0 : p0CVZp = zpCoupling(pID[0], "v");
428 0 : p2CAZp = zpCoupling(pID[2], "a");
429 0 : p2CVZp = zpCoupling(pID[2], "v");
430 : // Set the gamma/Z/Z' channels to include.
431 0 : if (abs(pID[4]) == 22) {
432 0 : includeGamma = true;
433 0 : } else if (abs(pID[4]) == 23) {
434 0 : int mode = settingsPtr->mode("WeakZ0:gmZmode");
435 0 : if (mode == 0) {includeGamma = true; includeZ = true;}
436 0 : else if (mode == 1) includeGamma = true;
437 0 : else if (mode == 2) includeZ = true;
438 0 : } else if (abs(pID[4]) == 32) {
439 0 : int mode = settingsPtr->mode("Zprime:gmZmode");
440 0 : if (mode == 0) {includeGamma = true; includeZ = true; includeZp = true;}
441 0 : else if (mode == 1) includeGamma = true;
442 0 : else if (mode == 2) includeZ = true;
443 0 : else if (mode == 3) includeZp = true;
444 0 : else if (mode == 4) {includeGamma = true; includeZ = true;}
445 0 : else if (mode == 5) {includeGamma = true; includeZp = true;}
446 0 : else if (mode == 6) {includeZ = true; includeZp = true;}
447 0 : }
448 : // Default behavior without settings pointer.
449 : } else {
450 0 : p0CAZp = p0CAZ;
451 0 : p0CVZp = p2CAZ;
452 0 : p2CAZp = p0CVZ;
453 0 : p2CVZp = p2CVZ;
454 0 : if (abs(pID[4]) == 22) includeGamma = true;
455 0 : else if (abs(pID[4]) == 23) includeZ = true;
456 0 : else if (abs(pID[4]) == 32) includeZp = true;
457 : }
458 :
459 0 : }
460 :
461 : //--------------------------------------------------------------------------
462 :
463 : // Initialize wave functions for the helicity matrix element.
464 :
465 : void HMETwoFermions2GammaZ2TwoFermions::initWaves(vector<HelicityParticle>& p)
466 : {
467 :
468 0 : vector< Wave4 > u4;
469 0 : u.clear();
470 0 : pMap.resize(4);
471 0 : setFermionLine(0, p[0], p[1]);
472 0 : setFermionLine(2, p[2], p[3]);
473 0 : u4.push_back(Wave4(p[2].p() + p[3].p()));
474 0 : u.push_back(u4);
475 : // Fermion line charge.
476 0 : p0Q = p[0].charge(); p2Q = p[2].charge();
477 : // Center of mass energy.
478 0 : s = max( 1., pow2(p[4].m()));
479 : // Check if incoming fermions are oriented along z-axis.
480 0 : zaxis = (p[0].pAbs() == fabs(p[0].pz())) &&
481 0 : (p[1].pAbs() == fabs(p[1].pz()));
482 :
483 0 : }
484 :
485 : //--------------------------------------------------------------------------
486 :
487 : // Return element for the helicity matrix element.
488 :
489 : complex HMETwoFermions2GammaZ2TwoFermions::calculateME(vector<int> h) {
490 :
491 0 : complex answer(0,0);
492 0 : if (includeGamma)
493 0 : answer += calculateGammaME(h);
494 0 : if (includeZ)
495 0 : answer += calculateZME(h, zM, zG, p0CAZ, p2CAZ, p0CVZ, p2CVZ);
496 0 : if (includeZp)
497 0 : answer += calculateZME(h, zpM, zpG, p0CAZp, p2CAZp, p0CVZp, p2CVZp);
498 0 : return answer;
499 :
500 0 : }
501 :
502 : //--------------------------------------------------------------------------
503 :
504 : // Return gamma element for the helicity matrix element.
505 :
506 : complex HMETwoFermions2GammaZ2TwoFermions::calculateGammaME(vector<int> h) {
507 :
508 0 : complex answer(0,0);
509 0 : for (int mu = 0; mu <= 3; mu++) {
510 0 : answer += (u[1][h[pMap[1]]] * gamma[mu] * u[0][h[pMap[0]]])
511 0 : * gamma[4](mu,mu) * (u[3][h[pMap[3]]] * gamma[mu] * u[2][h[pMap[2]]]);
512 : }
513 0 : return p0Q*p2Q * answer / s;
514 :
515 0 : }
516 :
517 : //--------------------------------------------------------------------------
518 :
519 : // Return Z/Z' element for helicity matrix element.
520 :
521 : complex HMETwoFermions2GammaZ2TwoFermions::calculateZME(
522 : vector<int> h, double m, double g, double p0CA, double p2CA, double p0CV,
523 : double p2CV) {
524 :
525 0 : complex answer(0,0);
526 : // Return zero if correct helicity conditions.
527 0 : if (h[0] == h[1] && zaxis) return answer;
528 0 : for (int mu = 0; mu <= 3; mu++) {
529 0 : for (int nu = 0; nu <= 3; nu++) {
530 0 : answer +=
531 0 : (u[1][h[pMap[1]]] * gamma[mu] * (p0CV - p0CA * gamma[5]) *
532 0 : u[0][h[pMap[0]]]) *
533 0 : (gamma[4](mu,nu) - gamma[4](mu,mu)*u[4][0](mu) *
534 0 : gamma[4](nu,nu) * u[4][0](nu) / (zM*zM)) *
535 0 : (u[3][h[pMap[3]]] * gamma[nu] * (p2CV - p2CA * gamma[5]) *
536 0 : u[2][h[pMap[2]]]);
537 : }
538 : }
539 0 : return answer / (16 * pow2(sin2W * cos2W) * (s - m*m + complex(0, s*g/m)));
540 :
541 0 : }
542 :
543 : //--------------------------------------------------------------------------
544 :
545 : // Return the Z' vector or axial coupling for a fermion.
546 :
547 : double HMETwoFermions2GammaZ2TwoFermions::zpCoupling(int id, string type) {
548 :
549 0 : if (!settingsPtr) return 0;
550 0 : id = abs(id);
551 0 : string name;
552 0 : if (id == 1) name = "d";
553 0 : else if (id == 2) name = "u";
554 0 : else if (id == 3) name = "s";
555 0 : else if (id == 4) name = "c";
556 0 : else if (id == 5) name = "b";
557 0 : else if (id == 6) name = "t";
558 0 : else if (id == 7) name = "b'";
559 0 : else if (id == 8) name = "t'";
560 0 : else if (id == 11) name = "e";
561 0 : else if (id == 12) name = "nue";
562 0 : else if (id == 13) name = "mu";
563 0 : else if (id == 14) name = "numu";
564 0 : else if (id == 15) name = "tau";
565 0 : else if (id == 16) name = "nutau";
566 0 : else return 0;
567 0 : return settingsPtr->parm("Zprime:" + type + name);
568 :
569 0 : }
570 :
571 : //==========================================================================
572 :
573 : // Helicity matrix element for X -> two fermions.
574 :
575 : // Base class for the W, photon, and Z -> two fermions helicity matrix
576 : // elements.
577 :
578 : //--------------------------------------------------------------------------
579 :
580 : // Initialize wave functions for the helicity matrix element.
581 :
582 : void HMEX2TwoFermions::initWaves(vector<HelicityParticle>& p) {
583 :
584 0 : u.clear();
585 0 : pMap.resize(4);
586 : // Initialize W wave function.
587 0 : vector< Wave4 > u1;
588 0 : pMap[1] = 1;
589 0 : for (int h = 0; h < p[pMap[1]].spinStates(); h++)
590 0 : u1.push_back(p[pMap[1]].wave(h));
591 0 : u.push_back(u1);
592 : // Initialize fermion wave functions.
593 0 : setFermionLine(2, p[2], p[3]);
594 :
595 0 : }
596 :
597 : //==========================================================================
598 :
599 : // Helicity matrix element for W/W' -> two fermions.
600 :
601 : // This matrix element is used when the production of the W is from an
602 : // unknown process.
603 :
604 : //--------------------------------------------------------------------------
605 :
606 : // Initialize the constants for the helicity matrix element.
607 :
608 : void HMEW2TwoFermions::initConstants() {
609 :
610 : // Set the constants for the W'.
611 0 : if (abs(pID[0]) == 34 && settingsPtr) {
612 0 : if (abs(pID[2]) < 11) {
613 0 : p2CA = settingsPtr->parm("Wprime:aq");
614 0 : p2CV = settingsPtr->parm("Wprime:vq");
615 0 : } else {
616 0 : p2CA = settingsPtr->parm("Wprime:al");
617 0 : p2CV = settingsPtr->parm("Wprime:vl");
618 : }
619 :
620 : // The default constants (SM W).
621 : } else {
622 0 : p2CA = -1; p2CV = 1;
623 : }
624 :
625 0 : }
626 :
627 : //--------------------------------------------------------------------------
628 :
629 : // Return element for helicity matrix element.
630 :
631 : complex HMEW2TwoFermions::calculateME(vector<int> h) {
632 :
633 0 : complex answer(0,0);
634 0 : for (int mu = 0; mu <= 3; mu++) {
635 0 : answer +=
636 0 : u[0][h[pMap[1]]](mu) * (u[2][h[pMap[3]]] * gamma[mu]
637 0 : * (p2CV + p2CA * gamma[5]) * u[1][h[pMap[2]]]);
638 : }
639 0 : return answer;
640 :
641 0 : }
642 :
643 : //==========================================================================
644 :
645 : // Helicity matrix element for photon -> two fermions.
646 :
647 : // This matrix element is used when the production of the photon is from an
648 : // unknown process.
649 :
650 : //--------------------------------------------------------------------------
651 :
652 : // Return element for helicity matrix element.
653 :
654 : complex HMEGamma2TwoFermions::calculateME(vector<int> h) {
655 :
656 0 : complex answer(0,0);
657 0 : for (int mu = 0; mu <= 3; mu++) {
658 0 : answer +=
659 0 : u[0][h[pMap[1]]](mu) * (u[2][h[pMap[3]]] * gamma[mu] * u[1][h[pMap[2]]]);
660 : }
661 0 : return answer;
662 0 : }
663 :
664 : //==========================================================================
665 :
666 : // Helicity matrix element for Z/Z' -> two fermions.
667 :
668 : // This matrix element is used when the production of the Z is from an
669 : // unknown process.
670 :
671 : // p2CA: axial coupling of particle 2 to the Z
672 : // p2CV: vector coupling of particle 2 to the Z
673 :
674 : //--------------------------------------------------------------------------
675 :
676 : // Initialize the constants for the helicity matrix element.
677 :
678 : void HMEZ2TwoFermions::initConstants() {
679 :
680 : // Set the vector and axial couplings to the fermions.
681 0 : p2CA = couplingsPtr->af(abs(pID[2]));
682 0 : p2CV = couplingsPtr->vf(abs(pID[2]));
683 0 : if (settingsPtr && abs(pID[0]) == 32) {
684 0 : p2CA = zpCoupling(abs(pID[2]), "a");
685 0 : p2CV = zpCoupling(abs(pID[2]), "v");
686 0 : }
687 :
688 0 : }
689 :
690 : //--------------------------------------------------------------------------
691 :
692 : // Return element for helicity matrix element.
693 :
694 : complex HMEZ2TwoFermions::calculateME(vector<int> h) {
695 :
696 0 : complex answer(0,0);
697 0 : for (int mu = 0; mu <= 3; mu++) {
698 0 : answer +=
699 0 : u[0][h[pMap[1]]](mu) * (u[2][h[pMap[3]]] * gamma[mu]
700 0 : * (p2CV - p2CA * gamma[5]) * u[1][h[pMap[2]]]);
701 : }
702 0 : return answer;
703 0 : }
704 :
705 : //--------------------------------------------------------------------------
706 :
707 : // Return the Z' vector or axial coupling for a fermion.
708 :
709 : double HMEZ2TwoFermions::zpCoupling(int id, string type) {
710 :
711 0 : if (!settingsPtr) return 0;
712 0 : id = abs(id);
713 0 : string name;
714 0 : if (id == 1) name = "d";
715 0 : else if (id == 2) name = "u";
716 0 : else if (id == 3) name = "s";
717 0 : else if (id == 4) name = "c";
718 0 : else if (id == 5) name = "b";
719 0 : else if (id == 6) name = "t";
720 0 : else if (id == 7) name = "b'";
721 0 : else if (id == 8) name = "t'";
722 0 : else if (id == 11) name = "e";
723 0 : else if (id == 12) name = "nue";
724 0 : else if (id == 13) name = "mu";
725 0 : else if (id == 14) name = "numu";
726 0 : else if (id == 15) name = "tau";
727 0 : else if (id == 16) name = "nutau";
728 0 : else return 0;
729 0 : return settingsPtr->parm("Zprime:" + type + name);
730 :
731 0 : }
732 :
733 : //==========================================================================
734 :
735 : // Helicity matrix element for the decay of a Higgs to two fermions.
736 :
737 : // All SM and MSSM Higgses couple to fermions with a vertex factor of
738 : // (pfCV - pfCA * gamma[5]) where pf indicates the type of fermion line. For
739 : // simplicity for the SM and MSSM CP even Higgses pfCV is set to one, and
740 : // pfCA to zero, as this matrix element is used only for calculating helicity
741 : // density matrices.
742 :
743 : // p2CA: in the SM/MSSM this coupling is zero for CP-even and for CP-odd:
744 : // -g_w * m_f / (2 * m_W)
745 : // * cot(beta) for A^0 u-type
746 : // * tan(beta) for A^0 d-type
747 : // p2CA: for the charged Higgs in the MSSM this coupling is given by:
748 : // +/- i * g / (sqrt(8.) * m_W) * (m_d * tan(beta) - m_u * cot(beta))
749 : // p2CV: in the SM/MSSM this coupling is zero for CP-odd and for CP-even:
750 : // i * g_w * m_f / (2 * m_W)
751 : // * -1 for the SM H
752 : // * -sin(alpha) / sin(beta) for H^0 u-type
753 : // * -cos(alpha) / cos(beta) for H^0 d-type
754 : // * -cos(alpha) / sin(beta) for h^0 u-type
755 : // * sin(alpha) / cos(beta) for h^0 d-type
756 : // p2CV: for the charged Higgs in the MSSM this coupling is given by:
757 : // i * g / (sqrt(8.) * m_W) * (m_d * tan(beta) + m_u * cot(beta))
758 :
759 : //--------------------------------------------------------------------------
760 :
761 : // Initialize the constant for the helicity matrix element.
762 :
763 : void HMEHiggs2TwoFermions::initConstants() {
764 :
765 : // Set the H4 constants.
766 0 : p2CA = 0; p2CV = 0;
767 0 : if (abs(pID[1]) == 37) {
768 0 : p2CA = pID[1] == 37 ? 1 : -1; p2CV = 1;
769 :
770 : // Neutral constants; settings available.
771 0 : } else if (settingsPtr) {
772 : int mode;
773 : double eta, phi;
774 : // Set the H1 mixing.
775 0 : if (abs(pID[1]) == 25) {
776 0 : mode = settingsPtr->mode("HiggsH1:parity");
777 0 : eta = settingsPtr->parm("HiggsH1:etaParity");
778 0 : phi = settingsPtr->parm("HiggsH1:phiParity");
779 0 : if (mode == 2) {p2CA = 1; p2CV = 0;}
780 0 : else if (mode == 3) {p2CA = eta; p2CV = complex(0, 1);}
781 0 : else if (mode == 4) {p2CA = cos(phi); p2CV = complex(0, 1) * sin(phi);}
782 0 : else {p2CA = 0; p2CV = complex(0, 1);}
783 : // Set the H2 mixing.
784 0 : } else if (abs(pID[1]) == 35) {
785 0 : mode = settingsPtr->mode("HiggsH2:parity");
786 0 : eta = settingsPtr->parm("HiggsH2:etaParity");
787 0 : phi = settingsPtr->parm("HiggsH2:phiParity");
788 0 : if (mode == 2) {p2CA = 1; p2CV = 0;}
789 0 : else if (mode == 3) {p2CA = eta; p2CV = complex(0, 1);}
790 0 : else if (mode == 4) {p2CA = cos(phi); p2CV = complex(0, 1) * sin(phi);}
791 0 : else {p2CA = 0; p2CV = complex(0, 1);}
792 : // Set the A3 mixing.
793 0 : } else if (abs(pID[1]) == 36) {
794 0 : mode = settingsPtr->mode("HiggsA3:parity");
795 0 : eta = settingsPtr->parm("HiggsA3:etaParity");
796 0 : phi = settingsPtr->parm("HiggsA3:phiParity");
797 0 : if (mode == 1) {p2CA = 0; p2CV = complex(0, 1);}
798 0 : else if (mode == 3) {p2CA = eta; p2CV = complex(0, 1);}
799 0 : else if (mode == 4) {p2CA = cos(phi); p2CV = complex(0, 1) * sin(phi);}
800 0 : else {p2CA = 1; p2CV = 0;}
801 : }
802 :
803 : // Neutral constants; default SM/MSSM.
804 0 : } else {
805 0 : if (abs(pID[1]) == 25) {p2CA = 0; p2CV = complex(0, 1);}
806 0 : else if (abs(pID[1]) == 35) {p2CA = 0; p2CV = complex(0, 1);}
807 0 : else if (abs(pID[1]) == 36) {p2CA = 1; p2CV = 0;}
808 : }
809 0 : }
810 :
811 : //--------------------------------------------------------------------------
812 :
813 : // Initialize wave functions for the helicity matrix element.
814 :
815 : void HMEHiggs2TwoFermions::initWaves(vector<HelicityParticle>& p) {
816 :
817 0 : u.clear();
818 0 : pMap.resize(4);
819 0 : setFermionLine(2, p[2], p[3]);
820 :
821 0 : }
822 :
823 : //--------------------------------------------------------------------------
824 :
825 : // Return element for the helicity matrix element.
826 :
827 : complex HMEHiggs2TwoFermions::calculateME(vector<int> h) {
828 :
829 0 : return (u[1][h[pMap[3]]] * (p2CV + p2CA * gamma[5]) * u[0][h[pMap[2]]]);
830 :
831 0 : }
832 :
833 : //==========================================================================
834 :
835 : // Base class for all tau decay matrix elements. This class derives from
836 : // the HelicityMatrixElement class and redefines some of the virtual functions.
837 :
838 : // One new method, initHadronicCurrent is defined which initializes the
839 : // hadronic current in the initWaves method. For each tau decay matrix element
840 : // the hadronic current method must be redefined accordingly, but initWaves
841 : // should not be redefined.
842 :
843 : //--------------------------------------------------------------------------
844 :
845 : // Initialize wave functions for the helicity matrix element.
846 : void HMETauDecay::initWaves(vector<HelicityParticle>& p) {
847 :
848 0 : u.clear();
849 0 : pMap.resize(p.size());
850 0 : setFermionLine(0, p[0], p[1]);
851 0 : initHadronicCurrent(p);
852 :
853 0 : }
854 :
855 : //--------------------------------------------------------------------------
856 :
857 : // Return element for the helicity matrix element.
858 : complex HMETauDecay::calculateME(vector<int> h) {
859 :
860 0 : complex answer(0,0);
861 0 : for (int mu = 0; mu <= 3; mu++) {
862 0 : answer +=
863 0 : (u[1][h[pMap[1]]] * gamma[mu] * (1 - gamma[5]) * u[0][h[pMap[0]]])
864 0 : * gamma[4](mu,mu) * u[2][0](mu);
865 : }
866 0 : return answer;
867 :
868 0 : }
869 :
870 : //--------------------------------------------------------------------------
871 :
872 : // Return the maximum decay weight for the helicity matrix element.
873 :
874 : double HMETauDecay::decayWeightMax(vector<HelicityParticle>& p) {
875 :
876 : // Determine the maximum on-diagonal element of rho.
877 0 : double on = real(p[0].rho[0][0]) > real(p[0].rho[1][1]) ?
878 0 : real(p[0].rho[0][0]) : real(p[0].rho[1][1]);
879 : // Determine the maximum off-diagonal element of rho.
880 0 : double off = fabs(real(p[0].rho[0][1])) + fabs(imag(p[0].rho[0][1]));
881 0 : return DECAYWEIGHTMAX * (on + off);
882 :
883 : }
884 :
885 : //--------------------------------------------------------------------------
886 :
887 : // Calculate complex resonance weights given a phase and amplitude vector.
888 :
889 : void HMETauDecay::calculateResonanceWeights(vector<double>& phase,
890 : vector<double>& amplitude, vector<complex>& weight) {
891 :
892 0 : for (unsigned int i = 0; i < phase.size(); i++)
893 0 : weight.push_back(amplitude[i] * (cos(phase[i]) +
894 0 : complex(0,1) * sin(phase[i])));
895 :
896 0 : }
897 :
898 : //==========================================================================
899 :
900 : // Tau decay matrix element for tau decay into a single scalar meson.
901 :
902 : // The maximum decay weight for this matrix element can be found analytically
903 : // to be 4 * m_tau^2 * (m_tau^2 - m_meson^2). However, because m_tau >> m_meson
904 : // for the relevant tau decay channels, this expression is approximated by
905 : // m_tau^4.
906 :
907 : //--------------------------------------------------------------------------
908 :
909 : // Initialize constants for the helicity matrix element.
910 :
911 : void HMETau2Meson::initConstants() {
912 :
913 0 : DECAYWEIGHTMAX = 4*pow4(pM[0]);
914 :
915 0 : }
916 :
917 : //--------------------------------------------------------------------------
918 :
919 : // Initialize the hadronic current for the helicity matrix element.
920 :
921 : void HMETau2Meson::initHadronicCurrent(vector<HelicityParticle>& p) {
922 :
923 0 : vector< Wave4 > u2;
924 0 : pMap[2] = 2;
925 0 : u2.push_back(Wave4(p[2].p()));
926 0 : u.push_back(u2);
927 :
928 0 : }
929 :
930 : //==========================================================================
931 :
932 : // Tau decay matrix element for tau decay into two leptons. Because there is
933 : // no hadronic current, but rather a leptonic current, the calculateME and
934 : // initWaves methods must be redefined.
935 :
936 : //--------------------------------------------------------------------------
937 :
938 : // Initialize constants for the helicity matrix element.
939 :
940 : void HMETau2TwoLeptons::initConstants() {
941 :
942 0 : DECAYWEIGHTMAX = 16*pow4(pM[0]);
943 :
944 0 : }
945 :
946 : //--------------------------------------------------------------------------
947 :
948 : // Initialize spinors for the helicity matrix element.
949 :
950 : void HMETau2TwoLeptons::initWaves(vector<HelicityParticle>& p) {
951 :
952 0 : u.clear();
953 0 : pMap.resize(4);
954 0 : setFermionLine(0,p[0],p[1]);
955 0 : setFermionLine(2,p[2],p[3]);
956 :
957 0 : }
958 :
959 : //--------------------------------------------------------------------------
960 :
961 : // Return element for the helicity matrix element.
962 :
963 : complex HMETau2TwoLeptons::calculateME(vector<int> h) {
964 :
965 0 : complex answer(0,0);
966 0 : for (int mu = 0; mu <= 3; mu++) {
967 0 : answer += (u[1][h[pMap[1]]] * gamma[mu] * (1 - gamma[5])
968 0 : * u[0][h[pMap[0]]]) * gamma[4](mu,mu) * (u[3][h[pMap[3]]]
969 0 : * gamma[mu] * (1 - gamma[5]) * u[2][h[pMap[2]]]);
970 : }
971 0 : return answer;
972 :
973 0 : }
974 :
975 : //==========================================================================
976 :
977 : // Tau decay matrix element for tau decay into two mesons through an
978 : // intermediate vector meson. This matrix element is used for pi^0 + pi^-
979 : // decays (rho resonances), K^0 + K^- decays (rho resonances), and eta + K^-
980 : // decays (K^* resonances). Note that for the rho resonances the pi^0 + pi^-
981 : // running width dominates while for the K^* resonances the pi^- + K^0 running
982 : // width dominates.
983 :
984 : // vecM: on-shell masses for the vector resonances
985 : // vecG: on-shell widths for the vector resonances
986 : // vecP: phases used to calculate vector resonance weights
987 : // vecA: amplitudes used to calculate vector resonance weights
988 : // vecW: vector resonance weights
989 :
990 : //--------------------------------------------------------------------------
991 :
992 : // Initialize constants for the helicity matrix element.
993 :
994 : void HMETau2TwoMesonsViaVector::initConstants() {
995 :
996 : // Clear the vectors from previous decays.
997 0 : vecM.clear(); vecG.clear(); vecP.clear(); vecA.clear(); vecW.clear();
998 :
999 : // Decay through K^* resonances (eta + K^- decay).
1000 0 : if (abs(pID[2]) == 221) {
1001 0 : DECAYWEIGHTMAX = 10;
1002 0 : pM[2] = particleDataPtr->m0(211); pM[3] = particleDataPtr->m0(311);
1003 0 : vecM.push_back(0.8921); vecM.push_back(1.700);
1004 0 : vecG.push_back(0.0513); vecG.push_back(0.235);
1005 0 : vecP.push_back(0); vecP.push_back(M_PI);
1006 0 : vecA.push_back(1); vecA.push_back(0.038);
1007 0 : }
1008 :
1009 : // Decay through rho resonances (pi^0 + pi^- and K^0 + K^- decays).
1010 : else {
1011 0 : if (abs(pID[2]) == 111) DECAYWEIGHTMAX = 800;
1012 0 : else if (abs(pID[2]) == 311) DECAYWEIGHTMAX = 6;
1013 0 : pM[2] = particleDataPtr->m0(111); pM[3] = particleDataPtr->m0(211);
1014 0 : vecM.push_back(0.7746); vecM.push_back(1.4080); vecM.push_back(1.700);
1015 0 : vecG.push_back(0.1490); vecG.push_back(0.5020); vecG.push_back(0.235);
1016 0 : vecP.push_back(0); vecP.push_back(M_PI); vecP.push_back(0);
1017 0 : vecA.push_back(1.0); vecA.push_back(0.167); vecA.push_back(0.050);
1018 : }
1019 0 : calculateResonanceWeights(vecP, vecA, vecW);
1020 :
1021 0 : }
1022 :
1023 : //--------------------------------------------------------------------------
1024 :
1025 : // Initialize the hadronic current for the helicity matrix element.
1026 :
1027 : void HMETau2TwoMesonsViaVector::initHadronicCurrent(
1028 : vector<HelicityParticle>& p) {
1029 :
1030 0 : vector< Wave4 > u2;
1031 0 : Wave4 u3(p[3].p() - p[2].p());
1032 0 : Wave4 u4(p[2].p() + p[3].p());
1033 0 : double s1 = m2(u3, u4);
1034 0 : double s2 = m2(u4);
1035 0 : complex sumBW = 0;
1036 0 : for (unsigned int i = 0; i < vecW.size(); i++)
1037 0 : sumBW += vecW[i] * pBreitWigner(pM[2], pM[3], s2, vecM[i], vecG[i]);
1038 0 : u2.push_back((u3 - s1 / s2 * u4) * sumBW);
1039 0 : u.push_back(u2);
1040 :
1041 0 : }
1042 :
1043 : //==========================================================================
1044 :
1045 : // Tau decay matrix element for tau decay into two mesons through both
1046 : // intermediate vector and scalar mesons.
1047 :
1048 : // scaC: scalar coupling constant
1049 : // scaM: on-shell masses for the scalar resonances
1050 : // scaG: on-shell widths for the scalar resonances
1051 : // scaP: phases used to calculate scalar resonance weights
1052 : // scaA: amplitudes used to calculate scalar resonance weights
1053 : // scaW: scalar resonance weights
1054 : // vecC: scalar coupling constant
1055 : // vecM: on-shell masses for the vector resonances
1056 : // vecG: on-shell widths for the vector resonances
1057 : // vecP: phases used to calculate vector resonance weights
1058 : // vecA: amplitudes used to calculate vector resonance weights
1059 : // vecW: vector resonance weights
1060 :
1061 : //--------------------------------------------------------------------------
1062 :
1063 : // Initialize constants for the helicity matrix element.
1064 :
1065 : void HMETau2TwoMesonsViaVectorScalar::initConstants() {
1066 :
1067 0 : DECAYWEIGHTMAX = 5400;
1068 : // Clear the vectors from previous decays.
1069 0 : scaM.clear(); scaG.clear(); scaP.clear(); scaA.clear(); scaW.clear();
1070 0 : vecM.clear(); vecG.clear(); vecP.clear(); vecA.clear(); vecW.clear();
1071 : // Scalar resonance parameters.
1072 0 : scaC = 0.465;
1073 0 : scaM.push_back(0.878);
1074 0 : scaG.push_back(0.499);
1075 0 : scaP.push_back(0);
1076 0 : scaA.push_back(1);
1077 0 : calculateResonanceWeights(scaP, scaA, scaW);
1078 : // Vector resonance parameters.
1079 0 : vecC = 1;
1080 0 : vecM.push_back(0.89547); vecM.push_back(1.414);
1081 0 : vecG.push_back(0.04619); vecG.push_back(0.232);
1082 0 : vecP.push_back(0); vecP.push_back(1.4399);
1083 0 : vecA.push_back(1); vecA.push_back(0.075);
1084 0 : calculateResonanceWeights(vecP, vecA, vecW);
1085 :
1086 0 : }
1087 :
1088 : //--------------------------------------------------------------------------
1089 :
1090 : // Initialize the hadronic current for the helicity matrix element.
1091 :
1092 : void HMETau2TwoMesonsViaVectorScalar::initHadronicCurrent(
1093 : vector<HelicityParticle>& p) {
1094 :
1095 0 : vector< Wave4 > u2;
1096 0 : Wave4 u3(p[3].p() - p[2].p());
1097 0 : Wave4 u4(p[2].p() + p[3].p());
1098 0 : double s1 = m2(u3,u4);
1099 0 : double s2 = m2(u4);
1100 0 : complex scaSumBW = 0; complex scaSumW = 0;
1101 0 : complex vecSumBW = 0; complex vecSumW = 0; complex vecSumBWM = 0;
1102 0 : for (unsigned int i = 0; i < scaW.size(); i++) {
1103 0 : scaSumBW += scaW[i] * sBreitWigner(pM[2], pM[3], s2, scaM[i], scaG[i]);
1104 0 : scaSumW += scaW[i];
1105 : }
1106 0 : for (unsigned int i = 0; i < vecW.size(); i++) {
1107 0 : vecSumBW += vecW[i] * pBreitWigner(pM[2], pM[3], s2, vecM[i], vecG[i]);
1108 0 : vecSumBWM += vecW[i] * pBreitWigner(pM[2], pM[3], s2, vecM[i], vecG[i]) /
1109 0 : pow2(vecM[i]);
1110 0 : vecSumW += vecW[i];
1111 : }
1112 0 : u2.push_back(vecC * (vecSumBW * u3 - s1 * vecSumBWM * u4) / vecSumW +
1113 0 : scaC * u4 * scaSumBW / scaSumW);
1114 0 : u.push_back(u2);
1115 :
1116 0 : }
1117 :
1118 : //==========================================================================
1119 :
1120 : // Tau decay matrix element for tau decay into three mesons. This matrix
1121 : // element provides a base class for all implemented three meson decays.
1122 :
1123 : // mode: three meson decay mode of the tau
1124 : // initMode(): initialize the decay mode
1125 : // initResonances(): initialize the resonance constants
1126 : // s1, s2, s3, s4: center-of-mass energies
1127 : // q, q2, q3, q4: summed and individual hadronic momentum four-vectors
1128 : // a1BW: stored value of a1BreitWigner for speed
1129 : // a1PhaseSpace(s): phase space factor for the a1
1130 : // a1BreitWigner(s): Breit-Wigner for the a1
1131 : // T(m0, m1, s, M, G, W): sum weighted p-wave Breit-Wigners
1132 : // T(s, M, G, W): sum weighted fixed width Breit-Wigners
1133 : // F1(), F2(), F3(), F4(): sub-current form factors
1134 :
1135 : //--------------------------------------------------------------------------
1136 :
1137 : // Initialize constants for the helicity matrix element.
1138 :
1139 : void HMETau2ThreeMesons::initConstants() {
1140 :
1141 0 : initMode();
1142 0 : initResonances();
1143 :
1144 0 : }
1145 :
1146 : //--------------------------------------------------------------------------
1147 :
1148 : // Initialize the hadronic current for the helicity matrix element.
1149 :
1150 : void HMETau2ThreeMesons::initHadronicCurrent(vector<HelicityParticle>& p) {
1151 :
1152 0 : vector< Wave4 > u2;
1153 :
1154 : // Initialize the momenta.
1155 0 : initMomenta(p);
1156 :
1157 : // Calculate the center of mass energies.
1158 0 : s1 = m2(q);
1159 0 : s2 = m2(q3 + q4);
1160 0 : s3 = m2(q2 + q4);
1161 0 : s4 = m2(q2 + q3);
1162 :
1163 : // Calculate the form factors.
1164 0 : a1BW = a1BreitWigner(s1);
1165 0 : complex f1 = F1();
1166 0 : complex f2 = F2();
1167 0 : complex f3 = F3();
1168 0 : complex f4 = F4();
1169 :
1170 : // Calculate the hadronic current.
1171 0 : Wave4 u3 = (f3 - f2) * q2 + (f1 - f3) * q3 + (f2 - f1) * q4;
1172 0 : u3 = u3 - (u3 * gamma[4] * q / s1) * q;
1173 0 : if (f4 != complex(0, 0))
1174 0 : u3 = u3 + complex(0, 1) * f4 * epsilon(q2, q3, q4);
1175 0 : u2.push_back(u3);
1176 0 : u.push_back(u2);
1177 :
1178 0 : }
1179 :
1180 : //--------------------------------------------------------------------------
1181 :
1182 : // Initialize the tau decay mode.
1183 :
1184 : void HMETau2ThreeMesons::initMode() {
1185 :
1186 0 : if (abs(pID[2]) == 111 && abs(pID[3]) == 111 && abs(pID[4]) == 211)
1187 0 : mode = Pi0Pi0Pim;
1188 0 : else if (abs(pID[2]) == 211 && abs(pID[3]) == 211 && abs(pID[4]) == 211)
1189 0 : mode = PimPimPip;
1190 0 : else if (abs(pID[2]) == 111 && abs(pID[3]) == 211 && abs(pID[4]) == 311)
1191 0 : mode = Pi0PimK0b;
1192 0 : else if (abs(pID[2]) == 211 && abs(pID[3]) == 211 && abs(pID[4]) == 321)
1193 0 : mode = PimPipKm;
1194 0 : else if (abs(pID[2]) == 111 && abs(pID[3]) == 211 && abs(pID[4]) == 221)
1195 0 : mode = Pi0PimEta;
1196 0 : else if (abs(pID[2]) == 211 && abs(pID[3]) == 321 && abs(pID[4]) == 321)
1197 0 : mode = PimKmKp;
1198 0 : else if (abs(pID[2]) == 111 && abs(pID[3]) == 311 && abs(pID[4]) == 321)
1199 0 : mode = Pi0K0Km;
1200 0 : else if (abs(pID[2]) == 130 && abs(pID[3]) == 211 && abs(pID[4]) == 310)
1201 0 : mode = KlPimKs;
1202 0 : else if (abs(pID[2]) == 111 && abs(pID[3]) == 111 && abs(pID[4]) == 321)
1203 0 : mode = Pi0Pi0Km;
1204 0 : else if (abs(pID[2]) == 130 && abs(pID[3]) == 130 && abs(pID[4]) == 211)
1205 0 : mode = KlKlPim;
1206 0 : else if (abs(pID[2]) == 211 && abs(pID[3]) == 310 && abs(pID[4]) == 310)
1207 0 : mode = PimKsKs;
1208 0 : else if (abs(pID[2]) == 211 && abs(pID[3]) == 311 && abs(pID[4]) == 311)
1209 0 : mode = PimK0bK0;
1210 : else
1211 0 : mode = Uknown;
1212 0 : }
1213 :
1214 : //--------------------------------------------------------------------------
1215 :
1216 : // Initialize the momenta for the helicity matrix element.
1217 :
1218 : void HMETau2ThreeMesons::initMomenta(vector<HelicityParticle>& p) {
1219 :
1220 0 : q = Wave4(p[2].p() + p[3].p() + p[4].p());
1221 : // pi-, pi-, pi+ decay and pi0, pi0, pi- decay.
1222 0 : if (mode == PimPimPip || mode == Pi0Pi0Pim) {
1223 0 : q2 = Wave4(p[2].p()); q3 = Wave4(p[3].p()); q4 = Wave4(p[4].p());
1224 : // K-, pi-, K+ decay.
1225 0 : } else if (mode == PimKmKp) {
1226 0 : q2 = Wave4(p[3].p()); q3 = Wave4(p[2].p()); q4 = Wave4(p[4].p());
1227 : // K0, pi-, Kbar0 decay.
1228 0 : } else if (mode == PimK0bK0) {
1229 0 : q2 = Wave4(p[3].p()); q3 = Wave4(p[2].p()); q4 = Wave4(p[4].p());
1230 : // K_S0, pi-, K_S0 decay.
1231 0 : } else if (mode == PimKsKs) {
1232 0 : q2 = Wave4(p[3].p()); q3 = Wave4(p[2].p()); q4 = Wave4(p[4].p());
1233 : // K_L0, pi-, K_L0 decay.
1234 0 : } else if (mode == KlKlPim) {
1235 0 : q2 = Wave4(p[2].p()); q3 = Wave4(p[4].p()); q4 = Wave4(p[3].p());
1236 : // K_S0, pi-, K_L0 decay.
1237 0 : } else if (mode == KlPimKs) {
1238 0 : q2 = Wave4(p[4].p()); q3 = Wave4(p[3].p()); q4 = Wave4(p[2].p());
1239 0 : } // K-, pi0, K0 decay.
1240 0 : else if (mode == Pi0K0Km) {
1241 0 : q2 = Wave4(p[4].p()); q3 = Wave4(p[2].p()); q4 = Wave4(p[3].p());
1242 0 : } // pi0, pi0, K- decay.
1243 0 : else if (mode == Pi0Pi0Km) {
1244 0 : q2 = Wave4(p[2].p()); q3 = Wave4(p[3].p()); q4 = Wave4(p[4].p());
1245 0 : } // K-, pi-, pi+ decay.
1246 0 : else if (mode == PimPipKm) {
1247 0 : q2 = Wave4(p[4].p()); q3 = Wave4(p[2].p()); q4 = Wave4(p[3].p());
1248 0 : } // pi-, Kbar0, pi0 decay.
1249 0 : else if (mode == Pi0PimK0b) {
1250 0 : q2 = Wave4(p[3].p()); q3 = Wave4(p[4].p()); q4 = Wave4(p[2].p());
1251 0 : } // pi-, pi0, eta decay.
1252 0 : else if (mode == Pi0PimEta) {
1253 0 : q2 = Wave4(p[3].p()); q3 = Wave4(p[2].p()); q4 = Wave4(p[4].p());
1254 0 : }
1255 :
1256 0 : }
1257 :
1258 : //--------------------------------------------------------------------------
1259 :
1260 : // Return the phase space factor for the a1. Implements equation 3.16 of Z.
1261 : // Phys. C48 (1990) 445-452.
1262 :
1263 : double HMETau2ThreeMesons::a1PhaseSpace(double s) {
1264 :
1265 : double piM = 0.13957; // Mass of the charged pion.
1266 : double rhoM = 0.773; // Mass of the rho.
1267 0 : if (s < pow2(3 * piM))
1268 0 : return 0;
1269 0 : else if (s < pow2(rhoM + piM)) {
1270 0 : double sum = (s - 9 * piM * piM);
1271 0 : return 4.1 * sum * sum * sum * (1 - 3.3 * sum + 5.8 * sum * sum);
1272 : }
1273 : else
1274 0 : return s * (1.623 + 10.38 / s - 9.32 / (s * s) + 0.65 / (s * s * s));
1275 :
1276 0 : }
1277 :
1278 : //--------------------------------------------------------------------------
1279 :
1280 : // Return the Breit-Wigner for the a1. Implements equation 3.18
1281 : // of Z. Phys. C48 (1990) 445-452.
1282 :
1283 : complex HMETau2ThreeMesons::a1BreitWigner(double s) {
1284 :
1285 0 : double a1M = 1.251; // Mass of the a1.
1286 0 : double a1G = 0.475; // Width of the a1.
1287 0 : return a1M * a1M / (a1M * a1M - s - complex(0,1) * a1M * a1G
1288 0 : * a1PhaseSpace(s) / a1PhaseSpace(a1M * a1M));
1289 :
1290 0 : }
1291 :
1292 : //--------------------------------------------------------------------------
1293 :
1294 : // Return summed weighted running p Breit-Wigner resonances.
1295 :
1296 : complex HMETau2ThreeMesons::T(double m0, double m1, double s,
1297 : vector<double> &M, vector<double> &G, vector<double> &W) {
1298 :
1299 0 : complex num(0, 0);
1300 0 : double den(0);
1301 0 : for (unsigned int i = 0; i < M.size(); i++) {
1302 0 : num += W[i] * pBreitWigner(m0, m1, s, M[i], G[i]);
1303 0 : den += W[i];
1304 : }
1305 0 : return num / den;
1306 :
1307 0 : }
1308 :
1309 : //--------------------------------------------------------------------------
1310 :
1311 : // Return summed weighted fixed width Breit-Wigner resonances.
1312 :
1313 : complex HMETau2ThreeMesons::T(double s, vector<double> &M,
1314 : vector<double> &G, vector<double> &W) {
1315 :
1316 0 : complex num(0, 0);
1317 0 : double den(0);
1318 0 : for (unsigned int i = 0; i < M.size(); i++) {
1319 0 : num += W[i] * breitWigner(s, M[i], G[i]);
1320 0 : den += W[i];
1321 : }
1322 0 : return num / den;
1323 :
1324 0 : }
1325 :
1326 : //==========================================================================
1327 :
1328 : // Tau decay matrix element for tau decay into three pions. This matrix element
1329 : // is taken from the Herwig++ implementation based on the CLEO fits.
1330 :
1331 : // rhoM: on-shell masses for the rho resonances
1332 : // rhoG: on-shell widths for the rho resonances
1333 : // rhoPp: p-wave phase for the rho coupling to the a1
1334 : // rhoAp: p-wave amplitude for the rho coupling to the a1
1335 : // rhoPd: d-wave phase for the rho coupling to the a1
1336 : // rhoAd: d-wave amplitude for the rho coupling to the a1
1337 : // f0M: f0 on-shell mass
1338 : // f0G: f0 on-shell width
1339 : // f0P: phase for the coupling of the f0 to the a1
1340 : // f0A: amplitude for the coupling of the f0 to the a1
1341 : // f2M: f2 on-shell mass
1342 : // f2G: f2 on-shell width
1343 : // f2P: phase for the coupling of the f2 to the a1
1344 : // f2P: amplitude for the coupling of the f2 to the a1
1345 : // sigM: sigma on-shell mass
1346 : // sigG: sigma on-shell width
1347 : // sigP: phase for the coupling of the sigma to the a1
1348 : // sigA: amplitude for the coupling of the sigma to the a1
1349 :
1350 : //--------------------------------------------------------------------------
1351 :
1352 : // Initialize resonance constants for the helicity matrix element.
1353 :
1354 : void HMETau2ThreePions::initResonances() {
1355 :
1356 : // Three charged pion decay.
1357 0 : if (mode == PimPimPip) DECAYWEIGHTMAX = 6000;
1358 :
1359 : // Two neutral and one charged pion decay.
1360 0 : else DECAYWEIGHTMAX = 3000;
1361 :
1362 : // Clear the vectors from previous decays.
1363 0 : rhoM.clear(); rhoG.clear();
1364 0 : rhoPp.clear(); rhoAp.clear(); rhoWp.clear();
1365 0 : rhoPd.clear(); rhoAd.clear(); rhoWd.clear();
1366 :
1367 : // Rho parameters.
1368 0 : rhoM.push_back(.7743); rhoM.push_back(1.370); rhoM.push_back(1.720);
1369 0 : rhoG.push_back(.1491); rhoG.push_back(.386); rhoG.push_back(.250);
1370 0 : rhoPp.push_back(0); rhoPp.push_back(3.11018); rhoPp.push_back(0);
1371 0 : rhoAp.push_back(1); rhoAp.push_back(0.12); rhoAp.push_back(0);
1372 0 : rhoPd.push_back(-0.471239); rhoPd.push_back(1.66504); rhoPd.push_back(0);
1373 0 : rhoAd.push_back(0.37); rhoAd.push_back(0.87); rhoAd.push_back(0);
1374 :
1375 : // Scalar and tensor parameters.
1376 0 : f0M = 1.186; f2M = 1.275; sigM = 0.860;
1377 0 : f0G = 0.350; f2G = 0.185; sigG = 0.880;
1378 0 : f0P = -1.69646; f2P = 1.75929; sigP = 0.722566;
1379 0 : f0A = 0.77; f2A = 0.71; sigA = 2.1;
1380 :
1381 : // Calculate the weights from the phases and amplitudes.
1382 0 : calculateResonanceWeights(rhoPp, rhoAp, rhoWp);
1383 0 : calculateResonanceWeights(rhoPd, rhoAd, rhoWd);
1384 0 : f0W = f0A * (cos(f0P) + complex(0,1) * sin(f0P));
1385 0 : f2W = f2A * (cos(f2P) + complex(0,1) * sin(f2P));
1386 0 : sigW = sigA * (cos(sigP) + complex(0,1) * sin(sigP));
1387 :
1388 0 : }
1389 :
1390 : //--------------------------------------------------------------------------
1391 :
1392 : // Return the first form factor.
1393 :
1394 : complex HMETau2ThreePions::F1() {
1395 :
1396 0 : complex answer(0,0);
1397 :
1398 : // Three charged pion decay.
1399 0 : if (mode == PimPimPip) {
1400 0 : for (unsigned int i = 0; i < rhoM.size(); i++) {
1401 0 : answer += - rhoWp[i] * pBreitWigner(pM[3], pM[4], s2, rhoM[i], rhoG[i])
1402 0 : - rhoWd[i] / 3.0 * pBreitWigner(pM[2], pM[4], s3, rhoM[i], rhoG[i])
1403 0 : * (s2 - s4);
1404 : }
1405 0 : answer += -2.0 / 3.0 * (sigW * sBreitWigner(pM[2], pM[4], s3, sigM, sigG)
1406 0 : + f0W * sBreitWigner(pM[2], pM[4], s3, f0M, f0G));
1407 0 : answer += f2W * (0.5 * (s4 - s3)
1408 0 : * dBreitWigner(pM[3], pM[4], s2, f2M, f2G)
1409 0 : - 1.0 / (18 * s3) * (4 * pow2(pM[2]) - s3)
1410 0 : * (s1 + s3 - pow2(pM[2]))
1411 0 : * dBreitWigner(pM[2], pM[4], s3, f2M, f2G));
1412 0 : }
1413 :
1414 : // Two neutral and one charged pion decay.
1415 : else {
1416 0 : for (unsigned int i = 0; i < rhoM.size(); i++) {
1417 0 : answer += rhoWp[i] * pBreitWigner(pM[3], pM[4], s2, rhoM[i], rhoG[i])
1418 0 : - rhoWd[i] / 3.0 * pBreitWigner(pM[2], pM[4], s3, rhoM[i], rhoG[i])
1419 0 : * (s4 - s2 - pow2(pM[4]) + pow2(pM[2]));
1420 : }
1421 0 : answer += 2.0 / 3.0 * (sigW * sBreitWigner(pM[2], pM[3], s4, sigM, sigG)
1422 0 : + f0W * sBreitWigner(pM[2], pM[3], s4, f0M, f0G));
1423 0 : answer += f2W / (18 * s4) * (s1 - pow2(pM[4]) + s4)
1424 0 : * (4 * pow2(pM[2]) - s4) * dBreitWigner(pM[2], pM[3], s4, f2M, f2G);
1425 : }
1426 0 : return a1BW * answer;
1427 :
1428 0 : }
1429 :
1430 : //--------------------------------------------------------------------------
1431 :
1432 : // Return the second form factor.
1433 :
1434 : complex HMETau2ThreePions::F2() {
1435 :
1436 0 : complex answer(0,0);
1437 :
1438 : // Three charged pion decay.
1439 0 : if (mode == PimPimPip) {
1440 0 : for (unsigned int i = 0; i < rhoM.size(); i++) {
1441 0 : answer += -rhoWp[i] * pBreitWigner(pM[2], pM[4], s3, rhoM[i], rhoG[i])
1442 0 : - rhoWd[i] / 3.0 * pBreitWigner(pM[3], pM[4], s2, rhoM[i], rhoG[i])
1443 0 : * (s3 - s4);
1444 : }
1445 0 : answer += -2.0 / 3.0 * (sigW * sBreitWigner(pM[3], pM[4], s2, sigM, sigG)
1446 0 : + f0W * sBreitWigner(pM[3], pM[4], s2, f0M, f0G));
1447 0 : answer += f2W * (0.5 * (s4 - s2)
1448 0 : * dBreitWigner(pM[2], pM[4], s3, f2M, f2G)
1449 0 : - 1.0 / (18 * s2) * (4 * pow2(pM[2]) - s2) * (s1 + s2 - pow2(pM[2]))
1450 0 : * dBreitWigner(pM[3], pM[4], s2, f2M, f2G));
1451 0 : }
1452 :
1453 : // Two neutral and one charged pion decay.
1454 : else {
1455 0 : for (unsigned int i = 0; i < rhoM.size(); i++) {
1456 0 : answer += -rhoWp[i] / 3.0 *
1457 0 : pBreitWigner(pM[2], pM[4], s3, rhoM[i], rhoG[i]) -
1458 0 : rhoWd[i] * pBreitWigner(pM[3], pM[4], s2, rhoM[i], rhoG[i]) *
1459 0 : (s4 - s3 - pow2(pM[4]) + pow2(pM[3]));
1460 : }
1461 0 : answer += 2.0 / 3.0 * (sigW * sBreitWigner(pM[2], pM[3], s4, sigM, sigG)
1462 0 : + f0W * sBreitWigner(pM[2], pM[3], s4, f0M, f0G));
1463 0 : answer += f2W / (18 * s4) * (s1 - pow2(pM[4]) + s4) *
1464 0 : (4 * pow2(pM[2]) - s4) * dBreitWigner(pM[2], pM[3], s4, f2M, f2G);
1465 : }
1466 0 : return -a1BW * answer;
1467 :
1468 0 : }
1469 :
1470 : //--------------------------------------------------------------------------
1471 :
1472 : // Return the third form factor.
1473 :
1474 : complex HMETau2ThreePions::F3() {
1475 :
1476 0 : complex answer(0,0);
1477 :
1478 : // Three charged pion decay.
1479 0 : if (mode == PimPimPip) {
1480 0 : for (unsigned int i = 0; i < rhoM.size(); i++) {
1481 0 : answer += -rhoWd[i] * (1.0 / 3.0 * (s3 - s4)
1482 0 : * pBreitWigner(pM[3], pM[4], s2, rhoM[i], rhoG[i]) - 1.0 / 3.0
1483 0 : * (s2 - s4) * pBreitWigner(pM[2], pM[4], s3, rhoM[i], rhoG[i]));
1484 : }
1485 0 : answer += -2.0 / 3.0 * (sigW * sBreitWigner(pM[3], pM[4], s2, sigM, sigG)
1486 0 : + f0W * sBreitWigner(pM[3], pM[4], s2, f0M, f0G));
1487 0 : answer += 2.0 / 3.0 * (sigW * sBreitWigner(pM[2], pM[4], s3, sigM, sigG)
1488 0 : + f0W * sBreitWigner(pM[2], pM[4], s3, f0M, f0G));
1489 0 : answer += f2W * (-1.0 / (18 * s2) * (4 * pow2(pM[2]) - s2)
1490 0 : * (s1 + s2 - pow2(pM[2])) * dBreitWigner(pM[3], pM[4], s2, f2M, f2G)
1491 0 : + 1.0 / (18 * s3) * (4 * pow2(pM[2]) - s3) * (s1 + s3 - pow2(pM[2]))
1492 0 : * dBreitWigner(pM[2], pM[4], s3, f2M, f2G));
1493 0 : }
1494 :
1495 : // Two neutral and one charged pion decay.
1496 : else {
1497 0 : for (unsigned int i = 0; i < rhoM.size(); i++) {
1498 0 : answer += rhoWd[i] * (-1.0 / 3.0 * (s4 - s3 - pow2(pM[4]) + pow2(pM[3]))
1499 0 : * pBreitWigner(pM[3], pM[4], s2, rhoM[i], rhoG[i])
1500 0 : + 1.0 / 3.0 * (s4 - s2 - pow2(pM[4]) + pow2(pM[2]))
1501 0 : * pBreitWigner(pM[2], pM[4], s3, rhoM[i], rhoG[i]));
1502 : }
1503 0 : answer += -f2W / 2.0 * (s2 - s3)
1504 0 : * dBreitWigner(pM[2], pM[3], s4, f2M, f2G);
1505 : }
1506 0 : return a1BW * answer;
1507 :
1508 0 : }
1509 :
1510 : //--------------------------------------------------------------------------
1511 :
1512 : // Return the running width for the a1 (multiplied by a factor of a1M).
1513 :
1514 : double HMETau2ThreePions::a1PhaseSpace(double s) {
1515 :
1516 : double picM = 0.1753; // (m_pi^- + m_pi^- + m_pi^+)^2
1517 : double pinM = 0.1676; // (m_pi^0 + m_pi^0 + m_pi^-)^2
1518 : double kM = 0.496; // K mass.
1519 : double ksM = 0.894; // K^* mass.
1520 : double picG = 0; // Width contribution from three charged pions.
1521 : double pinG = 0; // Width contribution from two neutral one charged.
1522 : double kG = 0; // Width contributions from s-wave K K^*.
1523 0 : double piW = pow2(0.2384)/1.0252088; // Overall weight factor.
1524 0 : double kW = pow2(4.7621); // K K^* width weight factor.
1525 :
1526 : // Three charged pion width contribution.
1527 0 : if (s < picM)
1528 0 : picG = 0;
1529 0 : else if (s < 0.823)
1530 0 : picG = 5.80900 * pow3(s - picM) * (1.0 - 3.00980 * (s - picM) +
1531 0 : 4.5792 * pow2(s - picM));
1532 : else
1533 0 : picG = -13.91400 + 27.67900 * s - 13.39300 * pow2(s) + 3.19240 * pow3(s)
1534 0 : - 0.10487 * pow4(s);
1535 :
1536 : // Two neutral and one charged pion width contribution.
1537 0 : if (s < pinM)
1538 0 : pinG = 0;
1539 0 : else if (s < 0.823)
1540 0 : pinG = 6.28450 * pow3(s - pinM) * (1.0 - 2.95950 * (s - pinM) +
1541 0 : 4.33550 * pow2(s - pinM));
1542 : else
1543 0 : pinG = -15.41100 + 32.08800 * s - 17.66600 * pow2(s) + 4.93550 * pow3(s)
1544 0 : - 0.37498 * pow4(s);
1545 :
1546 : // K and K^* width contribution.
1547 0 : if (s > pow2(ksM + kM))
1548 0 : kG = 0.5 * sqrt((s - pow2(ksM + kM)) * (s - pow2(ksM - kM))) / s;
1549 0 : return piW*(picG + pinG + kW*kG);
1550 :
1551 : }
1552 :
1553 : //--------------------------------------------------------------------------
1554 :
1555 : // Return the Breit-Wigner for the a1.
1556 :
1557 : complex HMETau2ThreePions::a1BreitWigner(double s) {
1558 :
1559 : double a1M = 1.331; // Mass of the a1.
1560 0 : return a1M*a1M/(a1M*a1M - s - complex(0,1)*a1PhaseSpace(s));
1561 :
1562 0 : }
1563 :
1564 : //==========================================================================
1565 :
1566 : // Tau decay matrix element for tau decay into three mesons with kaons.
1567 : // The form factors are taken from hep-ph/9503474.
1568 :
1569 : // rhoMa(v): on-shell masses for the axial (vector) rho resonances
1570 : // rhoGa(v): widths for the axial (vector) rho resonances
1571 : // rhoWa(v): weights for the axial (vector) rho resonances
1572 : // kstarXa(v): on-shell masses, widths, and weights for the K* resonances
1573 : // k1Xa(b): on-shell masses, width, and weight for the K1 resonances,
1574 : // the a constants are for K1 -> K* pi, K* -> K pi
1575 : // the b constants are for K1 -> rho K, rho -> pi pi
1576 : // omegaX: on-shell masses, width, and weights for the omega reosnances
1577 : // kM: kaon mass
1578 : // piM: charged pion mass
1579 : // piW: pion coupling, f_W
1580 :
1581 : //--------------------------------------------------------------------------
1582 :
1583 : // Initialize resonance constants for the helicity matrix element.
1584 :
1585 : void HMETau2ThreeMesonsWithKaons::initResonances() {
1586 :
1587 : // K-, pi-, K+ decay.
1588 0 : if (mode == PimKmKp) DECAYWEIGHTMAX = 130;
1589 : // K0, pi-, Kbar0 decay.
1590 0 : else if (mode == PimK0bK0) DECAYWEIGHTMAX = 115;
1591 : // K_S0, pi-, K_S0 and K_L0, pi-, K_L0 decay.
1592 0 : else if (mode == PimKsKs || mode == KlKlPim) DECAYWEIGHTMAX = 230;
1593 : // K_S0, pi-, K_L0 decay.
1594 0 : else if (mode == KlPimKs) DECAYWEIGHTMAX = 230;
1595 : // K-, pi0, K0 decay.
1596 0 : else if (mode == Pi0K0Km) DECAYWEIGHTMAX = 125;
1597 : // pi0, pi0, K- decay.
1598 0 : else if (mode == Pi0Pi0Km) DECAYWEIGHTMAX = 2.5e4;
1599 : // K-, pi-, pi+ decay.
1600 0 : else if (mode == PimPipKm) DECAYWEIGHTMAX = 1.8e4;
1601 : // pi-, Kbar0, pi0 decay.
1602 0 : else if (mode == Pi0PimK0b) DECAYWEIGHTMAX = 3.9e4;
1603 :
1604 : // Clear the vectors from previous decays.
1605 0 : rhoMa.clear(); rhoGa.clear(); rhoWa.clear();
1606 0 : rhoMv.clear(); rhoGv.clear(); rhoWv.clear();
1607 0 : kstarMa.clear(); kstarGa.clear(); kstarWa.clear();
1608 0 : kstarMv.clear(); kstarGv.clear(); kstarWv.clear();
1609 0 : k1Ma.clear(); k1Ga.clear(); k1Wa.clear();
1610 0 : k1Mb.clear(); k1Gb.clear(); k1Wb.clear();
1611 0 : omegaM.clear(); omegaG.clear(); omegaW.clear();
1612 :
1613 : // Rho parameters.
1614 0 : rhoMa.push_back(0.773); rhoGa.push_back(0.145); rhoWa.push_back(1);
1615 0 : rhoMa.push_back(1.370); rhoGa.push_back(0.510); rhoWa.push_back(-0.145);
1616 0 : rhoMv.push_back(0.773); rhoGv.push_back(0.145); rhoWv.push_back(1);
1617 0 : rhoMv.push_back(1.500); rhoGv.push_back(0.220); rhoWv.push_back(-6.5 / 26.0);
1618 0 : rhoMv.push_back(1.750); rhoGv.push_back(0.120); rhoWv.push_back(-1.0 / 26.0);
1619 :
1620 : // Kstar parameters.
1621 0 : kstarMa.push_back(0.892); kstarGa.push_back(0.050);
1622 0 : kstarMa.push_back(1.412); kstarGa.push_back(0.227);
1623 0 : kstarWa.push_back(1);
1624 0 : kstarWa.push_back(-0.135);
1625 0 : kstarMv.push_back(0.892); kstarGv.push_back(0.050);
1626 0 : kstarMv.push_back(1.412); kstarGv.push_back(0.227);
1627 0 : kstarMv.push_back(1.714); kstarGv.push_back(0.323);
1628 0 : kstarWv.push_back(1);
1629 0 : kstarWv.push_back(-6.5 / 26.0);
1630 0 : kstarWv.push_back(-1.0 / 26.0);
1631 :
1632 : // K1 parameters.
1633 0 : k1Ma.push_back(1.270); k1Ga.push_back(0.090); k1Wa.push_back(0.33);
1634 0 : k1Ma.push_back(1.402); k1Ga.push_back(0.174); k1Wa.push_back(1);
1635 0 : k1Mb.push_back(1.270); k1Gb.push_back(0.090); k1Wb.push_back(1);
1636 :
1637 : // Omega and phi parameters.
1638 0 : omegaM.push_back(0.782); omegaG.push_back(0.00843); omegaW.push_back(1);
1639 0 : omegaM.push_back(1.020); omegaG.push_back(0.00443); omegaW.push_back(0.05);
1640 :
1641 : // Kaon and pion parameters
1642 0 : kM = 0.49765; piM = 0.13957; piW = 0.0942;
1643 :
1644 0 : }
1645 :
1646 : //--------------------------------------------------------------------------
1647 :
1648 : // Return the first form factor.
1649 :
1650 : complex HMETau2ThreeMesonsWithKaons::F1() {
1651 :
1652 0 : complex answer;
1653 : // K-, pi-, K+ decay.
1654 0 : if (mode == PimKmKp)
1655 0 : answer = a1BW * T(piM, kM, s2, kstarMa, kstarGa, kstarWa) / 2.0;
1656 : // K0, pi-, Kbar0 decay.
1657 0 : else if (mode == PimK0bK0)
1658 0 : answer = a1BW * T(piM, kM, s2, kstarMa, kstarGa, kstarWa) / 2.0;
1659 : // K_S0, pi-, K_S0 decay and K_L0, pi-, K_L0 decay.
1660 0 : else if (mode == PimKsKs || mode == KlKlPim)
1661 0 : answer = -a1BW * (T(piM, kM, s2, kstarMa, kstarGa, kstarWa)
1662 0 : + T(piM, kM, s4, kstarMa, kstarGa, kstarWa)) / 2.0;
1663 : // K_S0, pi-, K_L0 decay.
1664 0 : else if (mode == KlPimKs)
1665 0 : answer = a1BW * (T(piM, kM, s2, kstarMa, kstarGa, kstarWa)
1666 0 : - T(piM, kM, s4, kstarMa, kstarGa, kstarWa)) / 2.0;
1667 : // K-, pi0, K0 decay.
1668 0 : else if (mode == Pi0K0Km)
1669 0 : answer = a1BW * (T(piM, kM, s2, kstarMa, kstarGa, kstarWa)
1670 0 : - T(piM, kM, s4, kstarMa, kstarGa, kstarWa)) / 2.0;
1671 : // pi0, pi0, K- decay.
1672 0 : else if (mode == Pi0Pi0Km)
1673 0 : answer = T(s1, k1Ma, k1Ga, k1Wa)
1674 0 : * T(piM, kM, s2, kstarMa, kstarGa, kstarWa);
1675 : // K-, pi-, pi+ decay.
1676 0 : else if (mode == PimPipKm)
1677 0 : answer = T(s1, k1Mb, k1Gb, k1Wb)
1678 0 : * T(piM, piM, s2, rhoMa, rhoGa, rhoWa);
1679 : // pi-, Kbar0, pi0 decay.
1680 0 : else if (mode == Pi0PimK0b)
1681 0 : answer = T(s1, k1Ma, k1Ga, k1Wa)
1682 0 : * (T(piM, kM, s2, kstarMa, kstarGa, kstarWa)
1683 0 : - T(piM, kM, s4, kstarMa, kstarGa, kstarWa));
1684 0 : return -1.0 / 3.0 * answer;
1685 0 : }
1686 :
1687 : //--------------------------------------------------------------------------
1688 :
1689 : // Return the second form factor.
1690 :
1691 : complex HMETau2ThreeMesonsWithKaons::F2() {
1692 :
1693 0 : complex answer;
1694 : // K-, pi-, K+ decay.
1695 0 : if (mode == PimKmKp)
1696 0 : answer = a1BW * T(piM, piM, s3, rhoMa, rhoGa, rhoWa) / 2.0;
1697 : // K0, pi-, Kbar0 decay.
1698 0 : else if (mode == PimK0bK0)
1699 0 : answer = a1BW * T(piM, piM, s3, rhoMa, rhoGa, rhoWa) / 2.0;
1700 : // K_S0, pi-, K_S0 decay and K_L0, pi-, K_L0 decay.
1701 0 : else if (mode == PimKsKs || mode == KlKlPim)
1702 0 : answer = a1BW * T(piM, kM, s4, kstarMa, kstarGa, kstarWa) / 2.0;
1703 : // K_S0, pi-, K_L0 decay.
1704 0 : else if (mode == KlPimKs)
1705 0 : answer = a1BW * (2.0 * T(piM, piM, s3, rhoMa, rhoGa, rhoWa)
1706 0 : + T(piM, kM, s4, kstarMa, kstarGa, kstarWa)) / 2.0;
1707 : // K-, pi0, K0 decay.
1708 0 : else if (mode == Pi0K0Km)
1709 0 : answer = a1BW * (2.0 * T(piM, piM, s3, rhoMa, rhoGa, rhoWa)
1710 0 : + T(piM, kM, s4, kstarMa, kstarGa, kstarWa)) / 2.0;
1711 : // pi0, pi0, K- decay.
1712 0 : else if (mode == Pi0Pi0Km)
1713 0 : answer = T(s1, k1Ma, k1Ga, k1Wa)
1714 0 : * T(piM, kM, s3, kstarMa, kstarGa, kstarWa);
1715 : // K-, pi-, pi+ decay.
1716 0 : else if (mode == PimPipKm)
1717 0 : answer = T(s1, k1Ma, k1Ga, k1Wa)
1718 0 : * T(piM, kM, s3, kstarMa, kstarGa, kstarWa);
1719 : // pi-, Kbar0, pi0 decay.
1720 0 : else if (mode == Pi0PimK0b)
1721 0 : answer = 2.0 * T(s1, k1Mb, k1Gb, k1Wb)
1722 0 : * T(piM, piM, s3, rhoMa, rhoGa, rhoWa)
1723 0 : + T(s1, k1Ma, k1Ga, k1Wa) * T(piM, kM, s4, kstarMa, kstarGa, kstarWa);
1724 0 : return 1.0 / 3.0 * answer;
1725 :
1726 0 : }
1727 :
1728 : //--------------------------------------------------------------------------
1729 :
1730 : // Return the fourth form factor.
1731 :
1732 : complex HMETau2ThreeMesonsWithKaons::F4() {
1733 :
1734 0 : complex answer;
1735 : // K-, pi-, K+ decay.
1736 0 : if (mode == PimKmKp)
1737 0 : answer = (sqrt(2.) - 1) * T(piM, piM, s1, rhoMv, rhoGv, rhoWv)
1738 0 : * (sqrt(2.) * T(s3, omegaM, omegaG, omegaW)
1739 0 : + T(piM, kM, s2, kstarMa, kstarGa, kstarWa));
1740 : // K0, pi-, Kbar0 decay.
1741 0 : else if (mode == PimK0bK0)
1742 0 : answer = -(sqrt(2.) - 1) * T(piM, piM, s1, rhoMv, rhoGv, rhoWv)
1743 0 : * (sqrt(2.) * T(s3, omegaM, omegaG, omegaW)
1744 0 : + T(piM, kM, s2, kstarMa, kstarGa, kstarWa));
1745 : // K_S0, pi-, K_S0 decay and K_L0, pi-, K_L0 decay.
1746 0 : else if (mode == PimKsKs || mode == KlKlPim)
1747 0 : answer = (sqrt(2.) - 1) * T(piM, piM, s1, rhoMv, rhoGv, rhoWv)
1748 0 : * (T(piM, kM, s2, kstarMa, kstarGa, kstarWa)
1749 0 : - T(piM, kM, s4, kstarMa, kstarGa, kstarWa));
1750 : // K_S0, pi-, K_L0 decay.
1751 0 : else if (mode == KlPimKs)
1752 0 : answer = -(sqrt(2.) - 1) * T(piM, piM, s1, rhoMv, rhoGv, rhoWv)
1753 0 : * (2 * sqrt(2.) * T(s3, omegaM, omegaG, omegaW)
1754 0 : + T(piM, kM, s2, kstarMa, kstarGa, kstarWa)
1755 0 : + T(piM, kM, s4, kstarMa, kstarGa, kstarWa));
1756 : // K-, pi0, K0 decay.
1757 0 : else if (mode == Pi0K0Km)
1758 0 : answer = -(sqrt(2.) - 1) * T(piM, piM, s1, rhoMv, rhoGv, rhoWv)
1759 0 : * (T(piM, kM, s4, kstarMa, kstarGa, kstarWa)
1760 0 : - T(piM, kM, s2, kstarMa, kstarGa, kstarWa));
1761 : // pi0, pi0, K- decay.
1762 0 : else if (mode == Pi0Pi0Km)
1763 0 : answer = T(piM, kM, s1, kstarMv, kstarGv, kstarWv)
1764 0 : * (T(piM, kM, s2, kstarMa, kstarGa, kstarWa)
1765 0 : - T(piM, kM, s3, kstarMa, kstarGa, kstarWa));
1766 : // K-, pi-, pi+ decay.
1767 0 : else if (mode == PimPipKm)
1768 0 : answer = -T(piM, kM, s1, kstarMv, kstarGv, kstarWv)
1769 0 : * (T(piM, piM, s2, rhoMa, rhoGa, rhoWa)
1770 0 : + T(piM, kM, s3, kstarMa, kstarGa, kstarWa));
1771 : // pi-, Kbar0, pi0 decay.
1772 0 : else if (mode == Pi0PimK0b)
1773 0 : answer = T(piM, kM, s1, kstarMv, kstarGv, kstarWv)
1774 0 : * (2.0 * T(piM, piM, s3, rhoMa, rhoGa, rhoWa)
1775 0 : + T(piM, kM, s2, kstarMa, kstarGa, kstarWa)
1776 0 : + T(piM, kM, s4, kstarMa, kstarGa, kstarWa));
1777 0 : return 1.0 / (8.0 * M_PI * M_PI * piW * piW) * answer;
1778 :
1779 0 : }
1780 :
1781 : //==========================================================================
1782 :
1783 : // Tau decay matrix element for tau decay into three mesons. The form
1784 : // factors are taken from Comp. Phys. Com. 76 (1993) 361-380.
1785 :
1786 : // rhoMa(v): on-shell masses for the axial (vector) rho resonances
1787 : // rhoGa(v): widths for the axial (vector) rho resonances
1788 : // rhoWa(v): weights for the axial (vector) rho resonances
1789 : // kstarX: on-shell masses, widths, and weights for the K* resonances
1790 : // k1X: on-shell masses, width, and weight for the K1 resonances
1791 : // kM: kaon mass
1792 : // piM: charged pion mass
1793 : // piW: pion coupling, f_W
1794 :
1795 : //--------------------------------------------------------------------------
1796 :
1797 : // Initialize resonances for the helicity matrix element.
1798 :
1799 : void HMETau2ThreeMesonsGeneric::initResonances() {
1800 :
1801 : // pi-, pi-, pi+ decay and pi0, pi0, pi- decay.
1802 0 : if (mode == PimPimPip || mode == Pi0Pi0Pim) DECAYWEIGHTMAX = 1.3e4;
1803 : // K-, pi-, K+ decay.
1804 0 : else if (mode == PimKmKp) DECAYWEIGHTMAX = 330;
1805 : // K0, pi-, Kbar0 decay.
1806 0 : else if (mode == PimK0bK0) DECAYWEIGHTMAX = 300;
1807 : // K-, pi0, K0 decay.
1808 0 : else if (mode == Pi0K0Km) DECAYWEIGHTMAX = 40;
1809 : // pi0, pi0, K- decay.
1810 0 : else if (mode == Pi0Pi0Km) DECAYWEIGHTMAX = 9.4e4;
1811 : // K-, pi-, pi+ decay.
1812 0 : else if (mode == PimPipKm) DECAYWEIGHTMAX = 9.0e3;
1813 : // pi-, Kbar0, pi0 decay.
1814 0 : else if (mode == Pi0PimK0b) DECAYWEIGHTMAX = 1.2e4;
1815 : // pi-, pi0, eta decay.
1816 0 : else if (mode == Pi0PimEta) DECAYWEIGHTMAX = 360;
1817 :
1818 : // Clear the vectors from previous decays.
1819 0 : rhoMa.clear(); rhoGa.clear(); rhoWa.clear();
1820 0 : rhoMv.clear(); rhoGv.clear(); rhoWv.clear();
1821 0 : kstarM.clear(); kstarG.clear(); kstarW.clear();
1822 0 : k1M.clear(); k1G.clear(); k1W.clear();
1823 :
1824 : // Rho parameters.
1825 0 : rhoMa.push_back(0.773); rhoGa.push_back(0.145); rhoWa.push_back(1);
1826 0 : rhoMa.push_back(1.370); rhoGa.push_back(0.510); rhoWa.push_back(-0.145);
1827 0 : rhoMv.push_back(0.773); rhoGv.push_back(0.145); rhoWv.push_back(-26);
1828 0 : rhoMv.push_back(1.5); rhoGv.push_back(0.220); rhoWv.push_back(6.5);
1829 0 : rhoMv.push_back(1.75); rhoGv.push_back(0.120); rhoWv.push_back(1);
1830 :
1831 : // Kaon parameters.
1832 0 : kstarM.push_back(0.892); kstarG.push_back(0.0513); kstarW.push_back(1);
1833 0 : k1M.push_back(1.402); k1G.push_back(0.174); k1W.push_back(1);
1834 :
1835 : // Kaon and pion parameters
1836 0 : kM = 0.49765; piM = 0.13957; piW = 0.0942;
1837 :
1838 0 : }
1839 :
1840 : //--------------------------------------------------------------------------
1841 :
1842 : // Return the first form factor.
1843 :
1844 : complex HMETau2ThreeMesonsGeneric::F1() {
1845 :
1846 0 : complex answer;
1847 : // pi-, pi-, pi+ decay and pi0, pi0, pi- decay.
1848 0 : if (mode == PimPimPip || mode == Pi0Pi0Pim)
1849 0 : answer = a1BW * T(piM, piM, s2, rhoMa, rhoGa, rhoWa);
1850 : // K-, pi-, K+ decay.
1851 0 : else if (mode == PimKmKp)
1852 0 : answer = -a1BW * T(piM, kM, s2, kstarM, kstarG, kstarW) / 3.0;
1853 : // K0, pi-, Kbar0 decay.
1854 0 : else if (mode == PimK0bK0)
1855 0 : answer = -a1BW * T(piM, kM, s2, kstarM, kstarG, kstarW) / 3.0;
1856 : // K-, pi0, K0 decay.
1857 0 : else if (mode == Pi0K0Km)
1858 0 : answer = 0;
1859 : // pi0, pi0, K- decay.
1860 0 : else if (mode == Pi0Pi0Km)
1861 0 : answer = T(s1, k1M, k1G, k1W) * T(piM, kM, s2, kstarM, kstarG, kstarW);
1862 : // K-, pi-, pi+ decay.
1863 0 : else if (mode == PimPipKm)
1864 0 : answer = -T(s1, k1M, k1G, k1W) * T(piM, piM, s2, rhoMa, rhoGa, rhoWa)
1865 0 : / 3.0;
1866 : // pi-, Kbar0, pi0 decay.
1867 0 : else if (mode == Pi0PimK0b)
1868 0 : answer = 0;
1869 : // pi-, pi0, eta decay.
1870 0 : else if (mode == Pi0PimEta)
1871 0 : answer = 0;
1872 0 : return answer;
1873 :
1874 : }
1875 :
1876 : //--------------------------------------------------------------------------
1877 :
1878 : // Return the second form factor.
1879 :
1880 : complex HMETau2ThreeMesonsGeneric::F2() {
1881 :
1882 0 : complex answer;
1883 : // pi-, pi-, pi+ decay and pi0, pi0, pi- decay.
1884 0 : if (mode == PimPimPip || mode == Pi0Pi0Pim)
1885 0 : answer = -a1BW * T(piM, piM, s3, rhoMa, rhoGa, rhoWa);
1886 : // K-, pi-, K+ decay.
1887 0 : else if (mode == PimKmKp)
1888 0 : answer = a1BW * T(piM, piM, s3, rhoMa, rhoGa, rhoWa) / 3.0;
1889 : // K0, pi-, Kbar0 decay.
1890 0 : else if (mode == PimK0bK0)
1891 0 : answer = a1BW * T(piM, piM, s3, rhoMa, rhoGa, rhoWa) / 3.0;
1892 : // K-, pi0, K0 decay.
1893 0 : else if (mode == Pi0K0Km)
1894 0 : answer = a1BW * T(piM, piM, s3, rhoMa, rhoGa, rhoWa);
1895 : // pi0, pi0, K- decay.
1896 0 : else if (mode == Pi0Pi0Km)
1897 0 : answer = -T(s1, k1M, k1G, k1W) * T(piM, kM, s3, kstarM, kstarG, kstarW);
1898 : // K-, pi-, pi+ decay.
1899 0 : else if (mode == PimPipKm)
1900 0 : answer = T(s1, k1M, k1G, k1W)
1901 0 : * T(piM, kM, s3, kstarM, kstarG, kstarW) / 3.0;
1902 : // pi-, Kbar0, pi0 decay.
1903 0 : else if (mode == Pi0PimK0b)
1904 0 : answer = T(s1, k1M, k1G, k1W) * T(piM, piM, s3, rhoMa, rhoGa, rhoWa);
1905 : // pi-, pi0, eta decay.
1906 0 : else if (mode == Pi0PimEta)
1907 0 : answer = 0;
1908 0 : return answer;
1909 :
1910 : }
1911 :
1912 : //--------------------------------------------------------------------------
1913 :
1914 : // Return the fourth form factor.
1915 :
1916 : complex HMETau2ThreeMesonsGeneric::F4() {
1917 :
1918 0 : complex answer;
1919 : // pi-, pi-, pi+ decay and pi0, pi0, pi- decay.
1920 0 : if (mode == PimPimPip || mode == Pi0Pi0Pim)
1921 0 : answer = 0;
1922 : // K-, pi-, K+ decay.
1923 0 : else if (mode == PimKmKp)
1924 0 : answer = T(piM, piM, s1, rhoMv, rhoGv, rhoWv)
1925 0 : * (T(piM, piM, s3, rhoMa, rhoGa, rhoWa)
1926 0 : - 0.2 * T(piM, kM, s2, kstarM, kstarG, kstarW)) * (1.25);
1927 : // K0, pi-, Kbar0 decay.
1928 0 : else if (mode == PimK0bK0)
1929 0 : answer = -T(piM, piM, s1, rhoMv, rhoGv, rhoWv)
1930 0 : * (T(piM, piM, s3, rhoMa, rhoGa, rhoWa)
1931 0 : - 0.2 * T(piM, kM, s2, kstarM, kstarG, kstarW)) * (1.25);
1932 : // K-, pi0, K0 decay.
1933 0 : else if (mode == Pi0K0Km)
1934 0 : answer = 0;
1935 : // pi0, pi0, K- decay.
1936 0 : else if (mode == Pi0Pi0Km)
1937 0 : answer = 0;
1938 : // K-, pi-, pi+ decay.
1939 0 : else if (mode == PimPipKm)
1940 0 : answer = -T(piM, kM, s1, kstarM, kstarG, kstarW)
1941 0 : * (T(piM, piM, s2, rhoMa, rhoGa, rhoWa)
1942 0 : - 0.2 * T(piM, kM, s3, kstarM, kstarG, kstarW)) * (1.25);
1943 : // pi-, Kbar0, pi0 decay.
1944 0 : else if (mode == Pi0PimK0b)
1945 0 : answer = 2.0 * T(piM, kM, s1, kstarM, kstarG, kstarW)
1946 0 : * (T(piM, piM, s3, rhoMa, rhoGa, rhoWa)
1947 0 : - 0.2 * T(piM, kM, s2, kstarM, kstarG, kstarW)) * (1.25);
1948 : // pi-, pi0, eta decay.
1949 0 : else if (mode == Pi0PimEta)
1950 0 : answer = T(piM, piM, s1, rhoMv, rhoGv, rhoWv)
1951 0 : * T(piM, piM, s4, rhoMa, rhoGa, rhoWa);
1952 0 : return 1.0 / (4.0 * M_PI * M_PI * piW * piW) * answer;
1953 :
1954 0 : }
1955 :
1956 : //==========================================================================
1957 :
1958 : // Tau decay matrix element for tau decay into two pions with a photons taken
1959 : // from Comp. Phys. Com. 76 (1993) 361-380. Because a photon is in the final
1960 : // state the matrix element is reimplented to handle the two spin states.
1961 :
1962 : // F(s, M, G, W): form factor for resonance
1963 : // rhoM: on-shell mass of the rho resonances
1964 : // rhoG: width of the rho resonances
1965 : // rhoW: weight of the rho resonances
1966 : // omegaX: on-shell mass, width, and weight of the omega resonances
1967 : // piM: charged pion mass
1968 :
1969 : //--------------------------------------------------------------------------
1970 :
1971 : // Initialize constants for the helicity matrix element.
1972 :
1973 : void HMETau2TwoPionsGamma::initConstants() {
1974 :
1975 0 : DECAYWEIGHTMAX = 4e4;
1976 :
1977 : // Clear the vectors from previous decays.
1978 0 : rhoM.clear(); rhoG.clear(); rhoW.clear();
1979 0 : omegaM.clear(); omegaG.clear(); omegaW.clear();
1980 :
1981 : // Set parameters.
1982 0 : rhoM.push_back(0.773); rhoG.push_back(0.145); rhoW.push_back(1);
1983 0 : rhoM.push_back(1.7); rhoG.push_back(0.26); rhoW.push_back(-0.1);
1984 0 : omegaM.push_back(0.782); omegaG.push_back(0.0085); omegaW.push_back(1);
1985 0 : piM = 0.13957;
1986 :
1987 0 : }
1988 :
1989 : //--------------------------------------------------------------------------
1990 :
1991 : // Initialize wave functions for the helicity matrix element.
1992 : void HMETau2TwoPionsGamma::initWaves(vector<HelicityParticle>& p) {
1993 :
1994 : // Calculate the hadronic current.
1995 0 : u.clear();
1996 0 : pMap.resize(p.size());
1997 0 : setFermionLine(0, p[0], p[1]);
1998 :
1999 : // Calculate the hadronic current.
2000 0 : vector< Wave4 > u2;
2001 0 : Wave4 q(p[2].p() + p[3].p() + p[4].p());
2002 0 : Wave4 q2(p[2].p()), q3(p[3].p()), q4(p[4].p());
2003 0 : double s1 = m2(q);
2004 0 : double s2 = m2(q3 + q2);
2005 0 : complex f = F(s1, rhoM, rhoG, rhoW) * F(0, rhoM, rhoG, rhoW)
2006 0 : * F(s2, omegaM, omegaG, omegaW);
2007 0 : double q4q2 = m2(q4, q2);
2008 0 : double q4q3 = m2(q4, q3);
2009 0 : double q3q2 = m2(q3, q2);
2010 0 : for (int h = 0; h < 2; h++) {
2011 0 : Wave4 e = p[2].wave(h);
2012 0 : complex q4e = q4*gamma[4]*e;
2013 0 : complex q3e = q3*gamma[4]*e;
2014 0 : u2.push_back(f * (e * (piM*piM*q4q2 - q3q2*(q4q3 - q4q2))
2015 0 : - q3 * (q3e*q4q2 - q4e*q3q2)
2016 0 : + q2 * (q3e*q4q3 - q4e*(piM*piM + q3q2))));
2017 0 : }
2018 0 : u.push_back(u2);
2019 :
2020 0 : }
2021 :
2022 : //--------------------------------------------------------------------------
2023 :
2024 : // Return element for the helicity matrix element.
2025 : complex HMETau2TwoPionsGamma::calculateME(vector<int> h) {
2026 :
2027 0 : complex answer(0,0);
2028 0 : for (int mu = 0; mu <= 3; mu++) {
2029 0 : answer +=
2030 0 : (u[1][h[pMap[1]]] * gamma[mu] * (1 - gamma[5]) * u[0][h[pMap[0]]])
2031 0 : * gamma[4](mu,mu) * u[2][h[2]](mu);
2032 : }
2033 0 : return answer;
2034 :
2035 0 : }
2036 :
2037 : //--------------------------------------------------------------------------
2038 :
2039 : // Return the form factor.
2040 : complex HMETau2TwoPionsGamma::F(double s, vector<double> M, vector<double> G,
2041 : vector<double> W) {
2042 :
2043 0 : complex answer(0, 0);
2044 0 : for (unsigned int i = 0; i < M.size(); i++)
2045 0 : answer += W[i] / (M[i]*M[i] - s - complex(0, 1) * M[i] * G[i]);
2046 0 : return answer;
2047 :
2048 : }
2049 :
2050 : //==========================================================================
2051 :
2052 : // Tau decay matrix element for tau decay into pions using the Novosibirsk
2053 : // model of Comp. Phys. Com. 174 (2006) 818-835.
2054 :
2055 : // G(i,s): G-factor for index 1, 2, or 3
2056 : // tX(q,q1,q2,q3,q4): combined resonance current
2057 : // a1D(s): a1 Breit-Wigner denominator
2058 : // rhoD(s): rho Breit-Wigner denominator
2059 : // sigD(s): sigma Breit-Wigner denominator
2060 : // omeD(s): omega Breit-Wigner denominator
2061 : // a1FormFactor(s): form factor for the a1 resonance
2062 : // rhoFormFactor1(2)(s): form factor for the rho resonance
2063 : // omeFormFactor(s): form factor for the omega resonance
2064 : // sigM: on-shell mass of the sigma resonance
2065 : // sigG: width of the sigma resonance
2066 : // sigA: amplitude of the sigma resonance
2067 : // sigP: phase of the sigma resonance
2068 : // sigW: weight of the sigma resonance (from amplitude and weight)
2069 : // omeX: mass, width, amplitude, phase, and weight of the omega resonance
2070 : // a1X: mass and width of the a1 resonance
2071 : // rhoX: mass and width of the rho resonance
2072 : // picM: charge pion mass
2073 : // pinM: neutral pion mass
2074 : // lambda2: a1 form factor cut-off
2075 :
2076 : //--------------------------------------------------------------------------
2077 :
2078 : // Initialize constants for the helicity matrix element.
2079 :
2080 : void HMETau2FourPions::initConstants() {
2081 :
2082 0 : if (abs(pID[3]) == 111) DECAYWEIGHTMAX = 5e8;
2083 0 : else DECAYWEIGHTMAX = 5e9;
2084 0 : pinM = particleDataPtr->m0(111);
2085 0 : picM = particleDataPtr->m0(211);
2086 0 : sigM = 0.8; omeM = 0.782; a1M = 1.23; rhoM = 0.7761;
2087 0 : sigG = 0.8; omeG = 0.00841; a1G = 0.45; rhoG = 0.1445;
2088 0 : sigP = 0.43585; omeP = 0.0;
2089 0 : sigA = 1.39987; omeA = 1.0;
2090 0 : sigW = sigA*(cos(sigP)+complex(0,1)*sin(sigP));
2091 0 : omeW = omeA*(cos(omeP)+complex(0,1)*sin(omeP));
2092 0 : lambda2 = 1.2;
2093 :
2094 0 : }
2095 :
2096 : //--------------------------------------------------------------------------
2097 :
2098 : // Initialize the hadronic current for the helicity matrix element.
2099 :
2100 : void HMETau2FourPions::initHadronicCurrent(vector<HelicityParticle>& p) {
2101 :
2102 0 : vector< Wave4 > u2;
2103 :
2104 : // Create pion momenta.
2105 0 : Wave4 q(p[2].p() + p[3].p() + p[4].p()+ p[5].p());
2106 0 : Wave4 q2(p[2].p()), q3(p[3].p()), q4(p[4].p()), q5(p[5].p());
2107 :
2108 : // Calculate the four pion system energy.
2109 0 : double s = m2(q);
2110 :
2111 : // Create the hadronic current for the 3 neutral pion channel.
2112 0 : if (abs(pID[3]) == 111)
2113 0 : u2.push_back(G(1,s)*(t1(q,q3,q4,q5,q2) + t1(q,q3,q2,q5,q4) +
2114 0 : t1(q,q4,q3,q5,q2) + t1(q,q4,q2,q5,q3) +
2115 0 : t1(q,q2,q3,q5,q4) + t1(q,q2,q4,q5,q3) +
2116 0 : t2(q,q3,q5,q4,q2) + t2(q,q4,q5,q3,q2) +
2117 0 : t2(q,q2,q5,q4,q3) - t2(q,q5,q3,q4,q2) -
2118 0 : t2(q,q5,q4,q3,q2) - t2(q,q5,q2,q4,q3)));
2119 :
2120 : // Create the hadronic current for the 3 charged pion channel.
2121 0 : else if (abs(pID[3]) == 211)
2122 0 : u2.push_back(G(2,s)*(t1(q,q3,q5,q4,q2) + t1(q,q4,q5,q3,q2) +
2123 0 : t1(q,q3,q4,q5,q2) + t1(q,q4,q3,q5,q2) +
2124 0 : t1(q,q2,q4,q3,q5) + t1(q,q2,q3,q4,q5) +
2125 0 : t2(q,q2,q4,q3,q5) + t2(q,q2,q3,q4,q5) -
2126 0 : t2(q,q3,q2,q4,q5) - t2(q,q4,q2,q3,q5)) +
2127 0 : G(3,s)*(t3(q,q3,q5,q4,q2) + t3(q,q4,q5,q3,q2) -
2128 0 : t3(q,q3,q4,q5,q2) - t3(q,q4,q3,q5,q2) -
2129 0 : t3(q,q3,q2,q4,q5) - t3(q,q4,q2,q3,q5)));
2130 0 : u.push_back(u2);
2131 :
2132 0 : }
2133 :
2134 : //--------------------------------------------------------------------------
2135 :
2136 : // Return the first t-vector.
2137 :
2138 : Wave4 HMETau2FourPions::t1(Wave4 &q, Wave4 &q1, Wave4 &q2,
2139 : Wave4 &q3, Wave4 &q4) {
2140 :
2141 0 : Wave4 a1Q(q2 + q3 + q4);
2142 0 : Wave4 rhoQ(q3 + q4);
2143 0 : double a1S = m2(a1Q);
2144 0 : double rhoS = m2(rhoQ);
2145 :
2146 : // Needed to match Herwig++.
2147 0 : double gM = sqrtpos(rhoM*rhoM - 4*picM*picM) * (rhoM*rhoM - 4*picM*picM)
2148 0 : / rhoM;
2149 0 : double dm = (rhoFormFactor1(0) - rhoFormFactor1(rhoM*rhoM) +
2150 0 : rhoM*rhoM * rhoFormFactor2(rhoM*rhoM)) / gM;
2151 0 : return - a1FormFactor(a1S) / (a1D(a1S) * rhoD(rhoS)) * pow2(a1M) *
2152 0 : (rhoM*rhoM + rhoM*rhoG*dm) *
2153 0 : (m2(q,a1Q) * (m2(q3,a1Q) * q4 - m2(q4,a1Q) * q3) +
2154 0 : (m2(q,q4) * m2(q1,q3) - m2(q,q3) * m2(q1,q4)) * a1Q);
2155 :
2156 0 : }
2157 :
2158 : //--------------------------------------------------------------------------
2159 :
2160 : // Return the second t-vector.
2161 :
2162 : Wave4 HMETau2FourPions::t2(Wave4 &q, Wave4 &/*q1*/, Wave4 &q2,
2163 : Wave4 &q3, Wave4 &q4) {
2164 :
2165 0 : Wave4 a1Q(q2 + q3 + q4);
2166 0 : Wave4 sigQ(q3 + q4);
2167 0 : double a1S = m2(a1Q);
2168 0 : double sigS = m2(sigQ);
2169 0 : return sigW * a1FormFactor(a1S) / (a1D(a1S) * sigD(sigS)) *
2170 0 : pow2(a1M) * pow2(sigM) * (m2(q,a1Q) * a1S * q2 - m2(q,q2) * a1S * a1Q);
2171 :
2172 0 : }
2173 :
2174 : //--------------------------------------------------------------------------
2175 :
2176 : // Return the third t-vector.
2177 :
2178 : Wave4 HMETau2FourPions::t3(Wave4 &q, Wave4 &q1, Wave4 &q2,
2179 : Wave4 &q3, Wave4 &q4) {
2180 0 : Wave4 omeQ(q2 + q3 + q4);
2181 0 : Wave4 rhoQ(q3 + q4);
2182 0 : double omeS = m2(omeQ);
2183 0 : double rhoS = m2(rhoQ);
2184 :
2185 : // Needed to match Herwig++.
2186 0 : double gM = sqrtpos(rhoM*rhoM - 4*picM*picM) * (rhoM*rhoM - 4*picM*picM)
2187 0 : / rhoM;
2188 0 : double dm = (rhoFormFactor1(0) - rhoFormFactor1(rhoM*rhoM) +
2189 0 : rhoM*rhoM * rhoFormFactor2(rhoM*rhoM)) / gM;
2190 0 : return omeW * omeFormFactor(omeS) / (omeD(omeS) * rhoD(rhoS)) *
2191 0 : pow2(omeM) * (rhoM*rhoM + rhoM*rhoG*dm) *
2192 0 : ((m2(q,q3) * m2(q1,q4) - m2(q,q4) * m2(q1,q3)) * q2 +
2193 0 : (m2(q,q4) * m2(q1,q2) - m2(q,q2) * m2(q1,q4)) * q3 +
2194 0 : (m2(q,q2) * m2(q1,q3) - m2(q,q3) * m2(q1,q2)) * q4);
2195 :
2196 0 : }
2197 :
2198 : //--------------------------------------------------------------------------
2199 :
2200 : // Return the D function for the a1(1260).
2201 :
2202 : complex HMETau2FourPions::a1D(double s) {
2203 :
2204 : // rG is defined as the running width.
2205 0 : double rG = 0;
2206 :
2207 : // The rho and pion cut off thresholds defined in the fit.
2208 : double piM = 0.16960;
2209 : double rM = 0.83425;
2210 :
2211 : // Fit of width below three pion mass threshold.
2212 0 : if (s < piM)
2213 0 : rG = 0;
2214 :
2215 : // Fit of width below pion and rho mass threshold.
2216 0 : else if (s < rM)
2217 0 : rG = 0.003052*pow3(s - piM)*(1.0 + 151.088*(s - piM) +
2218 0 : 174.495*pow2(s - piM));
2219 :
2220 : // Fit of width above pion and rho mass threshold.
2221 : else
2222 0 : rG = 2.60817 - 2.47790*s + 0.66539*pow2(s) - 0.0678183*pow3(s) +
2223 0 : 1.66577*(s-1.23701)/s;
2224 0 : return s - a1M*a1M + complex(0,1) * sqrtpos(s) * rG;
2225 :
2226 0 : }
2227 :
2228 : //--------------------------------------------------------------------------
2229 :
2230 : // Return the D function for the rho(770).
2231 :
2232 : complex HMETau2FourPions::rhoD(double s) {
2233 :
2234 0 : double gQ = sqrtpos(s - 4*picM*picM) * (s - 4*picM*picM) / sqrtpos(s);
2235 0 : double gM = sqrtpos(rhoM*rhoM - 4*picM*picM) * (rhoM*rhoM - 4*picM*picM)
2236 0 : / rhoM;
2237 0 : double dm = (rhoFormFactor1(s) - rhoFormFactor1(rhoM*rhoM) -
2238 0 : (s - rhoM*rhoM) * rhoFormFactor2(rhoM*rhoM)) / gM;
2239 :
2240 : // Ensure complex part is zero below available channel.
2241 0 : if (s < 4*picM*picM) gQ = 0;
2242 0 : return s - rhoM*rhoM - rhoM*rhoG*dm + complex(0,1)*rhoM*rhoG*(gQ/gM);
2243 :
2244 0 : }
2245 :
2246 : //--------------------------------------------------------------------------
2247 :
2248 : // Return the D function for the sigma(800) (just s-wave running width).
2249 :
2250 : complex HMETau2FourPions::sigD(double s) {
2251 :
2252 : // Sigma decay to two neutral pions for three neutral pion channel.
2253 0 : double piM = abs(pID[3]) == 111 ? pinM : picM;
2254 0 : double gQ = sqrtpos(1.0 - 4*piM*piM/s);
2255 0 : double gM = sqrtpos(1.0 - 4*piM*piM/(sigM*sigM));
2256 0 : return s - sigM*sigM + complex(0,1)*sigM*sigG*gQ/gM;
2257 :
2258 0 : }
2259 :
2260 : //--------------------------------------------------------------------------
2261 :
2262 : // Return the D function for the omega(782).
2263 :
2264 : complex HMETau2FourPions::omeD(double s) {
2265 :
2266 0 : double g = 0;
2267 0 : double q = sqrtpos(s);
2268 0 : double x = q - omeM;
2269 :
2270 : // Fit of width given in TAUOLA.
2271 0 : if (s < 1)
2272 0 : g = 1 + 17.560*x + 141.110*pow2(x) + 894.884*pow3(x) + 4977.35*pow4(x) +
2273 0 : 7610.66*pow5(x) - 42524.4*pow6(x);
2274 : else
2275 0 : g = -1333.26 + 4860*q - 6000.81*pow2(q) + 2504.97*pow3(q);
2276 0 : if (g < 0) g = 0;
2277 0 : return s - omeM*omeM + complex(0,1)*omeM*omeG*g;
2278 :
2279 0 : }
2280 :
2281 : //--------------------------------------------------------------------------
2282 :
2283 : // Return the form factor for the a1.
2284 :
2285 : double HMETau2FourPions::a1FormFactor(double s) {
2286 :
2287 0 : return pow2((1.0 + a1M*a1M/lambda2) / (1.0 + s/lambda2));
2288 :
2289 : }
2290 :
2291 : //--------------------------------------------------------------------------
2292 :
2293 : // Return the form factor for the rho(770) (equivalent to h(s) in TAUOLA).
2294 :
2295 : double HMETau2FourPions::rhoFormFactor1(double s) {
2296 :
2297 : double f = 0.;
2298 0 : if (s > 4. * picM * picM) {
2299 0 : double thr = sqrtpos(1 - 4. * picM * picM / s);
2300 0 : f = thr * log((1. + thr) / (1. - thr)) * (s - 4. * picM * picM) / M_PI;
2301 0 : } else if (s < 0.0000001) f = -8. * picM * picM / M_PI;
2302 0 : return f;
2303 :
2304 : }
2305 :
2306 : //--------------------------------------------------------------------------
2307 :
2308 : // Return the form factor for the rho(770) (equivalent to h(s) derivative).
2309 :
2310 : double HMETau2FourPions::rhoFormFactor2(double s) {
2311 :
2312 0 : double f = sqrtpos(1 - 4*picM*picM/s);
2313 0 : if (s > 4*picM*picM)
2314 0 : f = f / (M_PI * s) * (s*f + (2*picM*picM + s)*log((1 + f) / (1 - f)));
2315 : else
2316 : f = 0;
2317 0 : return f;
2318 :
2319 : }
2320 :
2321 : //--------------------------------------------------------------------------
2322 :
2323 : // Return the form factor for the omega(782).
2324 :
2325 : double HMETau2FourPions::omeFormFactor(double /*s*/) {
2326 :
2327 0 : return 1.0;
2328 :
2329 : }
2330 :
2331 : //--------------------------------------------------------------------------
2332 :
2333 : // Return the G-functions given in TAUOLA using a piece-wise fit.
2334 :
2335 : double HMETau2FourPions::G(int i, double s) {
2336 :
2337 : // Break points for the fits.
2338 : double s0(0), s1(0), s2(0), s3(0), s4(0), s5(0);
2339 :
2340 : // Parameters for the fits.
2341 : double a1(0), b1(0);
2342 : double a2(0), b2(0), c2(0), d2(0), e2(0);
2343 : double a3(0), b3(0), c3(0), d3(0), e3(0);
2344 : double a4(0), b4(0);
2345 : double a5(0), b5(0);
2346 :
2347 : // Three neutral pion parameters.
2348 0 : if (i == 1) {
2349 : s0 = 0.614403; s1 = 0.656264; s2 = 1.57896;
2350 : s3 = 3.08198; s4 = 3.12825; s5 = 3.17488;
2351 : a1 = -23383.7; b1 = 38059.2;
2352 : a2 = 230.368; b2 = -4.39368; c2 = 687.002;
2353 : d2 = -732.581; e2 = 207.087;
2354 : a3 = 1633.92; b3 = -2596.21; c3 = 1703.08;
2355 : d3 = -501.407; e3 = 54.5919;
2356 : a4 = -2982.44; b4 = 986.009;
2357 : a5 = 6948.99; b5 = -2188.74;
2358 0 : }
2359 :
2360 : // Three charged pion parameters.
2361 0 : else if (i == 2) {
2362 : s0 = 0.614403; s1 = 0.635161; s2 = 2.30794;
2363 : s3 = 3.08198; s4 = 3.12825; s5 = 3.17488;
2364 : a1 = -54171.5; b1 = 88169.3;
2365 : a2 = 454.638; b2 = -3.07152; c2 = -48.7086;
2366 : d2 = 81.9702; e2 = -24.0564;
2367 : a3 = -162.421; b3 = 308.977; c3 = -27.7887;
2368 : d3 = -48.5957; e3 = 10.6168;
2369 : a4 = -2650.29; b4 = 879.776;
2370 : a5 = 6936.99; b5 = -2184.97;
2371 0 : }
2372 :
2373 : // Omega mediated three charged pion parameters.
2374 0 : else if (i == 3) {
2375 : s0 = 0.81364; s1 = 0.861709; s2 = 1.92621;
2376 : s3 = 3.08198; s4 = 3.12825; s5 = 3.17488;
2377 : a1 = -84888.9; b1 = 104332;
2378 : a2 = 2698.15; b2 = -3.08302; c2 = 1936.11;
2379 : d2 = -1254.59; e2 = 201.291;
2380 : a3 = 7171.65; b3 = -6387.9; c3 = 3056.27;
2381 : d3 = -888.63; e3 = 108.632;
2382 : a4 = -5607.48; b4 = 1917.27;
2383 : a5 = 26573; b5 = -8369.76;
2384 0 : }
2385 :
2386 : // Return the appropriate fit.
2387 0 : if (s < s0)
2388 0 : return 0.0;
2389 0 : else if (s < s1)
2390 0 : return a1 + b1*s;
2391 0 : else if (s < s2)
2392 0 : return a2*pow(s,b2) + c2*pow2(s) + d2*pow3(s) + e2*pow4(s);
2393 0 : else if (s < s3)
2394 0 : return a3 + b3*s + c3*pow2(s) + d3*pow3(s) + e3*pow4(s);
2395 0 : else if (s < s4)
2396 0 : return a4 + b4*s;
2397 0 : else if (s < s5)
2398 0 : return a5 + b5*s;
2399 : else
2400 0 : return 0.0;
2401 :
2402 0 : }
2403 :
2404 : //==========================================================================
2405 :
2406 : // Tau decay matrix element for tau decay into five pions using the model given
2407 : // in hep-ph/0602162v1.
2408 :
2409 : // Ja(q,q1,q2,q3,q4,q5): current through rho and omega resonances
2410 : // Jb(q,q1,q2,q3,q4,q5): current through a1 and sigma resonances
2411 : // breitWigner(s,M,G): s-wave Breit-Wigner assuming massless products
2412 : // a1M: on-shell mass of the a1 resonance
2413 : // a1G: width of the a1 resonance
2414 : // rhoX: mass and width of the rho resonance
2415 : // omegaX: mass, width, and weight of the omega resonance
2416 : // sigmaX: mass, width, and weight of the sigma resonance
2417 :
2418 : //--------------------------------------------------------------------------
2419 :
2420 : // Initialize constants for the helicity matrix element.
2421 :
2422 : void HMETau2FivePions::initConstants() {
2423 :
2424 : // pi-, pi-, pi+, pi+, pi- decay.
2425 0 : if (abs(pID[2]) == 211 && abs(pID[3]) == 211 && abs(pID[4]) == 211 &&
2426 0 : abs(pID[5]) == 211 && abs(pID[6]) == 211)
2427 0 : DECAYWEIGHTMAX = 4e4;
2428 : // pi+, pi-, pi0, pi-, pi0 decay.
2429 0 : else if (abs(pID[2]) == 111 && abs(pID[3]) == 111 && abs(pID[4]) == 211 &&
2430 0 : abs(pID[5]) == 211 && abs(pID[6]) == 211)
2431 0 : DECAYWEIGHTMAX = 1e7;
2432 : // pi0, pi0, pi-, pi0, pi0 decay.
2433 0 : else if (abs(pID[2]) == 111 && abs(pID[3]) == 111 && abs(pID[4]) == 111 &&
2434 0 : abs(pID[5]) == 111 && abs(pID[6]) == 211)
2435 0 : DECAYWEIGHTMAX = 1e5;
2436 :
2437 : // Set resonances.
2438 0 : a1M = 1.260; a1G = 0.400;
2439 0 : rhoM = 0.776; rhoG = 0.150;
2440 0 : omegaM = 0.782; omegaG = 0.0085; omegaW = 11.5;
2441 0 : sigmaM = 0.800; sigmaG = 0.600; sigmaW = 1;
2442 :
2443 0 : }
2444 :
2445 : //--------------------------------------------------------------------------
2446 :
2447 : // Initialize the hadronic current for the helicity matrix element.
2448 :
2449 : void HMETau2FivePions::initHadronicCurrent(vector<HelicityParticle>& p) {
2450 :
2451 0 : vector< Wave4 > u2;
2452 :
2453 0 : Wave4 q(p[2].p() + p[3].p() + p[4].p() + p[5].p() + p[6].p());
2454 0 : Wave4 q2(p[2].p()), q3(p[3].p()), q4(p[4].p()), q5(p[5].p()), q6(p[6].p());
2455 : // pi-, pi-, pi+, pi+, pi- decay.
2456 0 : if (abs(pID[2]) == 211 && abs(pID[3]) == 211 && abs(pID[4]) == 211 &&
2457 0 : abs(pID[5]) == 211 && abs(pID[6]) == 211)
2458 0 : u2.push_back(Jb(q, q2, q3, q5, q6, q4) + Jb(q, q4, q3, q5, q6, q2)
2459 0 : + Jb(q, q2, q4, q5, q6, q3) + Jb(q, q2, q3, q6, q5, q4)
2460 0 : + Jb(q, q4, q3, q6, q5, q2) + Jb(q, q2, q4, q6, q5, q3));
2461 : // pi+, pi-, pi0, pi-, pi0 decay.
2462 0 : else if (abs(pID[2]) == 111 && abs(pID[3]) == 111 && abs(pID[4]) == 211 &&
2463 0 : abs(pID[5]) == 211 && abs(pID[6]) == 211)
2464 0 : u2.push_back(Ja(q, q6, q4, q2, q5, q3) + Ja(q, q6, q5, q2, q4, q3)
2465 0 : + Ja(q, q6, q4, q3, q5, q2) + Ja(q, q6, q5, q3, q4, q2)
2466 0 : + Jb(q, q4, q5, q6, q2, q3) + Jb(q, q2, q3, q4, q6, q5)
2467 0 : + Jb(q, q2, q3, q5, q6, q4));
2468 : // pi0, pi0, pi-, pi0, pi0 decay.
2469 0 : else if (abs(pID[2]) == 111 && abs(pID[3]) == 111 && abs(pID[4]) == 111 &&
2470 0 : abs(pID[5]) == 111 && abs(pID[6]) == 211)
2471 0 : u2.push_back(Jb(q, q2, q3, q6, q4, q5) + Jb(q, q5, q3, q6, q4, q2)
2472 0 : + Jb(q, q3, q4, q6, q2, q5) + Jb(q, q2, q4, q6, q3, q5)
2473 0 : + Jb(q, q2, q5, q6, q4, q3) + Jb(q, q4, q5, q6, q2, q3));
2474 :
2475 0 : u.push_back(u2);
2476 :
2477 0 : }
2478 :
2479 : //--------------------------------------------------------------------------
2480 :
2481 : // Return the omega-rho hadronic current.
2482 :
2483 : Wave4 HMETau2FivePions::Ja(Wave4 &q, Wave4 &q1, Wave4 &q2,
2484 : Wave4 &q3, Wave4 &q4, Wave4 &q5) {
2485 :
2486 0 : Wave4 j = epsilon(q1, q2, q3);
2487 0 : return omegaW * (breitWigner(m2(q), a1M, a1G)
2488 0 : * breitWigner(m2(q1 + q2 + q3), omegaM, omegaG)
2489 0 : * breitWigner(m2(q4 + q5), rhoM, rhoG)
2490 0 : * epsilon(q4 - q5, j, q)
2491 0 : * (breitWigner(m2(q2 + q3), rhoM, rhoG)
2492 0 : + breitWigner(m2(q1 + q3), rhoM, rhoG)
2493 0 : + breitWigner(m2(q1 + q2), rhoM, rhoG)));
2494 :
2495 0 : }
2496 :
2497 : //--------------------------------------------------------------------------
2498 :
2499 : // Return the a1-sigma hadronic current.
2500 :
2501 : Wave4 HMETau2FivePions::Jb(Wave4 &q, Wave4 &q1, Wave4 &q2,
2502 : Wave4 &q3, Wave4 &q4, Wave4 &q5) {
2503 :
2504 0 : double s = m2(q);
2505 0 : Wave4 a1Q = q1 + q2 + q3;
2506 0 : double a1S = m2(a1Q);
2507 0 : Wave4 j = (m2(q2, q1 - q3) / a1S * a1Q - q1 + q3)
2508 0 : * breitWigner(m2(q1 + q3), rhoM, rhoG)
2509 0 : + (m2(q1, q2 - q3) / a1S * a1Q - q2 + q3)
2510 0 : * breitWigner(m2(q2 + q3), rhoM, rhoG);
2511 0 : j = (j * gamma[4] * q / s) * q - j;
2512 0 : return sigmaW * (breitWigner(s, a1M, a1G) * breitWigner(a1S, a1M, a1G)
2513 0 : * breitWigner(m2(q4 + q5), sigmaM, sigmaG) * j);
2514 :
2515 0 : }
2516 :
2517 : complex HMETau2FivePions::breitWigner(double s, double M, double G) {
2518 :
2519 0 : return M * M / (M * M - s - complex(0, 1) * M * G);
2520 :
2521 : }
2522 :
2523 : //==========================================================================
2524 :
2525 : } // end namespace Pythia8
|
