a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
MPll.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2014 HEPfit Collaboration
3 *
4 *
5 * For the licensing terms see doc/COPYING.
6 */
7
8#include "StandardModel.h"
9#include "MPll.h"
10#include "std_make_vector.h"
11#include "gslpp_function_adapter.h"
12#include "F_1.h"
13#include "F_2.h"
14#include <gsl/gsl_sf.h>
15#include <boost/bind/bind.hpp>
16#include <limits>
17#include <TFitResult.h>
18using namespace boost::placeholders;
19
20MPll::MPll(const StandardModel& SM_i, QCD::meson meson_i, QCD::meson pseudoscalar_i, QCD::lepton lep_i)
21: mySM(SM_i), myF_1(new F_1()), myF_2(new F_2()),
22fplus_lat_cache(3, 0.),
23fT_lat_cache(3, 0.),
24f0_lat_cache(3, 0.),
25fplus_cache(2, 0.),
26fT_cache(2, 0.),
27k2_cache(2, 0.),
28SL_cache(2, 0.),
29N_cache(3, 0.),
30Ycache(2, 0.),
31H_V0cache(2, 0.),
32H_Scache(2, 0.),
33H_P_cache(4, 0.),
34Itree_cache(3, 0.),
35T_cache(5, 0.)
36{
37 lep = lep_i;
38 meson = meson_i;
39 pseudoscalar = pseudoscalar_i;
40 dispersion = true;
41 FixedWCbtos = false;
42 MPll_Lattice_flag = false;
43 MPll_GRvDV_flag = false;
44 MPll_DM_flag = false;
45 NeutrinoTree_flag = false;
46 mJ2 = 3.096 * 3.096;
47
48 I0_updated = 0;
49 I2_updated = 0;
50 I8_updated = 0;
51 Itree_updated = 0;
52
53 VL_updated = 0;
54 TL_updated = 0;
55 SL_updated = 0;
56
57 deltaTparpupdated = 0;
58 deltaTparmupdated = 0;
59
60 w_sigma = gsl_integration_cquad_workspace_alloc(100);
61 w_delta = gsl_integration_cquad_workspace_alloc(100);
62 w_sigmaTree = gsl_integration_cquad_workspace_alloc(100);
63 w_DTPPR = gsl_integration_cquad_workspace_alloc(100);
64
65 acc_Re_deltaC7_QCDF = gsl_interp_accel_alloc();
66 acc_Im_deltaC7_QCDF = gsl_interp_accel_alloc();
67 acc_Re_deltaC9_QCDF = gsl_interp_accel_alloc();
68 acc_Im_deltaC9_QCDF = gsl_interp_accel_alloc();
69
70 spline_Re_deltaC7_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
71 spline_Im_deltaC7_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
72 spline_Re_deltaC9_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
73 spline_Im_deltaC9_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
74}
75
76MPll::~MPll()
77{
78}
79
80std::vector<std::string> MPll::initializeMPllParameters()
81{
82 dispersion = mySM.getFlavour().getFlagUseDispersionRelation();
83 FixedWCbtos = mySM.getFlavour().getFlagFixedWCbtos();
84 MPll_Lattice_flag = mySM.getFlavour().getFlagMPll_FNALMILC();
85 MPll_GRvDV_flag = mySM.getFlavour().getFlagMPll_GRvDV();
86 MPll_DM_flag = mySM.getFlavour().getFlagMPll_DM();
87 NeutrinoTree_flag = mySM.getFlavour().getFlagNeutrinoTree();
88
89 if (pseudoscalar == StandardModel::K_P || pseudoscalar == StandardModel::K_0) {
90 if (MPll_Lattice_flag) mpllParameters = make_vector<std::string>()
91 << "b_0_fplus" << "b_1_fplus" << "b_2_fplus" << "m_fit_fplus_lat"
92 << "b_0_fT" << "b_1_fT" << "b_2_fT" << "m_fit_fT_lat"
93 << "b_0_f0" << "b_1_f0" << "b_2_f0" << "m_fit_f0_lat" ;
94 else if (MPll_GRvDV_flag) mpllParameters = make_vector<std::string>()
95 << "b_0_fplus" << "b_1_fplus" << "b_2_fplus" << "m_fit_fplus_lat"
96 << "b_0_fT" << "b_1_fT" << "b_2_fT" << "m_fit_fT_lat"
97 << "b_1_f0" << "b_2_f0" << "m_fit_f0_lat" ;
98 else if (MPll_DM_flag) mpllParameters = make_vector<std::string>()
99 << "b_0_fplus" << "b_1_fplus" << "b_2_fplus" << "m_fit_fplus_lat"
100 << "b_0_fT" << "b_1_fT" << "b_2_fT" << "m_fit_fT_lat"
101 << "b_1_f0" << "b_2_f0" ;
102 else mpllParameters = make_vector<std::string>()
103 << "r_1_fplus" << "r_2_fplus" << "m_fit2_fplus" << "r_1_fT" << "r_2_fT" << "m_fit2_fT" << "r_2_f0" << "m_fit2_f0";
104 } else {
105 std::stringstream out;
106 out << pseudoscalar;
107 throw std::runtime_error("MPll: pseudoscalar " + out.str() + " not implemented");
108 }
109
110 if (lep != QCD::NEUTRINO_1){
111 if (dispersion) {
112 mpllParameters.insert(mpllParameters.end(), { "r1_BK", "r2_BK", "deltaC9_BK", "phDC9_BK" });
113 } else {
114#if NFPOLARBASIS_MPLL
115 mpllParameters.insert(mpllParameters.end(), { "absh_0_MP", "argh_0_MP", "absh_1_MP", "argh_1_MP", "absh_2_MP", "argh_2_MP", "Delta_C7_U", "Delta_C9_U" });
116#else
117 mpllParameters.insert(mpllParameters.end(), { "reh_0_MP", "imh_0_MP", "reh_1_MP", "imh_1_MP", "reh_2_MP", "imh_2_MP", "Delta_C7_U", "Delta_C9_U" });
118#endif
119 }
120 }
121
122 if (FixedWCbtos)
123 if (lep != QCD::NEUTRINO_1) mpllParameters.insert(mpllParameters.end(), { "C7_SM", "C9_SM", "C10_SM" });
124 else mpllParameters.insert(mpllParameters.end(), { "CLnunu_SM" });
125 mySM.initializeMeson(meson);
126 mySM.initializeMeson(pseudoscalar);
127 return mpllParameters;
128}
129
130void MPll::updateParameters()
131{
132 if (!mySM.getFlavour().getUpdateFlag(meson, pseudoscalar, lep)) return;
133
134
135
136 GF = mySM.getGF();
137 ale = mySM.getAle();
138 if (lep == QCD::NEUTRINO_1){
139 Mlep = 0.;
140 }
141 else{
142 Mlep = mySM.getLeptons(lep).getMass();
143 }
144
145 MM = mySM.getMesons(meson).getMass();
146 MM3 = MM * MM * MM;
147 MP = mySM.getMesons(pseudoscalar).getMass();
148 Mb = mySM.getQuarks(QCD::BOTTOM).getMass(); // add the PS b mass
149 Mc = mySM.getQuarks(QCD::CHARM).getMass();
150 mb_pole = mySM.Mbar2Mp(Mb,QCD::BOTTOM); /* Conversion to pole mass*/
151 mc_pole = mySM.Mbar2Mp(Mc,QCD::CHARM); /* Conversion to pole mass*/
152 Ms = mySM.getQuarks(QCD::STRANGE).getMass();
153 MW = mySM.Mw();
154 lambda_t = mySM.getCKM().computelamt_s();
155 mu_b = mySM.getMub();
156 mu_h = sqrt(mu_b * .5); // From Beneke Neubert
157 width = mySM.getMesons(meson).computeWidth();
158 alpha_s_mub = mySM.Als(mu_b);
159
160 switch (pseudoscalar) {
161 case StandardModel::K_P:
162 case StandardModel::K_0:
163 if (MPll_Lattice_flag) {
164 b_0_fplus = mySM.getOptionalParameter("b_0_fplus");
165 b_1_fplus = mySM.getOptionalParameter("b_1_fplus");
166 b_2_fplus = mySM.getOptionalParameter("b_2_fplus");
167 m_fit2_fplus_lat = mySM.getOptionalParameter("m_fit_fplus_lat") * mySM.getOptionalParameter("m_fit_fplus_lat");
168 b_0_fT = mySM.getOptionalParameter("b_0_fT");
169 b_1_fT = mySM.getOptionalParameter("b_1_fT");
170 b_2_fT = mySM.getOptionalParameter("b_2_fT");
171 m_fit2_fT_lat = mySM.getOptionalParameter("m_fit_fT_lat") * mySM.getOptionalParameter("m_fit_fT_lat");
172 b_0_f0 = mySM.getOptionalParameter("b_0_f0");
173 b_1_f0 = mySM.getOptionalParameter("b_1_f0");
174 b_2_f0 = mySM.getOptionalParameter("b_2_f0");
175 m_fit2_f0_lat = mySM.getOptionalParameter("m_fit_f0_lat") * mySM.getOptionalParameter("m_fit_f0_lat");
176 } else if (MPll_GRvDV_flag) {
177 b_0_fplus = mySM.getOptionalParameter("b_0_fplus");
178 b_1_fplus = mySM.getOptionalParameter("b_1_fplus");
179 b_2_fplus = mySM.getOptionalParameter("b_2_fplus");
180 m_fit2_fplus_lat = mySM.getOptionalParameter("m_fit_fplus_lat") * mySM.getOptionalParameter("m_fit_fplus_lat");
181 b_0_fT = mySM.getOptionalParameter("b_0_fT");
182 b_1_fT = mySM.getOptionalParameter("b_1_fT");
183 b_2_fT = mySM.getOptionalParameter("b_2_fT");
184 m_fit2_fT_lat = mySM.getOptionalParameter("m_fit_fT_lat") * mySM.getOptionalParameter("m_fit_fT_lat");
185 b_0_f0 = b_0_fplus;
186 b_1_f0 = mySM.getOptionalParameter("b_1_f0");
187 b_2_f0 = mySM.getOptionalParameter("b_2_f0");
188 m_fit2_f0_lat = mySM.getOptionalParameter("m_fit_f0_lat") * mySM.getOptionalParameter("m_fit_f0_lat");
189 } else if (MPll_DM_flag) {
190 b_0_fplus = mySM.getOptionalParameter("b_0_fplus");
191 b_1_fplus = mySM.getOptionalParameter("b_1_fplus");
192 b_2_fplus = mySM.getOptionalParameter("b_2_fplus");
193 m_fit2_fplus_lat = mySM.getOptionalParameter("m_fit_fplus_lat") * mySM.getOptionalParameter("m_fit_fplus_lat");
194 b_0_fT = mySM.getOptionalParameter("b_0_fT");
195 b_1_fT = mySM.getOptionalParameter("b_1_fT");
196 b_2_fT = mySM.getOptionalParameter("b_2_fT");
197 m_fit2_fT_lat = mySM.getOptionalParameter("m_fit_fT_lat") * mySM.getOptionalParameter("m_fit_fT_lat");
198 b_1_f0 = mySM.getOptionalParameter("b_1_f0");
199 b_2_f0 = mySM.getOptionalParameter("b_2_f0");
200 b_0_f0 = fplus_DM(0.,b_0_fplus,b_1_fplus,b_2_fplus,m_fit2_fplus_lat)*phi0_DM(0) - b_1_f0*zeta_DM(0) - b_2_f0*zeta_DM(0)*zeta_DM(0);
201 } else {
202 r_1_fplus = mySM.getOptionalParameter("r_1_fplus");
203 r_2_fplus = mySM.getOptionalParameter("r_2_fplus");
204 m_fit2_fplus = mySM.getOptionalParameter("m_fit2_fplus");
205 r_1_fT = mySM.getOptionalParameter("r_1_fT");
206 r_2_fT = mySM.getOptionalParameter("r_2_fT");
207 m_fit2_fT = mySM.getOptionalParameter("m_fit2_fT");
208 r_2_f0 = mySM.getOptionalParameter("r_2_f0");
209 m_fit2_f0 = mySM.getOptionalParameter("m_fit2_f0");
210 }
211
212 if (pseudoscalar == StandardModel::K_P) spectator_charge = mySM.getQuarks(QCD::UP).getCharge();
213 else spectator_charge = mySM.getQuarks(QCD::DOWN).getCharge();
214
215 etaP = -1;
216 angmomP = 0.;
217
218 break;
219 default:
220 std::stringstream out;
221 out << pseudoscalar;
222 throw std::runtime_error("MPll: pseudoscalar " + out.str() + " not implemented");
223 }
224
225 if (lep != QCD::NEUTRINO_1){
226 if (!dispersion) {
227#if NFPOLARBASIS_MPLL
228 h_0 = gslpp::complex(mySM.getOptionalParameter("absh_0_MP"), mySM.getOptionalParameter("argh_0_MP"), true);
229 h_1 = gslpp::complex(mySM.getOptionalParameter("absh_1_MP"), mySM.getOptionalParameter("argh_1_MP"), true);
230 h_2 = gslpp::complex(mySM.getOptionalParameter("absh_2_MP"), mySM.getOptionalParameter("argh_2_MP"), true);
231
232 r_1 = 0.;
233 r_2 = 0.;
234 Delta_C9 = 0.;
235 exp_Phase = 0.;
236
237 Delta_C7_U = mySM.getOptionalParameter("Delta_C7_U");
238 Delta_C9_U = mySM.getOptionalParameter("Delta_C9_U");
239#else
240 h_0 = gslpp::complex(mySM.getOptionalParameter("reh_0_MP"), mySM.getOptionalParameter("imh_0_MP"), false);
241 h_1 = gslpp::complex(mySM.getOptionalParameter("reh_1_MP"), mySM.getOptionalParameter("imh_1_MP"), false);
242 h_2 = gslpp::complex(mySM.getOptionalParameter("reh_2_MP"), mySM.getOptionalParameter("imh_2_MP"), false);
243
244 r_1 = 0.;
245 r_2 = 0.;
246 Delta_C9 = 0.;
247 exp_Phase = 0.;
248
249 Delta_C7_U = mySM.getOptionalParameter("Delta_C7_U");
250 Delta_C9_U = mySM.getOptionalParameter("Delta_C9_U");
251#endif
252 } else {
253 h_0 = 0.;
254 h_1 = 0.;
255 h_2 = 0.;
256
257 r_1 = mySM.getOptionalParameter("r1_BK");
258 r_2 = mySM.getOptionalParameter("r2_BK");
259 Delta_C9 = mySM.getOptionalParameter("deltaC9_BK");
260 exp_Phase = exp(gslpp::complex::i() * mySM.getOptionalParameter("phDC9_BK"));
261
262 Delta_C7_U = 0.;
263 Delta_C9_U = 0.;
264 }
265 }
266
267 if (lep == QCD::NEUTRINO_1){
268 VusVub_abs2 = (mySM.getCKM().computelamu_s() * mySM.getCKM().computelamu_s().conjugate()).abs();
269 GF4 = GF * GF * GF * GF;
270 fM2 = mySM.getMesons(meson).getDecayconst() * mySM.getMesons(meson).getDecayconst();
271 fP2 = mySM.getMesons(pseudoscalar).getDecayconst() * mySM.getMesons(pseudoscalar).getDecayconst();
272 mtau = mySM.getLeptons(QCD::TAU).getMass();
273 mtau2 = mtau * mtau;
274 //from PDG 2024 tau lifetime: need SM prediction
275 Gammatau = HCUT / 0.2903;
276
277 allcoeff_nu = mySM.getFlavour().ComputeCoeffsnunu(QCD::NEUTRINO_1);
278 C_R_nunu_e = ((*(allcoeff_nu[LO]))(1) + (*(allcoeff_nu[NLO]))(1) + (*(allcoeff_nu[NLO_QED11]))(1));
279 if (FixedWCbtos) {
280 allcoeff_noSM_nu = mySM.getFlavour().ComputeCoeffsnunu(QCD::NEUTRINO_1,true); //check the mass scale, scheme fixed to NDR
281 C_L_nunu_e = mySM.getOptionalParameter("CLnunu_SM") + ((*(allcoeff_noSM_nu[LO]))(0) + (*(allcoeff_noSM_nu[NLO]))(0) + (*(allcoeff_noSM_nu[NLO_QED11]))(0));
282 } else
283 C_L_nunu_e = ((*(allcoeff_nu[LO]))(0) + (*(allcoeff_nu[NLO]))(0) + (*(allcoeff_nu[NLO_QED11]))(0));
284
285 allcoeff_nu = mySM.getFlavour().ComputeCoeffsnunu(QCD::NEUTRINO_2);
286 C_R_nunu_mu = ((*(allcoeff_nu[LO]))(1) + (*(allcoeff_nu[NLO]))(1) + (*(allcoeff_nu[NLO_QED11]))(1));
287 if (FixedWCbtos) {
288 allcoeff_noSM_nu = mySM.getFlavour().ComputeCoeffsnunu(QCD::NEUTRINO_2,true); //check the mass scale, scheme fixed to NDR
289 C_L_nunu_mu = mySM.getOptionalParameter("CLnunu_SM") + ((*(allcoeff_noSM_nu[LO]))(0) + (*(allcoeff_noSM_nu[NLO]))(0) + (*(allcoeff_noSM_nu[NLO_QED11]))(0));
290 } else
291 C_L_nunu_mu = ((*(allcoeff_nu[LO]))(0) + (*(allcoeff_nu[NLO]))(0) + (*(allcoeff_nu[NLO_QED11]))(0));
292
293 allcoeff_nu = mySM.getFlavour().ComputeCoeffsnunu(QCD::NEUTRINO_3);
294 C_R_nunu_tau = ((*(allcoeff_nu[LO]))(1) + (*(allcoeff_nu[NLO]))(1) + (*(allcoeff_nu[NLO_QED11]))(1));
295 if (FixedWCbtos) {
296 allcoeff_noSM_nu = mySM.getFlavour().ComputeCoeffsnunu(QCD::NEUTRINO_3,true); //check the mass scale, scheme fixed to NDR
297 C_L_nunu_tau = mySM.getOptionalParameter("CLnunu_SM") + ((*(allcoeff_noSM_nu[LO]))(0) + (*(allcoeff_noSM_nu[NLO]))(0) + (*(allcoeff_noSM_nu[NLO_QED11]))(0));
298 } else
299 C_L_nunu_tau = ((*(allcoeff_nu[LO]))(0) + (*(allcoeff_nu[NLO]))(0) + (*(allcoeff_nu[NLO_QED11]))(0));
300 }
301 else{
302 allcoeff = mySM.getFlavour().ComputeCoeffBMll(mu_b, lep); //check the mass scale, scheme fixed to NDR
303 allcoeffprime = mySM.getFlavour().ComputeCoeffprimeBMll(mu_b, lep); //check the mass scale, scheme fixed to NDR
304
305 C_1 = (*(allcoeff[LO]))(0) + (*(allcoeff[NLO]))(0);
306 C_1L_bar = (*(allcoeff[LO]))(0) / 2.;
307 C_2 = ((*(allcoeff[LO]))(1) + (*(allcoeff[NLO]))(1));
308 C_2L_bar = (*(allcoeff[LO]))(1) - (*(allcoeff[LO]))(0) / 6.;
309 C_3 = ((*(allcoeff[LO]))(2) + (*(allcoeff[NLO]))(2));
310 C_4 = (*(allcoeff[LO]))(3) + (*(allcoeff[NLO]))(3);
311 C_5 = ((*(allcoeff[LO]))(4) + (*(allcoeff[NLO]))(4));
312 C_6 = ((*(allcoeff[LO]))(5) + (*(allcoeff[NLO]))(5));
313 C_8 = ((*(allcoeff[LO]))(7) + (*(allcoeff[NLO]))(7));
314 C_8L = (*(allcoeff[LO]))(7);
315 C_S = MW / Mb * ((*(allcoeff[LO]))(10) + (*(allcoeff[NLO]))(10));
316 C_P = MW / Mb * ((*(allcoeff[LO]))(11) + (*(allcoeff[NLO]))(11));
317 C_9p = (*(allcoeffprime[LO]))(8) + (*(allcoeffprime[NLO]))(8);
318 C_10p = (*(allcoeffprime[LO]))(9) + (*(allcoeffprime[NLO]))(9);
319 C_Sp = MW / Mb * ((*(allcoeffprime[LO]))(10) + (*(allcoeffprime[NLO]))(10));
320 C_Pp = MW / Mb * ((*(allcoeffprime[LO]))(11) + (*(allcoeffprime[NLO]))(11));
321
322 if (FixedWCbtos) {
323 allcoeff_noSM = mySM.getFlavour().ComputeCoeffBMll(mu_b, lep, true); //check the mass scale, scheme fixed to NDR
324 C_7 = mySM.getOptionalParameter("C7_SM") + ((*(allcoeff_noSM[LO]))(6) + (*(allcoeff_noSM[NLO]))(6));
325 C_9 = mySM.getOptionalParameter("C9_SM") + ((*(allcoeff_noSM[LO]))(8) + (*(allcoeff_noSM[NLO]))(8));
326 C_10 = mySM.getOptionalParameter("C10_SM") + ((*(allcoeff_noSM[LO]))(9) + (*(allcoeff_noSM[NLO]))(9));
327 } else {
328 C_7 = ((*(allcoeff[LO]))(6) + (*(allcoeff[NLO]))(6));
329 C_9 = ((*(allcoeff[LO]))(8) + (*(allcoeff[NLO]))(8));
330 C_10 = ((*(allcoeff[LO]))(9) + (*(allcoeff[NLO]))(9));
331 }
332 C_7p = MsoMb * ((*(allcoeffprime[LO]))(6) + (*(allcoeffprime[NLO]))(6));
333 C_7p -= MsoMb * (C_7 + 1. / 3. * C_3 + 4 / 9 * C_4 + 20. / 3. * C_5 + 80. / 9. * C_6);
334
335 allcoeffh = mySM.getFlavour().ComputeCoeffBMll(mu_h, lep); //check the mass scale, scheme fixed to NDR
336
337 C_1Lh_bar = (*(allcoeffh[LO]))(0) / 2.;
338 C_2Lh_bar = (*(allcoeffh[LO]))(1) - (*(allcoeff[LO]))(0) / 6.;
339 C_8Lh = (*(allcoeffh[LO]))(7);
340
341 if (!dispersion) {
342 C_7 += Delta_C7_U;
343 C_9 += Delta_C9_U;
344 }
345 }
346
347 checkCache();
348
349 H_0_pre = 8. / 27. + 4. / 9. * gslpp::complex::i() * M_PI;
350 H_0_WC = (C_3 + 4. / 3. * C_4 + 16. * C_5 + 64. / 3. * C_6);
351 H_c_WC = (4. / 3. * C_1 + C_2 + 6. * C_3 + 60. * C_5);
352 H_b_WC = (7. * C_3 + 4. / 3. * C_4 + 76. * C_5 + 64. / 3. * C_6);
353 fournineth = 4. / 9.;
354 half = 1. / 2.;
355 twothird = 2. / 3.;
356 ihalfMPI = gslpp::complex::i() * M_PI / 2.;
357 Mc2 = Mc*Mc;
358 Mb2 = Mb*Mb;
359 mu_b2 = mu_b*mu_b;
360 logMc = log(Mc2 / mu_b2);
361 logMb = log(Mb2 / mu_b2);
362 fourMc2 = 4. * Mc2;
363 fourMb2 = 4. * Mb2;
364 Mlep2 = Mlep*Mlep;
365 NN = ((4. * GF * MM * ale * lambda_t) / (sqrt(2.)*4. * M_PI)).abs2();
366 MM2 = MM*MM;
367 MM4 = MM2*MM2;
368 MP2 = MP*MP;
369 MP4 = MP2*MP2;
370 MM2mMP2 = MM2 - MP2;
371 twoMP2 = 2. * MP2;
372 twoMM = 2. * MM;
373 twoMM2 = 2. * MM2;
374 twoMM2_MMpMP = twoMM2 * (MM + MP);
375 twoMM_MbpMs = twoMM * (Mb + Ms);
376 S_L_pre = -(MM2mMP2 / twoMM_MbpMs) * (1 + MsoMb) / (1 - MsoMb);
377 fourMM2 = 4. * MM2;
378 twoMboMM = 2 * Mb / MM;
379 sixteenM_PI2 = 16. * M_PI*M_PI;
380 ninetysixM_PI3MM3 = 96. * M_PI * M_PI * M_PI * MM * MM*MM;
381 MboMW = Mb / MW;
382 MboMM = Mb / MM;
383 MsoMb = Ms / Mb;
384 twoMlepMb = 2. * Mlep*Mb;
385 threeGegen0 = mySM.getMesons(pseudoscalar).getGegenalpha(0)*3.;
386 threeGegen1otwo = mySM.getMesons(pseudoscalar).getGegenalpha(1)*3. / 2.;
387 M_PI2osix = M_PI * M_PI / 6.;
388 twoMc2 = 2. * Mc2;
389 sixMMoMb = 6. * MM / Mb;
390 CF = 4. / 3.;
391 deltaT_0 = alpha_s_mub * MboMM / 4. / M_PI;
392 deltaT_1par = mySM.Als(mu_h) * CF / 4. * M_PI / 3. * mySM.getMesons(meson).getDecayconst() *
393 mySM.getMesons(pseudoscalar).getDecayconst() / MM * MboMM;
394
395 F87_0 = -32. / 9. * log(mu_b / Mb) + 8. / 27. * M_PI * M_PI - 44. / 9. - 8. / 9. * gslpp::complex::i() * M_PI;
396 F87_1 = (4. / 3. * M_PI * M_PI - 40. / 3.);
397 F87_2 = (32. / 9. * M_PI * M_PI - 316. / 9.);
398 F87_3 = (200. / 27. * M_PI * M_PI - 658. / 9.);
399
400 F89_0 = 104. / 9. - 32. / 27. * M_PI * M_PI;
401 F89_1 = 1184. / 27. - 40. / 9. * M_PI * M_PI;
402 F89_2 = (-32. / 3. * M_PI * M_PI + 14212. / 135.);
403 F89_3 = (-560. / 27. * M_PI * M_PI + 193444. / 945.);
404
405 F29_0 = (256. / 243. - 32. / 81. * gslpp::complex::i() * M_PI - 128. / 9. * log(Mc / Mb)) * log(mu_b / Mb) + 512. / 81. * log(mu_b / Mb) * log(mu_b / Mb) + 5.4082 - 1.0934 * gslpp::complex::i();
406 F29_L1 = 32. / 81. * log(mu_b / Mb) + (0.48576 + 0.31119 * gslpp::complex::i());
407 F29_1 = (-32. / 405. + 64. / 45. / Mc2 * Mb2) * log(mu_b / Mb) + (1.9061 + 0.80843 * gslpp::complex::i());
408 F29_2 = (-8. / 945. + 16. / 105. / Mc2 * Mb2 / Mc2 * Mb2) * log(mu_b / Mb) + (-1.8286 + 2.8428 * gslpp::complex::i());
409 F29_3 = (-32. / 25515. + 64. / 2835. / Mc2 * Mb2 / Mc2 * Mb2 / Mc2 * Mb2) * log(mu_b / Mb) + (-12.113 + 8.1251 * gslpp::complex::i());
410 F29_L1_1 = (0.21951 - 0.14852 * gslpp::complex::i());
411 F29_L1_2 = (0.13015 - 0.22155 * gslpp::complex::i());
412 F29_L1_3 = (-0.079692 - 0.31214 * gslpp::complex::i());
413
414 F27_0 = 416. / 81. * log(mu_b / Mb) + 3.8367 + 0.3531 * gslpp::complex::i();
415 F27_1 = (1.3098 + 0.60185 * gslpp::complex::i());
416 F27_2 = (0.13507 + 0.89014 * gslpp::complex::i());
417 F27_3 = (-1.0271 + 0.77168 * gslpp::complex::i());
418 F27_L1_1 = (-0.031936 - 0.10981 * gslpp::complex::i());
419 F27_L1_2 = (-0.14169 - 0.035553 * gslpp::complex::i());
420 F27_L1_3 = (-0.13592 + 0.093 * gslpp::complex::i());
421
422 std::map<std::pair<double, double>, unsigned int >::iterator it;
423
424 if (I0_updated == 0) for (it = sigma0Cached.begin(); it != sigma0Cached.end(); ++it) it->second = 0;
425 if (I2_updated == 0) for (it = sigma2Cached.begin(); it != sigma2Cached.end(); ++it) it->second = 0;
426 if (I8_updated == 0) for (it = sigma8Cached.begin(); it != sigma8Cached.end(); ++it) it->second = 0;
427
428 if (I0_updated == 0) for (it = delta0Cached.begin(); it != delta0Cached.end(); ++it) it->second = 0;
429 if (I2_updated == 0) for (it = delta2Cached.begin(); it != delta2Cached.end(); ++it) it->second = 0;
430
431 if (Itree_updated) for (it = sigmaTreeCached.begin(); it != sigmaTreeCached.end(); ++it) it->second = 0;
432
433 std::map<double, unsigned int >::iterator iti;
434 if (deltaTparpupdated == 0) for (iti = deltaTparpCached.begin(); iti != deltaTparpCached.end(); ++iti) iti->second = 0;
435 if (deltaTparmupdated == 0) for (iti = deltaTparmCached.begin(); iti != deltaTparmCached.end(); ++iti) iti->second = 0;
436
437 if (deltaTparpupdated * deltaTparmupdated == 0) for (it = I1Cached.begin(); it != I1Cached.end(); ++it) it->second = 0;
438
439#if SPLINE
440 if (mySM.getFlavour().getUpdateFlag(meson, pseudoscalar, lep)) spline_QCDF_func();
441#else
442 fit_DeltaC9_mumu();
443#endif
444
445 mySM.getFlavour().setUpdateFlag(meson, pseudoscalar, lep, false);
446
447 //std::cout << "fplus_DM(4)= " << f_plus(4) << std::endl;
448 //std::cout << "f0_DM(4)= " << f_0(4) << std::endl;
449 //std::cout << "fT_DM(4)= " << f_T(4) << std::endl;
450
451 return;
452
453}
454
455void MPll::checkCache()
456{
457
458 if (MM == k2_cache(0) && MP == k2_cache(1)) {
459 k2_updated = 1;
460 } else {
461 k2_updated = 0;
462 k2_cache(0) = MM;
463 k2_cache(1) = MP;
464 }
465
466 if (MPll_Lattice_flag || MPll_GRvDV_flag || MPll_DM_flag) {
467 if (b_0_fplus == fplus_lat_cache(0) && b_1_fplus == fplus_lat_cache(1) && b_2_fplus == fplus_lat_cache(2)) {
468 fplus_lat_updated = 1;
469 } else {
470 fplus_lat_updated = 0;
471 fplus_lat_cache(0) = b_0_fplus;
472 fplus_lat_cache(1) = b_1_fplus;
473 fplus_lat_cache(2) = b_2_fplus;
474 }
475
476 if (b_0_fT == fT_lat_cache(0) && b_1_fT == fT_lat_cache(1) && b_2_fT == fT_lat_cache(2)) {
477 fT_lat_updated = 1;
478 } else {
479 fT_lat_updated = 0;
480 fT_lat_cache(0) = b_0_fT;
481 fT_lat_cache(1) = b_1_fT;
482 fT_lat_cache(2) = b_2_fT;
483 }
484
485 if (b_0_f0 == f0_lat_cache(0) && b_1_f0 == f0_lat_cache(1) && b_2_f0 == f0_lat_cache(2)) {
486 f0_lat_updated = 1;
487 } else {
488 f0_lat_updated = 0;
489 f0_lat_cache(0) = b_0_f0;
490 f0_lat_cache(1) = b_1_f0;
491 f0_lat_cache(2) = b_2_f0;
492 }
493 } else {
494 if (r_1_fplus == fplus_cache(0) && r_2_fplus == fplus_cache(1)) {
495 fplus_updated = 1;
496 } else {
497 fplus_updated = 0;
498 fplus_cache(0) = r_1_fplus;
499 fplus_cache(1) = r_2_fplus;
500 }
501
502 if (r_1_fT == fT_cache(0) && r_2_fT == fT_cache(1)) {
503 fT_updated = 1;
504 } else {
505 fT_updated = 0;
506 fT_cache(0) = r_1_fT;
507 fT_cache(1) = r_2_fT;
508 }
509
510 if (r_2_f0 == f0_cache) {
511 f0_updated = 1;
512 } else {
513 f0_updated = 0;
514 f0_cache = r_2_f0;
515 }
516 }
517
518 if (Mlep == beta_cache) {
519 beta_updated = 1;
520 } else {
521 beta_updated = 0;
522 beta_cache = Mlep;
523 }
524
525 lambda_updated = k2_updated;
526 F_updated = lambda_updated * beta_updated;
527
528 if (MPll_Lattice_flag || MPll_GRvDV_flag || MPll_DM_flag) {
529 VL_updated = k2_updated * fplus_lat_updated;
530 TL_updated = k2_updated * fT_lat_updated;
531 } else {
532 VL_updated = k2_updated * fplus_updated;
533 TL_updated = k2_updated * fT_updated;
534 }
535
536 if (Mb == SL_cache(0) && Ms == SL_cache(1)) {
537 if (MPll_Lattice_flag || MPll_GRvDV_flag || MPll_DM_flag)
538 SL_updated = k2_updated * f0_lat_updated;
539 else
540 SL_updated = k2_updated * f0_updated;
541 } else {
542 SL_updated = 0;
543 SL_cache(0) = Mb;
544 SL_cache(1) = Ms;
545 }
546
547 if (GF == N_cache(0) && ale == N_cache(1) && MM == N_cache(2) && lambda_t == Nc_cache) {
548 N_updated = 1;
549 } else {
550 N_updated = 0;
551 N_cache(0) = GF;
552 N_cache(1) = ale;
553 N_cache(2) = MM;
554 Nc_cache = lambda_t;
555 }
556
557 if (C_1 == C_1_cache) {
558 C_1_updated = 1;
559 } else {
560 C_1_updated = 0;
561 C_1_cache = C_1;
562 }
563
564 if (C_2 == C_2_cache) {
565 C_2_updated = 1;
566 } else {
567 C_2_updated = 0;
568 C_2_cache = C_2;
569 }
570
571 if (C_3 == C_3_cache) {
572 C_3_updated = 1;
573 } else {
574 C_3_updated = 0;
575 C_3_cache = C_3;
576 }
577
578 if (C_4 == C_4_cache) {
579 C_4_updated = 1;
580 } else {
581 C_4_updated = 0;
582 C_4_cache = C_4;
583 }
584
585 if (C_5 == C_5_cache) {
586 C_5_updated = 1;
587 } else {
588 C_5_updated = 0;
589 C_5_cache = C_5;
590 }
591
592 if (C_6 == C_6_cache) {
593 C_6_updated = 1;
594 } else {
595 C_6_updated = 0;
596 C_6_cache = C_6;
597 }
598
599 if (C_7 == C_7_cache) {
600 C_7_updated = 1;
601 } else {
602 C_7_updated = 0;
603 C_7_cache = C_7;
604 }
605
606 if (C_9 == C_9_cache) {
607 C_9_updated = 1;
608 } else {
609 C_9_updated = 0;
610 C_9_cache = C_9;
611 }
612
613 if (C_10 == C_10_cache) {
614 C_10_updated = 1;
615 } else {
616 C_10_updated = 0;
617 C_10_cache = C_10;
618 }
619
620 if (C_S == C_S_cache) {
621 C_S_updated = 1;
622 } else {
623 C_S_updated = 0;
624 C_S_cache = C_S;
625 }
626
627 if (C_P == C_P_cache) {
628 C_P_updated = 1;
629 } else {
630 C_P_updated = 0;
631 C_P_cache = C_P;
632 }
633
634 if (C_7p == C_7p_cache) {
635 C_7p_updated = 1;
636 } else {
637 C_7p_updated = 0;
638 C_7p_cache = C_7p;
639 }
640
641 if (C_9p == C_9p_cache) {
642 C_9p_updated = 1;
643 } else {
644 C_9p_updated = 0;
645 C_9p_cache = C_9p;
646 }
647
648 if (C_10p == C_10p_cache) {
649 C_10p_updated = 1;
650 } else {
651 C_10p_updated = 0;
652 C_10p_cache = C_10p;
653 }
654
655 if (C_Sp == C_Sp_cache) {
656 C_Sp_updated = 1;
657 } else {
658 C_Sp_updated = 0;
659 C_Sp_cache = C_Sp;
660 }
661
662 if (C_Pp == C_Pp_cache) {
663 C_Pp_updated = 1;
664 } else {
665 C_Pp_updated = 0;
666 C_Pp_cache = C_Pp;
667 }
668
669 if (C_2Lh_bar == C_2Lh_cache) {
670 C_2Lh_updated = 1;
671 } else {
672 C_2Lh_updated = 0;
673 C_2Lh_cache = C_2Lh_bar;
674 }
675
676 if (C_8Lh == C_8Lh_cache) {
677 C_8Lh_updated = 1;
678 } else {
679 C_8Lh_updated = 0;
680 C_8Lh_cache = C_8Lh;
681 }
682
683 if (C_L_nunu_e == C_L_nunu_e_cache) {
684 C_L_nunu_e_updated = 1;
685 } else {
686 C_L_nunu_e_updated = 0;
687 C_L_nunu_e_cache = C_L_nunu_e;
688 }
689
690 if (C_L_nunu_mu == C_L_nunu_mu_cache) {
691 C_L_nunu_mu_updated = 1;
692 } else {
693 C_L_nunu_mu_updated = 0;
694 C_L_nunu_mu_cache = C_L_nunu_mu;
695 }
696
697 if (C_L_nunu_tau == C_L_nunu_tau_cache) {
698 C_L_nunu_tau_updated = 1;
699 } else {
700 C_L_nunu_tau_updated = 0;
701 C_L_nunu_tau_cache = C_L_nunu_tau;
702 }
703
704 if (C_R_nunu_e == C_R_nunu_e_cache) {
705 C_R_nunu_e_updated = 1;
706 } else {
707 C_R_nunu_e_updated = 0;
708 C_R_nunu_e_cache = C_R_nunu_e;
709 }
710
711 if (C_R_nunu_mu == C_R_nunu_mu_cache) {
712 C_R_nunu_mu_updated = 1;
713 } else {
714 C_R_nunu_mu_updated = 0;
715 C_R_nunu_mu_cache = C_R_nunu_mu;
716 }
717
718 if (C_R_nunu_tau == C_R_nunu_tau_cache) {
719 C_R_nunu_tau_updated = 1;
720 } else {
721 C_R_nunu_tau_updated = 0;
722 C_R_nunu_tau_cache = C_R_nunu_tau;
723 }
724
725 if (Mb == Ycache(0) && Mc == Ycache(1)) {
726 Yupdated = C_1_updated * C_2_updated * C_3_updated * C_4_updated * C_5_updated * C_6_updated;
727 } else {
728 Yupdated = 0;
729 Ycache(0) = Mb;
730 Ycache(1) = Mc;
731 }
732
733 if (lep == QCD::NEUTRINO_1){
734 H_V0updated = N_updated * VL_updated * C_L_nunu_e_updated * C_L_nunu_mu_updated * C_L_nunu_tau_updated * C_R_nunu_e_updated * C_R_nunu_mu_updated * C_R_nunu_tau_updated;
735 H_A0updated = H_V0updated;
736 } else {
737
738 if (!dispersion) {
739 if (MM == H_V0cache(0) && Mb == H_V0cache(1) && h_0 == H_V0Ccache[0] && h_1 == H_V0Ccache[1] && h_2 == H_V0Ccache[2]) {
740 H_V0updated = N_updated * C_9_updated * Yupdated * VL_updated * C_9p_updated * C_7_updated * TL_updated * C_7p_updated;
741 } else {
742 H_V0updated = 0;
743 H_V0cache(0) = MM;
744 H_V0cache(1) = Mb;
745 H_V0Ccache[0] = h_0;
746 H_V0Ccache[1] = h_1;
747 H_V0Ccache[2] = h_2;
748 }
749 } else {
750 if (MM == H_V0cache(0) && Mb == H_V0cache(1) && r_1 == H_V0Ccache_dispersion[0] && r_2 == H_V0Ccache_dispersion[1] && Delta_C9 == H_V0Ccache_dispersion[2] && exp_Phase == H_V0Ccache_dispersion[3]) {
751 H_V0updated = N_updated * C_9_updated * Yupdated * VL_updated * C_9p_updated * C_7_updated * TL_updated * C_7p_updated;
752 } else {
753 H_V0updated = 0;
754 H_V0cache(0) = MM;
755 H_V0cache(1) = Mb;
756 H_V0Ccache_dispersion[0] = r_1;
757 H_V0Ccache_dispersion[1] = r_2;
758 H_V0Ccache_dispersion[2] = Delta_C9;
759 H_V0Ccache_dispersion[3] = exp_Phase;
760 }
761 }
762
763 H_A0updated = N_updated * C_10_updated * VL_updated * C_10p_updated;
764 }
765
766 if (Mb == H_Scache(0) && MW == H_Scache(1)) {
767 H_Supdated = N_updated * C_S_updated * SL_updated * C_Sp_updated;
768 } else {
769 H_Supdated = 0;
770 H_Scache(0) = Mb;
771 H_Scache(1) = MW;
772 }
773
774 if (Mb == H_P_cache(0) && MW == H_P_cache(1) && Mlep == H_P_cache(2) && Ms == H_P_cache(3)) {
775 H_P_updated = N_updated * C_P_updated * SL_updated * C_Pp_updated * SL_updated * C_10_updated * C_10p_updated;
776 } else {
777 H_P_updated = 0;
778 H_P_cache(0) = Mb;
779 H_P_cache(1) = MW;
780 H_P_cache(2) = Mlep;
781 H_P_cache(3) = Ms;
782 }
783
784 if (MM == T_cache(0) && Mb == T_cache(1) && Mc == T_cache(2) &&
785 mySM.getMesons(pseudoscalar).getGegenalpha(0) == T_cache(3) && mySM.getMesons(pseudoscalar).getGegenalpha(1) == T_cache(4)) {
786 T_updated = 1;
787 } else {
788 T_updated = 0;
789 T_cache(0) = MM;
790 T_cache(1) = Mb;
791 T_cache(2) = Mc;
792 T_cache(3) = mySM.getMesons(pseudoscalar).getGegenalpha(0);
793 T_cache(4) = mySM.getMesons(pseudoscalar).getGegenalpha(1);
794 }
795
796 deltaTparpupdated = C_2Lh_updated * T_updated;
797 deltaTparmupdated = C_2Lh_updated * C_8Lh_updated * T_updated;
798
799 I0_updated = F_updated * H_V0updated * H_A0updated * H_P_updated * beta_updated * H_Supdated * deltaTparmupdated;
800 I2_updated = F_updated * beta_updated * H_V0updated * H_A0updated * deltaTparmupdated;
801 I8_updated = F_updated * beta_updated * H_Supdated * H_V0updated * deltaTparmupdated;
802
803 if (MM2 == Itree_cache(0) && mtau2 == Itree_cache(1) && MP2 == Itree_cache(2)) {
804 Itree_updated = 1;
805 } else {
806 Itree_updated = 0;
807 Itree_cache(0) = MM2;
808 Itree_cache(1) = mtau2;
809 Itree_cache(2) = MP2;
810 }
811
812}
813
814/*******************************************************************************
815 * Transverse Form Factors *
816 * ****************************************************************************/
817double MPll::LCSR_fit1(double q2, double r_1, double r_2, double m_fit2)
818{
819 return r_1 / (1. - q2 / m_fit2) + r_2 / pow((1. - q2 / m_fit2), 2.);
820
821}
822
823double MPll::LCSR_fit2(double q2, double r_2, double m_fit2)
824{
825 return r_2 / (1. - q2 / m_fit2);
826}
827
828double MPll::LCSR_fit3(double q2, double b_0, double b_1, double b_2, double m_fit2)
829{
830 return 1. / (1. - q2 / m_fit2) * (b_0 + b_1 * (zeta(q2) - zeta(0)) + b_2 * (zeta(q2) - zeta(0)) * (zeta(q2) - zeta(0)));
831}
832
833double MPll::zeta(double q2)
834{
835 double tp, t0;
836
837 tp = (MM + MP)*(MM + MP);
838 t0 = (MM + MP)*(sqrt(MM) - sqrt(MP))*(sqrt(MM) - sqrt(MP));
839
840 return (sqrt(tp - q2) - sqrt(tp - t0)) / (sqrt(tp - q2) + sqrt(tp - t0));
841}
842
843double MPll::zeta_DM(double q2)
844{
845 double tp, tm;
846
847 tp = (MM + MP)*(MM + MP);
848 tm = (MM - MP)*(MM - MP);
849
850 return (sqrt(tp - q2) - sqrt(tp - tm)) / (sqrt(tp - q2) + sqrt(tp - tm));
851}
852
853double MPll::phiplus_DM(double q2, double m_fit2)
854{
855 double z = zeta_DM(q2);
856 double z_M = zeta_DM(m_fit2);
857 double rP = MP/MM;
858 double Chi1minus = 0.000623174575;
859
860 return 16.*rP*rP/MM*sqrt(4./3./Chi1minus/M_PI) * (1. + z)*(1. + z)*sqrt(1. - z)/pow((1. + rP)*(1. - z)+2.*sqrt(rP)*(1. + z),5) * (z - z_M)/(1. - z_M*z);
861}
862
863double MPll::phi0_DM(double q2)
864{
865 double z = zeta_DM(q2);
866 double rP = MP/MM;
867 double Chi0plus = 0.0142;
868
869 return 2.*rP*(1 - rP*rP)*sqrt(4./Chi0plus/M_PI)* (1. - z*z)*sqrt(1. - z)/pow((1. + rP)*(1. - z)+2.*sqrt(rP)*(1. + z),4);
870}
871
872double MPll::phiT_DM(double q2, double m_fit2)
873{
874 double z = zeta_DM(q2);
875 double z_M = zeta_DM(m_fit2);
876 double rP = MP/MM;
877 double ChiTT = 0.0454644444;
878
879 return 16*rP*rP/MM/(1. + rP)*sqrt(4./3./ChiTT/M_PI) * (1. + z)*(1. + z)/sqrt(1. - z)/pow((1. + rP)*(1. - z)+2.*sqrt(rP)*(1. + z),4) * (z - z_M)/(1. - z_M*z);
880}
881
882double MPll::fplus_DM(double q2, double b_0, double b_1, double b_2, double m_fit2)
883{
884 double z = zeta_DM(q2);
885
886 return (b_0 + b_1*z + b_2*z*z) / phiplus_DM(q2, m_fit2);
887}
888
889double MPll::f0_DM(double q2, double b_0, double b_1, double b_2)
890{
891 double z = zeta_DM(q2);
892
893 return (b_0 + b_1*z + b_2*z*z) / phi0_DM(q2);
894}
895
896double MPll::fT_DM(double q2, double b_0, double b_1, double b_2, double m_fit2)
897{
898 double z = zeta_DM(q2);
899
900 return (b_0 + b_1*z + b_2*z*z) / phiT_DM(q2, m_fit2);
901}
902
903
904double MPll::LATTICE_fit1(double q2, double b_0, double b_1, double b_2, double m_fit2)
905{
906 double z2 = zeta(q2) * zeta(q2);
907 double z3 = zeta(q2) * z2;
908
909 return 1. / (1. - q2 / m_fit2) * (b_0 + b_1 * (zeta(q2) - 1. / 3. * z3) + b_2 * (z2 + 2. / 3. * z3));
910
911}
912
913double MPll::LATTICE_fit2(double q2, double b_0, double b_1, double b_2, double m_fit2)
914{
915 return 1. / (1. - q2 / m_fit2) * (b_0 + b_1 * zeta(q2) + b_2 * zeta(q2) * zeta(q2));
916}
917
918double MPll::f_plus(double q2)
919{
920 if (MPll_Lattice_flag)
921 return LATTICE_fit1(q2, b_0_fplus, b_1_fplus, b_2_fplus, m_fit2_fplus_lat);
922 else if (MPll_GRvDV_flag)
923 return LCSR_fit3(q2, b_0_fplus, b_1_fplus, b_2_fplus, m_fit2_fplus_lat);
924 else if (MPll_DM_flag)
925 return fplus_DM(q2, b_0_fplus, b_1_fplus, b_2_fplus, m_fit2_fplus_lat);
926 else
927 return LCSR_fit1(q2, r_1_fplus, r_2_fplus, m_fit2_fplus);
928}
929
930double MPll::f_T(double q2)
931{
932 if (MPll_Lattice_flag)
933 return LATTICE_fit1(q2, b_0_fT, b_1_fT, b_2_fT, m_fit2_fT_lat);
934 else if (MPll_GRvDV_flag)
935 return LCSR_fit3(q2, b_0_fT, b_1_fT, b_2_fT, m_fit2_fT_lat);
936 else if (MPll_DM_flag)
937 return fT_DM(q2, b_0_fT, b_1_fT, b_2_fT, m_fit2_fT_lat);
938 else
939 return LCSR_fit1(q2, r_1_fT, r_2_fT, m_fit2_fT);
940}
941
942double MPll::f_0(double q2)
943{
944 if (MPll_Lattice_flag)
945 return LATTICE_fit2(q2, b_0_f0, b_1_f0, b_2_f0, m_fit2_f0_lat);
946 else if (MPll_GRvDV_flag)
947 return LCSR_fit3(q2, b_0_f0, b_1_f0, b_2_f0, m_fit2_f0_lat);
948 else if (MPll_DM_flag)
949 return f0_DM(q2, b_0_f0, b_1_f0, b_2_f0);
950 else
951 return LCSR_fit2(q2, r_2_f0, m_fit2_f0);
952}
953
954gslpp::complex MPll::V_L(double q2)
955{
956 return /*gslpp::complex::i() */ sqrt(lambda(q2)) / (twoMM * sqrt(q2)) * f_plus(q2);
957}
958
959gslpp::complex MPll::T_L(double q2)
960{
961 return /*gslpp::complex::i() */ sqrt(lambda(q2) * q2) / (twoMM2_MMpMP) * f_T(q2);
962}
963
964double MPll::S_L(double q2)
965{
966 return S_L_pre * f_0(q2);
967}
968
969/*******************************************************************************
970 * QCDF *
971 * ****************************************************************************/
972
973gslpp::complex MPll::I1(double u, double q2)
974{
975 std::pair<double, double > uq2 = std::make_pair(u, q2);
976
977 if (I1Cached[uq2] == 0) {
978 ubar = 1. - u;
979 xp = .5 + sqrt(0.25 - (Mc2 - gslpp::complex::i()*1.e-10) / (ubar * MM2 + u * q2));
980 xm = .5 - sqrt(0.25 - (Mc2 - gslpp::complex::i()*1.e-10) / (ubar * MM2 + u * q2));
981 yp = .5 + sqrt(0.25 - (Mc2 - gslpp::complex::i()*1.e-10) / q2);
982 ym = .5 - sqrt(0.25 - (Mc2 - gslpp::complex::i()*1.e-10) / q2);
983 L1xp = log(1. - 1. / xp) * log(1. - xp) - M_PI2osix + dilog(xp / (xp - 1.));
984 L1xm = log(1. - 1. / xm) * log(1. - xm) - M_PI2osix + dilog(xm / (xm - 1.));
985 L1yp = log(1. - 1. / yp) * log(1. - yp) - M_PI2osix + dilog(yp / (yp - 1.));
986 L1ym = log(1. - 1. / ym) * log(1. - ym) - M_PI2osix + dilog(ym / (ym - 1.));
987
988 cacheI1[uq2] = 1. + twoMc2 / ubar / (MM2 - q2)*(L1xp + L1xm - L1yp - L1ym);
989 I1Cached[uq2] = 1;
990 }
991
992 return cacheI1[uq2];
993}
994
995gslpp::complex MPll::Tparplus(double u, double q2)
996{
997 Ee = (MM2 - q2) / twoMM;
998 ubar = 1. - u;
999 arg1 = (fourMc2 - gslpp::complex::i()*1.e-10) / (ubar * MM2 + u * q2) - 1.;
1000 B01 = -2. * sqrt(arg1) * arctan(1. / sqrt(arg1));
1001 arg1 = (fourMc2 - gslpp::complex::i()*1.e-10) / q2 - 1.;
1002 B00 = -2. * sqrt(arg1) * arctan(1. / sqrt(arg1));
1003
1004 gslpp::complex tpar = twoMM / Ee / ubar * I1(u, q2) + (ubar * MM2 + u * q2) / Ee / Ee / ubar / ubar * (B01 - B00);
1005 return -MM / Mb * mySM.getQuarks(QCD::CHARM).getCharge() * tpar*C_2Lh_bar;
1006}
1007
1008gslpp::complex MPll::Tparminus(double u, double q2)
1009{
1010 double ubar = 1. - u;
1011 return -spectator_charge * (8. * C_8Lh / (ubar + u * q2 / MM2)
1012 + sixMMoMb * H_c(ubar * MM2 + u * q2, mu_h * mu_h) * C_2Lh_bar);
1013}
1014
1015double MPll::Integrand_ReTparplus(double up)
1016{
1017 double u = up;
1018 return ((Tparplus(u, tmpq2)*6. * u * (1. - u)*
1019 (1 + threeGegen0 * (2. * u - 1)
1020 + threeGegen1otwo * ((10. * u - 5.)*(2. * u - 1.) - 1.))) / mySM.getMesons(meson).getLambdaM()).real();
1021}
1022
1023double MPll::Integrand_ImTparplus(double up)
1024{
1025 double u = up;
1026 return ((Tparplus(u, tmpq2)*6. * u * (1. - u)*
1027 (1 + threeGegen0 * (2. * u - 1)
1028 + threeGegen1otwo * ((10. * u - 5.)*(2. * u - 1.) - 1.))) / mySM.getMesons(meson).getLambdaM()).imag();
1029}
1030
1031double MPll::Integrand_ReTparminus(double up)
1032{
1033 double Lambdaplus = mySM.getMesons(meson).getLambdaM();
1034 gslpp::complex Lambdamin = exp(-tmpq2 / MM / Lambdaplus) / Lambdaplus * (-gsl_sf_expint_Ei(tmpq2 / MM / Lambdaplus) + gslpp::complex::i() * M_PI);
1035
1036 double u = up;
1037 return ((Tparminus(u, tmpq2)*6. * u * (1. - u)*
1038 (1 + threeGegen0 * (2. * u - 1)
1039 + threeGegen1otwo * ((10. * u - 5.)*(2. * u - 1.) - 1.))) / Lambdamin).real();
1040}
1041
1042double MPll::Integrand_ImTparminus(double up)
1043{
1044 double Lambdaplus = mySM.getMesons(meson).getLambdaM();
1045 gslpp::complex Lambdamin = exp(-tmpq2 / MM / Lambdaplus) / Lambdaplus * (-gsl_sf_expint_Ei(tmpq2 / MM / Lambdaplus) + gslpp::complex::i() * M_PI);
1046
1047 double u = up;
1048 return ((Tparminus(u, tmpq2)*6. * u * (1. - u)*
1049 (1 + threeGegen0 * (2. * u - 1)
1050 + threeGegen1otwo * ((10. * u - 5.)*(2. * u - 1.) - 1.))) / Lambdamin).imag();
1051}
1052
1053double MPll::Integrand_ImTpar_pm(double up)
1054{
1055 return Integrand_ImTparplus(up) + Integrand_ImTparminus(up);
1056}
1057
1058double MPll::Integrand_ReTpar_pm(double up)
1059{
1060 return Integrand_ReTparplus(up) + Integrand_ReTparminus(up);
1061}
1062
1063gslpp::complex MPll::F19(double q2)
1064{
1065 double s = q2 / Mb2;
1066 double s2 = s*s;
1067 double Ls = log(s);
1068 double Lc = log(Mc / Mb);
1069 double Lm = log(mu_b / Mb);
1070 gslpp::complex i = gslpp::complex::i();
1071 return (-1424. / 729. + 16. / 243. * i * M_PI + 64. / 27. * Lc)*Lm - 16. / 243. * Lm * Ls + (16. / 1215. - 32. / 135. / Mc2 * Mb2) * Lm * s
1072 + (4. / 2835. - 8. / 315. / Mc2 * Mb2 / Mc2 * Mb2) * Lm * s2 + (16. / 76545. - 32. / 8505 / Mc2 * Mb2 / Mc2 * Mb2 / Mc2 * Mb2) * Lm * s * s2 - 256. / 243. * Lm * Lm
1073 + (-11.65 + 0.18223 * i + (-24.709 - 0.13474 * i) * s + (-43.588 - 0.4738 * i) * s2 + (-86.22 - 1.3542 * i) * s * s2
1074 + (-0.080959 - 0.051864 * i + (-0.036585 + 0.024753 * i) * s + (-0.021692 + 0.036925 * i) * s2 + (0.013282 + 0.052023 * i) * s * s2) * Ls);
1075}
1076
1077gslpp::complex MPll::F27(double q2)
1078{
1079 double s = q2 / Mb2;
1080 double s2 = s*s;
1081 double Ls = log(s);
1082 gslpp::complex i = gslpp::complex::i();
1083 return F27_0 + F27_1 * s + F27_2 * s2 + F27_3 * s * s2 + F27_L1_1 * Ls * s + F27_L1_2 * Ls * s2 + F27_L1_3 * Ls * s * s2;
1084}
1085
1086gslpp::complex MPll::F29(double q2)
1087{
1088 double s = q2 / Mb2;
1089 double s2 = s*s;
1090 double Ls = log(s);
1091 gslpp::complex i = gslpp::complex::i();
1092 return F29_0 + F29_L1 * Ls + F29_1 * s + F29_2 * s2 + F29_3 * s * s2 + F29_L1_1 * Ls * s + F29_L1_2 * Ls * s2 + F29_L1_3 * Ls * s2 *s;
1093}
1094
1095gslpp::complex MPll::F87(double q2)
1096{
1097 double s = q2 / Mb2;
1098 double s2 = s*s;
1099 return F87_0 + F87_1 * s + F87_2 * s2 + F87_3 * s * s2 - 0.888889 * log(s)*(1. + s + s2 + s * s2);
1100}
1101
1102double MPll::F89(double q2)
1103{
1104 double s = q2 / Mb2;
1105 double s2 = s*s;
1106 return F89_0 + F89_1 * s + F89_2 * s2 + F89_3 * s * s2 + 1.77778 * log(s)*(1. + s + s2 + s * s2);
1107}
1108
1109gslpp::complex MPll::Cpar(double q2)
1110{
1111 return -(C_2L_bar * F27(q2) + C_8L * F87(q2) + MM / 2. / Mb *
1112 (C_2L_bar * F29(q2) + 2. * C_1L_bar * (F19(q2) + F29(q2) / 6.) + C_8L * F89(q2)));
1113}
1114
1115gslpp::complex MPll::deltaTpar(double q2)
1116{
1117 tmpq2 = q2;
1118
1119 //old_handler = gsl_set_error_handler_off();
1120
1121 if (deltaTparpCached[q2] == 0) {
1122
1123 DTPPR = convertToGslFunction(bind(&MPll::Integrand_ReTpar_pm, &(*this), _1));
1124 if (gsl_integration_cquad(&DTPPR, 0., 1., 1.e-2, 1.e-1, w_DTPPR, &avaDTPPR, &errDTPPR, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
1125 double ReTppint = avaDTPPR;
1126
1127 DTPPR = convertToGslFunction(bind(&MPll::Integrand_ImTpar_pm, &(*this), _1));
1128 if (gsl_integration_cquad(&DTPPR, 0., 1., 1.e-2, 1.e-1, w_DTPPR, &avaDTPPR, &errDTPPR, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
1129 double ImTppint = avaDTPPR;
1130
1131 cacheDeltaTparp[q2] = (ReTppint + gslpp::complex::i() * ImTppint);
1132 deltaTparpCached[q2] = 1;
1133 }
1134
1135 //gsl_set_error_handler(old_handler);
1136
1137 return deltaT_0 * Cpar(q2) + deltaT_1par * cacheDeltaTparp[q2] / f_plus(q2);
1138}
1139
1140double MPll::reDC9fit(double* x, double* p)
1141{
1142 return p[0] / x[0] + p[1] + p[2] * x[0] + p[3] * x[0] * x[0] + p[4] * x[0] * x[0] * x[0] + p[5] * x[0] * x[0] * x[0] * x[0] + p[6] * x[0] * x[0] * x[0] * x[0] * x[0];
1143}
1144
1145double MPll::imDC9fit(double* x, double* p)
1146{
1147 return p[0] / x[0] + p[1] + p[2] * x[0] + p[3] * x[0] * x[0] + p[4] * x[0] * x[0] * x[0] + p[5] * x[0] * x[0] * x[0] * x[0] + p[6] * x[0] * x[0] * x[0] * x[0] * x[0] + p[7] * x[0] * x[0] * x[0] * x[0] * x[0] * x[0];
1148
1149 //double thr = 4.*Mc2;
1150}
1151
1152void MPll::fit_DeltaC9_mumu()
1153{
1154 int dim = 0;
1155 for (double i = 0.1; i < MPllSWITCH; i += 0.4) {
1156 double q2tmp = i;
1157 myq2.push_back(q2tmp);
1158 ReDeltaC9.push_back((deltaTpar(q2tmp)).real());
1159 ImDeltaC9.push_back((deltaTpar(q2tmp)).imag());
1160 dim++;
1161 }
1162 for (double i = MPllSWITCH; i < 8.2; i += 0.4) {
1163 double q2tmp = i;
1164 myq2.push_back(q2tmp);
1165 ReDeltaC9.push_back(q2tmp * (deltaTpar(q2tmp)).real());
1166 ImDeltaC9.push_back(q2tmp * (deltaTpar(q2tmp)).imag());
1167 dim++;
1168 }
1169
1170 gr1 = TGraph(dim, myq2.data(), ReDeltaC9.data());
1171 gr2 = TGraph(dim, myq2.data(), ImDeltaC9.data());
1172
1173 reffit = TF1("reffit", this, &MPll::reDC9fit, 0, 8.1, 7);
1174 imffit = TF1("imffit", this, &MPll::imDC9fit, 0, 8.1, 8);
1175
1176 refres = gr1.Fit(&reffit, "SQN0+rob=0.99");
1177 imfres = gr2.Fit(&imffit, "SQN0+rob=0.99");
1178
1179 ReDeltaC9.clear();
1180 ImDeltaC9.clear();
1181 myq2.clear();
1182}
1183
1184gslpp::complex MPll::fDeltaC9(double q2)
1185{
1186 if (q2 < MPllSWITCH) return (reDC9fit(&q2, const_cast<double *> (refres->GetParams()))
1187 + gslpp::complex::i() * imDC9fit(&q2, const_cast<double *> (imfres->GetParams())));
1188 else return (reDC9fit(&q2, const_cast<double *> (refres->GetParams()))
1189 + gslpp::complex::i() * imDC9fit(&q2, const_cast<double *> (imfres->GetParams()))) / q2;
1190
1191}
1192
1193gslpp::complex MPll::DeltaC9(double q2)
1194{
1195 return deltaTpar(q2);
1196}
1197
1198gslpp::complex MPll::deltaC7_QCDF(double q2, bool spline)
1199{
1200#if SPLINE
1201 if (spline) return gsl_spline_eval(spline_Re_deltaC7_QCDF, q2, acc_Re_deltaC7_QCDF);
1202#endif
1203
1204 double muh = mu_b / mb_pole;
1205 double z = mc_pole * mc_pole / mb_pole / mb_pole;
1206 double sh = q2 / mb_pole / mb_pole;
1207 double sh2 = sh*sh;
1208
1209 //gslpp::complex A_Sdl = A_Seidel(q2, mb_pole*mb_pole); /* hep-ph/0403185v2.*/
1210 //gslpp::complex Fu_17 = -A_Sdl; /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1211 //gslpp::complex Fu_27 = 6. * A_Sdl; /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1212 gslpp::complex F_17 = myF_1->F_17re(muh, z, sh, 20) + gslpp::complex::i() * myF_1->F_17im(muh, z, sh, 20); /*q^2 = 0 gives nan. Independent of how small q^2 is. arXiv:0810.4077*/
1213 gslpp::complex F_27 = myF_2->F_27re(muh, z, sh, 20) + gslpp::complex::i() * myF_2->F_27im(muh, z, sh, 20); /*q^2 = 0 gives nan. Independent of how small q^2 is. arXiv:0810.4077*/
1214 gslpp::complex F_87 = F87_0 + F87_1 * sh + F87_2 * sh2 + F87_3 * sh * sh2 - 8. / 9. * log(sh) * (sh + sh2 + sh * sh2);
1215
1216 gslpp::complex delta = C_1 * F_17 + C_2 * F_27;
1217 gslpp::complex delta_t = C_8 * F_87 + delta;
1218 //gslpp::complex delta_u = delta + C_1 * Fu_17 + C_2 * Fu_27;
1219
1220 return -alpha_s_mub / (4. * M_PI) * (delta_t /*- lambda_u / lambda_t * delta_u */);
1221}
1222
1223gslpp::complex MPll::deltaC9_QCDF(double q2, bool spline)
1224{
1225#if SPLINE
1226 if (spline) return gsl_spline_eval(spline_Re_deltaC9_QCDF, q2, acc_Re_deltaC9_QCDF);
1227#endif
1228 double muh = mu_b / mb_pole;
1229 double z = mc_pole * mc_pole / mb_pole / mb_pole;
1230 double sh = q2 / mb_pole / mb_pole;
1231 double sh2 = sh*sh;
1232
1233 //gslpp::complex B_Sdl = B_Seidel(q2, mb_pole*mb_pole); /* hep-ph/0403185v2.*/
1234 //gslpp::complex C_Sdl = C_Seidel(q2); /* hep-ph/0403185v2.*/
1235 //gslpp::complex Fu_19 = -(B_Sdl + 4. * C_Sdl); /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1236 //gslpp::complex Fu_29 = -(-6. * B_Sdl + 3. * C_Sdl); /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1237 gslpp::complex F_19 = myF_1->F_19re(muh, z, sh, 20) + gslpp::complex::i() * myF_1->F_19im(muh, z, sh, 20); /*q^2 = 0 gives nan. Independent of how small q^2 is. arXiv:0810.4077*/
1238 gslpp::complex F_29 = myF_2->F_29re(muh, z, sh, 20) + gslpp::complex::i() * myF_2->F_29im(muh, z, sh, 20); /*q^2 = 0 gives nan. Independent of how small q^2 is. arXiv:0810.4077*/
1239 gslpp::complex F_89 = (F89_0 + F89_1 * sh + F89_2 * sh2 + F89_3 * sh * sh2 + 16. / 9. * log(sh) * (1. + sh + sh2 + sh * sh2));
1240
1241 gslpp::complex delta = C_1 * F_19 + C_2 * F_29;
1242 gslpp::complex delta_t = C_8 * F_89 + delta;
1243 //gslpp::complex delta_u = delta + C_1 * Fu_19 + C_2 * Fu_29;
1244
1245 return -alpha_s_mub / (4. * M_PI) * (delta_t /*- lambda_u / lambda_t * delta_u*/);
1246}
1247
1248void MPll::spline_QCDF_func()
1249{
1250 int dim_DC = GSL_INTERP_DIM_DC;
1251 double min = 0.001;
1252 double interval_DC = (9.9 - min) / ((double) dim_DC);
1253 double q2_spline_DC[dim_DC];
1254 double fq2_Re_deltaC7_QCDF[dim_DC], fq2_Im_deltaC7_QCDF[dim_DC], fq2_Re_deltaC9_QCDF[dim_DC], fq2_Im_deltaC9_QCDF[dim_DC];
1255
1256 for (int i = 0; i < dim_DC; i++) {
1257 q2_spline_DC[i] = min + (double) i*interval_DC;
1258 fq2_Re_deltaC7_QCDF[i] = deltaC7_QCDF(q2_spline_DC[i], false).real();
1259 fq2_Im_deltaC7_QCDF[i] = deltaC7_QCDF(q2_spline_DC[i], false).imag();
1260 fq2_Re_deltaC9_QCDF[i] = deltaC9_QCDF(q2_spline_DC[i], false).real();
1261 fq2_Im_deltaC9_QCDF[i] = deltaC9_QCDF(q2_spline_DC[i], false).imag();
1262
1263 }
1264
1265 gsl_spline_init(spline_Re_deltaC7_QCDF, q2_spline_DC, fq2_Re_deltaC7_QCDF, dim_DC);
1266 gsl_spline_init(spline_Im_deltaC7_QCDF, q2_spline_DC, fq2_Im_deltaC7_QCDF, dim_DC);
1267 gsl_spline_init(spline_Re_deltaC9_QCDF, q2_spline_DC, fq2_Re_deltaC9_QCDF, dim_DC);
1268 gsl_spline_init(spline_Im_deltaC9_QCDF, q2_spline_DC, fq2_Im_deltaC9_QCDF, dim_DC);
1269
1270
1271}
1272
1273/*******************************************************************************
1274 * Helicity amplitudes *
1275 * ****************************************************************************/
1276gslpp::complex MPll::H_c(double q2, double mu2)
1277{
1278 double x = fourMc2 / q2;
1279 gslpp::complex par;
1280
1281 if (x > 1.) par = sqrt(x - 1.) * atan(1. / sqrt(x - 1.));
1282 else par = sqrt(1. - x) * (log((1. + sqrt(1. - x)) / sqrt(x)) - ihalfMPI);
1283 return -fournineth * (log(Mc2 / mu2) - twothird - x) - fournineth * (2. + x) * par;
1284}
1285
1286gslpp::complex MPll::H_b(double q2, double mu2)
1287{
1288 double x = fourMb2 / q2;
1289 gslpp::complex par;
1290
1291 if (x > 1.) par = sqrt(x - 1.) * atan(1. / sqrt(x - 1.));
1292 else par = sqrt(1. - x) * (log((1. + sqrt(1. - x)) / sqrt(x)) - ihalfMPI);
1293
1294 return -fournineth * (log(Mb2 / mu2) - twothird - x) - fournineth * (2. + x) * par;
1295}
1296
1297gslpp::complex MPll::H_0(double q2)
1298{
1299 return (H_0_pre - fournineth * log(q2 / mu_b2));
1300}
1301
1302gslpp::complex MPll::Y(double q2)
1303{
1304 return -half * H_0(q2) * H_0_WC + H_c(q2, mu_b2) * H_c_WC - half * H_b(q2, mu_b2) * H_b_WC;
1305}
1306
1307gslpp::complex MPll::funct_g(double q2)
1308{
1309 if (q2 < 4. * Mc * Mc)
1310 return -8. / 9. * log(Mc / Mb) + 8. / 27. + 16. / 9. * Mc * Mc / q2 - 4. / 9. * (2. + 4. * Mc * Mc / q2) * (sqrt(4. * Mc * Mc / q2 - 1.) * atan(1. / sqrt(4. * Mc * Mc / q2 - 1.)));
1311 else
1312 return -8. / 9. * log(Mc / Mb) + 8. / 27. + 16. / 9. * Mc * Mc / q2 - 4. / 9. * (2. + 4. * Mc * Mc / q2) * (sqrt(1. - 4. * Mc * Mc / q2) * (log(1. + sqrt(1. - 4. * Mc * Mc / q2) / sqrt(4. * Mc * Mc / q2)) - gslpp::complex::i() * M_PI_2));
1313}
1314
1315gslpp::complex MPll::DeltaC9_KD(double q2)
1316{
1317 return ((Delta_C9 + r_1 * q2 / mJ2) / (1. - r_2 * q2 / mJ2) - (3. * (-0.267) + 1.117) * funct_g(q2))*exp_Phase;
1318 /* C_1 = -0.267 and C_2 = 1.117 in KMPW */
1319}
1320
1321gslpp::complex MPll::h_lambda(double q2)
1322{
1323 if (!dispersion) {
1324// if (q2 <= 1.) return 1.3e-4/MM2 * (1. + gslpp::complex::i()) / sqrt(2.);
1325// else {
1326// h_2 = (-sixteenM_PI2*4.9e-7/MM2 * (1. + gslpp::complex::i()) / sqrt(2.) - h_1 / MM2 * V_L(1.) - h_2 * sixteenM_PI2) / twoMboMM / T_L(1.);
1327// h_2 = 1.3e-4/MM2 * (1. + gslpp::complex::i()) / sqrt(2.) - h_1 * V_L(1.)/sixteenM_PI2/MM2;
1328 //else return 4.9e-7/MM2 + h_1 * (q2 * V_L(q2) - T_L(q2). * V_L(1)/T_L(1.)) / MM2 / sixteenM_PI2;
1329
1330 // return (twoMboMM * h_0 * T_L(q2) + h_1 * q2 / MM2 * V_L(q2)) / sixteenM_PI2 + h_2 * q2 * sqrt(q2);
1331 return (h_1 + h_2 * q2 ) * sqrt(q2) * sqrt(lambda(q2)) / MM3 / 2.;
1332// }
1333 } else {
1334 return -q2 / (MM2 * sixteenM_PI2) * V_L(q2) * DeltaC9_KD(q2);
1335 }
1336}
1337
1338gslpp::complex MPll::H_V(double q2)
1339{
1340 return -((C_9 + deltaC9_QCDF(q2, SPLINE) + Y(q2) /*+ fDeltaC9(q2)*/ - etaP * pow(-1, angmomP) * C_9p) * V_L(q2)
1341 + MM2 / q2 * (twoMboMM * (C_7 + deltaC7_QCDF(q2, SPLINE) - etaP * pow(-1, angmomP) * C_7p) * T_L(q2)
1342 - sixteenM_PI2 * h_lambda(q2)));
1343}
1344
1345gslpp::complex MPll::H_A(double q2)
1346{
1347 return (-C_10 + etaP * pow(-1, angmomP) * C_10p) *V_L(q2);
1348}
1349
1350gslpp::complex MPll::H_S(double q2)
1351{
1352 return MboMW * (C_S - etaP * pow(-1, angmomP) * C_Sp) * S_L(q2);
1353}
1354
1355gslpp::complex MPll::H_P(double q2)
1356{
1357 return ( MboMW * (C_P - etaP * pow(-1, angmomP) * C_Pp) + twoMlepMb / q2 * (C_10 * (1. + etaP * pow(-1, angmomP) * MsoMb) - C_10p * (etaP * pow(-1, angmomP) + MsoMb))) * S_L(q2);
1358}
1359
1360gslpp::complex MPll::H_nunu(double q2, QCD::lepton lep)
1361{
1362 if (lep == QCD::NEUTRINO_1)
1363 return -(C_L_nunu_e - etaP * pow(-1, angmomP) * C_R_nunu_e) * V_L(q2);
1364 else if (lep == QCD::NEUTRINO_2)
1365 return -(C_L_nunu_mu - etaP * pow(-1, angmomP) * C_R_nunu_mu) * V_L(q2);
1366 else if (lep == QCD::NEUTRINO_3)
1367 return -(C_L_nunu_tau - etaP * pow(-1, angmomP) * C_R_nunu_tau) * V_L(q2);
1368 else throw std::runtime_error("MPll::H_nunu: lepton not supported");
1369}
1370
1371/*******************************************************************************
1372 * Angular coefficients *
1373 * ****************************************************************************/
1374double MPll::k2(double q2)
1375{
1376 return (MM4 + q2 * q2 + MP4 - twoMP2 * q2 - twoMM2 * (q2 + MP2)) / fourMM2;
1377}
1378
1379double MPll::beta(double q2)
1380{
1381 return sqrt(1. - 4. * Mlep2 / q2);
1382}
1383
1384double MPll::beta2(double q2)
1385{
1386 return 1. - 4. * Mlep2 / q2;
1387}
1388
1389double MPll::lambda(double q2)
1390{
1391 return 4. * MM2 * k2(q2);
1392}
1393
1394double MPll::F(double q2)
1395{
1396 return sqrt(lambda(q2)) * beta(q2) * q2 / (ninetysixM_PI3MM3);
1397}
1398
1399double MPll::I_1c(double q2)
1400{
1401 if (lep == QCD::NEUTRINO_1) return F(q2)*H_nunu(q2,QCD::NEUTRINO_1).abs2() + F(q2)*H_nunu(q2,QCD::NEUTRINO_2).abs2() + F(q2)*H_nunu(q2,QCD::NEUTRINO_3).abs2();
1402
1403 else return F(q2)*((H_V(q2).abs2() + H_A(q2).abs2()) / 2. + H_P(q2).abs2() + 2. * Mlep2 / q2 * (H_V(q2).abs2()
1404 - H_A(q2).abs2()) + beta2(q2) * H_S(q2).abs2());
1405}
1406
1407double MPll::I_2c(double q2)
1408{
1409 if (lep == QCD::NEUTRINO_1) return - F(q2)*H_nunu(q2,QCD::NEUTRINO_1).abs2() - F(q2)*H_nunu(q2,QCD::NEUTRINO_2).abs2() - F(q2)*H_nunu(q2,QCD::NEUTRINO_3).abs2();
1410
1411 else return -F(q2) * beta2(q2) / 2. * (H_V(q2).abs2() + H_A(q2).abs2());
1412}
1413
1414double MPll::I_6c(double q2)
1415{
1416 return 4. * F(q2) * beta(q2) * Mlep / sqrt(q2)*(H_S(q2).conjugate() * H_V(q2)).real();
1417}
1418
1419double MPll::Delta(int i, double q2)
1420{
1421 return 0; /* FIX CPV */
1422 //return (I(i, q2,0) - I(i, q2,1))/2;
1423}
1424
1425double MPll::integrateSigma(int i, double q_min, double q_max)
1426{
1427 updateParameters();
1428
1429 std::pair<double, double > qbin = std::make_pair(q_min, q_max);
1430
1431 old_handler = gsl_set_error_handler_off();
1432
1433 switch (i) {
1434 case 0:
1435 if (sigma0Cached[qbin] == 0) {
1436 FS = convertToGslFunction(bind(&MPll::getSigma1c, &(*this), _1));
1437 if (gsl_integration_cquad(&FS, q_min, q_max, 1.e-2, 1.e-1, w_sigma, &avaSigma, &errSigma, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
1438 cacheSigma0[qbin] = NN*avaSigma;
1439 sigma0Cached[qbin] = 1;
1440 }
1441 return cacheSigma0[qbin];
1442 break;
1443 case 2:
1444 if (sigma2Cached[qbin] == 0) {
1445 FS = convertToGslFunction(bind(&MPll::getSigma2c, &(*this), _1));
1446 if (gsl_integration_cquad(&FS, q_min, q_max, 1.e-2, 1.e-1, w_sigma, &avaSigma, &errSigma, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
1447 cacheSigma2[qbin] = NN*avaSigma;
1448 sigma2Cached[qbin] = 1;
1449 }
1450 return cacheSigma2[qbin];
1451 break;
1452 case 8:
1453 if (sigma8Cached[qbin] == 0) {
1454 FS = convertToGslFunction(bind(&MPll::getSigma6c, &(*this), _1));
1455 if (gsl_integration_cquad(&FS, q_min, q_max, 1.e-2, 1.e-1, w_sigma, &avaSigma, &errSigma, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
1456 cacheSigma8[qbin] = NN*avaSigma;
1457 sigma8Cached[qbin] = 1;
1458 }
1459 return cacheSigma8[qbin];
1460 break;
1461 default:
1462 std::stringstream out;
1463 out << i;
1464 throw std::runtime_error("MPll::integrateSigma: index " + out.str() + " not implemented");
1465 }
1466
1467 gsl_set_error_handler(old_handler);
1468
1469}
1470
1471double MPll::getSigma(int i, double q_2)
1472{
1473 updateParameters();
1474
1475 switch (i) {
1476 case 0:
1477 return getSigma1c(q_2);
1478 break;
1479 case 2:
1480 return getSigma2c(q_2);
1481 break;
1482 case 8:
1483 return getSigma6c(q_2);
1484 break;
1485 default:
1486 std::stringstream out;
1487 out << i;
1488 throw std::runtime_error("MPll::getSigma: index " + out.str() + " not implemented");
1489 }
1490}
1491
1492double MPll::integrateDelta(int i, double q_min, double q_max)
1493{
1494 updateParameters();
1495
1496 std::pair<double, double > qbin = std::make_pair(q_min, q_max);
1497
1498 old_handler = gsl_set_error_handler_off();
1499
1500 switch (i) {
1501 case 0:
1502 if (delta0Cached[qbin] == 0) {
1503 FD = convertToGslFunction(bind(&MPll::getDelta1c, &(*this), _1));
1504 if (gsl_integration_cquad(&FD, q_min, q_max, 1.e-2, 1.e-1, w_delta, &avaDelta, &errDelta, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
1505 cacheDelta0[qbin] = NN*avaDelta;
1506 delta0Cached[qbin] = 1;
1507 }
1508 return cacheDelta0[qbin];
1509 break;
1510 case 2:
1511 if (delta2Cached[qbin] == 0) {
1512 FD = convertToGslFunction(bind(&MPll::getDelta2c, &(*this), _1));
1513 if (gsl_integration_cquad(&FD, q_min, q_max, 1.e-2, 1.e-1, w_delta, &avaDelta, &errDelta, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
1514 cacheDelta2[qbin] = NN*avaDelta;
1515 delta2Cached[qbin] = 1;
1516 }
1517 return cacheDelta2[qbin];
1518 break;
1519 default:
1520 std::stringstream out;
1521 out << i;
1522 throw std::runtime_error("MPll::integrateDelta: index " + out.str() + " not implemented");
1523 }
1524
1525 gsl_set_error_handler(old_handler);
1526
1527}
1528
1529double MPll::integrateSigmaTree(double q_min, double q_max)
1530{
1531 if (lep != QCD::NEUTRINO_1 or meson != QCD::B_P or !NeutrinoTree_flag) return 0.;
1532
1533 updateParameters();
1534
1535 //phase space limit where tree-level contribution is relevant (0908.1174)
1536 double q_cut = (mtau2 - MP2) * (MM2 - mtau2) / mtau2;
1537 if (q_max >= q_cut) {
1538 if (q_min == 0.) return getintegratedSigmaTree();
1539 q_max = q_cut;
1540 }
1541
1542 double prefactor = mySM.getMesons(meson).getLifetime() / HCUT * GF4 * VusVub_abs2 * fP2 * fM2 / (64. * M_PI2 * MM3 * Gammatau) * mtau2 * mtau;
1543
1544 std::pair<double, double > qbin = std::make_pair(q_min, q_max);
1545
1546 old_handler = gsl_set_error_handler_off();
1547
1548 if (sigmaTreeCached[qbin] == 0) {
1549 FD = convertToGslFunction(bind(&MPll::SigmaTree, &(*this), _1));
1550 if (gsl_integration_cquad(&FD, q_min, q_max, 1.e-2, 1.e-1, w_sigmaTree, &avaSigmaTree, &errSigmaTree, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
1551 cacheSigmaTree[qbin] = avaSigmaTree;
1552 sigmaTreeCached[qbin] = 1;
1553 }
1554 return prefactor * cacheSigmaTree[qbin];
1555
1556 gsl_set_error_handler(old_handler);
1557}
1558
1559double MPll::SigmaTree(double q2)
1560{
1561 return MM2 * (mtau2 - MP2) - mtau2 * (mtau2 + q2 - MP2);
1562}
1563
1564double MPll::getintegratedSigmaTree()
1565{
1566 return mySM.getMesons(meson).getLifetime() / HCUT * GF4 * VusVub_abs2 * fP2 * fM2 / (128. * M_PI2 * MM3 * Gammatau) * mtau * (mtau2 - MP2) * (mtau2 - MP2) * (MM2 - mtau2) * (MM2 - mtau2);
1567}
F_1.h
F_2.h
LO
@ LO
Definition: OrderScheme.h:34
NLO
@ NLO
Definition: OrderScheme.h:35
NLO_QED11
@ NLO_QED11
Definition: OrderScheme.h:59
CKM::computelamt_s
const gslpp::complex computelamt_s() const
The product of the CKM elements .
Definition: CKM.cpp:174
CKM::computelamu_s
const gslpp::complex computelamu_s() const
The product of the CKM elements .
Definition: CKM.cpp:184
F_1
Definition: F_1.h:15
F_2
Definition: F_2.h:15
Flavour::getFlagUseDispersionRelation
bool getFlagUseDispersionRelation() const
Definition: Flavour.h:343
Flavour::getFlagMPll_DM
bool getFlagMPll_DM() const
Definition: Flavour.h:375
Flavour::ComputeCoeffsnunu
gslpp::vector< gslpp::complex > ** ComputeCoeffsnunu(QCD::lepton lepton=QCD::NOLEPTON, bool noSM=false) const
Definition: Flavour.cpp:170
Flavour::getFlagFixedWCbtos
bool getFlagFixedWCbtos() const
Definition: Flavour.h:363
Flavour::getFlagNeutrinoTree
bool getFlagNeutrinoTree() const
Definition: Flavour.h:383
Flavour::setUpdateFlag
void setUpdateFlag(QCD::meson meson_i, QCD::meson meson_j, QCD::lepton lep_i, bool updated_i) const
sets the update flag for the initial and final state dependent object for .
Definition: Flavour.cpp:315
Flavour::getFlagMPll_FNALMILC
bool getFlagMPll_FNALMILC() const
Definition: Flavour.h:367
Flavour::ComputeCoeffprimeBMll
gslpp::vector< gslpp::complex > ** ComputeCoeffprimeBMll(double mu, QCD::lepton lepton, schemes scheme=NDR) const
Computes the chirality flipped Wilson coefficient for the process .
Definition: Flavour.cpp:205
Flavour::getUpdateFlag
bool getUpdateFlag(QCD::meson meson_i, QCD::meson meson_j, QCD::lepton lep_i) const
gets the update flag for the initial and final state dependent object for .
Definition: Flavour.cpp:340
Flavour::ComputeCoeffBMll
gslpp::vector< gslpp::complex > ** ComputeCoeffBMll(double mu, QCD::lepton lepton, bool noSM=false, schemes scheme=NDR) const
Computes the Wilson coefficient for the process .
Definition: Flavour.cpp:195
Flavour::getFlagMPll_GRvDV
bool getFlagMPll_GRvDV() const
Definition: Flavour.h:371
MPll::initializeMPllParameters
std::vector< std::string > initializeMPllParameters()
A method for initializing the parameters necessary for MPll.
Definition: MPll.cpp:80
MPll::Integrand_ReTpar_pm
double Integrand_ReTpar_pm(double up)
The sum of Integrand_ReTparplus() and Integrand_ReTparminus().
Definition: MPll.cpp:1058
MPll::getSigma1c
double getSigma1c(double q2)
The CP average .
Definition: MPll.h:1016
MPll::FixedWCbtos
bool FixedWCbtos
Definition: MPll.h:323
MPll::SigmaTree
double SigmaTree(double q2)
Definition: MPll.cpp:1559
MPll::MPll_GRvDV_flag
bool MPll_GRvDV_flag
Definition: MPll.h:325
MPll::H_A
gslpp::complex H_A(double q2)
The helicity amplitude .
Definition: MPll.cpp:1345
MPll::Mb
double Mb
Definition: MPll.h:335
MPll::h_lambda
gslpp::complex h_lambda(double q2)
The non-pertubative ccbar contributions to the helicity amplitudes.
Definition: MPll.cpp:1321
MPll::F89
double F89(double q2)
The correction from .
Definition: MPll.cpp:1102
MPll::mc_pole
double mc_pole
Definition: MPll.h:340
MPll::deltaC9_QCDF
gslpp::complex deltaC9_QCDF(double q2, bool spline=false)
QCDF Correction from various BFS papers (hep-ph/0403185, hep-ph/0412400) and Greub et....
Definition: MPll.cpp:1223
MPll::DeltaC9_KD
gslpp::complex DeltaC9_KD(double q2)
Definition: MPll.cpp:1315
MPll::H_P
gslpp::complex H_P(double q2)
The helicity amplitude .
Definition: MPll.cpp:1355
MPll::pseudoscalar
QCD::meson pseudoscalar
Definition: MPll.h:318
MPll::H_S
gslpp::complex H_S(double q2)
The helicity amplitude .
Definition: MPll.cpp:1350
MPll::F87
gslpp::complex F87(double q2)
The correction from .
Definition: MPll.cpp:1095
MPll::mpllParameters
std::vector< std::string > mpllParameters
Definition: MPll.h:319
MPll::Integrand_ImTparminus
double Integrand_ImTparminus(double up)
The imaginary part of the integral involving at fixed , according to .
Definition: MPll.cpp:1042
MPll::myF_2
std::unique_ptr< F_2 > myF_2
Definition: MPll.h:321
MPll::MPll_Lattice_flag
bool MPll_Lattice_flag
Definition: MPll.h:324
MPll::Mlep
double Mlep
Definition: MPll.h:332
MPll::mu_h
double mu_h
Definition: MPll.h:337
MPll::I1
gslpp::complex I1(double u, double q2)
The function from .
Definition: MPll.cpp:973
MPll::Integrand_ReTparminus
double Integrand_ReTparminus(double up)
The real part of the integral involving at fixed , according to .
Definition: MPll.cpp:1031
MPll::MPll_DM_flag
bool MPll_DM_flag
Definition: MPll.h:326
MPll::fit_DeltaC9_mumu
void fit_DeltaC9_mumu()
The fitting routine for the QCDF correction in the muon channel.
Definition: MPll.cpp:1152
MPll::Tparminus
gslpp::complex Tparminus(double u, double q2)
The function from .
Definition: MPll.cpp:1008
MPll::imDC9fit
double imDC9fit(double *x, double *p)
The fit function for the imaginary part of the QCDF correction .
Definition: MPll.cpp:1145
MPll::getDelta2c
double getDelta2c(double q2)
The CP asymmetry .
Definition: MPll.h:1056
MPll::meson
QCD::meson meson
Definition: MPll.h:317
MPll::width
double width
Definition: MPll.h:343
MPll::F19
gslpp::complex F19(double q2)
The correction from .
Definition: MPll.cpp:1063
MPll::integrateSigma
double integrateSigma(int i, double q_min, double q_max)
The integral of from to .
Definition: MPll.cpp:1425
MPll::Cpar
gslpp::complex Cpar(double q2)
The correction from .
Definition: MPll.cpp:1109
MPll::reDC9fit
double reDC9fit(double *x, double *p)
The fit function for the real part of the QCDF correction .
Definition: MPll.cpp:1140
MPll::Mc
double Mc
Definition: MPll.h:338
MPll::mySM
const StandardModel & mySM
Definition: MPll.h:315
MPll::V_L
gslpp::complex V_L(double q2)
The helicity form factor .
Definition: MPll.cpp:954
MPll::H_nunu
gslpp::complex H_nunu(double q2, QCD::lepton lep)
The helicity amplitude for the invisible decay.
Definition: MPll.cpp:1360
MPll::getSigma
double getSigma(int i, double q_2)
The value of from to .
Definition: MPll.cpp:1471
MPll::lep
QCD::lepton lep
Definition: MPll.h:316
MPll::~MPll
virtual ~MPll()
Destructor.
Definition: MPll.cpp:76
MPll::funct_g
gslpp::complex funct_g(double q2)
Definition: MPll.cpp:1307
MPll::mJ2
double mJ2
Definition: MPll.h:328
MPll::GF
double GF
Definition: MPll.h:330
MPll::Integrand_ImTpar_pm
double Integrand_ImTpar_pm(double up)
The sum of Integrand_ImTparplus() and Integrand_ImTparminus().
Definition: MPll.cpp:1053
MPll::Integrand_ReTparplus
double Integrand_ReTparplus(double up)
The real part of the integral involving at fixed , according to .
Definition: MPll.cpp:1015
MPll::T_L
gslpp::complex T_L(double q2)
The helicity form factor .
Definition: MPll.cpp:959
MPll::alpha_s_mub
double alpha_s_mub
Definition: MPll.h:344
MPll::getSigma2c
double getSigma2c(double q2)
The CP average .
Definition: MPll.h:1026
MPll::Tparplus
gslpp::complex Tparplus(double u, double q2)
The function from .
Definition: MPll.cpp:995
MPll::MM
double MM
Definition: MPll.h:333
MPll::myF_1
std::unique_ptr< F_1 > myF_1
Definition: MPll.h:320
MPll::NeutrinoTree_flag
bool NeutrinoTree_flag
Definition: MPll.h:327
MPll::integrateSigmaTree
double integrateSigmaTree(double q_min, double q_max)
The integral of from to (arxiv/2301.06990)
Definition: MPll.cpp:1529
MPll::getintegratedSigmaTree
double getintegratedSigmaTree()
The integral of from 0 to .
Definition: MPll.cpp:1564
MPll::Ms
double Ms
Definition: MPll.h:341
MPll::getDelta1c
double getDelta1c(double q2)
The CP asymmetry .
Definition: MPll.h:1046
MPll::MP
double MP
Definition: MPll.h:334
MPll::H_V
gslpp::complex H_V(double q2)
The helicity amplitude .
Definition: MPll.cpp:1338
MPll::DeltaC9
gslpp::complex DeltaC9(double q2)
The total QCDF correction computed integrating over .
Definition: MPll.cpp:1193
MPll::deltaC7_QCDF
gslpp::complex deltaC7_QCDF(double q2, bool spline=false)
QCDF Correction from various BFS papers (hep-ph/0403185, hep-ph/0412400) and Greub et....
Definition: MPll.cpp:1198
MPll::dispersion
bool dispersion
Definition: MPll.h:322
MPll::getSigma6c
double getSigma6c(double q2)
The CP average .
Definition: MPll.h:1036
MPll::mb_pole
double mb_pole
Definition: MPll.h:339
MPll::F27
gslpp::complex F27(double q2)
The correction from .
Definition: MPll.cpp:1077
MPll::ale
double ale
Definition: MPll.h:331
MPll::mu_b
double mu_b
Definition: MPll.h:336
MPll::fDeltaC9
gslpp::complex fDeltaC9(double q2)
The total QCDF correction computed fitting over .
Definition: MPll.cpp:1184
MPll::S_L
double S_L(double q2)
The helicity form factor .
Definition: MPll.cpp:964
MPll::spectator_charge
double spectator_charge
Definition: MPll.h:342
MPll::MPll
MPll(const StandardModel &SM_i, QCD::meson meson_i, QCD::meson pseudoscalar_i, QCD::lepton lep_i)
Constructor.
Definition: MPll.cpp:20
MPll::deltaTpar
gslpp::complex deltaTpar(double q2)
The total correction from .
Definition: MPll.cpp:1115
MPll::integrateDelta
double integrateDelta(int i, double q_min, double q_max)
The integral of from to .
Definition: MPll.cpp:1492
MPll::F29
gslpp::complex F29(double q2)
The correction from .
Definition: MPll.cpp:1086
MPll::spline_QCDF_func
void spline_QCDF_func()
Definition: MPll.cpp:1248
MPll::Integrand_ImTparplus
double Integrand_ImTparplus(double up)
The imaginary part of the integral involving at fixed , according to .
Definition: MPll.cpp:1023
Meson::getLambdaM
const double & getLambdaM() const
Definition: Meson.h:402
Meson::computeWidth
double computeWidth() const
A method to compute the width of the meson from its lifetime.
Definition: Meson.cpp:521
Meson::getDecayconst
const double & getDecayconst() const
A get method for the decay constant of the meson.
Definition: Meson.h:360
Meson::getGegenalpha
const double & getGegenalpha(int i) const
A get method to get the Gegenbaur coefficient.
Definition: Meson.h:394
Meson::getLifetime
double getLifetime() const
A get method for the lifetime of the meson.
Definition: Meson.h:351
Particle::getMass
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
Particle::getCharge
double getCharge() const
A get method to access the particle charge.
Definition: Particle.h:97
QCD::meson
meson
An enum type for mesons.
Definition: QCD.h:336
QCD::K_P
@ K_P
Definition: QCD.h:340
QCD::B_P
@ B_P
Definition: QCD.h:345
QCD::K_0
@ K_0
Definition: QCD.h:339
QCD::getOptionalParameter
const double getOptionalParameter(std::string name) const
A method to get parameters that are specific to only one set of observables.
Definition: QCD.h:450
QCD::getMesons
const Meson & getMesons(const QCD::meson m) const
A get method to access a meson as an object of the type Meson.
Definition: QCD.h:526
QCD::UP
@ UP
Definition: QCD.h:324
QCD::BOTTOM
@ BOTTOM
Definition: QCD.h:329
QCD::DOWN
@ DOWN
Definition: QCD.h:325
QCD::STRANGE
@ STRANGE
Definition: QCD.h:327
QCD::CHARM
@ CHARM
Definition: QCD.h:326
QCD::lepton
lepton
An enum type for leptons.
Definition: QCD.h:310
QCD::NEUTRINO_2
@ NEUTRINO_2
Definition: QCD.h:313
QCD::NEUTRINO_1
@ NEUTRINO_1
Definition: QCD.h:311
QCD::NEUTRINO_3
@ NEUTRINO_3
Definition: QCD.h:315
QCD::TAU
@ TAU
Definition: QCD.h:316
QCD::getQuarks
const Particle & getQuarks(const QCD::quark q) const
A get method to access a quark as an object of the type Particle.
Definition: QCD.h:536
QCD::initializeMeson
void initializeMeson(QCD::meson meson_i) const
A method to initialize a meson.
Definition: QCD.cpp:280
QCD::getMub
const double getMub() const
A get method to access the threshold between five- and four-flavour theory in GeV.
Definition: QCD.h:582
QCD::Mbar2Mp
const double Mbar2Mp(const double mbar, const quark q, const orders order=FULLNNLO) const
Converts the mass to the pole mass.
Definition: QCD.cpp:1552
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:521
StandardModel::getLeptons
const Particle & getLeptons(const QCD::lepton p) const
A get method to retrieve the member object of a lepton.
Definition: StandardModel.h:756
StandardModel::getCKM
const CKM & getCKM() const
A get method to retrieve the member object of type CKM.
Definition: StandardModel.h:952
StandardModel::Mw
virtual const double Mw() const
The SM prediction for the -boson mass in the on-shell scheme, .
Definition: StandardModel/src/StandardModel.cpp:1025
StandardModel::getFlavour
const Flavour & getFlavour() const
Definition: StandardModel.h:1065
StandardModel::getGF
const double getGF() const
A get method to retrieve the Fermi constant .
Definition: StandardModel.h:792
StandardModel::Als
const double Als(const double mu, const orders order, const bool Nf_thr, const bool qed_flag) const
The running QCD coupling in the scheme including QED corrections.
Definition: StandardModel.h:1120
StandardModel::getAle
const double getAle() const
A get method to retrieve the fine-structure constant .
Definition: StandardModel.h:801
h_0
A class for the correction in .
Definition: MVllObservables.h:1991
s
Test Observable.