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 logMc = log(Mc2 / mu_b2);
360 logMb = log(Mb2 / mu_b2);
361 mu_b2 = mu_b*mu_b;
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
427 if (I0_updated == 0) for (it = delta0Cached.begin(); it != delta0Cached.end(); ++it) it->second = 0;
428 if (I2_updated == 0) for (it = delta2Cached.begin(); it != delta2Cached.end(); ++it) it->second = 0;
429
430 if (Itree_updated) for (it = sigmaTreeCached.begin(); it != sigmaTreeCached.end(); ++it) it->second = 0;
431
432 std::map<double, unsigned int >::iterator iti;
433 if (deltaTparpupdated == 0) for (iti = deltaTparpCached.begin(); iti != deltaTparpCached.end(); ++iti) iti->second = 0;
434 if (deltaTparmupdated == 0) for (iti = deltaTparmCached.begin(); iti != deltaTparmCached.end(); ++iti) iti->second = 0;
435
436 if (deltaTparpupdated * deltaTparmupdated == 0) for (it = I1Cached.begin(); it != I1Cached.end(); ++it) it->second = 0;
437
438#if SPLINE
439 if (mySM.getFlavour().getUpdateFlag(meson, pseudoscalar, lep)) spline_QCDF_func();
440#else
441 fit_DeltaC9_mumu();
442#endif
443
444 mySM.getFlavour().setUpdateFlag(meson, pseudoscalar, lep, false);
445
446 //std::cout << "fplus_DM(4)= " << f_plus(4) << std::endl;
447 //std::cout << "f0_DM(4)= " << f_0(4) << std::endl;
448 //std::cout << "fT_DM(4)= " << f_T(4) << std::endl;
449
450 return;
451
452}
453
454void MPll::checkCache()
455{
456
457 if (MM == k2_cache(0) && MP == k2_cache(1)) {
458 k2_updated = 1;
459 } else {
460 k2_updated = 0;
461 k2_cache(0) = MM;
462 k2_cache(1) = MP;
463 }
464
465 if (MPll_Lattice_flag || MPll_GRvDV_flag || MPll_DM_flag) {
466 if (b_0_fplus == fplus_lat_cache(0) && b_1_fplus == fplus_lat_cache(1) && b_2_fplus == fplus_lat_cache(2)) {
467 fplus_lat_updated = 1;
468 } else {
469 fplus_lat_updated = 0;
470 fplus_lat_cache(0) = b_0_fplus;
471 fplus_lat_cache(1) = b_1_fplus;
472 fplus_lat_cache(2) = b_2_fplus;
473 }
474
475 if (b_0_fT == fT_lat_cache(0) && b_1_fT == fT_lat_cache(1) && b_2_fT == fT_lat_cache(2)) {
476 fT_lat_updated = 1;
477 } else {
478 fT_lat_updated = 0;
479 fT_lat_cache(0) = b_0_fT;
480 fT_lat_cache(1) = b_1_fT;
481 fT_lat_cache(2) = b_2_fT;
482 }
483
484 if (b_0_f0 == f0_lat_cache(0) && b_1_f0 == f0_lat_cache(1) && b_2_f0 == f0_lat_cache(2)) {
485 f0_lat_updated = 1;
486 } else {
487 f0_lat_updated = 0;
488 f0_lat_cache(0) = b_0_f0;
489 f0_lat_cache(1) = b_1_f0;
490 f0_lat_cache(2) = b_2_f0;
491 }
492 } else {
493 if (r_1_fplus == fplus_cache(0) && r_2_fplus == fplus_cache(1)) {
494 fplus_updated = 1;
495 } else {
496 fplus_updated = 0;
497 fplus_cache(0) = r_1_fplus;
498 fplus_cache(1) = r_2_fplus;
499 }
500
501 if (r_1_fT == fT_cache(0) && r_2_fT == fT_cache(1)) {
502 fT_updated = 1;
503 } else {
504 fT_updated = 0;
505 fT_cache(0) = r_1_fT;
506 fT_cache(1) = r_2_fT;
507 }
508
509 if (r_2_f0 == f0_cache) {
510 f0_updated = 1;
511 } else {
512 f0_updated = 0;
513 f0_cache = r_2_f0;
514 }
515 }
516
517 if (Mlep == beta_cache) {
518 beta_updated = 1;
519 } else {
520 beta_updated = 0;
521 beta_cache = Mlep;
522 }
523
524 lambda_updated = k2_updated;
525 F_updated = lambda_updated * beta_updated;
526
527 if (MPll_Lattice_flag || MPll_GRvDV_flag || MPll_DM_flag) {
528 VL_updated = k2_updated * fplus_lat_updated;
529 TL_updated = k2_updated * fT_lat_updated;
530 } else {
531 VL_updated = k2_updated * fplus_updated;
532 TL_updated = k2_updated * fT_updated;
533 }
534
535 if (Mb == SL_cache(0) && Ms == SL_cache(1)) {
536 if (MPll_Lattice_flag || MPll_GRvDV_flag || MPll_DM_flag)
537 SL_updated = k2_updated * f0_lat_updated;
538 else
539 SL_updated = k2_updated * f0_updated;
540 } else {
541 SL_updated = 0;
542 SL_cache(0) = Mb;
543 SL_cache(1) = Ms;
544 }
545
546 if (GF == N_cache(0) && ale == N_cache(1) && MM == N_cache(2) && lambda_t == Nc_cache) {
547 N_updated = 1;
548 } else {
549 N_updated = 0;
550 N_cache(0) = GF;
551 N_cache(1) = ale;
552 N_cache(2) = MM;
553 Nc_cache = lambda_t;
554 }
555
556 if (C_1 == C_1_cache) {
557 C_1_updated = 1;
558 } else {
559 C_1_updated = 0;
560 C_1_cache = C_1;
561 }
562
563 if (C_2 == C_2_cache) {
564 C_2_updated = 1;
565 } else {
566 C_2_updated = 0;
567 C_2_cache = C_2;
568 }
569
570 if (C_3 == C_3_cache) {
571 C_3_updated = 1;
572 } else {
573 C_3_updated = 0;
574 C_3_cache = C_3;
575 }
576
577 if (C_4 == C_4_cache) {
578 C_4_updated = 1;
579 } else {
580 C_4_updated = 0;
581 C_4_cache = C_4;
582 }
583
584 if (C_5 == C_5_cache) {
585 C_5_updated = 1;
586 } else {
587 C_5_updated = 0;
588 C_5_cache = C_5;
589 }
590
591 if (C_6 == C_6_cache) {
592 C_6_updated = 1;
593 } else {
594 C_6_updated = 0;
595 C_6_cache = C_6;
596 }
597
598 if (C_7 == C_7_cache) {
599 C_7_updated = 1;
600 } else {
601 C_7_updated = 0;
602 C_7_cache = C_7;
603 }
604
605 if (C_9 == C_9_cache) {
606 C_9_updated = 1;
607 } else {
608 C_9_updated = 0;
609 C_9_cache = C_9;
610 }
611
612 if (C_10 == C_10_cache) {
613 C_10_updated = 1;
614 } else {
615 C_10_updated = 0;
616 C_10_cache = C_10;
617 }
618
619 if (C_S == C_S_cache) {
620 C_S_updated = 1;
621 } else {
622 C_S_updated = 0;
623 C_S_cache = C_S;
624 }
625
626 if (C_P == C_P_cache) {
627 C_P_updated = 1;
628 } else {
629 C_P_updated = 0;
630 C_P_cache = C_P;
631 }
632
633 if (C_7p == C_7p_cache) {
634 C_7p_updated = 1;
635 } else {
636 C_7p_updated = 0;
637 C_7p_cache = C_7p;
638 }
639
640 if (C_9p == C_9p_cache) {
641 C_9p_updated = 1;
642 } else {
643 C_9p_updated = 0;
644 C_9p_cache = C_9p;
645 }
646
647 if (C_10p == C_10p_cache) {
648 C_10p_updated = 1;
649 } else {
650 C_10p_updated = 0;
651 C_10p_cache = C_10p;
652 }
653
654 if (C_Sp == C_Sp_cache) {
655 C_Sp_updated = 1;
656 } else {
657 C_Sp_updated = 0;
658 C_Sp_cache = C_Sp;
659 }
660
661 if (C_Pp == C_Pp_cache) {
662 C_Pp_updated = 1;
663 } else {
664 C_Pp_updated = 0;
665 C_Pp_cache = C_Pp;
666 }
667
668 if (C_2Lh_bar == C_2Lh_cache) {
669 C_2Lh_updated = 1;
670 } else {
671 C_2Lh_updated = 0;
672 C_2Lh_cache = C_2Lh_bar;
673 }
674
675 if (C_8Lh == C_8Lh_cache) {
676 C_8Lh_updated = 1;
677 } else {
678 C_8Lh_updated = 0;
679 C_8Lh_cache = C_8Lh;
680 }
681
682 if (C_L_nunu_e == C_L_nunu_e_cache) {
683 C_L_nunu_e_updated = 1;
684 } else {
685 C_L_nunu_e_updated = 0;
686 C_L_nunu_e_cache = C_L_nunu_e;
687 }
688
689 if (C_L_nunu_mu == C_L_nunu_mu_cache) {
690 C_L_nunu_mu_updated = 1;
691 } else {
692 C_L_nunu_mu_updated = 0;
693 C_L_nunu_mu_cache = C_L_nunu_mu;
694 }
695
696 if (C_L_nunu_tau == C_L_nunu_tau_cache) {
697 C_L_nunu_tau_updated = 1;
698 } else {
699 C_L_nunu_tau_updated = 0;
700 C_L_nunu_tau_cache = C_L_nunu_tau;
701 }
702
703 if (C_R_nunu_e == C_R_nunu_e_cache) {
704 C_R_nunu_e_updated = 1;
705 } else {
706 C_R_nunu_e_updated = 0;
707 C_R_nunu_e_cache = C_R_nunu_e;
708 }
709
710 if (C_R_nunu_mu == C_R_nunu_mu_cache) {
711 C_R_nunu_mu_updated = 1;
712 } else {
713 C_R_nunu_mu_updated = 0;
714 C_R_nunu_mu_cache = C_R_nunu_mu;
715 }
716
717 if (C_R_nunu_tau == C_R_nunu_tau_cache) {
718 C_R_nunu_tau_updated = 1;
719 } else {
720 C_R_nunu_tau_updated = 0;
721 C_R_nunu_tau_cache = C_R_nunu_tau;
722 }
723
724 if (Mb == Ycache(0) && Mc == Ycache(1)) {
725 Yupdated = C_1_updated * C_2_updated * C_3_updated * C_4_updated * C_5_updated * C_6_updated;
726 } else {
727 Yupdated = 0;
728 Ycache(0) = Mb;
729 Ycache(1) = Mc;
730 }
731
732 if (lep == QCD::NEUTRINO_1){
733 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;
734 H_A0updated = H_V0updated;
735 } else {
736
737 if (!dispersion) {
738 if (MM == H_V0cache(0) && Mb == H_V0cache(1) && h_0 == H_V0Ccache[0] && h_1 == H_V0Ccache[1] && h_2 == H_V0Ccache[2]) {
739 H_V0updated = N_updated * C_9_updated * Yupdated * VL_updated * C_9p_updated * C_7_updated * TL_updated * C_7p_updated;
740 } else {
741 H_V0updated = 0;
742 H_V0cache(0) = MM;
743 H_V0cache(1) = Mb;
744 H_V0Ccache[0] = h_0;
745 H_V0Ccache[1] = h_1;
746 H_V0Ccache[2] = h_2;
747 }
748 } else {
749 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]) {
750 H_V0updated = N_updated * C_9_updated * Yupdated * VL_updated * C_9p_updated * C_7_updated * TL_updated * C_7p_updated;
751 } else {
752 H_V0updated = 0;
753 H_V0cache(0) = MM;
754 H_V0cache(1) = Mb;
755 H_V0Ccache_dispersion[0] = r_1;
756 H_V0Ccache_dispersion[1] = r_2;
757 H_V0Ccache_dispersion[2] = Delta_C9;
758 H_V0Ccache_dispersion[3] = exp_Phase;
759 }
760 }
761
762 H_A0updated = N_updated * C_10_updated * VL_updated * C_10p_updated;
763 }
764
765 if (Mb == H_Scache(0) && MW == H_Scache(1)) {
766 H_Supdated = N_updated * C_S_updated * SL_updated * C_Sp_updated;
767 } else {
768 H_Supdated = 0;
769 H_Scache(0) = Mb;
770 H_Scache(1) = MW;
771 }
772
773 if (Mb == H_P_cache(0) && MW == H_P_cache(1) && Mlep == H_P_cache(2) && Ms == H_P_cache(3)) {
774 H_P_updated = N_updated * C_P_updated * SL_updated * C_Pp_updated * SL_updated * C_10_updated * C_10p_updated;
775 } else {
776 H_P_updated = 0;
777 H_P_cache(0) = Mb;
778 H_P_cache(1) = MW;
779 H_P_cache(2) = Mlep;
780 H_P_cache(3) = Ms;
781 }
782
783 if (MM == T_cache(0) && Mb == T_cache(1) && Mc == T_cache(2) &&
784 mySM.getMesons(pseudoscalar).getGegenalpha(0) == T_cache(3) && mySM.getMesons(pseudoscalar).getGegenalpha(1) == T_cache(4)) {
785 T_updated = 1;
786 } else {
787 T_updated = 0;
788 T_cache(0) = MM;
789 T_cache(1) = Mb;
790 T_cache(2) = Mc;
791 T_cache(3) = mySM.getMesons(pseudoscalar).getGegenalpha(0);
792 T_cache(4) = mySM.getMesons(pseudoscalar).getGegenalpha(1);
793 }
794
795 deltaTparpupdated = C_2Lh_updated * T_updated;
796 deltaTparmupdated = C_2Lh_updated * C_8Lh_updated * T_updated;
797
798 I0_updated = F_updated * H_V0updated * H_A0updated * H_P_updated * beta_updated * H_Supdated * deltaTparmupdated;
799 I2_updated = F_updated * beta_updated * H_V0updated * H_A0updated * deltaTparmupdated;
800 I8_updated = F_updated * beta_updated * H_Supdated * H_V0updated * deltaTparmupdated;
801
802 if (MM2 == Itree_cache(0) && mtau2 == Itree_cache(1) && MP2 == Itree_cache(2)) {
803 Itree_updated = 1;
804 } else {
805 Itree_updated = 0;
806 Itree_cache(0) = MM2;
807 Itree_cache(1) = mtau2;
808 Itree_cache(2) = MP2;
809 }
810
811}
812
813/*******************************************************************************
814 * Transverse Form Factors *
815 * ****************************************************************************/
816double MPll::LCSR_fit1(double q2, double r_1, double r_2, double m_fit2)
817{
818 return r_1 / (1. - q2 / m_fit2) + r_2 / pow((1. - q2 / m_fit2), 2.);
819
820}
821
822double MPll::LCSR_fit2(double q2, double r_2, double m_fit2)
823{
824 return r_2 / (1. - q2 / m_fit2);
825}
826
827double MPll::LCSR_fit3(double q2, double b_0, double b_1, double b_2, double m_fit2)
828{
829 return 1. / (1. - q2 / m_fit2) * (b_0 + b_1 * (zeta(q2) - zeta(0)) + b_2 * (zeta(q2) - zeta(0)) * (zeta(q2) - zeta(0)));
830}
831
832double MPll::zeta(double q2)
833{
834 double tp, t0;
835
836 tp = (MM + MP)*(MM + MP);
837 t0 = (MM + MP)*(sqrt(MM) - sqrt(MP))*(sqrt(MM) - sqrt(MP));
838
839 return (sqrt(tp - q2) - sqrt(tp - t0)) / (sqrt(tp - q2) + sqrt(tp - t0));
840}
841
842double MPll::zeta_DM(double q2)
843{
844 double tp, tm;
845
846 tp = (MM + MP)*(MM + MP);
847 tm = (MM - MP)*(MM - MP);
848
849 return (sqrt(tp - q2) - sqrt(tp - tm)) / (sqrt(tp - q2) + sqrt(tp - tm));
850}
851
852double MPll::phiplus_DM(double q2, double m_fit2)
853{
854 double z = zeta_DM(q2);
855 double z_M = zeta_DM(m_fit2);
856 double rP = MP/MM;
857 double Chi1minus = 0.000623174575;
858
859 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);
860}
861
862double MPll::phi0_DM(double q2)
863{
864 double z = zeta_DM(q2);
865 double rP = MP/MM;
866 double Chi0plus = 0.0142;
867
868 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);
869}
870
871double MPll::phiT_DM(double q2, double m_fit2)
872{
873 double z = zeta_DM(q2);
874 double z_M = zeta_DM(m_fit2);
875 double rP = MP/MM;
876 double ChiTT = 0.0454644444;
877
878 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);
879}
880
881double MPll::fplus_DM(double q2, double b_0, double b_1, double b_2, double m_fit2)
882{
883 double z = zeta_DM(q2);
884
885 return (b_0 + b_1*z + b_2*z*z) / phiplus_DM(q2, m_fit2);
886}
887
888double MPll::f0_DM(double q2, double b_0, double b_1, double b_2)
889{
890 double z = zeta_DM(q2);
891
892 return (b_0 + b_1*z + b_2*z*z) / phi0_DM(q2);
893}
894
895double MPll::fT_DM(double q2, double b_0, double b_1, double b_2, double m_fit2)
896{
897 double z = zeta_DM(q2);
898
899 return (b_0 + b_1*z + b_2*z*z) / phiT_DM(q2, m_fit2);
900}
901
902
903double MPll::LATTICE_fit1(double q2, double b_0, double b_1, double b_2, double m_fit2)
904{
905 double z2 = zeta(q2) * zeta(q2);
906 double z3 = zeta(q2) * z2;
907
908 return 1. / (1. - q2 / m_fit2) * (b_0 + b_1 * (zeta(q2) - 1. / 3. * z3) + b_2 * (z2 + 2. / 3. * z3));
909
910}
911
912double MPll::LATTICE_fit2(double q2, double b_0, double b_1, double b_2, double m_fit2)
913{
914 return 1. / (1. - q2 / m_fit2) * (b_0 + b_1 * zeta(q2) + b_2 * zeta(q2) * zeta(q2));
915}
916
917double MPll::f_plus(double q2)
918{
919 if (MPll_Lattice_flag)
920 return LATTICE_fit1(q2, b_0_fplus, b_1_fplus, b_2_fplus, m_fit2_fplus_lat);
921 else if (MPll_GRvDV_flag)
922 return LCSR_fit3(q2, b_0_fplus, b_1_fplus, b_2_fplus, m_fit2_fplus_lat);
923 else if (MPll_DM_flag)
924 return fplus_DM(q2, b_0_fplus, b_1_fplus, b_2_fplus, m_fit2_fplus_lat);
925 else
926 return LCSR_fit1(q2, r_1_fplus, r_2_fplus, m_fit2_fplus);
927}
928
929double MPll::f_T(double q2)
930{
931 if (MPll_Lattice_flag)
932 return LATTICE_fit1(q2, b_0_fT, b_1_fT, b_2_fT, m_fit2_fT_lat);
933 else if (MPll_GRvDV_flag)
934 return LCSR_fit3(q2, b_0_fT, b_1_fT, b_2_fT, m_fit2_fT_lat);
935 else if (MPll_DM_flag)
936 return fT_DM(q2, b_0_fT, b_1_fT, b_2_fT, m_fit2_fT_lat);
937 else
938 return LCSR_fit1(q2, r_1_fT, r_2_fT, m_fit2_fT);
939}
940
941double MPll::f_0(double q2)
942{
943 if (MPll_Lattice_flag)
944 return LATTICE_fit2(q2, b_0_f0, b_1_f0, b_2_f0, m_fit2_f0_lat);
945 else if (MPll_GRvDV_flag)
946 return LCSR_fit3(q2, b_0_f0, b_1_f0, b_2_f0, m_fit2_f0_lat);
947 else if (MPll_DM_flag)
948 return f0_DM(q2, b_0_f0, b_1_f0, b_2_f0);
949 else
950 return LCSR_fit2(q2, r_2_f0, m_fit2_f0);
951}
952
953gslpp::complex MPll::V_L(double q2)
954{
955 return /*gslpp::complex::i() */ sqrt(lambda(q2)) / (twoMM * sqrt(q2)) * f_plus(q2);
956}
957
958gslpp::complex MPll::T_L(double q2)
959{
960 return /*gslpp::complex::i() */ sqrt(lambda(q2) * q2) / (twoMM2_MMpMP) * f_T(q2);
961}
962
963double MPll::S_L(double q2)
964{
965 return S_L_pre * f_0(q2);
966}
967
968/*******************************************************************************
969 * QCDF *
970 * ****************************************************************************/
971
972gslpp::complex MPll::I1(double u, double q2)
973{
974 std::pair<double, double > uq2 = std::make_pair(u, q2);
975
976 if (I1Cached[uq2] == 0) {
977 ubar = 1. - u;
978 xp = .5 + sqrt(0.25 - (Mc2 - gslpp::complex::i()*1.e-10) / (ubar * MM2 + u * q2));
979 xm = .5 - sqrt(0.25 - (Mc2 - gslpp::complex::i()*1.e-10) / (ubar * MM2 + u * q2));
980 yp = .5 + sqrt(0.25 - (Mc2 - gslpp::complex::i()*1.e-10) / q2);
981 ym = .5 - sqrt(0.25 - (Mc2 - gslpp::complex::i()*1.e-10) / q2);
982 L1xp = log(1. - 1. / xp) * log(1. - xp) - M_PI2osix + dilog(xp / (xp - 1.));
983 L1xm = log(1. - 1. / xm) * log(1. - xm) - M_PI2osix + dilog(xm / (xm - 1.));
984 L1yp = log(1. - 1. / yp) * log(1. - yp) - M_PI2osix + dilog(yp / (yp - 1.));
985 L1ym = log(1. - 1. / ym) * log(1. - ym) - M_PI2osix + dilog(ym / (ym - 1.));
986
987 cacheI1[uq2] = 1. + twoMc2 / ubar / (MM2 - q2)*(L1xp + L1xm - L1yp - L1ym);
988 I1Cached[uq2] = 1;
989 }
990
991 return cacheI1[uq2];
992}
993
994gslpp::complex MPll::Tparplus(double u, double q2)
995{
996 Ee = (MM2 - q2) / twoMM;
997 ubar = 1. - u;
998 arg1 = (fourMc2 - gslpp::complex::i()*1.e-10) / (ubar * MM2 + u * q2) - 1.;
999 B01 = -2. * sqrt(arg1) * arctan(1. / sqrt(arg1));
1000 arg1 = (fourMc2 - gslpp::complex::i()*1.e-10) / q2 - 1.;
1001 B00 = -2. * sqrt(arg1) * arctan(1. / sqrt(arg1));
1002
1003 gslpp::complex tpar = twoMM / Ee / ubar * I1(u, q2) + (ubar * MM2 + u * q2) / Ee / Ee / ubar / ubar * (B01 - B00);
1004 return -MM / Mb * mySM.getQuarks(QCD::CHARM).getCharge() * tpar*C_2Lh_bar;
1005}
1006
1007gslpp::complex MPll::Tparminus(double u, double q2)
1008{
1009 double ubar = 1. - u;
1010 return -spectator_charge * (8. * C_8Lh / (ubar + u * q2 / MM2)
1011 + sixMMoMb * H_c(ubar * MM2 + u * q2, mu_h * mu_h) * C_2Lh_bar);
1012}
1013
1014double MPll::Integrand_ReTparplus(double up)
1015{
1016 double u = up;
1017 return ((Tparplus(u, tmpq2)*6. * u * (1. - u)*
1018 (1 + threeGegen0 * (2. * u - 1)
1019 + threeGegen1otwo * ((10. * u - 5.)*(2. * u - 1.) - 1.))) / mySM.getMesons(meson).getLambdaM()).real();
1020}
1021
1022double MPll::Integrand_ImTparplus(double up)
1023{
1024 double u = up;
1025 return ((Tparplus(u, tmpq2)*6. * u * (1. - u)*
1026 (1 + threeGegen0 * (2. * u - 1)
1027 + threeGegen1otwo * ((10. * u - 5.)*(2. * u - 1.) - 1.))) / mySM.getMesons(meson).getLambdaM()).imag();
1028}
1029
1030double MPll::Integrand_ReTparminus(double up)
1031{
1032 double Lambdaplus = mySM.getMesons(meson).getLambdaM();
1033 gslpp::complex Lambdamin = exp(-tmpq2 / MM / Lambdaplus) / Lambdaplus * (-gsl_sf_expint_Ei(tmpq2 / MM / Lambdaplus) + gslpp::complex::i() * M_PI);
1034
1035 double u = up;
1036 return ((Tparminus(u, tmpq2)*6. * u * (1. - u)*
1037 (1 + threeGegen0 * (2. * u - 1)
1038 + threeGegen1otwo * ((10. * u - 5.)*(2. * u - 1.) - 1.))) / Lambdamin).real();
1039}
1040
1041double MPll::Integrand_ImTparminus(double up)
1042{
1043 double Lambdaplus = mySM.getMesons(meson).getLambdaM();
1044 gslpp::complex Lambdamin = exp(-tmpq2 / MM / Lambdaplus) / Lambdaplus * (-gsl_sf_expint_Ei(tmpq2 / MM / Lambdaplus) + gslpp::complex::i() * M_PI);
1045
1046 double u = up;
1047 return ((Tparminus(u, tmpq2)*6. * u * (1. - u)*
1048 (1 + threeGegen0 * (2. * u - 1)
1049 + threeGegen1otwo * ((10. * u - 5.)*(2. * u - 1.) - 1.))) / Lambdamin).imag();
1050}
1051
1052double MPll::Integrand_ImTpar_pm(double up)
1053{
1054 return Integrand_ImTparplus(up) + Integrand_ImTparminus(up);
1055}
1056
1057double MPll::Integrand_ReTpar_pm(double up)
1058{
1059 return Integrand_ReTparplus(up) + Integrand_ReTparminus(up);
1060}
1061
1062gslpp::complex MPll::F19(double q2)
1063{
1064 double s = q2 / Mb2;
1065 double s2 = s*s;
1066 double Ls = log(s);
1067 double Lc = log(Mc / Mb);
1068 double Lm = log(mu_b / Mb);
1069 gslpp::complex i = gslpp::complex::i();
1070 return (-1424. / 729. + 16. / 243. * i * M_PI + 64. / 27. * Lc)*Lm - 16. / 243. * Lm * Ls + (16. / 1215. - 32. / 135. / Mc2 * Mb2) * Lm * s
1071 + (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
1072 + (-11.65 + 0.18223 * i + (-24.709 - 0.13474 * i) * s + (-43.588 - 0.4738 * i) * s2 + (-86.22 - 1.3542 * i) * s * s2
1073 + (-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);
1074}
1075
1076gslpp::complex MPll::F27(double q2)
1077{
1078 double s = q2 / Mb2;
1079 double s2 = s*s;
1080 double Ls = log(s);
1081 gslpp::complex i = gslpp::complex::i();
1082 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;
1083}
1084
1085gslpp::complex MPll::F29(double q2)
1086{
1087 double s = q2 / Mb2;
1088 double s2 = s*s;
1089 double Ls = log(s);
1090 gslpp::complex i = gslpp::complex::i();
1091 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;
1092}
1093
1094gslpp::complex MPll::F87(double q2)
1095{
1096 double s = q2 / Mb2;
1097 double s2 = s*s;
1098 return F87_0 + F87_1 * s + F87_2 * s2 + F87_3 * s * s2 - 0.888889 * log(s)*(1. + s + s2 + s * s2);
1099}
1100
1101double MPll::F89(double q2)
1102{
1103 double s = q2 / Mb2;
1104 double s2 = s*s;
1105 return F89_0 + F89_1 * s + F89_2 * s2 + F89_3 * s * s2 + 1.77778 * log(s)*(1. + s + s2 + s * s2);
1106}
1107
1108gslpp::complex MPll::Cpar(double q2)
1109{
1110 return -(C_2L_bar * F27(q2) + C_8L * F87(q2) + MM / 2. / Mb *
1111 (C_2L_bar * F29(q2) + 2. * C_1L_bar * (F19(q2) + F29(q2) / 6.) + C_8L * F89(q2)));
1112}
1113
1114gslpp::complex MPll::deltaTpar(double q2)
1115{
1116 tmpq2 = q2;
1117
1118 //old_handler = gsl_set_error_handler_off();
1119
1120 if (deltaTparpCached[q2] == 0) {
1121
1122 DTPPR = convertToGslFunction(bind(&MPll::Integrand_ReTpar_pm, &(*this), _1));
1123 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();
1124 double ReTppint = avaDTPPR;
1125
1126 DTPPR = convertToGslFunction(bind(&MPll::Integrand_ImTpar_pm, &(*this), _1));
1127 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();
1128 double ImTppint = avaDTPPR;
1129
1130 cacheDeltaTparp[q2] = (ReTppint + gslpp::complex::i() * ImTppint);
1131 deltaTparpCached[q2] = 1;
1132 }
1133
1134 //gsl_set_error_handler(old_handler);
1135
1136 return deltaT_0 * Cpar(q2) + deltaT_1par * cacheDeltaTparp[q2] / f_plus(q2);
1137}
1138
1139double MPll::reDC9fit(double* x, double* p)
1140{
1141 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];
1142}
1143
1144double MPll::imDC9fit(double* x, double* p)
1145{
1146 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];
1147
1148 //double thr = 4.*Mc2;
1149}
1150
1151void MPll::fit_DeltaC9_mumu()
1152{
1153 int dim = 0;
1154 for (double i = 0.1; i < MPllSWITCH; i += 0.4) {
1155 double q2tmp = i;
1156 myq2.push_back(q2tmp);
1157 ReDeltaC9.push_back((deltaTpar(q2tmp)).real());
1158 ImDeltaC9.push_back((deltaTpar(q2tmp)).imag());
1159 dim++;
1160 }
1161 for (double i = MPllSWITCH; i < 8.2; i += 0.4) {
1162 double q2tmp = i;
1163 myq2.push_back(q2tmp);
1164 ReDeltaC9.push_back(q2tmp * (deltaTpar(q2tmp)).real());
1165 ImDeltaC9.push_back(q2tmp * (deltaTpar(q2tmp)).imag());
1166 dim++;
1167 }
1168
1169 gr1 = TGraph(dim, myq2.data(), ReDeltaC9.data());
1170 gr2 = TGraph(dim, myq2.data(), ImDeltaC9.data());
1171
1172 reffit = TF1("reffit", this, &MPll::reDC9fit, 0, 8.1, 7);
1173 imffit = TF1("imffit", this, &MPll::imDC9fit, 0, 8.1, 8);
1174
1175 refres = gr1.Fit(&reffit, "SQN0+rob=0.99");
1176 imfres = gr2.Fit(&imffit, "SQN0+rob=0.99");
1177
1178 ReDeltaC9.clear();
1179 ImDeltaC9.clear();
1180 myq2.clear();
1181}
1182
1183gslpp::complex MPll::fDeltaC9(double q2)
1184{
1185 if (q2 < MPllSWITCH) return (reDC9fit(&q2, const_cast<double *> (refres->GetParams()))
1186 + gslpp::complex::i() * imDC9fit(&q2, const_cast<double *> (imfres->GetParams())));
1187 else return (reDC9fit(&q2, const_cast<double *> (refres->GetParams()))
1188 + gslpp::complex::i() * imDC9fit(&q2, const_cast<double *> (imfres->GetParams()))) / q2;
1189
1190}
1191
1192gslpp::complex MPll::DeltaC9(double q2)
1193{
1194 return deltaTpar(q2);
1195}
1196
1197gslpp::complex MPll::deltaC7_QCDF(double q2, bool spline)
1198{
1199#if SPLINE
1200 if (spline) return gsl_spline_eval(spline_Re_deltaC7_QCDF, q2, acc_Re_deltaC7_QCDF);
1201#endif
1202
1203 double muh = mu_b / mb_pole;
1204 double z = mc_pole * mc_pole / mb_pole / mb_pole;
1205 double sh = q2 / mb_pole / mb_pole;
1206 double sh2 = sh*sh;
1207
1208 //gslpp::complex A_Sdl = A_Seidel(q2, mb_pole*mb_pole); /* hep-ph/0403185v2.*/
1209 //gslpp::complex Fu_17 = -A_Sdl; /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1210 //gslpp::complex Fu_27 = 6. * A_Sdl; /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1211 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*/
1212 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*/
1213 gslpp::complex F_87 = F87_0 + F87_1 * sh + F87_2 * sh2 + F87_3 * sh * sh2 - 8. / 9. * log(sh) * (sh + sh2 + sh * sh2);
1214
1215 gslpp::complex delta = C_1 * F_17 + C_2 * F_27;
1216 gslpp::complex delta_t = C_8 * F_87 + delta;
1217 //gslpp::complex delta_u = delta + C_1 * Fu_17 + C_2 * Fu_27;
1218
1219 return -alpha_s_mub / (4. * M_PI) * (delta_t /*- lambda_u / lambda_t * delta_u */);
1220}
1221
1222gslpp::complex MPll::deltaC9_QCDF(double q2, bool spline)
1223{
1224#if SPLINE
1225 if (spline) return gsl_spline_eval(spline_Re_deltaC9_QCDF, q2, acc_Re_deltaC9_QCDF);
1226#endif
1227 double muh = mu_b / mb_pole;
1228 double z = mc_pole * mc_pole / mb_pole / mb_pole;
1229 double sh = q2 / mb_pole / mb_pole;
1230 double sh2 = sh*sh;
1231
1232 //gslpp::complex B_Sdl = B_Seidel(q2, mb_pole*mb_pole); /* hep-ph/0403185v2.*/
1233 //gslpp::complex C_Sdl = C_Seidel(q2); /* hep-ph/0403185v2.*/
1234 //gslpp::complex Fu_19 = -(B_Sdl + 4. * C_Sdl); /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1235 //gslpp::complex Fu_29 = -(-6. * B_Sdl + 3. * C_Sdl); /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1236 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*/
1237 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*/
1238 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));
1239
1240 gslpp::complex delta = C_1 * F_19 + C_2 * F_29;
1241 gslpp::complex delta_t = C_8 * F_89 + delta;
1242 //gslpp::complex delta_u = delta + C_1 * Fu_19 + C_2 * Fu_29;
1243
1244 return -alpha_s_mub / (4. * M_PI) * (delta_t /*- lambda_u / lambda_t * delta_u*/);
1245}
1246
1247void MPll::spline_QCDF_func()
1248{
1249 int dim_DC = GSL_INTERP_DIM_DC;
1250 double min = 0.001;
1251 double interval_DC = (9.9 - min) / ((double) dim_DC);
1252 double q2_spline_DC[dim_DC];
1253 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];
1254
1255 for (int i = 0; i < dim_DC; i++) {
1256 q2_spline_DC[i] = min + (double) i*interval_DC;
1257 fq2_Re_deltaC7_QCDF[i] = deltaC7_QCDF(q2_spline_DC[i], false).real();
1258 fq2_Im_deltaC7_QCDF[i] = deltaC7_QCDF(q2_spline_DC[i], false).imag();
1259 fq2_Re_deltaC9_QCDF[i] = deltaC9_QCDF(q2_spline_DC[i], false).real();
1260 fq2_Im_deltaC9_QCDF[i] = deltaC9_QCDF(q2_spline_DC[i], false).imag();
1261
1262 }
1263
1264 gsl_spline_init(spline_Re_deltaC7_QCDF, q2_spline_DC, fq2_Re_deltaC7_QCDF, dim_DC);
1265 gsl_spline_init(spline_Im_deltaC7_QCDF, q2_spline_DC, fq2_Im_deltaC7_QCDF, dim_DC);
1266 gsl_spline_init(spline_Re_deltaC9_QCDF, q2_spline_DC, fq2_Re_deltaC9_QCDF, dim_DC);
1267 gsl_spline_init(spline_Im_deltaC9_QCDF, q2_spline_DC, fq2_Im_deltaC9_QCDF, dim_DC);
1268
1269
1270}
1271
1272/*******************************************************************************
1273 * Helicity amplitudes *
1274 * ****************************************************************************/
1275gslpp::complex MPll::H_c(double q2, double mu2)
1276{
1277 double x = fourMc2 / q2;
1278 gslpp::complex par;
1279
1280 if (x > 1.) par = sqrt(x - 1.) * atan(1. / sqrt(x - 1.));
1281 else par = sqrt(1. - x) * (log((1. + sqrt(1. - x)) / sqrt(x)) - ihalfMPI);
1282 return -fournineth * (log(Mc2 / mu2) - twothird - x) - fournineth * (2. + x) * par;
1283}
1284
1285gslpp::complex MPll::H_b(double q2, double mu2)
1286{
1287 double x = fourMb2 / q2;
1288 gslpp::complex par;
1289
1290 if (x > 1.) par = sqrt(x - 1.) * atan(1. / sqrt(x - 1.));
1291 else par = sqrt(1. - x) * (log((1. + sqrt(1. - x)) / sqrt(x)) - ihalfMPI);
1292
1293 return -fournineth * (log(Mb2 / mu2) - twothird - x) - fournineth * (2. + x) * par;
1294}
1295
1296gslpp::complex MPll::H_0(double q2)
1297{
1298 return (H_0_pre - fournineth * log(q2 / mu_b2));
1299}
1300
1301gslpp::complex MPll::Y(double q2)
1302{
1303 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;
1304}
1305
1306gslpp::complex MPll::funct_g(double q2)
1307{
1308 if (q2 < 4. * Mc * Mc)
1309 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.)));
1310 else
1311 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));
1312}
1313
1314gslpp::complex MPll::DeltaC9_KD(double q2)
1315{
1316 return ((Delta_C9 + r_1 * q2 / mJ2) / (1. - r_2 * q2 / mJ2) - (3. * (-0.267) + 1.117) * funct_g(q2))*exp_Phase;
1317 /* C_1 = -0.267 and C_2 = 1.117 in KMPW */
1318}
1319
1320gslpp::complex MPll::h_lambda(double q2)
1321{
1322 if (!dispersion) {
1323// if (q2 <= 1.) return 1.3e-4/MM2 * (1. + gslpp::complex::i()) / sqrt(2.);
1324// else {
1325// 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.);
1326// h_2 = 1.3e-4/MM2 * (1. + gslpp::complex::i()) / sqrt(2.) - h_1 * V_L(1.)/sixteenM_PI2/MM2;
1327 //else return 4.9e-7/MM2 + h_1 * (q2 * V_L(q2) - T_L(q2). * V_L(1)/T_L(1.)) / MM2 / sixteenM_PI2;
1328
1329 //return (twoMboMM * h_0 * T_L(q2) + h_1 * q2 / MM2 * V_L(q2)) / sixteenM_PI2 + h_2 * q2 * sqrt(q2);
1330 return (h_1 + h_2 * q2 ) * sqrt(q2) * sqrt(lambda(q2)) / MM3;
1331// }
1332 } else {
1333 return -q2 / (MM2 * sixteenM_PI2) * V_L(q2) * DeltaC9_KD(q2);
1334 }
1335}
1336
1337gslpp::complex MPll::H_V(double q2)
1338{
1339 return -((C_9 + deltaC9_QCDF(q2, SPLINE) + Y(q2) /*+ fDeltaC9(q2)*/ - etaP * pow(-1, angmomP) * C_9p) * V_L(q2)
1340 + MM2 / q2 * (twoMboMM * (C_7 + deltaC7_QCDF(q2, SPLINE) - etaP * pow(-1, angmomP) * C_7p) * T_L(q2)
1341 - sixteenM_PI2 * h_lambda(q2)));
1342}
1343
1344gslpp::complex MPll::H_A(double q2)
1345{
1346 return (-C_10 + etaP * pow(-1, angmomP) * C_10p) *V_L(q2);
1347}
1348
1349gslpp::complex MPll::H_S(double q2)
1350{
1351 return MboMW * (C_S - etaP * pow(-1, angmomP) * C_Sp) * S_L(q2);
1352}
1353
1354gslpp::complex MPll::H_P(double q2)
1355{
1356 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);
1357}
1358
1359gslpp::complex MPll::H_nunu(double q2, QCD::lepton lep)
1360{
1361 if (lep == QCD::NEUTRINO_1)
1362 return -(C_L_nunu_e - etaP * pow(-1, angmomP) * C_R_nunu_e) * V_L(q2);
1363 else if (lep == QCD::NEUTRINO_2)
1364 return -(C_L_nunu_mu - etaP * pow(-1, angmomP) * C_R_nunu_mu) * V_L(q2);
1365 else if (lep == QCD::NEUTRINO_3)
1366 return -(C_L_nunu_tau - etaP * pow(-1, angmomP) * C_R_nunu_tau) * V_L(q2);
1367 else throw std::runtime_error("MPll::H_nunu: lepton not supported");
1368}
1369
1370/*******************************************************************************
1371 * Angular coefficients *
1372 * ****************************************************************************/
1373double MPll::k2(double q2)
1374{
1375 return (MM4 + q2 * q2 + MP4 - twoMP2 * q2 - twoMM2 * (q2 + MP2)) / fourMM2;
1376}
1377
1378double MPll::beta(double q2)
1379{
1380 return sqrt(1. - 4. * Mlep2 / q2);
1381}
1382
1383double MPll::beta2(double q2)
1384{
1385 return 1. - 4. * Mlep2 / q2;
1386}
1387
1388double MPll::lambda(double q2)
1389{
1390 return 4. * MM2 * k2(q2);
1391}
1392
1393double MPll::F(double q2)
1394{
1395 return sqrt(lambda(q2)) * beta(q2) * q2 / (ninetysixM_PI3MM3);
1396}
1397
1398double MPll::I_1c(double q2)
1399{
1400 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();
1401
1402 else return F(q2)*((H_V(q2).abs2() + H_A(q2).abs2()) / 2. + H_P(q2).abs2() + 2. * Mlep2 / q2 * (H_V(q2).abs2()
1403 - H_A(q2).abs2()) + beta2(q2) * H_S(q2).abs2());
1404}
1405
1406double MPll::I_2c(double q2)
1407{
1408 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();
1409
1410 else return -F(q2) * beta2(q2) / 2. * (H_V(q2).abs2() + H_A(q2).abs2());
1411}
1412
1413double MPll::I_6c(double q2)
1414{
1415 return 4. * F(q2) * beta(q2) * Mlep / sqrt(q2)*(H_S(q2).conjugate() * H_V(q2)).real();
1416}
1417
1418double MPll::Delta(int i, double q2)
1419{
1420 return 0; /* FIX CPV */
1421 //return (I(i, q2,0) - I(i, q2,1))/2;
1422}
1423
1424double MPll::integrateSigma(int i, double q_min, double q_max)
1425{
1426 updateParameters();
1427
1428 std::pair<double, double > qbin = std::make_pair(q_min, q_max);
1429
1430 old_handler = gsl_set_error_handler_off();
1431
1432 switch (i) {
1433 case 0:
1434 if (sigma0Cached[qbin] == 0) {
1435 FS = convertToGslFunction(bind(&MPll::getSigma1c, &(*this), _1));
1436 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();
1437 cacheSigma0[qbin] = NN*avaSigma;
1438 sigma0Cached[qbin] = 1;
1439 }
1440 return cacheSigma0[qbin];
1441 break;
1442 case 2:
1443 if (sigma2Cached[qbin] == 0) {
1444 FS = convertToGslFunction(bind(&MPll::getSigma2c, &(*this), _1));
1445 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();
1446 cacheSigma2[qbin] = NN*avaSigma;
1447 sigma2Cached[qbin] = 1;
1448 }
1449 return cacheSigma2[qbin];
1450 break;
1451 default:
1452 std::stringstream out;
1453 out << i;
1454 throw std::runtime_error("MPll::integrateSigma: index " + out.str() + " not implemented");
1455 }
1456
1457 gsl_set_error_handler(old_handler);
1458
1459}
1460
1461double MPll::getSigma(int i, double q_2)
1462{
1463 updateParameters();
1464
1465 switch (i) {
1466 case 0:
1467 return getSigma1c(q_2);
1468 break;
1469 case 2:
1470 return getSigma2c(q_2);
1471 break;
1472 case 8:
1473 return getSigma6c(q_2);
1474 break;
1475 default:
1476 std::stringstream out;
1477 out << i;
1478 throw std::runtime_error("MPll::getSigma: index " + out.str() + " not implemented");
1479 }
1480}
1481
1482double MPll::integrateDelta(int i, double q_min, double q_max)
1483{
1484 updateParameters();
1485
1486 std::pair<double, double > qbin = std::make_pair(q_min, q_max);
1487
1488 old_handler = gsl_set_error_handler_off();
1489
1490 switch (i) {
1491 case 0:
1492 if (delta0Cached[qbin] == 0) {
1493 FD = convertToGslFunction(bind(&MPll::getDelta1c, &(*this), _1));
1494 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();
1495 cacheDelta0[qbin] = NN*avaDelta;
1496 delta0Cached[qbin] = 1;
1497 }
1498 return cacheDelta0[qbin];
1499 break;
1500 case 2:
1501 if (delta2Cached[qbin] == 0) {
1502 FD = convertToGslFunction(bind(&MPll::getDelta2c, &(*this), _1));
1503 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();
1504 cacheDelta2[qbin] = NN*avaDelta;
1505 delta2Cached[qbin] = 1;
1506 }
1507 return cacheDelta2[qbin];
1508 break;
1509 default:
1510 std::stringstream out;
1511 out << i;
1512 throw std::runtime_error("MPll::integrateDelta: index " + out.str() + " not implemented");
1513 }
1514
1515 gsl_set_error_handler(old_handler);
1516
1517}
1518
1519double MPll::integrateSigmaTree(double q_min, double q_max)
1520{
1521 if (lep != QCD::NEUTRINO_1 or meson != QCD::B_P or !NeutrinoTree_flag) return 0.;
1522
1523 updateParameters();
1524
1525 //phase space limit where tree-level contribution is relevant (0908.1174)
1526 double q_cut = (mtau2 - MP2) * (MM2 - mtau2) / mtau2;
1527 if (q_max >= q_cut) {
1528 if (q_min == 0.) return getintegratedSigmaTree();
1529 q_max = q_cut;
1530 }
1531
1532 double prefactor = mySM.getMesons(meson).getLifetime() / HCUT * GF4 * VusVub_abs2 * fP2 * fM2 / (64. * M_PI2 * MM3 * Gammatau) * mtau2 * mtau;
1533
1534 std::pair<double, double > qbin = std::make_pair(q_min, q_max);
1535
1536 old_handler = gsl_set_error_handler_off();
1537
1538 if (sigmaTreeCached[qbin] == 0) {
1539 FD = convertToGslFunction(bind(&MPll::SigmaTree, &(*this), _1));
1540 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();
1541 cacheSigmaTree[qbin] = avaSigmaTree;
1542 sigmaTreeCached[qbin] = 1;
1543 }
1544 return prefactor * cacheSigmaTree[qbin];
1545
1546 gsl_set_error_handler(old_handler);
1547}
1548
1549double MPll::SigmaTree(double q2)
1550{
1551 return MM2 * (mtau2 - MP2) - mtau2 * (mtau2 + q2 - MP2);
1552}
1553
1554double MPll::getintegratedSigmaTree()
1555{
1556 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);
1557}
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:1057
MPll::getSigma1c
double getSigma1c(double q2)
The CP average .
Definition: MPll.h:1014
MPll::FixedWCbtos
bool FixedWCbtos
Definition: MPll.h:323
MPll::SigmaTree
double SigmaTree(double q2)
Definition: MPll.cpp:1549
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:1344
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:1320
MPll::F89
double F89(double q2)
The correction from .
Definition: MPll.cpp:1101
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:1222
MPll::DeltaC9_KD
gslpp::complex DeltaC9_KD(double q2)
Definition: MPll.cpp:1314
MPll::H_P
gslpp::complex H_P(double q2)
The helicity amplitude .
Definition: MPll.cpp:1354
MPll::pseudoscalar
QCD::meson pseudoscalar
Definition: MPll.h:318
MPll::H_S
gslpp::complex H_S(double q2)
The helicity amplitude .
Definition: MPll.cpp:1349
MPll::F87
gslpp::complex F87(double q2)
The correction from .
Definition: MPll.cpp:1094
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:1041
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:972
MPll::Integrand_ReTparminus
double Integrand_ReTparminus(double up)
The real part of the integral involving at fixed , according to .
Definition: MPll.cpp:1030
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:1151
MPll::Tparminus
gslpp::complex Tparminus(double u, double q2)
The function from .
Definition: MPll.cpp:1007
MPll::imDC9fit
double imDC9fit(double *x, double *p)
The fit function for the imaginary part of the QCDF correction .
Definition: MPll.cpp:1144
MPll::getDelta2c
double getDelta2c(double q2)
The CP asymmetry .
Definition: MPll.h:1054
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:1062
MPll::integrateSigma
double integrateSigma(int i, double q_min, double q_max)
The integral of from to .
Definition: MPll.cpp:1424
MPll::Cpar
gslpp::complex Cpar(double q2)
The correction from .
Definition: MPll.cpp:1108
MPll::reDC9fit
double reDC9fit(double *x, double *p)
The fit function for the real part of the QCDF correction .
Definition: MPll.cpp:1139
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:953
MPll::H_nunu
gslpp::complex H_nunu(double q2, QCD::lepton lep)
The helicity amplitude for the invisible decay.
Definition: MPll.cpp:1359
MPll::getSigma
double getSigma(int i, double q_2)
The value of from to .
Definition: MPll.cpp:1461
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:1306
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:1052
MPll::Integrand_ReTparplus
double Integrand_ReTparplus(double up)
The real part of the integral involving at fixed , according to .
Definition: MPll.cpp:1014
MPll::T_L
gslpp::complex T_L(double q2)
The helicity form factor .
Definition: MPll.cpp:958
MPll::alpha_s_mub
double alpha_s_mub
Definition: MPll.h:344
MPll::getSigma2c
double getSigma2c(double q2)
The CP average .
Definition: MPll.h:1024
MPll::Tparplus
gslpp::complex Tparplus(double u, double q2)
The function from .
Definition: MPll.cpp:994
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:1519
MPll::getintegratedSigmaTree
double getintegratedSigmaTree()
The integral of from 0 to .
Definition: MPll.cpp:1554
MPll::Ms
double Ms
Definition: MPll.h:341
MPll::getDelta1c
double getDelta1c(double q2)
The CP asymmetry .
Definition: MPll.h:1044
MPll::MP
double MP
Definition: MPll.h:334
MPll::H_V
gslpp::complex H_V(double q2)
The helicity amplitude .
Definition: MPll.cpp:1337
MPll::DeltaC9
gslpp::complex DeltaC9(double q2)
The total QCDF correction computed integrating over .
Definition: MPll.cpp:1192
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:1197
MPll::dispersion
bool dispersion
Definition: MPll.h:322
MPll::getSigma6c
double getSigma6c(double q2)
The CP average .
Definition: MPll.h:1034
MPll::mb_pole
double mb_pole
Definition: MPll.h:339
MPll::F27
gslpp::complex F27(double q2)
The correction from .
Definition: MPll.cpp:1076
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:1183
MPll::S_L
double S_L(double q2)
The helicity form factor .
Definition: MPll.cpp:963
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:1114
MPll::integrateDelta
double integrateDelta(int i, double q_min, double q_max)
The integral of from to .
Definition: MPll.cpp:1482
MPll::F29
gslpp::complex F29(double q2)
The correction from .
Definition: MPll.cpp:1085
MPll::spline_QCDF_func
void spline_QCDF_func()
Definition: MPll.cpp:1247
MPll::Integrand_ImTparplus
double Integrand_ImTparplus(double up)
The imaginary part of the integral involving at fixed , according to .
Definition: MPll.cpp:1022
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.