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