| File: | generators/kkmc/tauola/IBSD_RHN.cc |
| Warning: | line 58, column 16 Value stored to 'Mtau_3' during its initialization is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 1 | #include <math.h> |
| 2 | #include <cstdlib> |
| 3 | #include <cstdio> |
| 4 | #include <iostream> |
| 5 | #include <cmath> |
| 6 | #include <complex> |
| 7 | |
| 8 | extern "C" void pigamma_right_(const double& Mtau, const double& Mpi, const double& m_rho, const double& Gamma_rho, const double& m_a1, const double& Gamma_a1, |
| 9 | const double& CV_RL, const double& CV_RR, const double& CS_LL, const double& CS_LR, const double& MnuR, |
| 10 | const double *ptau, const double *pnu, const double *ppi, const double *k, double &omega, double *hj); |
| 11 | |
| 12 | void pigamma_right_(const double& Mtau, const double& Mpi, const double& m_rho, const double& Gamma_rho, const double& m_a1, const double& Gamma_a1, |
| 13 | const double& CV_RL, const double& CV_RR, const double& CS_LL, const double& CS_LR, const double& MnuR, |
| 14 | const double* ptau, const double* pnu, const double* ppi, const double* k, double& omega, double* hj) |
| 15 | { |
| 16 | // Physical constants |
| 17 | // const double Mtau = 1.777; |
| 18 | // const double Mpi = 0.139568; |
| 19 | // const double m_rho = 0.7749; // PRD 78, 2008, 072006 |
| 20 | // const double Gamma_rho = 0.1486; |
| 21 | // const double m_a1 = 1.251; |
| 22 | // const double Gamma_a1 = 0.599; |
| 23 | const double M_Borel = 3.35; // arXiv:2010.00549 [hep-ph] |
| 24 | const double M_u = 2.16e-3; // Mass of up quark in GeV |
| 25 | const double M_d = 4.67e-3; // Mass of down quark in GeV |
| 26 | const double f_pi = 0.092; |
| 27 | |
| 28 | // Wilson coefficients (real as requested) |
| 29 | // const double CV_RL = 0.0; |
| 30 | // const double CV_RR = 0.0; |
| 31 | // const double CS_LL = 0.0; |
| 32 | // const double CS_LR = 0.0; |
| 33 | const std::complex<double> I(0.0, 1.0); |
| 34 | |
| 35 | // RH neutrino mass (set to zero by default; update if needed) |
| 36 | // const double MnuR = 0.0; |
| 37 | const double MnuR_sq = MnuR * MnuR; |
| 38 | |
| 39 | // Lorentz products with (+,-,-,-) |
| 40 | const double PpiPtau = ptau[3] * ppi[3] - ptau[2] * ppi[2] - ptau[1] * ppi[1] - ptau[0] * ppi[0]; |
| 41 | const double PnuPpi = pnu[3] * ppi[3] - pnu[2] * ppi[2] - pnu[1] * ppi[1] - pnu[0] * ppi[0]; |
| 42 | const double PnuPtau = ptau[3] * pnu[3] - ptau[2] * pnu[2] - ptau[1] * pnu[1] - ptau[0] * pnu[0]; |
| 43 | const double PkPtau = k[3] * ptau[3] - k[2] * ptau[2] - k[1] * ptau[1] - k[0] * ptau[0]; |
| 44 | const double PkPpi = k[3] * ppi[3] - k[2] * ppi[2] - k[1] * ppi[1] - k[0] * ppi[0]; |
| 45 | const double PkPnu = k[3] * pnu[3] - k[2] * pnu[2] - k[1] * pnu[1] - k[0] * pnu[0]; |
| 46 | |
| 47 | const double t = (k[3] + ppi[3]) * (k[3] + ppi[3]) - (k[2] + ppi[2]) * (k[2] + ppi[2]) - (k[1] + ppi[1]) * (k[1] + ppi[1]) - (k[0] + ppi[0]) * (k[0] + ppi[0]); |
| 48 | |
| 49 | // Eq. (50), (51) |
| 50 | const std::complex<double> FV_pi = (f_pi / (m_rho * m_rho - t - I * m_rho * Gamma_rho)) * std::exp((m_rho * m_rho) / (M_Borel * M_Borel)); |
| 51 | const std::complex<double> FA_pi = (f_pi / (m_a1 * m_a1 - t - I * m_a1 * Gamma_a1)) * std::exp((m_a1 * m_a1) / (M_Borel * M_Borel)); |
| 52 | |
| 53 | const double Mtau_sq = Mtau * Mtau; |
| 54 | const double Mpi_sq = Mpi * Mpi; |
| 55 | const double fpi_sq = f_pi * f_pi; |
| 56 | const double PkPtau_sq = PkPtau * PkPtau; |
| 57 | const double PkPpi_sq = PkPpi * PkPpi; |
| 58 | const double Mtau_3 = Mtau_sq * Mtau; |
Value stored to 'Mtau_3' during its initialization is never read | |
| 59 | |
| 60 | // Real-coefficient combinations used in Eqs. (33)-(44) |
| 61 | const double mass_coeff = Mpi_sq / (Mtau * (M_u + M_d)); |
| 62 | const double CIB_R = (CV_RL - CV_RR) + mass_coeff * (CS_LL - CS_LR); |
| 63 | const double kappa_plus = Mtau_sq * (CIB_R * CIB_R) + MnuR_sq * ((CV_RR - CV_RL) * (CV_RR - CV_RL)); |
| 64 | const double kappa_minus = Mtau_sq * (CIB_R * CIB_R) - MnuR_sq * ((CV_RR - CV_RL) * (CV_RR - CV_RL)); |
| 65 | const double CV_sum = CV_RL + CV_RR; |
| 66 | const double CV_diff = CV_RL - CV_RR; |
| 67 | |
| 68 | double p1[3], p2[3], p3[3]; |
| 69 | double h1[3], h2[3], h3[3], h4[3], h5[3], h6[3]; |
| 70 | |
| 71 | // Eq. (33) |
| 72 | const double omega_1 = fpi_sq * ( |
| 73 | ((2.0 * MnuR_sq * Mtau_sq * CV_diff * CIB_R) - kappa_plus * PnuPtau) * (Mpi_sq / PkPpi_sq) |
| 74 | + (kappa_plus * (PkPtau + Mtau_sq) * (PkPnu - PnuPtau) + 2.0 * MnuR_sq * Mtau_sq * Mtau_sq * CV_diff * CIB_R) / PkPtau_sq |
| 75 | + (kappa_plus * (2.0 * PnuPtau * PpiPtau + PkPtau * PnuPpi - PkPnu * PpiPtau) - 4.0 * MnuR_sq * Mtau_sq * CV_diff * CIB_R * PpiPtau) / (PkPpi * PkPtau) |
| 76 | ); |
| 77 | |
| 78 | // Eq. (34) |
| 79 | for (int j = 0; j < 3; ++j) { |
| 80 | p1[j] = pnu[j] * (Mpi_sq / PkPpi_sq + Mtau_sq / PkPtau_sq - 2.0 * PpiPtau / (PkPpi * PkPtau)); |
| 81 | p2[j] = k[j] * ((PnuPtau - PkPnu) / PkPtau_sq - PnuPpi / (PkPpi * PkPtau)); |
| 82 | p3[j] = ppi[j] * (PkPnu / (PkPpi * PkPtau)); |
| 83 | h1[j] = (fpi_sq * Mtau * kappa_minus / omega_1) * (p1[j] + p2[j] + p3[j]); |
| 84 | } |
| 85 | |
| 86 | const double common_sd = (PkPnu * PkPpi * PpiPtau) - (Mpi_sq * PkPnu * PkPtau) + (PkPpi * PkPtau * PnuPpi); |
| 87 | |
| 88 | // Eq. (35)-(36) |
| 89 | const double omega_2 = 2.0 * (CV_sum * CV_sum) * std::norm(FV_pi) * common_sd; |
| 90 | for (int j = 0; j < 3; ++j) { |
| 91 | p1[j] = k[j] * (Mpi_sq * PkPnu - PkPpi * PnuPpi); |
| 92 | p2[j] = ppi[j] * (PkPnu * PkPpi); |
| 93 | h2[j] = (-2.0 * (CV_sum * CV_sum) * std::norm(FV_pi) * Mtau / omega_2) * (p1[j] - p2[j]); |
| 94 | } |
| 95 | |
| 96 | // Eq. (37)-(38) |
| 97 | const double omega_3 = 2.0 * (CV_diff * CV_diff) * std::norm(FA_pi) * common_sd; |
| 98 | for (int j = 0; j < 3; ++j) { |
| 99 | p1[j] = k[j] * (Mpi_sq * PkPnu - PkPpi * PnuPpi); |
| 100 | p2[j] = ppi[j] * (PkPnu * PkPpi); |
| 101 | h3[j] = (-2.0 * (CV_diff * CV_diff) * std::norm(FA_pi) * Mtau / omega_3) * (p1[j] - p2[j]); |
| 102 | } |
| 103 | |
| 104 | // Eq. (39)-(40) |
| 105 | const double re_va = CV_sum * CV_diff * std::real(FV_pi * std::conj(FA_pi)); |
| 106 | const double omega_4 = -4.0 * re_va * ((PkPnu * PkPpi * PpiPtau) - (PkPpi * PkPtau * PnuPpi)); |
| 107 | for (int j = 0; j < 3; ++j) { |
| 108 | p1[j] = k[j] * (PkPpi * PnuPpi); |
| 109 | p2[j] = ppi[j] * (PkPnu * PkPpi); |
| 110 | h4[j] = (4.0 * re_va * Mtau / omega_4) * (p1[j] - p2[j]); |
| 111 | } |
| 112 | |
| 113 | // Eq. (41)-(42) |
| 114 | const double re_iv = CV_sum * std::real(std::conj(FV_pi)); |
| 115 | const double omega_5 = (-2.0 * re_iv * f_pi / PkPtau) * ((Mtau_sq * CIB_R * PkPnu) - (MnuR_sq * CV_diff * PkPtau)) * PkPpi; |
| 116 | for (int j = 0; j < 3; ++j) { |
| 117 | p1[j] = k[j] * (CIB_R * ((PpiPtau * PkPnu) - (PpiPtau * PnuPtau) + (Mtau_sq * PnuPpi)) + MnuR_sq * CV_diff * PkPpi); |
| 118 | p2[j] = ppi[j] * (CIB_R * ((PkPnu * PkPtau) + Mtau_sq * PkPnu - PkPtau * PnuPtau)); |
| 119 | h5[j] = (2.0 * re_iv * f_pi * Mtau / (PkPtau * omega_5)) * (p1[j] - p2[j]); |
| 120 | } |
| 121 | |
| 122 | // Eq. (43)-(44) |
| 123 | const double re_ia = CV_diff * std::real(std::conj(FA_pi)); |
| 124 | const double omega_6 = (2.0 * f_pi * re_ia / (PkPpi * PkPtau)) * ( |
| 125 | PkPpi_sq * (MnuR_sq * CV_diff * (PkPtau + Mtau_sq) + Mtau_sq * CIB_R * (PkPnu - PnuPtau)) |
| 126 | + Mpi_sq * PkPtau * (MnuR_sq * CV_diff * PkPtau - Mtau_sq * CIB_R * PkPnu) |
| 127 | + PkPpi * PkPtau * (Mtau_sq * CIB_R * PnuPpi - 2.0 * MnuR_sq * CV_diff * PpiPtau) |
| 128 | + Mtau_sq * CIB_R * PkPpi * PkPnu * PpiPtau |
| 129 | ); |
| 130 | for (int j = 0; j < 3; ++j) { |
| 131 | p1[j] = k[j] * ( |
| 132 | Mpi_sq * PkPtau * (MnuR_sq * CV_diff - CIB_R * PnuPtau) |
| 133 | - PkPpi * PpiPtau * (CIB_R * (PkPnu - PnuPtau) + MnuR_sq * CV_diff) |
| 134 | + MnuR_sq * CV_diff * PkPpi_sq |
| 135 | ); |
| 136 | p2[j] = pnu[j] * CIB_R * (Mpi_sq * PkPtau_sq - 2.0 * PkPtau * PkPpi * PpiPtau + Mtau_sq * PkPpi_sq); |
| 137 | p3[j] = ppi[j] * PkPtau * PkPpi * (CIB_R * (PkPnu + PnuPtau) - CV_diff * MnuR_sq); |
| 138 | h6[j] = (2.0 * f_pi * re_ia * Mtau / (PkPpi * PkPtau * omega_6)) * (p1[j] + p2[j] + p3[j]); |
| 139 | } |
| 140 | |
| 141 | omega = omega_1 + omega_2 + omega_3 + omega_4 + omega_5 + omega_6; |
| 142 | for (int j = 0; j < 3; ++j) { |
| 143 | hj[j] = (omega_1 * h1[j] + omega_2 * h2[j] + omega_3 * h3[j] + omega_4 * h4[j] + omega_5 * h5[j] + omega_6 * h6[j]) / omega; |
| 144 | } |
| 145 | } |