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