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