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