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