PYTHIA8/pythia8175/src/HelicityMatrixElements.cxx

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