| File: | generators/kkmc/tauola/IBSD_LHN.cc |
| Warning: | line 59, column 16 Value stored to 'PnuPpi_sq' 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_left_(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_LL, const double& CV_LR, const double& CS_RL, const double& CS_RR, |
| 10 | const double *ptau, const double *pnu, const double *ppi, const double *k, double &omega, double *hj); |
| 11 | |
| 12 | |
| 13 | void pigamma_left_(const double& Mtau, const double& Mpi,const double& m_rho, const double& Gamma_rho,const double& m_a1,const double& Gamma_a1, |
| 14 | const double& CV_LL, const double& CV_LR, const double& CS_RL, const double& CS_RR, |
| 15 | const double* ptau, const double* pnu, const double* ppi, const double* k, double& omega, double* hj) |
| 16 | { |
| 17 | |
| 18 | // Physical constants |
| 19 | // const double Mtau = 1.777; |
| 20 | // const double Mpi = 0.139568; |
| 21 | // const double m_rho = 0.7749; // PRD 78, 2008, 072006 |
| 22 | // const double Gamma_rho = 0.1486; |
| 23 | // const double m_a1 = 1.251; |
| 24 | // const double Gamma_a1 = 0.599; |
| 25 | const double M_Borel = 3.35; // arXiv:2010.00549 [hep-ph] |
| 26 | const double M_u = 2.16e-3; // Mass of up quark in GeV |
| 27 | const double M_d = 4.67e-3; // Mass of down quark in GeV |
| 28 | const double f_pi = 0.092; |
| 29 | |
| 30 | // Wilson coefficients (real as requested) |
| 31 | // const double CV_LL = 0.0; |
| 32 | // const double CV_LR = 0.0; |
| 33 | // const double CS_RL = 0.0; |
| 34 | // const double CS_RR = 0.0; |
| 35 | const std::complex<double> I(0.0, 1.0); |
| 36 | |
| 37 | // Compute scalar products using Minkowski metric |
| 38 | const double PpiPtau = ptau[3] * ppi[3] - ptau[2] * ppi[2] - ptau[1] * ppi[1] - ptau[0] * ppi[0]; |
| 39 | const double PnuPpi = pnu[3] * ppi[3] - pnu[2] * ppi[2] - pnu[1] * ppi[1] - pnu[0] * ppi[0]; |
| 40 | const double PnuPtau = ptau[3] * pnu[3] - ptau[2] * pnu[2] - ptau[1] * pnu[1] - ptau[0] * pnu[0]; |
| 41 | const double PkPtau = k[3] * ptau[3] - k[2] * ptau[2] - k[1] * ptau[1] - k[0] * ptau[0]; |
| 42 | const double PkPpi = k[3] * ppi[3] - k[2] * ppi[2] - k[1] * ppi[1] - k[0] * ppi[0]; |
| 43 | const double PkPnu = k[3] * pnu[3] - k[2] * pnu[2] - k[1] * pnu[1] - k[0] * pnu[0]; |
| 44 | |
| 45 | 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]); |
| 46 | |
| 47 | // F_V^π(t) |
| 48 | 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)); |
| 49 | // F_A^π(t) |
| 50 | 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)); |
| 51 | |
| 52 | const double Mtau_sq = Mtau * Mtau; |
| 53 | const double Mpi_sq = Mpi * Mpi; |
| 54 | const double fpi_sq = (f_pi) * (f_pi); |
| 55 | const double PpiPtau_sq = PpiPtau * PpiPtau; |
| 56 | const double PkPtau_sq = PkPtau * PkPtau; |
| 57 | const double PkPpi_sq = PkPpi * PkPpi; |
| 58 | const double PkPnu_sq = PkPnu * PkPnu; |
| 59 | const double PnuPpi_sq = PnuPpi * PnuPpi; |
Value stored to 'PnuPpi_sq' during its initialization is never read | |
| 60 | const double PnuPtau_sq = PnuPtau * PnuPtau; |
| 61 | const double Mtau_3 = Mtau_sq * Mtau; |
| 62 | |
| 63 | // Compute decay amplitude squared |A|^2 |
| 64 | // ----------------------------------------------------------------------------------------------------- |
| 65 | // |
| 66 | // |A_L|^2 = |A_L^IB|^2 + |A_L^V|^2 + |A_L^A|^2 |
| 67 | // + 2Re(A_L^V A_L^A*) + 2Re(A_L^IB A_L^V*) + 2Re(A_L^IB A_L^A*) |
| 68 | // |
| 69 | // ----------------------------------------------------------------------------------------------------- |
| 70 | |
| 71 | // common combination of wilson coefficients and form factors |
| 72 | |
| 73 | |
| 74 | const double mass_coeff = Mpi_sq / (Mtau * (M_u + M_d)); |
| 75 | const double CL_sum = 1.0 + CV_LL + CV_LR; |
| 76 | const double CL_diff = 1.0 + CV_LL - CV_LR; |
| 77 | const double CL_IB = CL_diff + mass_coeff * (CS_RL - CS_RR); |
| 78 | |
| 79 | // Declare variables |
| 80 | double p1[3], p2[3], p3[3]; |
| 81 | double h1[3], h2[3], h3[3], h4[3], h5[3], h6[3]; |
| 82 | |
| 83 | // Compute omega, h_j for |A_L^IB|^2 |
| 84 | |
| 85 | |
| 86 | // \omega = |CL_IB|^2 f^2_\pi m_\tau^2 \left\{-m_\pi^2\frac{p_\nu \cdot p_\tau}{(k\cdot p_\pi)^2} |
| 87 | // + \frac{(k\cdot p_\tau + m_\tau^2)(k\cdot p_\nu - p_\nu\cdot p_\tau) }{(k\cdot p_\tau)^2} + \frac{2 \,(p_\nu\cdot p_\tau) (p_\pi\cdot p_\tau) + (k\cdot p_\tau) (p_\nu \cdot p_\pi) - (k\cdot p_\nu) (p_\pi \cdot p_\tau)}{(k\cdot p_\pi) (k\cdot p_\tau)}\right\}, |
| 88 | // |
| 89 | double omega_1 = (CL_IB * CL_IB) * fpi_sq * Mtau_sq * (-Mpi_sq * (PnuPtau / PkPpi_sq) + ((PkPtau + Mtau_sq) * (PkPnu - PnuPtau) / PkPtau_sq) + |
| 90 | ((2.0 * PnuPtau * PpiPtau + PkPtau * PnuPpi - PkPnu * PpiPtau) / (PkPpi * PkPtau))); |
| 91 | |
| 92 | // Compute h for each spatial component j = 1,2,3 |
| 93 | for (int j = 0; j < 3; j++) { |
| 94 | // p_\nu^j \left(\frac{m_\pi^2}{(k\cdot p_\pi)^2} + \frac{m_\tau^2}{(k\cdot p_\tau)^2} - |
| 95 | // \frac{2 \,p_\pi\cdot p_\tau}{(k\cdot p_\pi) (k\cdot p_\tau)}\right) |
| 96 | p1[j] = pnu[j] * (Mpi_sq / PkPpi_sq + Mtau_sq / PkPtau_sq - 2.0 * PpiPtau / (PkPpi * PkPtau)); |
| 97 | // k^j \left(\frac{p_\nu\cdot p_\tau - k\cdot p_\nu}{(k\cdot p_\tau)^2} |
| 98 | // - \frac{p_\nu \cdot p_\pi}{(k\cdot p_\pi )(k\cdot p_\tau)}\right) |
| 99 | p2[j] = k[j] * ((PnuPtau - PkPnu) / PkPtau_sq - PnuPpi / (PkPpi * PkPtau)); |
| 100 | // p_\pi^j \frac{k\cdot p_\nu}{(k\cdot p_\pi)(k\cdot p_\tau)} |
| 101 | p3[j] = ppi[j] * PkPnu / (PkPpi * PkPtau); |
| 102 | |
| 103 | h1[j] = (-(CL_IB * CL_IB) * fpi_sq * Mtau_3 / omega_1) * (p1[j] + p2[j] + p3[j]); |
| 104 | } |
| 105 | |
| 106 | // compute omega, h_j for |A_L^V|^2 |
| 107 | |
| 108 | // \omega = 2 |(1+C^V_{LL} + C^V_{LR})F_V^{(\pi)}|^2 \left[(k\cdot p_\nu)(k \cdot p_\pi)(p_\pi\cdot p_\tau) |
| 109 | // - m_\pi^2 (k\cdot p_\nu)(k \cdot p_\tau) + (k\cdot p_\pi)(k\cdot p_\tau)(p_\nu\cdot p_\pi) \right], |
| 110 | |
| 111 | const double common_sd = (PkPnu * PkPpi * PpiPtau) - (Mpi_sq * PkPnu * PkPtau) + (PkPpi * PkPtau * PnuPpi); |
| 112 | double omega_2 = 2.0 * (CL_sum * CL_sum) * std::norm(FV_pi) * common_sd; |
| 113 | |
| 114 | // Compute h for each spatial component j = 1,2,3 |
| 115 | for (int j = 0; j < 3; j++) { |
| 116 | // k^j \left\{m_\pi^2 (k\cdot p_\nu) - (k\cdot p_\pi)(p_\nu\cdot p_\pi)\right\} |
| 117 | p1[j] = k[j] * (Mpi_sq * PkPnu - PkPpi * PnuPpi); // j-1 ----> J |
| 118 | // p_\pi^j (k\cdot p_\nu)(k\cdot p_\pi) |
| 119 | p2[j] = ppi[j] * PkPnu * PkPpi; |
| 120 | // 2 |(1+C^V_{LL} + C^V_{LR})F_V^{(\pi)}|^2 \,\frac{m_\tau}{\omega}[p1[j] - p2[j]] |
| 121 | |
| 122 | h2[j] = (2.0 * (CL_sum * CL_sum) * std::norm(FV_pi) * Mtau / omega_2) * (p1[j] - p2[j]); |
| 123 | } |
| 124 | |
| 125 | // compute omega, h_j for |A_L^A|^2 |
| 126 | |
| 127 | // \omega = 2 |(1+C^V_{LL} - C^V_{LR})F_A^{(\pi)}|^2 \left[(k\cdot p_\nu)(k \cdot p_\pi)(p_\pi\cdot p_\tau) - |
| 128 | // m_\pi^2 (k\cdot p_\nu)(k \cdot p_\tau) + (k\cdot p_\pi)(k\cdot p_\tau)(p_\nu\cdot p_\pi) \right], |
| 129 | |
| 130 | double omega_3 = 2.0 * (CL_diff * CL_diff) * std::norm(FA_pi) * common_sd; |
| 131 | |
| 132 | // Compute h for each spatial component j = 1,2,3 |
| 133 | |
| 134 | for (int j = 0; j < 3; j++) { |
| 135 | // k^j \left\{m_\pi^2 (k\cdot p_\nu) - (k\cdot p_\pi)(p_\nu\cdot p_\pi)\right\} |
| 136 | p1[j] = k[j] * (Mpi_sq * PkPnu - PkPpi * PnuPpi); |
| 137 | // p_\pi^j (k\cdot p_\nu)(k\cdot p_\pi) |
| 138 | p2[j] = ppi[j] * PkPnu * PkPpi; |
| 139 | // h^j = 2 |(1+C^V_{LL} - C^V_{LR})F_A^{(\pi)}|^2 \,\frac{m_\tau}{\omega} [p1[j]- p2[j]] |
| 140 | |
| 141 | h3[j] = (2.0 * (CL_diff * CL_diff) * std::norm(FA_pi) * Mtau / omega_3) * (p1[j] - p2[j]); |
| 142 | } |
| 143 | |
| 144 | // compute interference terms 2Re(A_L^V A_L^A*) |
| 145 | |
| 146 | // \omega = 4 \operatorname{Re}[F_V^{(\pi)}F_A^{(\pi)\ast} (1+C^V_{LL} + C^V_{LR})(1+C^V_{LL} - C^V_{LR})^\ast] |
| 147 | // \left[(k\cdot p_\nu)(k \cdot p_\pi)(p_\pi\cdot p_\tau) - (k\cdot p_\pi)(k\cdot p_\tau)(p_\nu\cdot p_\pi) \right], |
| 148 | |
| 149 | const double re_va = CL_sum * CL_diff * std::real(FV_pi * std::conj(FA_pi)); |
| 150 | const double common_va = (PkPnu * PkPpi * PpiPtau) - (PkPpi * PkPtau * PnuPpi); |
| 151 | double omega_4 = 4.0 * re_va * common_va; |
| 152 | |
| 153 | // Compute each spatial component j = 1,2,3 |
| 154 | for (int j = 0; j < 3; j++) { |
| 155 | // k^j (k\cdot p_\pi)(p_\nu\cdot p_\pi) |
| 156 | p1[j] = k[j] * PkPpi * PnuPpi; |
| 157 | // p_\pi^j (k\cdot p_\nu)(k\cdot p_\pi) |
| 158 | p2[j] = ppi[j] * PkPnu * PkPpi; |
| 159 | |
| 160 | // 4 \operatorname{Re}[F_V^{(\pi)}F_A^{(\pi)\ast} (1+C^V_{LL} + C^V_{LR})(1+C^V_{LL} - C^V_{LR})^\ast]\frac{m_\tau}{\omega} |
| 161 | // [p1[j] - p2[j]] |
| 162 | |
| 163 | h4[j] = (4.0 * re_va * Mtau / omega_4) * (p1[j] - p2[j]); |
| 164 | } |
| 165 | |
| 166 | // compute interference terms 2Re(A_L^IB A_L^V*) |
| 167 | |
| 168 | // \omega = 2\operatorname{Re} [F_V^{(\pi)} (1+C^V_{LL} - C^V_{LR})(1+C^V_{LL} + C^V_{LR})^\ast] f_\pi m_\tau^2 |
| 169 | // \frac{(k\cdot p_\nu)(k\cdot p_\pi)}{k\cdot p_\tau}, |
| 170 | |
| 171 | const double re_iv = CL_diff * CL_IB * std::real(FV_pi); |
| 172 | double omega_5 = 2.0 * re_iv * f_pi * Mtau_sq * (PkPnu * PkPpi / PkPtau); |
| 173 | |
| 174 | // Compute each spatial component j = 1,2,3 |
| 175 | |
| 176 | for (int j = 0; j < 3; j++) { |
| 177 | // k^j\left\{(p_\pi\cdot p_\tau)(k\cdot p_\nu) - (p_\pi\cdot p_\tau)(p_\nu\cdot p_\tau) + m_\tau^2\, p_\nu\cdot p_\pi\right\} |
| 178 | p1[j] = k[j] * ((PpiPtau * PkPnu) - (PpiPtau * PnuPtau) + (Mtau_sq * PnuPpi)); |
| 179 | // p_\pi^j \left\{(k\cdot p_\nu)(k\cdot p_\tau) + m_\tau^2\, k\cdot p_\nu - (k\cdot p_\tau) (p_\nu\cdot p_\tau)\right\} |
| 180 | p2[j] = ppi[j] * ((PkPnu * PkPtau) + Mtau_sq * PkPnu - PkPtau * PnuPtau); |
| 181 | |
| 182 | // h^j = \frac{2\operatorname{Re} [F_V^{(\pi)} (1+C^V_{LL} - C^V_{LR})(1+C^V_{LL} + C^V_{LR})^\ast]f_\pi}{k\cdot p_\tau}\frac{m_\tau}{\omega}[ p1[j] - p2[j]] |
| 183 | |
| 184 | h5[j] = ((2.0 * re_iv * f_pi * Mtau) / (PkPtau * omega_5)) * (p1[j] - p2[j]); |
| 185 | } |
| 186 | |
| 187 | // compute interference terms 2Re(A_L^IB A_L^A*) |
| 188 | |
| 189 | // \omega =\frac{ 2 f_\pi m_\tau^2 |1+C^V_{LL} - C^V_{LR}|^2 \operatorname{Re}(F_A^{(\pi)}) }{(k\cdot p_\pi)(k\cdot p_\tau)} |
| 190 | // \left[ (k\cdot p_\nu)\left\{(k\cdot p_\pi)^2 + (k\cdot p_\pi)(p_\pi \cdot p_\tau) - m_\pi^2\, k\cdot p_\tau\right\} |
| 191 | // + (k\cdot p_\pi)(k\cdot p_\tau)(p_\nu\cdot p_\pi) - (k\cdot p_\pi)^2 (p_\nu \cdot p_\tau)\right], |
| 192 | |
| 193 | const double re_ia = CL_diff * CL_IB * std::real(FA_pi); |
| 194 | double omega_6 = ((2.0 * f_pi * Mtau_sq * re_ia) / (PkPpi * PkPtau)) * |
| 195 | ((PkPnu * ( PkPpi_sq + (PkPpi * PpiPtau) - (Mpi_sq * PkPtau))) + (PkPpi * PkPtau * PnuPpi) - ( PkPpi_sq * PnuPtau)); |
| 196 | |
| 197 | // Compute each spatial component j = 1,2,3 |
| 198 | |
| 199 | for (int j = 0; j < 3; j++) { |
| 200 | // k^j\left\{-m_\pi^2\, (k\cdot p_\tau) (p_\nu\cdot p_\tau) + (k\cdot p_\pi)(p_\pi\cdot p_\tau)(p_\nu \cdot p_\tau - k\cdot p_\nu)\right\} |
| 201 | p1[j] = k[j] * (-Mpi_sq * PkPtau * PnuPtau + PkPpi * PpiPtau * (PnuPtau - PkPnu)); |
| 202 | // p_\nu^j \left\{m_\pi^2 (k\cdot p_\tau)^2 - 2(k\cdot p_\tau)(k\cdot p_\pi)(p_\pi\cdot p_\tau) + m_\tau^2\, (k\cdot p_\pi)(k\cdot p_\pi) \right\} |
| 203 | p2[j] = pnu[j] * (Mpi_sq * PkPtau_sq - 2.0 * PkPtau * PkPpi * PpiPtau + Mtau_sq * PkPpi_sq); |
| 204 | // p_\pi^j \left\{(k\cdot p_\tau)(k\cdot p_\pi)(k\cdot p_\nu + p_\nu \cdot p_\tau)\right\} |
| 205 | p3[j] = ppi[j] * (PkPtau * PkPpi * (PkPnu + PnuPtau)); |
| 206 | |
| 207 | // h^j = -\frac{2 f_\pi \operatorname{Re}(F_A^{(\pi)} CL_diff CL_IB^*)}{(k\cdot p_\pi)(k\cdot p_\tau)}\frac{m_\tau}{\omega} [p1[j] + p2[j] + p3[j]] |
| 208 | |
| 209 | h6[j] = ((-2.0 * f_pi * re_ia * Mtau) / (PkPpi * PkPtau * omega_6)) * (p1[j] + p2[j] + p3[j]); |
| 210 | } |
| 211 | |
| 212 | // Sum total amplitude |
| 213 | omega = omega_1 + omega_2 + omega_3 + omega_4 + omega_5 + omega_6; |
| 214 | |
| 215 | for (int j = 0; j < 3 ; j++ ){ |
| 216 | 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; |
| 217 | } |
| 218 | |
| 219 | |
| 220 | return; |
| 221 | } |