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) 2012 Philip Ilten, Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Function definitions (not found in the header) for physics classes
7 // used in tau decays.
8 
9 #include "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].spinType(); i++) {
57  for (int j = 0; j < p[0].spinType(); 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].spinType(); h1[i]++) {
86  for (h2[i] = 0; h2[i] < p[i].spinType(); 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].spinType(); i++) {
107  for (int j = 0; j < p[idx].spinType(); 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].spinType(); h1[i]++) {
137  for (h2[i] = 0; h2[i] < p[i].spinType(); 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].spinType(); h1[i]++) {
189  for (h2[i] = 0; h2[i] < p[i].spinType(); 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.spinType(); h++) u0.push_back(p0.wave(h));
247  for (int h = 0; h < p1.spinType(); 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.spinType(); h++) u1.push_back(p0.waveBar(h));
253  for (int h = 0; h < p1.spinType(); h++) u0.push_back(p1.wave(h));
254  }
255  u.push_back(u0); u.push_back(u1);
256 
257 }
258 
259 //--------------------------------------------------------------------------
260 
261 // Return an s-wave BreitWigner.
262 
263 complex HelicityMatrixElement::sBreitWigner(double m0, double m1, double s,
264  double M, double G) {
265 
266  double gs = sqrtpos((s - pow2(m0+m1)) * (s - pow2(m0-m1))) / (2*sqrtpos(s));
267  double gM = sqrtpos((M*M - pow2(m0+m1)) * (M*M - pow2(m0-m1))) / (2*M);
268  return M*M / (M*M - s - complex(0,1)*G*M*M/sqrtpos(s)*(gs/gM));
269 
270 }
271 
272 //--------------------------------------------------------------------------
273 
274 // Return a p-wave BreitWigner.
275 
276 complex HelicityMatrixElement::pBreitWigner(double m0, double m1, double s,
277  double M, double G) {
278 
279  double gs = sqrtpos((s - pow2(m0+m1)) * (s - pow2(m0-m1))) / (2*sqrtpos(s));
280  double gM = sqrtpos((M*M - pow2(m0+m1)) * (M*M - pow2(m0-m1))) / (2*M);
281  return M*M / (M*M - s - complex(0,1)*G*M*M/sqrtpos(s)*pow3(gs/gM));
282 
283 }
284 
285 //--------------------------------------------------------------------------
286 
287 // Return a d-wave BreitWigner.
288 
289 complex HelicityMatrixElement::dBreitWigner(double m0, double m1, double s,
290  double M, double G) {
291 
292  double gs = sqrtpos((s - pow2(m0+m1)) * (s - pow2(m0-m1))) / (2*sqrtpos(s));
293  double gM = sqrtpos((M*M - pow2(m0+m1)) * (M*M - pow2(m0-m1))) / (2*M);
294  return M*M / (M*M - s - complex(0,1)*G*M*M/sqrtpos(s)*pow5(gs/gM));
295 
296 }
297 
298 //==========================================================================
299 
300 // Helicity matrix element for two fermions -> W -> two fermions. This matrix
301 // element handles s-channel hard processes in addition to t-channel, assuming
302 // the first two particles are a fermion line and the second two particles
303 // are a fermion line. This matrix element is not scaled with respect to W
304 // propagator energy as currently this matrix element is used only for
305 // calculating helicity density matrices.
306 
307 //--------------------------------------------------------------------------
308 
309 // Initialize spinors for the helicity matrix element.
310 
311 void HMETwoFermions2W2TwoFermions::initWaves(vector<HelicityParticle>& p) {
312 
313  u.clear();
314  pMap.resize(4);
315  setFermionLine(0,p[0],p[1]);
316  setFermionLine(2,p[2],p[3]);
317 
318 }
319 
320 //--------------------------------------------------------------------------
321 
322  // Return element for the helicity matrix element.
323 
324 complex HMETwoFermions2W2TwoFermions::calculateME(vector<int> h) {
325 
326  complex answer(0,0);
327  for (int mu = 0; mu <= 3; mu++) {
328  answer += (u[1][h[pMap[1]]] * gamma[mu] * (1 - gamma[5])
329  * u[0][h[pMap[0]]]) * gamma[4](mu,mu) * (u[3][h[pMap[3]]]
330  * gamma[mu] * (1 - gamma[5]) * u[2][h[pMap[2]]]);
331  }
332  return answer;
333 
334 }
335 
336 //==========================================================================
337 
338 // Helicity matrix element for two fermions -> photon -> two fermions. This
339 // matrix element can be combined with the Z matrix element to provide full
340 // interference effects.
341 
342 // p0Q: charge of the incoming fermion line
343 // p2Q: charge of the outgoing fermion line
344 // s: center of mass energy
345 
346 //--------------------------------------------------------------------------
347 
348 // Initialize wave functions for the helicity matrix element.
349 
350 void HMETwoFermions2Gamma2TwoFermions::initWaves(
351  vector<HelicityParticle>& p) {
352 
353  u.clear();
354  pMap.resize(4);
355  setFermionLine(0, p[0], p[1]);
356  setFermionLine(2, p[2], p[3]);
357  s = pow2(p[4].m());
358  p0Q = p[0].charge(); p2Q = p[2].charge();
359 
360 }
361 
362 //--------------------------------------------------------------------------
363 
364 // Return element for the helicity matrix element.
365 
366 
367 complex HMETwoFermions2Gamma2TwoFermions::calculateME(vector<int> h) {
368 
369  complex answer(0,0);
370  for (int mu = 0; mu <= 3; mu++) {
371  answer += (u[1][h[pMap[1]]] * gamma[mu] * u[0][h[pMap[0]]])
372  * gamma[4](mu,mu) * (u[3][h[pMap[3]]] * gamma[mu] * u[2][h[pMap[2]]]);
373  }
374  return p0Q*p2Q * answer / s;
375 
376 }
377 
378 //==========================================================================
379 
380 // Helicity matrix element for two fermions -> Z -> two fermions. This matrix
381 // element can be combined with the photon matrix element to provide full
382 // interference effects.
383 
384 // Note that there is a double contraction in the Z matrix element, which can
385 // be very time consuming. If the two incoming fermions are oriented along
386 // the z-axis, their helicities must be opposite for a non-zero matrix element
387 // term. Consequently, this check is made to help speed up the matrix element.
388 
389 // sin2W: sine of the Weinberg angle
390 // cos2W: cosine of the Weinberg angle
391 // zM: on-shell mass of the Z
392 // zG: on-shell width of the Z
393 // p0CA: axial coupling of particle 0 to the Z
394 // p2CA: axial coupling of particle 2 to the Z
395 // p0CV: vector coupling of particle 0 to the Z
396 // p2CV: vector coupling of particle 2 to the Z
397 // zaxis: true if the incoming fermions are oriented along the z-axis
398 
399 //--------------------------------------------------------------------------
400 
401 // Initialize the constant for the helicity matrix element.
402 
403 void HMETwoFermions2Z2TwoFermions::initConstants() {
404 
405  // Set the Weinberg angle.
406  sin2W = couplingsPtr->sin2thetaW();
407  cos2W = couplingsPtr->cos2thetaW();
408  // Set the on-shell Z mass and width.
409  zG = particleDataPtr->mWidth(23);
410  zM = particleDataPtr->m0(23);
411  // Set the vector and axial couplings to the fermions.
412  p0CA = couplingsPtr->af(abs(pID[0]));
413  p2CA = couplingsPtr->af(abs(pID[2]));
414  p0CV = couplingsPtr->vf(abs(pID[0]));
415  p2CV = couplingsPtr->vf(abs(pID[2]));
416 
417 }
418 
419 //--------------------------------------------------------------------------
420 
421 // Initialize wave functions for the helicity matrix element.
422 
423 void HMETwoFermions2Z2TwoFermions::initWaves(vector<HelicityParticle>& p) {
424 
425  vector< Wave4 > u4;
426  u.clear();
427  pMap.resize(4);
428  setFermionLine(0, p[0], p[1]);
429  setFermionLine(2, p[2], p[3]);
430  u4.push_back(Wave4(p[2].p() + p[3].p()));
431  u.push_back(u4);
432  // Center of mass energy.
433  s = pow2(p[4].m());
434  // Check if incoming fermions are oriented along z-axis.
435  zaxis = (p[0].pAbs() == fabs(p[0].pz())) &&
436  (p[1].pAbs() == fabs(p[1].pz()));
437 
438 }
439 
440 //--------------------------------------------------------------------------
441 
442 // Return element for helicity matrix element.
443 
444 complex HMETwoFermions2Z2TwoFermions::calculateME(vector<int> h) {
445 
446  complex answer(0,0);
447  // Return zero if correct helicity conditions.
448  if (h[0] == h[1] && zaxis) return answer;
449  for (int mu = 0; mu <= 3; mu++) {
450  for (int nu = 0; nu <= 3; nu++) {
451  answer +=
452  (u[1][h[pMap[1]]] * gamma[mu] * (p0CV - p0CA * gamma[5]) *
453  u[0][h[pMap[0]]]) *
454  (gamma[4](mu,nu) - gamma[4](mu,mu)*u[4][0](mu) *
455  gamma[4](nu,nu) * u[4][0](nu) / (zM*zM)) *
456  (u[3][h[pMap[3]]] * gamma[nu] * (p2CV - p2CA * gamma[5]) *
457  u[2][h[pMap[2]]]);
458  }
459  }
460  return answer / (16 * pow2(sin2W * cos2W) *
461  (s - zM*zM + complex(0, s*zG/zM)));
462 
463 }
464 
465 //==========================================================================
466 
467 // Helicity matrix element for two fermions -> photon/Z -> two fermions. Full
468 // interference is obtained by combining the photon and Z helicity matrix
469 // elements.
470 
471 // In general the initPointers and initChannel methods should not be
472 // redeclared.
473 
474 //--------------------------------------------------------------------------
475 
476 // Initialize the matrix element.
477 
478 void HMETwoFermions2GammaZ2TwoFermions::initPointers(
479  ParticleData* particleDataPtrIn, Couplings* couplingsPtrIn) {
480 
481  zHME.initPointers(particleDataPtrIn, couplingsPtrIn);
482  gHME.initPointers(particleDataPtrIn, couplingsPtrIn);
483 
484 }
485 
486 //--------------------------------------------------------------------------
487 
488 // Initialize the channel for the helicity matrix element.
489 
490  HelicityMatrixElement* HMETwoFermions2GammaZ2TwoFermions::initChannel(
491  vector<HelicityParticle>& p) {
492 
493  zHME.initChannel(p);
494  zHME.initChannel(p);
495  return this;
496 
497 }
498 
499 //--------------------------------------------------------------------------
500 
501 // Initialize wave functions for the helicity matrix element.
502 
503 void HMETwoFermions2GammaZ2TwoFermions::initWaves(
504  vector<HelicityParticle>& p) {
505 
506  zHME.initWaves(p);
507  gHME.initWaves(p);
508 
509 }
510 
511 //--------------------------------------------------------------------------
512 
513 // Return element for the helicity matrix element.
514 
515 complex HMETwoFermions2GammaZ2TwoFermions::calculateME(vector<int> h) {
516 
517  return zHME.calculateME(h) + gHME.calculateME(h);
518 
519 }
520 
521 //==========================================================================
522 
523 // Helicity matrix element for the decay of a CP even Higgs to two fermions.
524 // All SM and MSSM Higgses couple to fermions with a vertex factor of
525 // (pfCV - pfCA * gamma[5]) where pf indicates the type of fermion line. For
526 // simplicity for the SM and MSSM CP even Higgses pfCV is set to one, and
527 // pfCA to zero, as this matrix element is used only for calculating helicity
528 // density matrices.
529 
530 // p2CA: in the SM and MSSM this coupling is zero
531 // p2CV: in the SM and MSSM this coupling is given by:
532 // i * g_w * m_f / (2 * m_W)
533 // * -1 for the SM H
534 // * -sin(alpha) / sin(beta) for H^0 u-type
535 // * -cos(alpha) / cos(beta) for H^0 d-type
536 // * -cos(alpha) / sin(beta) for h^0 u-type
537 // * sin(alpha) / cos(beta) for h^0 d-type
538 
539 //--------------------------------------------------------------------------
540 
541 // Initialize wave functions for the helicity matrix element.
542 
543 void HMEHiggsEven2TwoFermions::initWaves(vector<HelicityParticle>& p) {
544 
545  u.clear();
546  pMap.resize(4);
547  p2CA = 0; p2CV = 1;
548  setFermionLine(2, p[2], p[3]);
549 
550 }
551 
552 //--------------------------------------------------------------------------
553 
554 // Return element for the helicity matrix element.
555 
556 complex HMEHiggsEven2TwoFermions::calculateME(vector<int> h) {
557 
558  return (u[1][h[pMap[3]]] * (p2CV - p2CA * gamma[5]) * u[0][h[pMap[2]]]);
559 
560 }
561 
562 //==========================================================================
563 
564 // Helicity matrix element for the decay of a CP odd Higgs to two fermions.
565 // See HMEHiggsEven2TwoFermions for more details. For the MSSM CP odd Higgs
566 // pfCA is set to one and pfCV is set to zero.
567 
568 // p2CA: in the MSSM this coupling is given by:
569 // -g_w * m_f / (2 * m_W)
570 // * cot(beta) for A^0 u-type
571 // * tan(beta) for A^0 d-type
572 // p2CV: in the MSSM this coupling is zero
573 
574 //--------------------------------------------------------------------------
575 
576 // Initialize wave functions for the helicity matrix element.
577 
578 void HMEHiggsOdd2TwoFermions::initWaves(vector<HelicityParticle>& p) {
579 
580  u.clear();
581  pMap.resize(4);
582  p2CA = 1; p2CV = 0;
583  setFermionLine(2, p[2], p[3]);
584 
585 }
586 
587 //--------------------------------------------------------------------------
588 
589 // Return element for the helicity matrix element.
590 
591 complex HMEHiggsOdd2TwoFermions::calculateME(vector<int> h) {
592 
593  return (u[1][h[pMap[3]]] * (p2CV - p2CA * gamma[5]) * u[0][h[pMap[2]]]);
594 
595 }
596 
597 //==========================================================================
598 
599 // Helicity matrix element for the decay of a charged Higgs to two fermions.
600 // See HMEHiggsEven2TwoFermions for more details. For the MSSM charged Higgs
601 // pfCA is set to +/- one given an H^+/- and pfCV is set to one.
602 
603 // p2CA: in the MSSM this coupling is given by:
604 // i * g / (sqrt(8) * m_W) * (m_d * tan(beta) + m_u * cot(beta))
605 // p2CV: in the MSSM this coupling is given by:
606 // +/- i * g / (sqrt(8) * m_W) * (m_d * tan(beta) - m_u * cot(beta))
607 
608 //--------------------------------------------------------------------------
609 
610 // Initialize wave functions for the helicity matrix element.
611 
612 void HMEHiggsCharged2TwoFermions::initWaves(vector<HelicityParticle>& p) {
613 
614  u.clear();
615  pMap.resize(4);
616  p2CV = 1;
617  if (pID[3] == 15 || pID[3] == -16) p2CA = 1;
618  else p2CA = -1;
619  setFermionLine(2, p[2], p[3]);
620 
621 }
622 
623 //--------------------------------------------------------------------------
624 
625 // Return element for the helicity matrix element.
626 
627 complex HMEHiggsCharged2TwoFermions::calculateME(vector<int> h) {
628 
629  return (u[1][h[pMap[3]]] * (p2CV - p2CA * gamma[5]) * u[0][h[pMap[2]]]);
630 
631 }
632 
633 //==========================================================================
634 
635 // Helicity matrix element which provides an unpolarized helicity
636 // density matrix. This matrix element is used for unkown hard processes.
637 
638 // Note that calculateRho is redefined for this special case, but that in
639 // general calculateRho should not be redefined.
640 
641 //--------------------------------------------------------------------------
642 
643 // Calculate a particle's helicity density matrix.
644 
645 void HMEUnpolarized::calculateRho(unsigned int idx,
646  vector<HelicityParticle>& p) {
647 
648  for (int i = 0; i < p[idx].spinType(); i++ ) {
649  for (int j = 1; j < p[idx].spinType(); j++) {
650  if (i == j) p[idx].rho[i][j] = 1.0 /
651  static_cast<double>(p[idx].spinType());
652  else p[idx].rho[i][j] = 0;
653  }
654  }
655 
656 }
657 
658 //==========================================================================
659 
660 // Base class for all tau decay matrix elements. This class derives from
661 // the HelicityMatrixElement class and redefines some of the virtual functions.
662 
663 // One new method, initHadronicCurrent is defined which initializes the
664 // hadronic current in the initWaves method. For each tau decay matrix element
665 // the hadronic current method must be redefined accordingly, but initWaves
666 // should not be redefined.
667 
668 //--------------------------------------------------------------------------
669 
670 // Initialize wave functions for the helicity matrix element.
671 void HMETauDecay::initWaves(vector<HelicityParticle>& p) {
672 
673  u.clear();
674  pMap.resize(p.size());
675  setFermionLine(0, p[0], p[1]);
676  initHadronicCurrent(p);
677 
678 }
679 
680 //--------------------------------------------------------------------------
681 
682 // Return element for the helicity matrix element.
683 complex HMETauDecay::calculateME(vector<int> h) {
684 
685  complex answer(0,0);
686  for (int mu = 0; mu <= 3; mu++) {
687  answer +=
688  (u[1][h[pMap[1]]] * gamma[mu] * (1 - gamma[5]) * u[0][h[pMap[0]]])
689  * gamma[4](mu,mu) * u[2][0](mu);
690  }
691  return answer;
692 
693 }
694 
695 //--------------------------------------------------------------------------
696 
697 // Return the maximum decay weight for the helicity matrix element.
698 
699 double HMETauDecay::decayWeightMax(vector<HelicityParticle>& p) {
700 
701  // Determine the maximum on-diagonal element of rho.
702  double on = real(p[0].rho[0][0]) > real(p[0].rho[1][1]) ?
703  real(p[0].rho[0][0]) : real(p[0].rho[1][1]);
704  // Determine the maximum off-diagonal element of rho.
705  double off = fabs(real(p[0].rho[0][1])) + fabs(imag(p[0].rho[0][1]));
706  return DECAYWEIGHTMAX * (on + off);
707 
708 }
709 
710 //--------------------------------------------------------------------------
711 
712 // Calculate complex resonance weights given a phase and amplitude vector.
713 
714 void HMETauDecay::calculateResonanceWeights(vector<double>& phase,
715  vector<double>& amplitude, vector<complex>& weight) {
716 
717  for (unsigned int i = 0; i < phase.size(); i++)
718  weight.push_back(amplitude[i] * (cos(phase[i]) +
719  complex(0,1) * sin(phase[i])));
720 
721 }
722 
723 //==========================================================================
724 
725 // Tau decay matrix element for tau decay into a single scalar meson.
726 
727 // The maximum decay weight for this matrix element can be found analytically
728 // to be 4 * m_tau^2 * (m_tau^2 - m_meson^2). However, because m_tau >> m_meson
729 // for the relevant tau decay channels, this expression is approximated by
730 // m_tau^4.
731 
732 //--------------------------------------------------------------------------
733 
734 // Initialize constants for the helicity matrix element.
735 
736 void HMETau2Meson::initConstants() {
737 
738  DECAYWEIGHTMAX = 4*pow4(pM[0]);
739 
740 }
741 
742 //--------------------------------------------------------------------------
743 
744 // Initialize the hadronic current for the helicity matrix element.
745 
746 void HMETau2Meson::initHadronicCurrent(vector<HelicityParticle>& p) {
747 
748  vector< Wave4 > u2;
749  pMap[2] = 2;
750  u2.push_back(Wave4(p[2].p()));
751  u.push_back(u2);
752 
753 }
754 
755 //==========================================================================
756 
757 // Tau decay matrix element for tau decay into two leptons. Because there is
758 // no hadronic current, but rather a leptonic current, the calculateME and
759 // initWaves methods must be redefined.
760 
761 //--------------------------------------------------------------------------
762 
763 // Initialize constants for the helicity matrix element.
764 
765 void HMETau2TwoLeptons::initConstants() {
766 
767  DECAYWEIGHTMAX = 16*pow4(pM[0]);
768 
769 }
770 
771 //--------------------------------------------------------------------------
772 
773 // Initialize spinors for the helicity matrix element.
774 
775 void HMETau2TwoLeptons::initWaves(vector<HelicityParticle>& p) {
776 
777  u.clear();
778  pMap.resize(4);
779  setFermionLine(0,p[0],p[1]);
780  setFermionLine(2,p[2],p[3]);
781 
782 }
783 
784 //--------------------------------------------------------------------------
785 
786 // Return element for the helicity matrix element.
787 
788 complex HMETau2TwoLeptons::calculateME(vector<int> h) {
789 
790  complex answer(0,0);
791  for (int mu = 0; mu <= 3; mu++) {
792  answer += (u[1][h[pMap[1]]] * gamma[mu] * (1 - gamma[5])
793  * u[0][h[pMap[0]]]) * gamma[4](mu,mu) * (u[3][h[pMap[3]]]
794  * gamma[mu] * (1 - gamma[5]) * u[2][h[pMap[2]]]);
795  }
796  return answer;
797 
798 }
799 
800 //==========================================================================
801 
802 // Tau decay matrix element for tau decay into two mesons through an
803 // intermediate vector meson. This matrix element is used for pi^0 + pi^-
804 // decays (rho resonances), K^0 + K^- decays (rho resonances), and eta + K^-
805 // decays (K^* resonances). Note that for the rho resonances the pi^0 + pi^-
806 // running width dominates while for the K^* resonances the pi^- + K^0 running
807 // width dominates.
808 
809 // vecM: on-shell masses for the vector resonances
810 // vecG: on-shell widths for the vector resonances
811 // vecP: phases used to calculate vector resonance weights
812 // vecA: amplitudes used to calculate vector resonance weights
813 // vecW: vector resonance weights
814 
815 //--------------------------------------------------------------------------
816 
817 // Initialize constants for the helicity matrix element.
818 
819 void HMETau2TwoMesonsViaVector::initConstants() {
820 
821  // Clear the vectors from previous decays.
822  vecM.clear(); vecG.clear(); vecP.clear(); vecA.clear(); vecW.clear();
823 
824  // Decay through K^* resonances (eta + K^- decay).
825  if (abs(pID[2]) == 221) {
826  DECAYWEIGHTMAX = 10;
827  pM[2] = particleDataPtr->m0(211); pM[3] = particleDataPtr->m0(311);
828  vecM.push_back(0.8921); vecM.push_back(1.700);
829  vecG.push_back(0.0513); vecG.push_back(0.235);
830  vecP.push_back(0); vecP.push_back(M_PI);
831  vecA.push_back(1); vecA.push_back(0.038);
832  }
833 
834  // Decay through rho resonances (pi^0 + pi^- and K^0 + K^- decays).
835  else {
836  if (abs(pID[2]) == 111) DECAYWEIGHTMAX = 800;
837  else if (abs(pID[2]) == 311) DECAYWEIGHTMAX = 6;
838  pM[2] = particleDataPtr->m0(111); pM[3] = particleDataPtr->m0(211);
839  vecM.push_back(0.7746); vecM.push_back(1.4080); vecM.push_back(1.700);
840  vecG.push_back(0.1490); vecG.push_back(0.5020); vecG.push_back(0.235);
841  vecP.push_back(0); vecP.push_back(M_PI); vecP.push_back(0);
842  vecA.push_back(1.0); vecA.push_back(0.167); vecA.push_back(0.050);
843  }
844  calculateResonanceWeights(vecP, vecA, vecW);
845 
846 }
847 
848 //--------------------------------------------------------------------------
849 
850 // Initialize the hadronic current for the helicity matrix element.
851 
852 void HMETau2TwoMesonsViaVector::initHadronicCurrent(
853  vector<HelicityParticle>& p) {
854 
855  vector< Wave4 > u2;
856  Wave4 u3(p[3].p() - p[2].p());
857  Wave4 u4(p[2].p() + p[3].p());
858  double s1 = m2(u3, u4);
859  double s2 = m2(u4, u4);
860  complex sumBW = 0;
861  for (unsigned int i = 0; i < vecW.size(); i++)
862  sumBW += vecW[i] * pBreitWigner(pM[2], pM[3], s2, vecM[i], vecG[i]);
863  u2.push_back((u3 - s1 / s2 * u4) * sumBW);
864  u.push_back(u2);
865 
866 }
867 
868 //==========================================================================
869 
870 // Tau decay matrix element for tau decay into two mesons through both
871 // intermediate vector and scalar mesons.
872 
873 // scaC: scalar coupling constant
874 // scaM: on-shell masses for the scalar resonances
875 // scaG: on-shell widths for the scalar resonances
876 // scaP: phases used to calculate scalar resonance weights
877 // scaA: amplitudes used to calculate scalar resonance weights
878 // scaW: scalar resonance weights
879 // vecC: scalar coupling constant
880 // vecM: on-shell masses for the vector resonances
881 // vecG: on-shell widths for the vector resonances
882 // vecP: phases used to calculate vector resonance weights
883 // vecA: amplitudes used to calculate vector resonance weights
884 // vecW: vector resonance weights
885 
886 //--------------------------------------------------------------------------
887 
888 // Initialize constants for the helicity matrix element.
889 
890 void HMETau2TwoMesonsViaVectorScalar::initConstants() {
891 
892  DECAYWEIGHTMAX = 5400;
893  // Clear the vectors from previous decays.
894  scaM.clear(); scaG.clear(); scaP.clear(); scaA.clear(); scaW.clear();
895  vecM.clear(); vecG.clear(); vecP.clear(); vecA.clear(); vecW.clear();
896  // Scalar resonance parameters.
897  scaC = 0.465;
898  scaM.push_back(0.878);
899  scaG.push_back(0.499);
900  scaP.push_back(0);
901  scaA.push_back(1);
902  calculateResonanceWeights(scaP, scaA, scaW);
903  // Vector resonance parameters.
904  vecC = 1;
905  vecM.push_back(0.89547); vecM.push_back(1.414);
906  vecG.push_back(0.04619); vecG.push_back(0.232);
907  vecP.push_back(0); vecP.push_back(1.4399);
908  vecA.push_back(1); vecA.push_back(0.075);
909  calculateResonanceWeights(vecP, vecA, vecW);
910 
911 }
912 
913 //--------------------------------------------------------------------------
914 
915 // Initialize the hadronic current for the helicity matrix element.
916 
917 void HMETau2TwoMesonsViaVectorScalar::initHadronicCurrent(
918  vector<HelicityParticle>& p) {
919 
920  vector< Wave4 > u2;
921  Wave4 u3(p[3].p() - p[2].p());
922  Wave4 u4(p[2].p() + p[3].p());
923  double s1 = m2(u3,u4);
924  double s2 = m2(u4,u4);
925  complex scaSumBW = 0; complex scaSumW = 0;
926  complex vecSumBW = 0; complex vecSumW = 0; complex vecSumBWM = 0;
927  for (unsigned int i = 0; i < scaW.size(); i++) {
928  scaSumBW += scaW[i] * sBreitWigner(pM[2], pM[3], s2, scaM[i], scaG[i]);
929  scaSumW += scaW[i];
930  }
931  for (unsigned int i = 0; i < vecW.size(); i++) {
932  vecSumBW += vecW[i] * pBreitWigner(pM[2], pM[3], s2, vecM[i], vecG[i]);
933  vecSumBWM += vecW[i] * pBreitWigner(pM[2], pM[3], s2, vecM[i], vecG[i]) /
934  pow2(vecM[i]);
935  vecSumW += vecW[i];
936  }
937  u2.push_back(vecC * (vecSumBW * u3 - s1 * vecSumBWM * u4) / vecSumW +
938  scaC * u4 * scaSumBW / scaSumW);
939  u.push_back(u2);
940 
941 }
942 
943 //==========================================================================
944 
945 // Tau decay matrix element for tau decay into three pions. This matrix element
946 // is taken from the Herwig++ implementation based on the CLEO fits.
947 
948 // F1(): return the first form factor (both neutral and charged channels)
949 // F2(): return the second form factor (both neutral and charged channels)
950 // F3(): return the third form factor (both neutral and charged channels)
951 // a1Width(s): running width for the a1 (more documentation in method)
952 // a1BreitWigner(s): Breit-Wigner for the a1, using a1Width(s)
953 
954 // rhoM: on-shell masses for the rho resonances
955 // rhoG: on-shell widths for the rho resonances
956 // rhoPp: p-wave phase for the rho coupling to the a1
957 // rhoAp: p-wave amplitude for the rho coupling to the a1
958 // rhoPd: d-wave phase for the rho coupling to the a1
959 // rhoAd: d-wave amplitude for the rho coupling to the a1
960 // f0M: f0 on-shell mass
961 // f0G: f0 on-shell width
962 // f0P: phase for the coupling of the f0 to the a1
963 // f0A: amplitude for the coupling of the f0 to the a1
964 // f2M: f2 on-shell mass
965 // f2G: f2 on-shell width
966 // f2P: phase for the coupling of the f2 to the a1
967 // f2P: amplitude for the coupling of the f2 to the a1
968 // sigM: sigma on-shell mass
969 // sigG: sigma on-shell width
970 // sigP: phase for the coupling of the sigma to the a1
971 // sigA: amplitude for the coupling of the sigma to the a1
972 
973 //--------------------------------------------------------------------------
974 
975 // Initialize constants for the helicity matrix element.
976 
977 void HMETau2ThreePions::initConstants() {
978 
979  // Three charged pion decay.
980  if (abs(pID[2] + pID[3] + pID[4]) == 211) DECAYWEIGHTMAX = 6000;
981 
982  // Two neutral and one charged pion decay.
983  else DECAYWEIGHTMAX = 3000;
984 
985  // Clear the vectors from previous decays.
986  rhoM.clear(); rhoG.clear();
987  rhoPp.clear(); rhoAp.clear(); rhoWp.clear();
988  rhoPd.clear(); rhoAd.clear(); rhoWd.clear();
989 
990  // Rho parameters.
991  rhoM.push_back(.7743); rhoM.push_back(1.370); rhoM.push_back(1.720);
992  rhoG.push_back(.1491); rhoG.push_back(.386); rhoG.push_back(.250);
993  rhoPp.push_back(0); rhoPp.push_back(3.11018); rhoPp.push_back(0);
994  rhoAp.push_back(1); rhoAp.push_back(0.12); rhoAp.push_back(0);
995  rhoPd.push_back(-0.471239); rhoPd.push_back(1.66504); rhoPd.push_back(0);
996  rhoAd.push_back(3.7e-07); rhoAd.push_back(8.7e-07); rhoAd.push_back(0);
997 
998  // Scalar and tensor parameters.
999  f0M = 1.186; f2M = 1.275; sigM = 0.860;
1000  f0G = 0.350; f2G = 0.185; sigG = 0.880;
1001  f0P = -1.69646; f2P = 1.75929; sigP = 0.722566;
1002  f0A = 0.77; f2A = 7.1e-07; sigA = 2.1;
1003 
1004  // Calculate the weights from the phases and amplitudes.
1005  calculateResonanceWeights(rhoPp, rhoAp, rhoWp);
1006  calculateResonanceWeights(rhoPd, rhoAd, rhoWd);
1007  f0W = f0A * (cos(f0P) + complex(0,1) * sin(f0P));
1008  f2W = f2A * (cos(f2P) + complex(0,1) * sin(f2P));
1009  sigW = sigA * (cos(sigP) + complex(0,1) * sin(sigP));
1010 
1011 }
1012 
1013 //--------------------------------------------------------------------------
1014 
1015 // Initialize the hadronic current for the helicity matrix element.
1016 
1017 void HMETau2ThreePions::initHadronicCurrent(vector<HelicityParticle>& p) {
1018 
1019  vector< Wave4 > u2;
1020 
1021  // Calculate the center of mass energies.
1022  s1 = (p[2].p() + p[3].p() + p[4].p()) * (p[2].p() + p[3].p() + p[4].p());
1023  s2 = (p[3].p() + p[4].p()) * (p[3].p() + p[4].p());
1024  s3 = (p[2].p() + p[4].p()) * (p[2].p() + p[4].p());
1025  s4 = (p[2].p() + p[3].p()) * (p[2].p() + p[3].p());
1026 
1027  // Calculate the form factors.
1028  complex f1 = F1();
1029  complex f2 = F2();
1030  complex f3 = F3();
1031 
1032  // Calculate the a1 Breit-Wigner.
1033  complex a1BW = a1BreitWigner(s1);
1034  Wave4 u3(p[2].p() + p[3].p() + p[4].p());
1035  Wave4 u4 = (a1BW * ((f2 - f1) * Wave4(p[4].p()) +
1036  (f1 - f3) * Wave4(p[3].p()) +
1037  (f3 - f2) * Wave4(p[2].p())));
1038  u2.push_back(u4 - (u4 * gamma[4] * u3 / s1) * u3);
1039  u.push_back(u2);
1040 
1041 }
1042 
1043 //--------------------------------------------------------------------------
1044 
1045 // Return the first form factor.
1046 
1047 complex HMETau2ThreePions::F1() {
1048 
1049  complex answer(0,0);
1050 
1051  // Three charged pion decay.
1052  if (abs(pID[2] + pID[3] + pID[4]) == 211) {
1053  for (unsigned int i = 0; i < rhoM.size(); i++) {
1054  answer += - rhoWp[i] * pBreitWigner(pM[3], pM[4], s2, rhoM[i], rhoG[i])
1055  - rhoWd[i] / 3.0 * pBreitWigner(pM[2], pM[4], s3, rhoM[i], rhoG[i])
1056  * (s2 - s4);
1057  }
1058  answer += -2.0 / 3.0 * (sigW * sBreitWigner(pM[2], pM[4], s3, sigM, sigG)
1059  + f0W * sBreitWigner(pM[2], pM[4], s3, f0M, f0G));
1060  answer += f2W * (0.5 * (s4 - s3)
1061  * dBreitWigner(pM[3], pM[4], s2, f2M, f2G)
1062  - 1.0 / (18 * s3) * (4 * pow2(pM[2]) - s3)
1063  * (s1 + s3 - pow2(pM[2]))
1064  * dBreitWigner(pM[2], pM[4], s3, f2M, f2G));
1065  }
1066 
1067  // Two neutral and one charged pion decay.
1068  else {
1069  for (unsigned int i = 0; i < rhoM.size(); i++) {
1070  answer += rhoWp[i] * pBreitWigner(pM[3], pM[4], s2, rhoM[i], rhoG[i])
1071  - rhoWd[i] / 3.0 * pBreitWigner(pM[2], pM[4], s3, rhoM[i], rhoG[i])
1072  * (s4 - s2 - pow2(pM[4]) + pow2(pM[2]));
1073  }
1074  answer += 2.0 / 3.0 * (sigW * sBreitWigner(pM[2], pM[3], s4, sigM, sigG)
1075  + f0W * sBreitWigner(pM[2], pM[3], s4, f0M, f0G));
1076  answer += f2W / (18 * s4) * (s1 - pow2(pM[4]) + s4)
1077  * (4 * pow2(pM[2]) - s4) * dBreitWigner(pM[2], pM[3], s4, f2M, f2G);
1078  }
1079  return answer;
1080 
1081 }
1082 
1083 //--------------------------------------------------------------------------
1084 
1085 // Return the second form factor.
1086 
1087 complex HMETau2ThreePions::F2() {
1088 
1089  complex answer(0,0);
1090 
1091  // Three charged pion decay.
1092  if (abs(pID[2] + pID[3] + pID[4]) == 211) {
1093  for (unsigned int i = 0; i < rhoM.size(); i++) {
1094  answer += -rhoWp[i] * pBreitWigner(pM[2], pM[4], s3, rhoM[i], rhoG[i])
1095  - rhoWd[i] / 3.0 * pBreitWigner(pM[3], pM[4], s2, rhoM[i], rhoG[i])
1096  * (s3 - s4);
1097  }
1098  answer += -2.0 / 3.0 * (sigW * sBreitWigner(pM[3], pM[4], s2, sigM, sigG)
1099  + f0W * sBreitWigner(pM[3], pM[4], s2, f0M, f0G));
1100  answer += f2W * (0.5 * (s4 - s2)
1101  * dBreitWigner(pM[2], pM[4], s3, f2M, f2G)
1102  - 1.0 / (18 * s2) * (4 * pow2(pM[2]) - s2) * (s1 + s2 - pow2(pM[2]))
1103  * dBreitWigner(pM[3], pM[4], s2, f2M, f2G));
1104  }
1105 
1106  // Two neutral and one charged pion decay.
1107  else {
1108  for (unsigned int i = 0; i < rhoM.size(); i++) {
1109  answer += -rhoWp[i] / 3.0 *
1110  pBreitWigner(pM[2], pM[4], s3, rhoM[i], rhoG[i]) -
1111  rhoWd[i] * pBreitWigner(pM[3], pM[4], s2, rhoM[i], rhoG[i]) *
1112  (s4 - s3 - pow2(pM[4]) + pow2(pM[3]));
1113  }
1114  answer += 2.0 / 3.0 * (sigW * sBreitWigner(pM[2], pM[3], s4, sigM, sigG)
1115  + f0W * sBreitWigner(pM[2], pM[3], s4, f0M, f0G));
1116  answer += f2W / (18 * s4) * (s1 - pow2(pM[4]) + s4) *
1117  (4 * pow2(pM[2]) - s4) * dBreitWigner(pM[2], pM[3], s4, f2M, f2G);
1118  }
1119  return -answer;
1120 
1121 }
1122 
1123 //--------------------------------------------------------------------------
1124 
1125 // Return the third form factor.
1126 
1127 complex HMETau2ThreePions::F3() {
1128 
1129  complex answer(0,0);
1130 
1131  // Three charged pion decay.
1132  if (abs(pID[2] + pID[3] + pID[4]) == 211) {
1133  for (unsigned int i = 0; i < rhoM.size(); i++) {
1134  answer += -rhoWd[i] * (1.0 / 3.0 * (s3 - s4) *
1135  pBreitWigner(pM[3], pM[4], s2, rhoM[i], rhoG[i])
1136  - 1.0 / 3.0 * (s2 - s4) *
1137  pBreitWigner(pM[2], pM[4], s3, rhoM[i],
1138  rhoG[i]));
1139  }
1140  answer += -2.0 / 3.0 * (sigW * sBreitWigner(pM[3], pM[4], s2, sigM, sigG)
1141  + f0W * sBreitWigner(pM[3], pM[4], s2, f0M, f0G));
1142  answer += 2.0 / 3.0 * (sigW * sBreitWigner(pM[2], pM[4], s3, sigM, sigG)
1143  + f0W * sBreitWigner(pM[2], pM[4], s3, f0M, f0G));
1144  answer += f2W * (-1.0 / (18 * s2) * (4 * pow2(pM[2]) - s2) *
1145  (s1 + s2 - pow2(pM[2])) *
1146  dBreitWigner(pM[3], pM[4], s2, f2M, f2G) +
1147  1.0 / (18 * s3) * (4 * pow2(pM[2]) - s3) *
1148  (s1 + s3 - pow2(pM[2])) *
1149  dBreitWigner(pM[2], pM[4], s3, f2M, f2G));
1150  }
1151 
1152  // Two neutral and one charged pion decay.
1153  else {
1154  for (unsigned int i = 0; i < rhoM.size(); i++) {
1155  answer += rhoWd[i] * (-1.0 / 3.0 *
1156  (s4 - s3 - pow2(pM[4]) + pow2(pM[3])) *
1157  pBreitWigner(pM[3], pM[4], s2, rhoM[i], rhoG[i]) +
1158  1.0 / 3.0 * (s4 - s2 - pow2(pM[4]) + pow2(pM[2]))
1159  * pBreitWigner(pM[2], pM[4], s3, rhoM[i],
1160  rhoG[i]));
1161  }
1162  answer += -f2W / 2.0 * (s2 - s3) *
1163  dBreitWigner(pM[2], pM[3], s4, f2M, f2G);
1164  }
1165  return answer;
1166 
1167 }
1168 
1169 //--------------------------------------------------------------------------
1170 
1171 // Return the running width for the a1 (multiplied by a factor of a1M).
1172 
1173 double HMETau2ThreePions::a1Width(double s) {
1174 
1175  double picM = 0.1753; // (m_pi^- + m_pi^- + m_pi^+)^2
1176  double pinM = 0.1676; // (m_pi^0 + m_pi^0 + m_pi^-)^2
1177  double kM = 0.496; // K mass.
1178  double ksM = 0.894; // K^* mass.
1179  double picG = 0; // Width contribution from three charged pions.
1180  double pinG = 0; // Width contribution from two neutral one charged.
1181  double kG = 0; // Width contributions from s-wave K K^*.
1182  double piW = pow2(0.2384)/1.0252088; // Overall weight factor.
1183  double kW = pow2(4.7621); // K K^* width weight factor.
1184 
1185  // Three charged pion width contribution.
1186  if (s < picM)
1187  picG = 0;
1188  else if (s < 0.823)
1189  picG = 5.80900 * pow3(s - picM) * (1.0 - 3.00980 * (s - picM) +
1190  4.5792 * pow2(s - picM));
1191  else
1192  picG = -13.91400 + 27.67900 * s - 13.39300 * pow2(s) + 3.19240 * pow3(s)
1193  - 0.10487 * pow4(s);
1194 
1195  // Two neutral and one charged pion width contribution.
1196  if (s < pinM)
1197  pinG = 0;
1198  else if (s < 0.823)
1199  pinG = 6.28450 * pow3(s - pinM) * (1.0 - 2.95950 * (s - pinM) +
1200  4.33550 * pow2(s - pinM));
1201  else
1202  pinG = -15.41100 + 32.08800 * s - 17.66600 * pow2(s) + 4.93550 * pow3(s)
1203  - 0.37498 * pow4(s);
1204 
1205  // K and K^* width contribution.
1206  if (s > pow2(ksM + kM))
1207  kG = 0.5 * sqrt((s - pow2(ksM + kM)) * (s - pow2(ksM - kM))) / s;
1208  return piW*(picG + pinG + kW*kG);
1209 
1210 }
1211 
1212 //--------------------------------------------------------------------------
1213 
1214 // Return the Breit-Wigner for the a1.
1215 
1216 complex HMETau2ThreePions::a1BreitWigner(double s) {
1217 
1218  double a1M = 1.331; // Mass of the a1.
1219  return a1M*a1M/(a1M*a1M - s - complex(0,1)*a1Width(s));
1220 
1221 }
1222 
1223 //==========================================================================
1224 
1225 // Tau decay matrix element for tau decay into two mesons through both
1226 // intermediate vector and scalar mesons.
1227 
1228 // More detailed documentation will be written.
1229 
1230 //--------------------------------------------------------------------------
1231 
1232 // Initialize constants for the helicity matrix element.
1233 
1234 void HMETau2FourPions::initConstants() {
1235 
1236  if (abs(pID[3]) == 111) DECAYWEIGHTMAX = 5e8;
1237  else DECAYWEIGHTMAX = 5e9;
1238  pinM = particleDataPtr->m0(111);
1239  picM = particleDataPtr->m0(211);
1240  sigM = 0.8; omeM = 0.782; a1M = 1.23; rhoM = 0.7761;
1241  sigG = 0.8; omeG = 0.00841; a1G = 0.45; rhoG = 0.1445;
1242  sigP = 0.43585; omeP = 0.0;
1243  sigA = 1.39987; omeA = 1.0;
1244  sigW = sigA*(cos(sigP)+complex(0,1)*sin(sigP));
1245  omeW = omeA*(cos(omeP)+complex(0,1)*sin(omeP));
1246  lambda2 = 1.2;
1247 
1248 }
1249 
1250 //--------------------------------------------------------------------------
1251 
1252 // Initialize the hadronic current for the helicity matrix element.
1253 
1254 void HMETau2FourPions::initHadronicCurrent(vector<HelicityParticle>& p) {
1255 
1256  vector< Wave4 > u2;
1257 
1258  // Push back pion momenta (temporarily use hadronic current vector).
1259  // Store neutrino momenta as well to keep consistant notation.
1260  u2.push_back(Wave4(p[2].p() + p[3].p() + p[4].p()+ p[5].p()));
1261  u2.push_back(Wave4(p[1].p()));
1262  u2.push_back(Wave4(p[2].p()));
1263  u2.push_back(Wave4(p[3].p()));
1264  u2.push_back(Wave4(p[4].p()));
1265  u2.push_back(Wave4(p[5].p()));
1266  u.push_back(u2); u2.clear();
1267 
1268  // Calculate the four pion system energy.
1269  double s = m2(u2[0],u2[0]);
1270 
1271  // Create the hadronic current for the 3 neutral pion channel.
1272  if (abs(pID[3]) == 111)
1273  u2.push_back(G(1,s) * (t1(3,4,5,2) + t1(3,2,5,4) + t1(4,3,5,2) +
1274  t1(4,2,5,3) + t1(2,3,5,4) + t1(2,4,5,3) +
1275  t2(3,5,4,2) + t2(4,5,3,2) + t2(2,5,4,3) -
1276  t2(5,3,4,2) - t2(5,4,3,2) - t2(5,2,4,3)));
1277 
1278  // Create the hadronic current for the 3 charged pion channel.
1279  else if (abs(pID[3]) == 211)
1280  u2.push_back(G(2,s) * (t1(3,5,4,2) + t1(4,5,3,2) + t1(3,4,5,2) +
1281  t1(4,3,5,2) + t1(2,4,3,5) + t1(2,3,4,5) +
1282  t2(2,4,3,5) + t2(2,3,4,5) -
1283  t2(3,2,4,5) - t2(4,2,3,5)) +
1284  G(3,s) * (t3(3,5,4,2) + t3(4,5,3,2) - t3(3,4,5,2) -
1285  t3(4,3,5,2) - t3(3,2,4,5) - t3(4,2,3,5)));
1286 
1287  // Push back the hadronic current (remove temporary momenta first).
1288  u.pop_back(); u.push_back(u2);
1289 
1290 }
1291 
1292 //--------------------------------------------------------------------------
1293 
1294 // Return the first t-vector.
1295 
1296 Wave4 HMETau2FourPions::t1(int q1, int q2, int q3, int q4) {
1297 
1298  Wave4 a1Q(u[2][q2] + u[2][q3] + u[2][q4]);
1299  Wave4 rhoQ(u[2][q3] + u[2][q4]);
1300  double a1S = m2(a1Q, a1Q);
1301  double rhoS = m2(rhoQ, rhoQ);
1302 
1303  // Needed to match Herwig++.
1304  double gM = sqrtpos(rhoM*rhoM - 4*picM*picM) * (rhoM*rhoM - 4*picM*picM)
1305  / rhoM;
1306  double dm = (rhoFormFactor1(0) - rhoFormFactor1(rhoM*rhoM) +
1307  rhoM*rhoM * rhoFormFactor2(rhoM*rhoM)) / gM;
1308  return - a1FormFactor(a1S) / (a1D(a1S) * rhoD(rhoS)) *
1309  pow2(a1M) * (rhoM*rhoM + rhoM*rhoG*dm) *
1310  (m2(u[2][0],a1Q) * (m2(u[2][q3],a1Q) * u[2][q4] -
1311  m2(u[2][q4],a1Q) * u[2][q3]) +
1312  (m2(u[2][0],u[2][q4]) * m2(u[2][q1],u[2][q3]) -
1313  m2(u[2][0],u[2][q3]) * m2(u[2][q1],u[2][q4])) * a1Q);
1314 
1315 }
1316 
1317 //--------------------------------------------------------------------------
1318 
1319 // Return the second t-vector.
1320 
1321 Wave4 HMETau2FourPions::t2(int /*q1*/, int q2, int q3, int q4) {
1322 
1323  Wave4 a1Q(u[2][q2] + u[2][q3] + u[2][q4]);
1324  Wave4 sigQ(u[2][q3] + u[2][q4]);
1325  double a1S = m2(a1Q, a1Q);
1326  double sigS = m2(sigQ, sigQ);
1327  return sigW * a1FormFactor(a1S) / (a1D(a1S) * sigD(sigS)) *
1328  pow2(a1M) * pow2(sigM) *
1329  (m2(u[2][0],a1Q) * a1S * u[2][q2] - m2(u[2][0],u[2][q2]) * a1S * a1Q);
1330 
1331 }
1332 
1333 //--------------------------------------------------------------------------
1334 
1335 // Return the third t-vector.
1336 
1337 Wave4 HMETau2FourPions::t3(int q1, int q2, int q3, int q4) {
1338  Wave4 omeQ(u[2][q2] + u[2][q3] + u[2][q4]);
1339  Wave4 rhoQ(u[2][q3] + u[2][q4]);
1340  double omeS = m2(omeQ, omeQ);
1341  double rhoS = m2(rhoQ, rhoQ);
1342 
1343  // Needed to match Herwig++.
1344  double gM = sqrtpos(rhoM*rhoM - 4*picM*picM) * (rhoM*rhoM - 4*picM*picM)
1345  / rhoM;
1346  double dm = (rhoFormFactor1(0) - rhoFormFactor1(rhoM*rhoM) +
1347  rhoM*rhoM * rhoFormFactor2(rhoM*rhoM)) / gM;
1348  return omeW * omeFormFactor(omeS) / (omeD(omeS) * rhoD(rhoS)) *
1349  pow2(omeM) * (rhoM*rhoM + rhoM*rhoG*dm) *
1350  ((m2(u[2][0],u[2][q3]) * m2(u[2][q1],u[2][q4]) -
1351  m2(u[2][0],u[2][q4]) * m2(u[2][q1],u[2][q3])) * u[2][q2] +
1352  (m2(u[2][0],u[2][q4]) * m2(u[2][q1],u[2][q2]) -
1353  m2(u[2][0],u[2][q2]) * m2(u[2][q1],u[2][q4])) * u[2][q3] +
1354  (m2(u[2][0],u[2][q2]) * m2(u[2][q1],u[2][q3]) -
1355  m2(u[2][0],u[2][q3]) * m2(u[2][q1],u[2][q2])) * u[2][q4]);
1356 
1357 }
1358 
1359 //--------------------------------------------------------------------------
1360 
1361 // Return the D function for the a1(1260).
1362 
1363 complex HMETau2FourPions::a1D(double s) {
1364 
1365  // rG is defined as the running width.
1366  double rG = 0;
1367 
1368  // The rho and pion cut off thresholds defined in the fit.
1369  double piM = 0.16960;
1370  double rM = 0.83425;
1371 
1372  // Fit of width below three pion mass threshold.
1373  if (s < piM)
1374  rG = 0;
1375 
1376  // Fit of width below pion and rho mass threshold.
1377  else if (s < rM)
1378  rG = 0.003052*pow3(s - piM)*(1.0 + 151.088*(s - piM) +
1379  174.495*pow2(s - piM));
1380 
1381  // Fit of width above pion and rho mass threshold.
1382  else
1383  rG = 2.60817 - 2.47790*s + 0.66539*pow2(s) - 0.0678183*pow3(s) +
1384  1.66577*(s-1.23701)/s;
1385  return s - a1M*a1M + complex(0,1) * sqrtpos(s) * rG;
1386 
1387 }
1388 
1389 //--------------------------------------------------------------------------
1390 
1391 // Return the D function for the rho(770).
1392 
1393 complex HMETau2FourPions::rhoD(double s) {
1394 
1395  double gQ = sqrtpos(s - 4*picM*picM) * (s - 4*picM*picM) / sqrtpos(s);
1396  double gM = sqrtpos(rhoM*rhoM - 4*picM*picM) * (rhoM*rhoM - 4*picM*picM)
1397  / rhoM;
1398  double dm = (rhoFormFactor1(s) - rhoFormFactor1(rhoM*rhoM) -
1399  (s - rhoM*rhoM) * rhoFormFactor2(rhoM*rhoM)) / gM;
1400 
1401  // Ensure complex part is zero below available channel.
1402  if (s < 4*picM*picM) gQ = 0;
1403  return s - rhoM*rhoM - rhoM*rhoG*dm + complex(0,1)*rhoM*rhoG*(gQ/gM);
1404 
1405 }
1406 
1407 //--------------------------------------------------------------------------
1408 
1409 // Return the D function for the sigma(800) (just s-wave running width).
1410 
1411 complex HMETau2FourPions::sigD(double s) {
1412 
1413  // Sigma decay to two neutral pions for three neutral pion channel.
1414  double piM = abs(pID[3]) == 111 ? pinM : picM;
1415  double gQ = sqrtpos(1.0 - 4*piM*piM/s);
1416  double gM = sqrtpos(1.0 - 4*piM*piM/(sigM*sigM));
1417  return s - sigM*sigM + complex(0,1)*sigM*sigG*gQ/gM;
1418 
1419 }
1420 
1421 //--------------------------------------------------------------------------
1422 
1423 // Return the D function for the omega(782).
1424 
1425 complex HMETau2FourPions::omeD(double s) {
1426 
1427  double g = 0;
1428  double q = sqrtpos(s);
1429  double x = q - omeM;
1430 
1431  // Fit of width given in TAUOLA.
1432  if (s < 1)
1433  g = 1 + 17.560*x + 141.110*pow2(x) + 894.884*pow3(x) + 4977.35*pow4(x) +
1434  7610.66*pow5(x) - 42524.4*pow6(x);
1435  else
1436  g = -1333.26 + 4860*q - 6000.81*pow2(q) + 2504.97*pow3(q);
1437  if (g < 0) g = 0;
1438  return s - omeM*omeM + complex(0,1)*omeM*omeG*g;
1439 
1440 }
1441 
1442 //--------------------------------------------------------------------------
1443 
1444 // Return the form factor for the a1.
1445 
1446 double HMETau2FourPions::a1FormFactor(double s) {
1447 
1448  return pow2((1.0 + a1M*a1M/lambda2) / (1.0 + s/lambda2));
1449 
1450 }
1451 
1452 //--------------------------------------------------------------------------
1453 
1454 // Return the form factor for the rho(770) (equivalent to h(s) in TAUOLA).
1455 
1456 double HMETau2FourPions::rhoFormFactor1(double s) {
1457 
1458  double f = sqrtpos(1 - 4*picM*picM/s);
1459  if (s > 4*picM*picM)
1460  f = f * log((1 + f) / (1 - f)) * (s - 4*picM*picM) / M_PI;
1461  else if (s < 0.0000001)
1462  f = -8 * picM*picM / M_PI;
1463  else
1464  f = 0;
1465  return f;
1466 
1467 }
1468 
1469 //--------------------------------------------------------------------------
1470 
1471 // Return the form factor for the rho(770) (equivalent to h(s) derivative).
1472 
1473 double HMETau2FourPions::rhoFormFactor2(double s) {
1474 
1475  double f = sqrtpos(1 - 4*picM*picM/s);
1476  if (s > 4*picM*picM)
1477  f = f / (M_PI * s) * (s*f + (2*picM*picM + s)*log((1 + f) / (1 - f)));
1478  else
1479  f = 0;
1480  return f;
1481 
1482 }
1483 
1484 //--------------------------------------------------------------------------
1485 
1486 // Return the form factor for the omega(782).
1487 
1488 double HMETau2FourPions::omeFormFactor(double /*s*/) {
1489 
1490  return 1.0;
1491 
1492 }
1493 
1494 //--------------------------------------------------------------------------
1495 
1496 // Return the G-functions given in TAUOLA using a piece-wise fit.
1497 
1498 double HMETau2FourPions::G(int i, double s) {
1499 
1500  // Break points for the fits.
1501  double s0(0), s1(0), s2(0), s3(0), s4(0), s5(0);
1502 
1503  // Parameters for the fits.
1504  double a1(0), b1(0);
1505  double a2(0), b2(0), c2(0), d2(0), e2(0);
1506  double a3(0), b3(0), c3(0), d3(0), e3(0);
1507  double a4(0), b4(0);
1508  double a5(0), b5(0);
1509 
1510  // Three neutral pion parameters.
1511  if (i == 1) {
1512  s0 = 0.614403; s1 = 0.656264; s2 = 1.57896;
1513  s3 = 3.08198; s4 = 3.12825; s5 = 3.17488;
1514  a1 = -23383.7; b1 = 38059.2;
1515  a2 = 230.368; b2 = -4.39368; c2 = 687.002;
1516  d2 = -732.581; e2 = 207.087;
1517  a3 = 1633.92; b3 = -2596.21; c3 = 1703.08;
1518  d3 = -501.407; e3 = 54.5919;
1519  a4 = -2982.44; b4 = 986.009;
1520  a5 = 6948.99; b5 = -2188.74;
1521  }
1522 
1523  // Three charged pion parameters.
1524  else if (i == 2) {
1525  s0 = 0.614403; s1 = 0.635161; s2 = 2.30794;
1526  s3 = 3.08198; s4 = 3.12825; s5 = 3.17488;
1527  a1 = -54171.5; b1 = 88169.3;
1528  a2 = 454.638; b2 = -3.07152; c2 = -48.7086;
1529  d2 = 81.9702; e2 = -24.0564;
1530  a3 = -162.421; b3 = 308.977; c3 = -27.7887;
1531  d3 = -48.5957; e3 = 10.6168;
1532  a4 = -2650.29; b4 = 879.776;
1533  a5 = 6936.99; b5 = -2184.97;
1534  }
1535 
1536  // Omega mediated three charged pion parameters.
1537  else if (i == 3) {
1538  s0 = 0.81364; s1 = 0.861709; s2 = 1.92621;
1539  s3 = 3.08198; s4 = 3.12825; s5 = 3.17488;
1540  a1 = -84888.9; b1 = 104332;
1541  a2 = 2698.15; b2 = -3.08302; c2 = 1936.11;
1542  d2 = -1254.59; e2 = 201.291;
1543  a3 = 7171.65; b3 = -6387.9; c3 = 3056.27;
1544  d3 = -888.63; e3 = 108.632;
1545  a4 = -5607.48; b4 = 1917.27;
1546  a5 = 26573; b5 = -8369.76;
1547  }
1548 
1549  // Return the appropriate fit.
1550  if (s < s0)
1551  return 0.0;
1552  else if (s < s1)
1553  return a1 + b1*s;
1554  else if (s < s2)
1555  return a2*pow(s,b2) + c2*pow2(s) + d2*pow3(s) + e2*pow4(s);
1556  else if (s < s3)
1557  return a3 + b3*s + c3*pow2(s) + d3*pow3(s) + e3*pow4(s);
1558  else if (s < s4)
1559  return a4 + b4*s;
1560  else if (s < s5)
1561  return a5 + b5*s;
1562  else
1563  return 0.0;
1564 
1565 }
1566 
1567 //==========================================================================
1568 
1569 } // end namespace Pythia8
Generated by   doxygen 1.8.5