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.hpp>
16 #include <limits>
17 #include <TFitResult.h>
18 
19 MPll::MPll(const StandardModel& SM_i, QCD::meson meson_i, QCD::meson pseudoscalar_i, QCD::lepton lep_i)
20 : mySM(SM_i), myF_1(new F_1()), myF_2(new F_2()),
21 #if LATTICE
22 fplus_lat_cache(3, 0.),
23 fT_lat_cache(3, 0.),
24 f0_lat_cache(3, 0.),
25 #else
26 fplus_cache(2, 0.),
27 fT_cache(2, 0.),
28 f0_cache(2, 0.),
29 #endif
30 k2_cache(2, 0.),
31 SL_cache(2, 0.),
32 N_cache(3, 0.),
33 Ycache(2, 0.),
34 H_V0cache(2, 0.),
35 H_Scache(2, 0.),
36 H_P_cache(4, 0.),
37 T_cache(5, 0.)
38 {
39  lep = lep_i;
40  meson = meson_i;
41  pseudoscalar = pseudoscalar_i;
42  dispersion = true;
43  FixedWCbtos = false;
44  mJ2 = 3.096 * 3.096;
45 
46  I0_updated = 0;
47  I2_updated = 0;
48  I8_updated = 0;
49 
50  VL_updated = 0;
51  TL_updated = 0;
52  SL_updated = 0;
53 
54  deltaTparpupdated = 0;
55  deltaTparmupdated = 0;
56 
57  w_sigma = gsl_integration_cquad_workspace_alloc(100);
58  w_delta = gsl_integration_cquad_workspace_alloc(100);
59  w_DTPPR = gsl_integration_cquad_workspace_alloc(100);
60 
61  acc_Re_deltaC7_QCDF = gsl_interp_accel_alloc();
62  acc_Im_deltaC7_QCDF = gsl_interp_accel_alloc();
63  acc_Re_deltaC9_QCDF = gsl_interp_accel_alloc();
64  acc_Im_deltaC9_QCDF = gsl_interp_accel_alloc();
65 
66  spline_Re_deltaC7_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
67  spline_Im_deltaC7_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
68  spline_Re_deltaC9_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
69  spline_Im_deltaC9_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
70 }
71 
72 MPll::~MPll()
73 {
74 }
75 
76 std::vector<std::string> MPll::initializeMPllParameters()
77 {
78  dispersion = mySM.getFlavour().getFlagUseDispersionRelation();
79  FixedWCbtos = mySM.getFlavour().getFlagFixedWCbtos();
80 
81 #if NFPOLARBASIS_MPLL
82  if (pseudoscalar == StandardModel::K_P || pseudoscalar == StandardModel::K_0) mpllParameters = make_vector<std::string>()
83 #if LATTICE
84  << "b_0_fplus" << "b_1_fplus" << "b_2_fplus" << "m_fit_fplus_lat"
85  << "b_0_fT" << "b_1_fT" << "b_2_fT" << "m_fit_fT_lat"
86  << "b_0_f0" << "b_1_f0" << "b_2_f0" << "m_fit_f0_lat"
87 #else
88  << "r_1_fplus" << "r_2_fplus" << "m_fit2_fplus" << "r_1_fT" << "r_2_fT" << "m_fit2_fT" << "r_2_f0" << "m_fit2_f0"
89 #endif
90  << "absh_0_MP" << "argh_0_MP" << "absh_1_MP" << "argh_1_MP" << "absh_2_MP" << "argh_2_MP";
91 #else
92  if (pseudoscalar == StandardModel::K_P || pseudoscalar == StandardModel::K_0) mpllParameters = make_vector<std::string>()
93 #if LATTICE
94  << "b_0_fplus" << "b_1_fplus" << "b_2_fplus" << "m_fit_fplus_lat"
95  << "b_0_fT" << "b_1_fT" << "b_2_fT" << "m_fit_fT_lat"
96  << "b_0_f0" << "b_1_f0" << "b_2_f0" << "m_fit_f0_lat"
97 #else
98  << "r_1_fplus" << "r_2_fplus" << "m_fit2_fplus" << "r_1_fT" << "r_2_fT" << "m_fit2_fT" << "r_2_f0" << "m_fit2_f0"
99 #endif
100  << "reh_0_MP" << "imh_0_MP" << "reh_1_MP" << "imh_1_MP" << "reh_2_MP" << "imh_2_MP";
101 #endif
102  else {
103  std::stringstream out;
104  out << pseudoscalar;
105  throw std::runtime_error("MPll: pseudoscalar " + out.str() + " not implemented");
106  }
107 
108  if (dispersion) {
109  mpllParameters.clear();
110  if (pseudoscalar == StandardModel::K_P || pseudoscalar == StandardModel::K_0) mpllParameters = make_vector<std::string>()
111 #if LATTICE
112  << "b_0_fplus" << "b_1_fplus" << "b_2_fplus" << "m_fit_fplus_lat"
113  << "b_0_fT" << "b_1_fT" << "b_2_fT" << "m_fit_fT_lat"
114  << "b_0_f0" << "b_1_f0" << "b_2_f0" << "m_fit_f0_lat"
115 #else
116  << "r_1_fplus" << "r_2_fplus" << "m_fit2_fplus" << "r_1_fT" << "r_2_fT" << "m_fit2_fT" << "r_2_f0" << "m_fit2_f0"
117 #endif
118  << "r1_BK" << "r2_BK" << "deltaC9_BK" << "phDC9_BK";
119  }
120 
121  if (FixedWCbtos) mpllParameters.insert(mpllParameters.end(), { "C7_SM", "C9_SM", "C10_SM" });
122  mySM.initializeMeson(meson);
123  mySM.initializeMeson(pseudoscalar);
124  return mpllParameters;
125 }
126 
127 void MPll::updateParameters()
128 {
129  if (!mySM.getFlavour().getUpdateFlag(meson, pseudoscalar, lep)) return;
130 
131 
132 
133  GF = mySM.getGF();
134  ale = mySM.getAle();
135  Mlep = mySM.getLeptons(lep).getMass();
136  MM = mySM.getMesons(meson).getMass();
137  MP = mySM.getMesons(pseudoscalar).getMass();
138  Mb = mySM.getQuarks(QCD::BOTTOM).getMass(); // add the PS b mass
139  Mc = mySM.getQuarks(QCD::CHARM).getMass();
140  mb_pole = mySM.Mbar2Mp(Mb); /* Conversion to pole mass*/
141  mc_pole = mySM.Mbar2Mp(Mc,FULLNLO); /* Conversion to pole mass*/
142  Ms = mySM.getQuarks(QCD::STRANGE).getMass();
143  MW = mySM.Mw();
144  lambda_t = mySM.getCKM().computelamt_s();
145  mu_b = mySM.getMub();
146  mu_h = sqrt(mu_b * .5); // From Beneke Neubert
147  width = mySM.getMesons(meson).computeWidth();
148  alpha_s_mub = mySM.Als(mu_b);
149 
150  switch (pseudoscalar) {
151  case StandardModel::K_P:
152  case StandardModel::K_0:
153 #if LATTICE
154  b_0_fplus = mySM.getOptionalParameter("b_0_fplus");
155  b_1_fplus = mySM.getOptionalParameter("b_1_fplus");
156  b_2_fplus = mySM.getOptionalParameter("b_2_fplus");
157  m_fit2_fplus_lat = mySM.getOptionalParameter("m_fit_fplus_lat") * mySM.getOptionalParameter("m_fit_fplus_lat");
158  b_0_fT = mySM.getOptionalParameter("b_0_fT");
159  b_1_fT = mySM.getOptionalParameter("b_1_fT");
160  b_2_fT = mySM.getOptionalParameter("b_2_fT");
161  m_fit2_fT_lat = mySM.getOptionalParameter("m_fit_fT_lat") * mySM.getOptionalParameter("m_fit_fT_lat");
162  b_0_f0 = mySM.getOptionalParameter("b_0_f0");
163  b_1_f0 = mySM.getOptionalParameter("b_1_f0");
164  b_2_f0 = mySM.getOptionalParameter("b_2_f0");
165  m_fit2_f0_lat = mySM.getOptionalParameter("m_fit_f0_lat") * mySM.getOptionalParameter("m_fit_f0_lat");
166 #else
167  r_1_fplus = mySM.getOptionalParameter("r_1_fplus");
168  r_2_fplus = mySM.getOptionalParameter("r_2_fplus");
169  m_fit2_fplus = mySM.getOptionalParameter("m_fit2_fplus");
170  r_1_fT = mySM.getOptionalParameter("r_1_fT");
171  r_2_fT = mySM.getOptionalParameter("r_2_fT");
172  m_fit2_fT = mySM.getOptionalParameter("m_fit2_fT");
173  r_2_f0 = mySM.getOptionalParameter("r_2_f0");
174  m_fit2_f0 = mySM.getOptionalParameter("m_fit2_f0");
175 #endif
176 
177  if (pseudoscalar == StandardModel::K_P) spectator_charge = mySM.getQuarks(QCD::UP).getCharge();
178  else spectator_charge = mySM.getQuarks(QCD::DOWN).getCharge();
179 
180  etaP = -1;
181  angmomP = 0.;
182 
183  break;
184  default:
185  std::stringstream out;
186  out << pseudoscalar;
187  throw std::runtime_error("MPll: pseudoscalar " + out.str() + " not implemented");
188  }
189 
190  if (!dispersion) {
191 #if NFPOLARBASIS_MPLL
192  h_0 = gslpp::complex(mySM.getOptionalParameter("absh_0_MP"), mySM.getOptionalParameter("argh_0_MP"), true);
193  h_1 = gslpp::complex(mySM.getOptionalParameter("absh_1_MP"), mySM.getOptionalParameter("argh_1_MP"), true);
194  h_2 = gslpp::complex(mySM.getOptionalParameter("absh_2_MP"), mySM.getOptionalParameter("argh_2_MP"), true);
195 
196  r_1 = 0.;
197  r_2 = 0.;
198  Delta_C9 = 0.;
199  exp_Phase = 0.;
200 #else
201  h_0 = gslpp::complex(mySM.getOptionalParameter("reh_0_MP"), mySM.getOptionalParameter("imh_0_MP"), false);
202  h_1 = gslpp::complex(mySM.getOptionalParameter("reh_1_MP"), mySM.getOptionalParameter("imh_1_MP"), false);
203  h_2 = gslpp::complex(mySM.getOptionalParameter("reh_2_MP"), mySM.getOptionalParameter("imh_2_MP"), false);
204 
205  r_1 = 0.;
206  r_2 = 0.;
207  Delta_C9 = 0.;
208  exp_Phase = 0.;
209 #endif
210  } else {
211  h_0 = 0.;
212  h_1 = 0.;
213  h_2 = 0.;
214 
215  r_1 = mySM.getOptionalParameter("r1_BK");
216  r_2 = mySM.getOptionalParameter("r2_BK");
217  Delta_C9 = mySM.getOptionalParameter("deltaC9_BK");
218  exp_Phase = exp(gslpp::complex::i() * mySM.getOptionalParameter("phDC9_BK"));
219  }
220 
221  allcoeff = mySM.getFlavour().ComputeCoeffBMll(mu_b, lep); //check the mass scale, scheme fixed to NDR
222  allcoeffprime = mySM.getFlavour().ComputeCoeffprimeBMll(mu_b, lep); //check the mass scale, scheme fixed to NDR
223 
224  C_1 = (*(allcoeff[LO]))(0) + (*(allcoeff[NLO]))(0);
225  C_1L_bar = (*(allcoeff[LO]))(0) / 2.;
226  C_2 = ((*(allcoeff[LO]))(1) + (*(allcoeff[NLO]))(1));
227  C_2L_bar = (*(allcoeff[LO]))(1) - (*(allcoeff[LO]))(0) / 6.;
228  C_3 = ((*(allcoeff[LO]))(2) + (*(allcoeff[NLO]))(2));
229  C_4 = (*(allcoeff[LO]))(3) + (*(allcoeff[NLO]))(3);
230  C_5 = ((*(allcoeff[LO]))(4) + (*(allcoeff[NLO]))(4));
231  C_6 = ((*(allcoeff[LO]))(5) + (*(allcoeff[NLO]))(5));
232  C_8 = ((*(allcoeff[LO]))(7) + (*(allcoeff[NLO]))(7));
233  C_8L = (*(allcoeff[LO]))(7);
234  C_S = MW / Mb * ((*(allcoeff[LO]))(10) + (*(allcoeff[NLO]))(10));
235  C_P = MW / Mb * ((*(allcoeff[LO]))(11) + (*(allcoeff[NLO]))(11));
236  C_9p = (*(allcoeffprime[LO]))(8) + (*(allcoeffprime[NLO]))(8);
237  C_10p = (*(allcoeffprime[LO]))(9) + (*(allcoeffprime[NLO]))(9);
238  C_Sp = MW / Mb * ((*(allcoeffprime[LO]))(10) + (*(allcoeffprime[NLO]))(10));
239  C_Pp = MW / Mb * ((*(allcoeffprime[LO]))(11) + (*(allcoeffprime[NLO]))(11));
240 
241  if (FixedWCbtos) {
242  allcoeff_noSM = mySM.getFlavour().ComputeCoeffBMll(mu_b, lep, true); //check the mass scale, scheme fixed to NDR
243  C_7 = mySM.getOptionalParameter("C7_SM") + ((*(allcoeff_noSM[LO]))(6) + (*(allcoeff_noSM[NLO]))(6));
244  C_9 = mySM.getOptionalParameter("C9_SM") + ((*(allcoeff_noSM[LO]))(8) + (*(allcoeff_noSM[NLO]))(8));
245  C_10 = mySM.getOptionalParameter("C10_SM") + ((*(allcoeff_noSM[LO]))(9) + (*(allcoeff_noSM[NLO]))(9));
246  } else {
247  C_7 = ((*(allcoeff[LO]))(6) + (*(allcoeff[NLO]))(6));
248  C_9 = ((*(allcoeff[LO]))(8) + (*(allcoeff[NLO]))(8));
249  C_10 = ((*(allcoeff[LO]))(9) + (*(allcoeff[NLO]))(9));
250  }
251  C_7p = MsoMb * ((*(allcoeffprime[LO]))(6) + (*(allcoeffprime[NLO]))(6));
252  C_7p -= MsoMb * (C_7 + 1. / 3. * C_3 + 4 / 9 * C_4 + 20. / 3. * C_5 + 80. / 9. * C_6);
253 
254  allcoeffh = mySM.getFlavour().ComputeCoeffBMll(mu_h, lep); //check the mass scale, scheme fixed to NDR
255 
256  C_1Lh_bar = (*(allcoeffh[LO]))(0) / 2.;
257  C_2Lh_bar = (*(allcoeffh[LO]))(1) - (*(allcoeff[LO]))(0) / 6.;
258  C_8Lh = (*(allcoeffh[LO]))(7);
259 
260  checkCache();
261 
262  H_0_pre = 8. / 27. + 4. / 9. * gslpp::complex::i() * M_PI;
263  H_0_WC = (C_3 + 4. / 3. * C_4 + 16. * C_5 + 64. / 3. * C_6);
264  H_c_WC = (4. / 3. * C_1 + C_2 + 6. * C_3 + 60. * C_5);
265  H_b_WC = (7. * C_3 + 4. / 3. * C_4 + 76. * C_5 + 64. / 3. * C_6);
266  fournineth = 4. / 9.;
267  half = 1. / 2.;
268  twothird = 2. / 3.;
269  ihalfMPI = gslpp::complex::i() * M_PI / 2.;
270  Mc2 = Mc*Mc;
271  Mb2 = Mb*Mb;
272  logMc = log(Mc2 / mu_b2);
273  logMb = log(Mb2 / mu_b2);
274  mu_b2 = mu_b*mu_b;
275  fourMc2 = 4. * Mc2;
276  fourMb2 = 4. * Mb2;
277  Mlep2 = Mlep*Mlep;
278  NN = ((4. * GF * MM * ale * lambda_t) / (sqrt(2.)*4. * M_PI)).abs2();
279  MM2 = MM*MM;
280  MM4 = MM2*MM2;
281  MP2 = MP*MP;
282  MP4 = MP2*MP2;
283  MM2mMP2 = MM2 - MP2;
284  twoMP2 = 2. * MP2;
285  twoMM = 2. * MM;
286  twoMM2 = 2. * MM2;
287  twoMM2_MMpMP = twoMM2 * (MM + MP);
288  twoMM_MbpMs = twoMM * (Mb + Ms);
289  S_L_pre = -(MM2mMP2 / twoMM_MbpMs) * (1 + MsoMb) / (1 - MsoMb);
290  fourMM2 = 4. * MM2;
291  twoMboMM = 2 * Mb / MM;
292  sixteenM_PI2 = 16. * M_PI*M_PI;
293  ninetysixM_PI3MM3 = 96. * M_PI * M_PI * M_PI * MM * MM*MM;
294  MboMW = Mb / MW;
295  MboMM = Mb / MM;
296  MsoMb = Ms / Mb;
297  twoMlepMb = 2. * Mlep*Mb;
298  threeGegen0 = mySM.getMesons(pseudoscalar).getGegenalpha(0)*3.;
299  threeGegen1otwo = mySM.getMesons(pseudoscalar).getGegenalpha(1)*3. / 2.;
300  M_PI2osix = M_PI * M_PI / 6.;
301  twoMc2 = 2. * Mc2;
302  sixMMoMb = 6. * MM / Mb;
303  CF = 4. / 3.;
304  deltaT_0 = alpha_s_mub * MboMM / 4. / M_PI;
305  deltaT_1par = mySM.Als(mu_h) * CF / 4. * M_PI / 3. * mySM.getMesons(meson).getDecayconst() *
306  mySM.getMesons(pseudoscalar).getDecayconst() / MM * MboMM;
307 
308  F87_0 = -32. / 9. * log(mu_b / Mb) + 8. / 27. * M_PI * M_PI - 44. / 9. - 8. / 9. * gslpp::complex::i() * M_PI;
309  F87_1 = (4. / 3. * M_PI * M_PI - 40. / 3.);
310  F87_2 = (32. / 9. * M_PI * M_PI - 316. / 9.);
311  F87_3 = (200. / 27. * M_PI * M_PI - 658. / 9.);
312 
313  F89_0 = 104. / 9. - 32. / 27. * M_PI * M_PI;
314  F89_1 = 1184. / 27. - 40. / 9. * M_PI * M_PI;
315  F89_2 = (-32. / 3. * M_PI * M_PI + 14212. / 135.);
316  F89_3 = (-560. / 27. * M_PI * M_PI + 193444. / 945.);
317 
318  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();
319  F29_L1 = 32. / 81. * log(mu_b / Mb) + (0.48576 + 0.31119 * gslpp::complex::i());
320  F29_1 = (-32. / 405. + 64. / 45. / Mc2 * Mb2) * log(mu_b / Mb) + (1.9061 + 0.80843 * gslpp::complex::i());
321  F29_2 = (-8. / 945. + 16. / 105. / Mc2 * Mb2 / Mc2 * Mb2) * log(mu_b / Mb) + (-1.8286 + 2.8428 * gslpp::complex::i());
322  F29_3 = (-32. / 25515. + 64. / 2835. / Mc2 * Mb2 / Mc2 * Mb2 / Mc2 * Mb2) * log(mu_b / Mb) + (-12.113 + 8.1251 * gslpp::complex::i());
323  F29_L1_1 = (0.21951 - 0.14852 * gslpp::complex::i());
324  F29_L1_2 = (0.13015 - 0.22155 * gslpp::complex::i());
325  F29_L1_3 = (-0.079692 - 0.31214 * gslpp::complex::i());
326 
327  F27_0 = 416. / 81. * log(mu_b / Mb) + 3.8367 + 0.3531 * gslpp::complex::i();
328  F27_1 = (1.3098 + 0.60185 * gslpp::complex::i());
329  F27_2 = (0.13507 + 0.89014 * gslpp::complex::i());
330  F27_3 = (-1.0271 + 0.77168 * gslpp::complex::i());
331  F27_L1_1 = (-0.031936 - 0.10981 * gslpp::complex::i());
332  F27_L1_2 = (-0.14169 - 0.035553 * gslpp::complex::i());
333  F27_L1_3 = (-0.13592 + 0.093 * gslpp::complex::i());
334 
335  std::map<std::pair<double, double>, unsigned int >::iterator it;
336 
337  if (I0_updated == 0) for (it = sigma0Cached.begin(); it != sigma0Cached.end(); ++it) it->second = 0;
338  if (I2_updated == 0) for (it = sigma2Cached.begin(); it != sigma2Cached.end(); ++it) it->second = 0;
339 
340  if (I0_updated == 0) for (it = delta0Cached.begin(); it != delta0Cached.end(); ++it) it->second = 0;
341  if (I2_updated == 0) for (it = delta2Cached.begin(); it != delta2Cached.end(); ++it) it->second = 0;
342 
343  std::map<double, unsigned int >::iterator iti;
344  if (deltaTparpupdated == 0) for (iti = deltaTparpCached.begin(); iti != deltaTparpCached.end(); ++iti) iti->second = 0;
345  if (deltaTparmupdated == 0) for (iti = deltaTparmCached.begin(); iti != deltaTparmCached.end(); ++iti) iti->second = 0;
346 
347  if (deltaTparpupdated * deltaTparmupdated == 0) for (it = I1Cached.begin(); it != I1Cached.end(); ++it) it->second = 0;
348 
349 #if SPLINE
350  if (mySM.getFlavour().getUpdateFlag(meson, pseudoscalar, lep)) spline_QCDF_func();
351 #else
352  fit_DeltaC9_mumu();
353 #endif
354 
355  mySM.getFlavour().setUpdateFlag(meson, pseudoscalar, lep, false);
356 
357  return;
358 
359 }
360 
361 void MPll::checkCache()
362 {
363 
364  if (MM == k2_cache(0) && MP == k2_cache(1)) {
365  k2_updated = 1;
366  } else {
367  k2_updated = 0;
368  k2_cache(0) = MM;
369  k2_cache(1) = MP;
370  }
371 
372 #if LATTICE
373  if (b_0_fplus == fplus_lat_cache(0) && b_1_fplus == fplus_lat_cache(1) && b_2_fplus == fplus_lat_cache(2)) {
374  fplus_lat_updated = 1;
375  } else {
376  fplus_lat_updated = 0;
377  fplus_lat_cache(0) = b_0_fplus;
378  fplus_lat_cache(1) = b_1_fplus;
379  fplus_lat_cache(2) = b_2_fplus;
380  }
381 
382  if (b_0_fT == fT_lat_cache(0) && b_1_fT == fT_lat_cache(1) && b_2_fT == fT_lat_cache(2)) {
383  fT_lat_updated = 1;
384  } else {
385  fT_lat_updated = 0;
386  fT_lat_cache(0) = b_0_fT;
387  fT_lat_cache(1) = b_1_fT;
388  fT_lat_cache(2) = b_2_fT;
389  }
390 
391  if (b_0_f0 == f0_lat_cache(0) && b_1_f0 == f0_lat_cache(1) && b_2_f0 == f0_lat_cache(2)) {
392  f0_lat_updated = 1;
393  } else {
394  f0_lat_updated = 0;
395  f0_lat_cache(0) = b_0_f0;
396  f0_lat_cache(1) = b_1_f0;
397  f0_lat_cache(2) = b_2_f0;
398  }
399 #else
400  if (r_1_fplus == fplus_cache(0) && r_2_fplus == fplus_cache(1)) {
401  fplus_updated = 1;
402  } else {
403  fplus_updated = 0;
404  fplus_cache(0) = r_1_fplus;
405  fplus_cache(1) = r_2_fplus;
406  }
407 
408  if (r_1_fT == fT_cache(0) && r_2_fT == fT_cache(1)) {
409  fT_updated = 1;
410  } else {
411  fT_updated = 0;
412  fT_cache(0) = r_1_fT;
413  fT_cache(1) = r_2_fT;
414  }
415 
416  if (r_2_f0 == f0_cache) {
417  f0_updated = 1;
418  } else {
419  f0_updated = 0;
420  f0_cache = r_2_f0;
421  }
422 #endif
423  if (Mlep == beta_cache) {
424  beta_updated = 1;
425  } else {
426  beta_updated = 0;
427  beta_cache = Mlep;
428  }
429 
430  lambda_updated = k2_updated;
431  F_updated = lambda_updated * beta_updated;
432 
433 #if LATTICE
434  VL_updated = k2_updated * fplus_lat_updated;
435 #else
436  VL_updated = k2_updated * fplus_updated;
437 #endif
438 
439 #if LATTICE
440  TL_updated = k2_updated * fT_lat_updated;
441 #else
442  TL_updated = k2_updated * fT_updated;
443 #endif
444 
445  if (Mb == SL_cache(0) && Ms == SL_cache(1)) {
446 #if LATTICE
447  SL_updated = k2_updated * f0_lat_updated;
448 #else
449  SL_updated = k2_updated * f0_updated;
450 #endif
451  } else {
452  SL_updated = 0;
453  SL_cache(0) = Mb;
454  SL_cache(1) = Ms;
455  }
456 
457  if (GF == N_cache(0) && ale == N_cache(1) && MM == N_cache(2) && lambda_t == Nc_cache) {
458  N_updated = 1;
459  } else {
460  N_updated = 0;
461  N_cache(0) = GF;
462  N_cache(1) = ale;
463  N_cache(2) = MM;
464  Nc_cache = lambda_t;
465  }
466 
467  if (C_1 == C_1_cache) {
468  C_1_updated = 1;
469  } else {
470  C_1_updated = 0;
471  C_1_cache = C_1;
472  }
473 
474  if (C_2 == C_2_cache) {
475  C_2_updated = 1;
476  } else {
477  C_2_updated = 0;
478  C_2_cache = C_2;
479  }
480 
481  if (C_3 == C_3_cache) {
482  C_3_updated = 1;
483  } else {
484  C_3_updated = 0;
485  C_3_cache = C_3;
486  }
487 
488  if (C_4 == C_4_cache) {
489  C_4_updated = 1;
490  } else {
491  C_4_updated = 0;
492  C_4_cache = C_4;
493  }
494 
495  if (C_5 == C_5_cache) {
496  C_5_updated = 1;
497  } else {
498  C_5_updated = 0;
499  C_5_cache = C_5;
500  }
501 
502  if (C_6 == C_6_cache) {
503  C_6_updated = 1;
504  } else {
505  C_6_updated = 0;
506  C_6_cache = C_6;
507  }
508 
509  if (C_7 == C_7_cache) {
510  C_7_updated = 1;
511  } else {
512  C_7_updated = 0;
513  C_7_cache = C_7;
514  }
515 
516  if (C_9 == C_9_cache) {
517  C_9_updated = 1;
518  } else {
519  C_9_updated = 0;
520  C_9_cache = C_9;
521  }
522 
523  if (C_10 == C_10_cache) {
524  C_10_updated = 1;
525  } else {
526  C_10_updated = 0;
527  C_10_cache = C_10;
528  }
529 
530  if (C_S == C_S_cache) {
531  C_S_updated = 1;
532  } else {
533  C_S_updated = 0;
534  C_S_cache = C_S;
535  }
536 
537  if (C_P == C_P_cache) {
538  C_P_updated = 1;
539  } else {
540  C_P_updated = 0;
541  C_P_cache = C_P;
542  }
543 
544  if (C_7p == C_7p_cache) {
545  C_7p_updated = 1;
546  } else {
547  C_7p_updated = 0;
548  C_7p_cache = C_7p;
549  }
550 
551  if (C_9p == C_9p_cache) {
552  C_9p_updated = 1;
553  } else {
554  C_9p_updated = 0;
555  C_9p_cache = C_9p;
556  }
557 
558  if (C_10p == C_10p_cache) {
559  C_10p_updated = 1;
560  } else {
561  C_10p_updated = 0;
562  C_10p_cache = C_10p;
563  }
564 
565  if (C_Sp == C_Sp_cache) {
566  C_Sp_updated = 1;
567  } else {
568  C_Sp_updated = 0;
569  C_Sp_cache = C_Sp;
570  }
571 
572  if (C_Pp == C_Pp_cache) {
573  C_Pp_updated = 1;
574  } else {
575  C_Pp_updated = 0;
576  C_Pp_cache = C_Pp;
577  }
578 
579  if (C_2Lh_bar == C_2Lh_cache) {
580  C_2Lh_updated = 1;
581  } else {
582  C_2Lh_updated = 0;
583  C_2Lh_cache = C_2Lh_bar;
584  }
585 
586  if (C_8Lh == C_8Lh_cache) {
587  C_8Lh_updated = 1;
588  } else {
589  C_8Lh_updated = 0;
590  C_8Lh_cache = C_8Lh;
591  }
592 
593  if (Mb == Ycache(0) && Mc == Ycache(1)) {
594  Yupdated = C_1_updated * C_2_updated * C_3_updated * C_4_updated * C_5_updated * C_6_updated;
595  } else {
596  Yupdated = 0;
597  Ycache(0) = Mb;
598  Ycache(1) = Mc;
599  }
600 
601  if (!dispersion) {
602  if (MM == H_V0cache(0) && Mb == H_V0cache(1) && h_0 == H_V0Ccache[0] && h_1 == H_V0Ccache[1] && h_2 == H_V0Ccache[2]) {
603  H_V0updated = N_updated * C_9_updated * Yupdated * VL_updated * C_9p_updated * C_7_updated * TL_updated * C_7p_updated;
604  } else {
605  H_V0updated = 0;
606  H_V0cache(0) = MM;
607  H_V0cache(1) = Mb;
608  H_V0Ccache[0] = h_0;
609  H_V0Ccache[1] = h_1;
610  H_V0Ccache[2] = h_2;
611  }
612  } else {
613  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]) {
614  H_V0updated = N_updated * C_9_updated * Yupdated * VL_updated * C_9p_updated * C_7_updated * TL_updated * C_7p_updated;
615  } else {
616  H_V0updated = 0;
617  H_V0cache(0) = MM;
618  H_V0cache(1) = Mb;
619  H_V0Ccache_dispersion[0] = r_1;
620  H_V0Ccache_dispersion[1] = r_2;
621  H_V0Ccache_dispersion[2] = Delta_C9;
622  H_V0Ccache_dispersion[3] = exp_Phase;
623  }
624  }
625 
626  H_A0updated = N_updated * C_10_updated * VL_updated * C_10p_updated;
627 
628  if (Mb == H_Scache(0) && MW == H_Scache(1)) {
629  H_Supdated = N_updated * C_S_updated * SL_updated * C_Sp_updated;
630  } else {
631  H_Supdated = 0;
632  H_Scache(0) = Mb;
633  H_Scache(1) = MW;
634  }
635 
636  if (Mb == H_P_cache(0) && MW == H_P_cache(1) && Mlep == H_P_cache(2) && Ms == H_P_cache(3)) {
637  H_P_updated = N_updated * C_P_updated * SL_updated * C_Pp_updated * SL_updated * C_10_updated * C_10p_updated;
638  } else {
639  H_P_updated = 0;
640  H_P_cache(0) = Mb;
641  H_P_cache(1) = MW;
642  H_P_cache(2) = Mlep;
643  H_P_cache(3) = Ms;
644  }
645 
646  if (MM == T_cache(0) && Mb == T_cache(1) && Mc == T_cache(2) &&
647  mySM.getMesons(pseudoscalar).getGegenalpha(0) == T_cache(3) && mySM.getMesons(pseudoscalar).getGegenalpha(1) == T_cache(4)) {
648  T_updated = 1;
649  } else {
650  T_updated = 0;
651  T_cache(0) = MM;
652  T_cache(1) = Mb;
653  T_cache(2) = Mc;
654  T_cache(3) = mySM.getMesons(pseudoscalar).getGegenalpha(0);
655  T_cache(4) = mySM.getMesons(pseudoscalar).getGegenalpha(1);
656  }
657 
658  deltaTparpupdated = C_2Lh_updated * T_updated;
659  deltaTparmupdated = C_2Lh_updated * C_8Lh_updated * T_updated;
660 
661  I0_updated = F_updated * H_V0updated * H_A0updated * H_P_updated * beta_updated * H_Supdated * deltaTparmupdated;
662  I2_updated = F_updated * beta_updated * H_V0updated * H_A0updated * deltaTparmupdated;
663  I8_updated = F_updated * beta_updated * H_Supdated * H_V0updated * deltaTparmupdated;
664 
665 }
666 
667 /*******************************************************************************
668  * Transverse Form Factors *
669  * ****************************************************************************/
670 double MPll::LCSR_fit1(double q2, double r_1, double r_2, double m_fit2)
671 {
672  return r_1 / (1. - q2 / m_fit2) + r_2 / pow((1. - q2 / m_fit2), 2.);
673 
674 }
675 
676 double MPll::LCSR_fit2(double q2, double r_2, double m_fit2)
677 {
678  return r_2 / (1. - q2 / m_fit2);
679 }
680 
681 double MPll::zeta(double q2)
682 {
683  double tp, t0;
684 
685  tp = (MM + MP)*(MM + MP);
686  t0 = (MM + MP)*(sqrt(MM) - sqrt(MP))*(sqrt(MM) - sqrt(MP));
687 
688  return (sqrt(tp - q2) - sqrt(tp - t0)) / (sqrt(tp - q2) + sqrt(tp - t0));
689 }
690 
691 double MPll::LATTICE_fit1(double q2, double b_0, double b_1, double b_2, double m_fit2)
692 {
693  double z2 = zeta(q2) * zeta(q2);
694  double z3 = zeta(q2) * z2;
695 
696  return 1. / (1. - q2 / m_fit2) * (b_0 + b_1 * (zeta(q2) - 1. / 3. * z3) + b_2 * (z2 + 2. / 3. * z3));
697 
698 }
699 
700 double MPll::LATTICE_fit2(double q2, double b_0, double b_1, double b_2, double m_fit2)
701 {
702  return 1. / (1. - q2 / m_fit2) * (b_0 + b_1 * zeta(q2) + b_2 * zeta(q2) * zeta(q2));
703 }
704 
705 double MPll::f_plus(double q2)
706 {
707 #if LATTICE
708  return LATTICE_fit1(q2, b_0_fplus, b_1_fplus, b_2_fplus, m_fit2_fplus_lat);
709 #else
710  return LCSR_fit1(q2, r_1_fplus, r_2_fplus, m_fit2_fplus);
711 #endif
712 }
713 
714 double MPll::f_T(double q2)
715 {
716 #if LATTICE
717  return LATTICE_fit1(q2, b_0_fT, b_1_fT, b_2_fT, m_fit2_fT_lat);
718 #else
719  return LCSR_fit1(q2, r_1_fT, r_2_fT, m_fit2_fT);
720 #endif
721 }
722 
723 double MPll::f_0(double q2)
724 {
725 #if LATTICE
726  return LATTICE_fit2(q2, b_0_f0, b_1_f0, b_2_f0, m_fit2_f0_lat);
727 #else
728  return LCSR_fit2(q2, r_2_f0, m_fit2_f0);
729 #endif
730 }
731 
732 gslpp::complex MPll::V_L(double q2)
733 {
734  return gslpp::complex::i() * sqrt(lambda(q2)) / (twoMM * sqrt(q2)) * f_plus(q2);
735 }
736 
737 gslpp::complex MPll::T_L(double q2)
738 {
739  return gslpp::complex::i() * sqrt(lambda(q2) * q2) / (twoMM2_MMpMP) * f_T(q2);
740 }
741 
742 double MPll::S_L(double q2)
743 {
744  return S_L_pre * f_0(q2);
745 }
746 
747 /*******************************************************************************
748  * QCDF *
749  * ****************************************************************************/
750 
751 gslpp::complex MPll::I1(double u, double q2)
752 {
753  std::pair<double, double > uq2 = std::make_pair(u, q2);
754 
755  if (I1Cached[uq2] == 0) {
756  ubar = 1. - u;
757  xp = .5 + sqrt(0.25 - (Mc2 - gslpp::complex::i()*1.e-10) / (ubar * MM2 + u * q2));
758  xm = .5 - sqrt(0.25 - (Mc2 - gslpp::complex::i()*1.e-10) / (ubar * MM2 + u * q2));
759  yp = .5 + sqrt(0.25 - (Mc2 - gslpp::complex::i()*1.e-10) / q2);
760  ym = .5 - sqrt(0.25 - (Mc2 - gslpp::complex::i()*1.e-10) / q2);
761  L1xp = log(1. - 1. / xp) * log(1. - xp) - M_PI2osix + dilog(xp / (xp - 1.));
762  L1xm = log(1. - 1. / xm) * log(1. - xm) - M_PI2osix + dilog(xm / (xm - 1.));
763  L1yp = log(1. - 1. / yp) * log(1. - yp) - M_PI2osix + dilog(yp / (yp - 1.));
764  L1ym = log(1. - 1. / ym) * log(1. - ym) - M_PI2osix + dilog(ym / (ym - 1.));
765 
766  cacheI1[uq2] = 1. + twoMc2 / ubar / (MM2 - q2)*(L1xp + L1xm - L1yp - L1ym);
767  I1Cached[uq2] = 1;
768  }
769 
770  return cacheI1[uq2];
771 }
772 
773 gslpp::complex MPll::Tparplus(double u, double q2)
774 {
775  Ee = (MM2 - q2) / twoMM;
776  ubar = 1. - u;
777  arg1 = (fourMc2 - gslpp::complex::i()*1.e-10) / (ubar * MM2 + u * q2) - 1.;
778  B01 = -2. * sqrt(arg1) * arctan(1. / sqrt(arg1));
779  arg1 = (fourMc2 - gslpp::complex::i()*1.e-10) / q2 - 1.;
780  B00 = -2. * sqrt(arg1) * arctan(1. / sqrt(arg1));
781 
782  gslpp::complex tpar = twoMM / Ee / ubar * I1(u, q2) + (ubar * MM2 + u * q2) / Ee / Ee / ubar / ubar * (B01 - B00);
783  return -MM / Mb * mySM.getQuarks(QCD::CHARM).getCharge() * tpar*C_2Lh_bar;
784 }
785 
786 gslpp::complex MPll::Tparminus(double u, double q2)
787 {
788  double ubar = 1. - u;
789  return -spectator_charge * (8. * C_8Lh / (ubar + u * q2 / MM2)
790  + sixMMoMb * H_c(ubar * MM2 + u * q2, mu_h * mu_h) * C_2Lh_bar);
791 }
792 
793 double MPll::Integrand_ReTparplus(double up)
794 {
795  double u = up;
796  return ((Tparplus(u, tmpq2)*6. * u * (1. - u)*
797  (1 + threeGegen0 * (2. * u - 1)
798  + threeGegen1otwo * ((10. * u - 5.)*(2. * u - 1.) - 1.))) / mySM.getMesons(meson).getLambdaM()).real();
799 }
800 
801 double MPll::Integrand_ImTparplus(double up)
802 {
803  double u = up;
804  return ((Tparplus(u, tmpq2)*6. * u * (1. - u)*
805  (1 + threeGegen0 * (2. * u - 1)
806  + threeGegen1otwo * ((10. * u - 5.)*(2. * u - 1.) - 1.))) / mySM.getMesons(meson).getLambdaM()).imag();
807 }
808 
809 double MPll::Integrand_ReTparminus(double up)
810 {
811  double Lambdaplus = mySM.getMesons(meson).getLambdaM();
812  gslpp::complex Lambdamin = exp(-tmpq2 / MM / Lambdaplus) / Lambdaplus * (-gsl_sf_expint_Ei(tmpq2 / MM / Lambdaplus) + gslpp::complex::i() * M_PI);
813 
814  double u = up;
815  return ((Tparminus(u, tmpq2)*6. * u * (1. - u)*
816  (1 + threeGegen0 * (2. * u - 1)
817  + threeGegen1otwo * ((10. * u - 5.)*(2. * u - 1.) - 1.))) / Lambdamin).real();
818 }
819 
820 double MPll::Integrand_ImTparminus(double up)
821 {
822  double Lambdaplus = mySM.getMesons(meson).getLambdaM();
823  gslpp::complex Lambdamin = exp(-tmpq2 / MM / Lambdaplus) / Lambdaplus * (-gsl_sf_expint_Ei(tmpq2 / MM / Lambdaplus) + gslpp::complex::i() * M_PI);
824 
825  double u = up;
826  return ((Tparminus(u, tmpq2)*6. * u * (1. - u)*
827  (1 + threeGegen0 * (2. * u - 1)
828  + threeGegen1otwo * ((10. * u - 5.)*(2. * u - 1.) - 1.))) / Lambdamin).imag();
829 }
830 
831 double MPll::Integrand_ImTpar_pm(double up)
832 {
833  return Integrand_ImTparplus(up) + Integrand_ImTparminus(up);
834 }
835 
836 double MPll::Integrand_ReTpar_pm(double up)
837 {
838  return Integrand_ReTparplus(up) + Integrand_ReTparminus(up);
839 }
840 
841 gslpp::complex MPll::F19(double q2)
842 {
843  double s = q2 / Mb2;
844  double s2 = s*s;
845  double Ls = log(s);
846  double Lc = log(Mc / Mb);
847  double Lm = log(mu_b / Mb);
848  gslpp::complex i = gslpp::complex::i();
849  return (-1424. / 729. + 16. / 243. * i * M_PI + 64. / 27. * Lc)*Lm - 16. / 243. * Lm * Ls + (16. / 1215. - 32. / 135. / Mc2 * Mb2) * Lm * s
850  + (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
851  + (-11.65 + 0.18223 * i + (-24.709 - 0.13474 * i) * s + (-43.588 - 0.4738 * i) * s2 + (-86.22 - 1.3542 * i) * s * s2
852  + (-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);
853 }
854 
855 gslpp::complex MPll::F27(double q2)
856 {
857  double s = q2 / Mb2;
858  double s2 = s*s;
859  double Ls = log(s);
860  gslpp::complex i = gslpp::complex::i();
861  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;
862 }
863 
864 gslpp::complex MPll::F29(double q2)
865 {
866  double s = q2 / Mb2;
867  double s2 = s*s;
868  double Ls = log(s);
869  gslpp::complex i = gslpp::complex::i();
870  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;
871 }
872 
873 gslpp::complex MPll::F87(double q2)
874 {
875  double s = q2 / Mb2;
876  double s2 = s*s;
877  return F87_0 + F87_1 * s + F87_2 * s2 + F87_3 * s * s2 - 0.888889 * log(s)*(1. + s + s2 + s * s2);
878 }
879 
880 double MPll::F89(double q2)
881 {
882  double s = q2 / Mb2;
883  double s2 = s*s;
884  return F89_0 + F89_1 * s + F89_2 * s2 + F89_3 * s * s2 + 1.77778 * log(s)*(1. + s + s2 + s * s2);
885 }
886 
887 gslpp::complex MPll::Cpar(double q2)
888 {
889  return -(C_2L_bar * F27(q2) + C_8L * F87(q2) + MM / 2. / Mb *
890  (C_2L_bar * F29(q2) + 2. * C_1L_bar * (F19(q2) + F29(q2) / 6.) + C_8L * F89(q2)));
891 }
892 
893 gslpp::complex MPll::deltaTpar(double q2)
894 {
895  tmpq2 = q2;
896 
897  //old_handler = gsl_set_error_handler_off();
898 
899  if (deltaTparpCached[q2] == 0) {
900 
901  DTPPR = convertToGslFunction(boost::bind(&MPll::Integrand_ReTpar_pm, &(*this), _1));
902  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();
903  double ReTppint = avaDTPPR;
904 
905  DTPPR = convertToGslFunction(boost::bind(&MPll::Integrand_ImTpar_pm, &(*this), _1));
906  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();
907  double ImTppint = avaDTPPR;
908 
909  cacheDeltaTparp[q2] = (ReTppint + gslpp::complex::i() * ImTppint);
910  deltaTparpCached[q2] = 1;
911  }
912 
913  //gsl_set_error_handler(old_handler);
914 
915  return deltaT_0 * Cpar(q2) + deltaT_1par * cacheDeltaTparp[q2] / f_plus(q2);
916 }
917 
918 double MPll::reDC9fit(double* x, double* p)
919 {
920  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];
921 }
922 
923 double MPll::imDC9fit(double* x, double* p)
924 {
925  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];
926 
927  //double thr = 4.*Mc2;
928 }
929 
930 void MPll::fit_DeltaC9_mumu()
931 {
932  int dim = 0;
933  for (double i = 0.1; i < MPllSWITCH; i += 0.4) {
934  double q2tmp = i;
935  myq2.push_back(q2tmp);
936  ReDeltaC9.push_back((deltaTpar(q2tmp)).real());
937  ImDeltaC9.push_back((deltaTpar(q2tmp)).imag());
938  dim++;
939  }
940  for (double i = MPllSWITCH; i < 8.2; i += 0.4) {
941  double q2tmp = i;
942  myq2.push_back(q2tmp);
943  ReDeltaC9.push_back(q2tmp * (deltaTpar(q2tmp)).real());
944  ImDeltaC9.push_back(q2tmp * (deltaTpar(q2tmp)).imag());
945  dim++;
946  }
947 
948  gr1 = TGraph(dim, myq2.data(), ReDeltaC9.data());
949  gr2 = TGraph(dim, myq2.data(), ImDeltaC9.data());
950 
951  reffit = TF1("reffit", this, &MPll::reDC9fit, 0, 8.1, 7, "MPll", "reDC9fit");
952  imffit = TF1("imffit", this, &MPll::imDC9fit, 0, 8.1, 8, "MPll", "imDC9fit");
953 
954  refres = gr1.Fit(&reffit, "SQN0+rob=0.99");
955  imfres = gr2.Fit(&imffit, "SQN0+rob=0.99");
956 
957  ReDeltaC9.clear();
958  ImDeltaC9.clear();
959  myq2.clear();
960 }
961 
962 gslpp::complex MPll::fDeltaC9(double q2)
963 {
964  if (q2 < MPllSWITCH) return (reDC9fit(&q2, const_cast<double *> (refres->GetParams()))
965  + gslpp::complex::i() * imDC9fit(&q2, const_cast<double *> (imfres->GetParams())));
966  else return (reDC9fit(&q2, const_cast<double *> (refres->GetParams()))
967  + gslpp::complex::i() * imDC9fit(&q2, const_cast<double *> (imfres->GetParams()))) / q2;
968 
969 }
970 
971 gslpp::complex MPll::DeltaC9(double q2)
972 {
973  return deltaTpar(q2);
974 }
975 
976 gslpp::complex MPll::deltaC7_QCDF(double q2, bool spline)
977 {
978 #if SPLINE
979  if (spline) return gsl_spline_eval(spline_Re_deltaC7_QCDF, q2, acc_Re_deltaC7_QCDF);
980 #endif
981 
982  double muh = mu_b / mb_pole;
983  double z = mc_pole * mc_pole / mb_pole / mb_pole;
984  double sh = q2 / mb_pole / mb_pole;
985  double sh2 = sh*sh;
986 
987  //gslpp::complex A_Sdl = A_Seidel(q2, mb_pole*mb_pole); /* hep-ph/0403185v2.*/
988  //gslpp::complex Fu_17 = -A_Sdl; /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
989  //gslpp::complex Fu_27 = 6. * A_Sdl; /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
990  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*/
991  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*/
992  gslpp::complex F_87 = F87_0 + F87_1 * sh + F87_2 * sh2 + F87_3 * sh * sh2 - 8. / 9. * log(sh) * (sh + sh2 + sh * sh2);
993 
994  gslpp::complex delta = C_1 * F_17 + C_2 * F_27;
995  gslpp::complex delta_t = C_8 * F_87 + delta;
996  //gslpp::complex delta_u = delta + C_1 * Fu_17 + C_2 * Fu_27;
997 
998  return -alpha_s_mub / (4. * M_PI) * (delta_t /*- lambda_u / lambda_t * delta_u */);
999 }
1000 
1001 gslpp::complex MPll::deltaC9_QCDF(double q2, bool spline)
1002 {
1003 #if SPLINE
1004  if (spline) return gsl_spline_eval(spline_Re_deltaC9_QCDF, q2, acc_Re_deltaC9_QCDF);
1005 #endif
1006  double muh = mu_b / mb_pole;
1007  double z = mc_pole * mc_pole / mb_pole / mb_pole;
1008  double sh = q2 / mb_pole / mb_pole;
1009  double sh2 = sh*sh;
1010 
1011  //gslpp::complex B_Sdl = B_Seidel(q2, mb_pole*mb_pole); /* hep-ph/0403185v2.*/
1012  //gslpp::complex C_Sdl = C_Seidel(q2); /* hep-ph/0403185v2.*/
1013  //gslpp::complex Fu_19 = -(B_Sdl + 4. * C_Sdl); /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1014  //gslpp::complex Fu_29 = -(-6. * B_Sdl + 3. * C_Sdl); /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1015  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*/
1016  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*/
1017  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));
1018 
1019  gslpp::complex delta = C_1 * F_19 + C_2 * F_29;
1020  gslpp::complex delta_t = C_8 * F_89 + delta;
1021  //gslpp::complex delta_u = delta + C_1 * Fu_19 + C_2 * Fu_29;
1022 
1023  return -alpha_s_mub / (4. * M_PI) * (delta_t /*- lambda_u / lambda_t * delta_u*/);
1024 }
1025 
1026 void MPll::spline_QCDF_func()
1027 {
1028  int dim_DC = GSL_INTERP_DIM_DC;
1029  double min = 0.001;
1030  double interval_DC = (9.9 - min) / ((double) dim_DC);
1031  double q2_spline_DC[dim_DC];
1032  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];
1033 
1034  for (int i = 0; i < dim_DC; i++) {
1035  q2_spline_DC[i] = min + (double) i*interval_DC;
1036  fq2_Re_deltaC7_QCDF[i] = deltaC7_QCDF(q2_spline_DC[i], false).real();
1037  fq2_Im_deltaC7_QCDF[i] = deltaC7_QCDF(q2_spline_DC[i], false).imag();
1038  fq2_Re_deltaC9_QCDF[i] = deltaC9_QCDF(q2_spline_DC[i], false).real();
1039  fq2_Im_deltaC9_QCDF[i] = deltaC9_QCDF(q2_spline_DC[i], false).imag();
1040 
1041  }
1042 
1043  gsl_spline_init(spline_Re_deltaC7_QCDF, q2_spline_DC, fq2_Re_deltaC7_QCDF, dim_DC);
1044  gsl_spline_init(spline_Im_deltaC7_QCDF, q2_spline_DC, fq2_Im_deltaC7_QCDF, dim_DC);
1045  gsl_spline_init(spline_Re_deltaC9_QCDF, q2_spline_DC, fq2_Re_deltaC9_QCDF, dim_DC);
1046  gsl_spline_init(spline_Im_deltaC9_QCDF, q2_spline_DC, fq2_Im_deltaC9_QCDF, dim_DC);
1047 
1048 
1049 }
1050 
1051 /*******************************************************************************
1052  * Helicity amplitudes *
1053  * ****************************************************************************/
1054 gslpp::complex MPll::H_c(double q2, double mu2)
1055 {
1056  double x = fourMc2 / q2;
1057  gslpp::complex par;
1058 
1059  if (x > 1.) par = sqrt(x - 1.) * atan(1. / sqrt(x - 1.));
1060  else par = sqrt(1. - x) * (log((1. + sqrt(1. - x)) / sqrt(x)) - ihalfMPI);
1061  return -fournineth * (log(Mc2 / mu2) - twothird - x) - fournineth * (2. + x) * par;
1062 }
1063 
1064 gslpp::complex MPll::H_b(double q2, double mu2)
1065 {
1066  double x = fourMb2 / q2;
1067  gslpp::complex par;
1068 
1069  if (x > 1.) par = sqrt(x - 1.) * atan(1. / sqrt(x - 1.));
1070  else par = sqrt(1. - x) * (log((1. + sqrt(1. - x)) / sqrt(x)) - ihalfMPI);
1071 
1072  return -fournineth * (log(Mb2 / mu2) - twothird - x) - fournineth * (2. + x) * par;
1073 }
1074 
1075 gslpp::complex MPll::H_0(double q2)
1076 {
1077  return (H_0_pre - fournineth * log(q2 / mu_b2));
1078 }
1079 
1080 gslpp::complex MPll::Y(double q2)
1081 {
1082  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;
1083 }
1084 
1085 gslpp::complex MPll::funct_g(double q2)
1086 {
1087  if (q2 < 4. * Mc * Mc)
1088  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.)));
1089  else
1090  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));
1091 }
1092 
1093 gslpp::complex MPll::DeltaC9_KD(double q2)
1094 {
1095  return ((Delta_C9 + r_1 * q2 / mJ2) / (1. - r_2 * q2 / mJ2) - (3. * (-0.267) + 1.117) * funct_g(q2))*exp_Phase;
1096  /* C_1 = -0.267 and C_2 = 1.117 in KMPW */
1097 }
1098 
1099 gslpp::complex MPll::h_lambda(double q2)
1100 {
1101  if (!dispersion) {
1102  return (twoMboMM * h_0 * T_L(q2) + h_1 * q2 / MM2 * V_L(q2)) / sixteenM_PI2 + h_2 * q2 * q2;
1103  } else {
1104  return -q2 / (MM2 * sixteenM_PI2) * V_L(q2) * DeltaC9_KD(q2);
1105  }
1106 }
1107 
1108 gslpp::complex MPll::H_V(double q2)
1109 {
1110  return -((C_9 + deltaC9_QCDF(q2, SPLINE) + Y(q2) /*+ fDeltaC9(q2)*/ - etaP * pow(-1, angmomP) * C_9p) * V_L(q2) + MM2 / q2 * (twoMboMM * (C_7 + deltaC7_QCDF(q2, SPLINE) - etaP * pow(-1, angmomP) * C_7p) * T_L(q2) - sixteenM_PI2 * h_lambda(q2)));
1111 }
1112 
1113 gslpp::complex MPll::H_A(double q2)
1114 {
1115  return (-C_10 + etaP * pow(-1, angmomP) * C_10p) *V_L(q2);
1116 }
1117 
1118 gslpp::complex MPll::H_S(double q2)
1119 {
1120  return MboMW * (C_S - etaP * pow(-1, angmomP) * C_Sp) * S_L(q2);
1121 }
1122 
1123 gslpp::complex MPll::H_P(double q2)
1124 {
1125  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);
1126 }
1127 
1128 /*******************************************************************************
1129  * Angular coefficients *
1130  * ****************************************************************************/
1131 double MPll::k2(double q2)
1132 {
1133  return (MM4 + q2 * q2 + MP4 - twoMP2 * q2 - twoMM2 * (q2 + MP2)) / fourMM2;
1134 }
1135 
1136 double MPll::beta(double q2)
1137 {
1138  return sqrt(1. - 4. * Mlep2 / q2);
1139 }
1140 
1141 double MPll::beta2(double q2)
1142 {
1143  return 1. - 4. * Mlep2 / q2;
1144 }
1145 
1146 double MPll::lambda(double q2)
1147 {
1148  return 4. * MM2 * k2(q2);
1149 }
1150 
1151 double MPll::F(double q2)
1152 {
1153  return sqrt(lambda(q2)) * beta(q2) * q2 / (ninetysixM_PI3MM3);
1154 }
1155 
1156 double MPll::I_1c(double q2)
1157 {
1158  return F(q2)*((H_V(q2).abs2() + H_A(q2).abs2()) / 2. + H_P(q2).abs2() + 2. * Mlep2 / q2 * (H_V(q2).abs2()
1159  - H_A(q2).abs2()) + beta2(q2) * H_S(q2).abs2());
1160 }
1161 
1162 double MPll::I_2c(double q2)
1163 {
1164  return -F(q2) * beta2(q2) / 2. * (H_V(q2).abs2() + H_A(q2).abs2());
1165 }
1166 
1167 double MPll::I_6c(double q2)
1168 {
1169  return 4. * F(q2) * beta(q2) * Mlep / sqrt(q2)*(H_S(q2).conjugate() * H_V(q2)).real();
1170 }
1171 
1172 double MPll::Delta(int i, double q2)
1173 {
1174  return 0; /* FIX CPV */
1175  //return (I(i, q2,0) - I(i, q2,1))/2;
1176 }
1177 
1178 double MPll::integrateSigma(int i, double q_min, double q_max)
1179 {
1180  updateParameters();
1181 
1182  std::pair<double, double > qbin = std::make_pair(q_min, q_max);
1183 
1184  old_handler = gsl_set_error_handler_off();
1185 
1186  switch (i) {
1187  case 0:
1188  if (sigma0Cached[qbin] == 0) {
1189  FS = convertToGslFunction(boost::bind(&MPll::getSigma1c, &(*this), _1));
1190  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();
1191  cacheSigma0[qbin] = NN*avaSigma;
1192  sigma0Cached[qbin] = 1;
1193  }
1194  return cacheSigma0[qbin];
1195  break;
1196  case 2:
1197  if (sigma2Cached[qbin] == 0) {
1198  FS = convertToGslFunction(boost::bind(&MPll::getSigma2c, &(*this), _1));
1199  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();
1200  cacheSigma2[qbin] = NN*avaSigma;
1201  sigma2Cached[qbin] = 1;
1202  }
1203  return cacheSigma2[qbin];
1204  break;
1205  default:
1206  std::stringstream out;
1207  out << i;
1208  throw std::runtime_error("MPll::integrateSigma: index " + out.str() + " not implemented");
1209  }
1210 
1211  gsl_set_error_handler(old_handler);
1212 
1213 }
1214 
1215 double MPll::getSigma(int i, double q_2)
1216 {
1217  updateParameters();
1218 
1219  switch (i) {
1220  case 0:
1221  return getSigma1c(q_2);
1222  break;
1223  case 2:
1224  return getSigma2c(q_2);
1225  break;
1226  case 8:
1227  return getSigma6c(q_2);
1228  break;
1229  default:
1230  std::stringstream out;
1231  out << i;
1232  throw std::runtime_error("MPll::getSigma: index " + out.str() + " not implemented");
1233  }
1234 }
1235 
1236 double MPll::integrateDelta(int i, double q_min, double q_max)
1237 {
1238  updateParameters();
1239 
1240  std::pair<double, double > qbin = std::make_pair(q_min, q_max);
1241 
1242  old_handler = gsl_set_error_handler_off();
1243 
1244  switch (i) {
1245  case 0:
1246  if (delta0Cached[qbin] == 0) {
1247  FD = convertToGslFunction(boost::bind(&MPll::getDelta1c, &(*this), _1));
1248  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();
1249  cacheDelta0[qbin] = NN*avaDelta;
1250  delta0Cached[qbin] = 1;
1251  }
1252  return cacheDelta0[qbin];
1253  break;
1254  case 2:
1255  if (delta2Cached[qbin] == 0) {
1256  FD = convertToGslFunction(boost::bind(&MPll::getDelta2c, &(*this), _1));
1257  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();
1258  cacheDelta2[qbin] = NN*avaDelta;
1259  delta2Cached[qbin] = 1;
1260  }
1261  return cacheDelta2[qbin];
1262  break;
1263  default:
1264  std::stringstream out;
1265  out << i;
1266  throw std::runtime_error("MPll::integrateDelta: index " + out.str() + " not implemented");
1267  }
1268 
1269  gsl_set_error_handler(old_handler);
1270 
1271 }
MPll::dispersion
bool dispersion
Definition: MPll.h:282
gslpp_special_functions::zeta
double zeta(int i)
Definition: gslpp_special_functions.cpp:20
MPll::H_P
gslpp::complex H_P(double q2)
The helicity amplitude .
Definition: MPll.cpp:1120
std_make_vector.h
MPll::H_V
gslpp::complex H_V(double q2)
The helicity amplitude .
Definition: MPll.cpp:1105
QCD::BOTTOM
Definition: QCD.h:329
MPll::mc_pole
double mc_pole
Definition: MPll.h:296
MPll::MP
double MP
Definition: MPll.h:290
gslpp::arctan
complex arctan(const complex &z)
Definition: gslpp_complex.cpp:492
MPll::pseudoscalar
QCD::meson pseudoscalar
Definition: MPll.h:278
MPll.h
gslpp::dilog
complex dilog(const complex &z)
Definition: gslpp_complex.cpp:370
F_2.h
MPll::meson
QCD::meson meson
Definition: MPll.h:277
MPll::~MPll
virtual ~MPll()
Destructor.
Definition: MPll.cpp:72
CKM::computelamt_s
gslpp::complex computelamt_s() const
The product of the CKM elements .
Definition: CKM.cpp:136
QCD::initializeMeson
void initializeMeson(QCD::meson meson_i) const
A method to initialize a meson.
Definition: QCD.cpp:236
make_vector
Definition: std_make_vector.h:15
LO
Definition: OrderScheme.h:33
QCD::UP
Definition: QCD.h:324
MPll::H_A
gslpp::complex H_A(double q2)
The helicity amplitude .
Definition: MPll.cpp:1110
StandardModel.h
F_1.h
QCD::CHARM
Definition: QCD.h:326
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:218
MPll::Ms
double Ms
Definition: MPll.h:297
MPll::ale
double ale
Definition: MPll.h:287
gslpp::complex
A class for defining operations on and functions of complex numbers.
Definition: gslpp_complex.h:35
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:243
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
h_0
A class for the correction in .
Definition: MVllObservables.h:1587
Meson::getLambdaM
const double & getLambdaM() const
Definition: Meson.h:402
MPll::myF_1
std::unique_ptr< F_1 > myF_1
Definition: MPll.h:280
MPll::spectator_charge
double spectator_charge
Definition: MPll.h:298
gslpp::complex::abs2
double abs2() const
Definition: gslpp_complex.cpp:86
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:474
MPll::mpllParameters
std::vector< std::string > mpllParameters
Definition: MPll.h:279
QCD::Mbar2Mp
double Mbar2Mp(const double mbar, const orders order=FULLNNLO) const
Converts the mass to the pole mass.
Definition: QCD.cpp:1221
Meson::computeWidth
double computeWidth() const
A method to compute the width of the meson from its lifetime.
Definition: Meson.cpp:479
MPll::integrateSigma
double integrateSigma(int i, double q_min, double q_max)
The integral of from to .
Definition: MPll.cpp:1174
gslpp::complex::conjugate
complex conjugate() const
Definition: gslpp_complex.cpp:288
Particle::getMass
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
QCD::meson
meson
An enum type for mesons.
Definition: QCD.h:336
MPll::integrateDelta
double integrateDelta(int i, double q_min, double q_max)
The integral of from to .
Definition: MPll.cpp:1232
F_2
Definition: F_2.h:15
MPll::lep
QCD::lepton lep
Definition: MPll.h:276
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
StandardModel::getGF
double getGF() const
A get method to retrieve the Fermi constant .
Definition: StandardModel.h:736
gslpp::sqrt
complex sqrt(const complex &z)
Definition: gslpp_complex.cpp:385
MPll::MPll
MPll(const StandardModel &SM_i, QCD::meson meson_i, QCD::meson pseudoscalar_i, QCD::lepton lep_i)
Constructor.
Definition: MPll.cpp:19
QCD::getMesons
Meson getMesons(const QCD::meson m) const
A get method to access a meson as an object of the type Meson.
Definition: QCD.h:524
gslpp::complex::i
static const complex & i()
Definition: gslpp_complex.cpp:154
Particle::getCharge
double getCharge() const
A get method to access the particle charge.
Definition: Particle.h:97
Meson::getDecayconst
const double & getDecayconst() const
A get method for the decay constant of the meson.
Definition: Meson.h:360
MPll::Mc
double Mc
Definition: MPll.h:294
StandardModel::Als
double Als(double mu, orders order=FULLNLO, bool qed_flag=false, bool Nf_thr=true) const
The running QCD coupling in the scheme including QED corrections.
Definition: StandardModel.cpp:576
gslpp_function_adapter.h
QCD::getQuarks
Particle getQuarks(const QCD::quark q) const
A get method to access a quark as an object of the type Particle.
Definition: QCD.h:534
MPll::H_S
gslpp::complex H_S(double q2)
The helicity amplitude .
Definition: MPll.cpp:1115
MPll::FixedWCbtos
bool FixedWCbtos
Definition: MPll.h:283
MPll::Mb
double Mb
Definition: MPll.h:291
MPll::MM
double MM
Definition: MPll.h:289
QCD::K_P
Definition: QCD.h:340
StandardModel::getFlavour
const Flavour & getFlavour() const
Definition: StandardModel.h:1006
MPll::GF
double GF
Definition: MPll.h:286
C_4
Definition: Doxygen/examples-src/myModel/src/myObservables.h:76
MPll::mu_b
double mu_b
Definition: MPll.h:292
MPll::Mlep
double Mlep
Definition: MPll.h:288
F_1
Definition: F_1.h:15
MPll::mJ2
double mJ2
Definition: MPll.h:284
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:113
QCD::getMub
double getMub() const
A get method to access the threshold between five- and four-flavour theory in GeV.
Definition: QCD.h:570
QCD::K_0
Definition: QCD.h:339
QCD::STRANGE
Definition: QCD.h:327
gslpp::complex::real
const double & real() const
Definition: gslpp_complex.cpp:53
MPll::mu_h
double mu_h
Definition: MPll.h:293
MPll::getSigma
double getSigma(int i, double q_2)
The value of from to .
Definition: MPll.cpp:1211
MPll::mySM
const StandardModel & mySM
Definition: MPll.h:275
QCD::getOptionalParameter
double getOptionalParameter(std::string name) const
A method to get parameters that are specific to only one set of observables.
Definition: QCD.h:448
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:118
MPll::h_lambda
gslpp::complex h_lambda(double q2)
The non-pertubative ccbar contributions to the helicity amplitudes.
Definition: MPll.cpp:1096
NLO
Definition: OrderScheme.h:34
MPll::width
double width
Definition: MPll.h:299
StandardModel::Mw
virtual double Mw() const
The SM prediction for the -boson mass in the on-shell scheme, .
Definition: StandardModel.cpp:944
convertToGslFunction
gsl_function convertToGslFunction(const F &f)
Definition: gslpp_function_adapter.h:24
gslpp::exp
complex exp(const complex &z)
Definition: gslpp_complex.cpp:333
MPll::funct_g
gslpp::complex funct_g(double q2)
Definition: MPll.cpp:1082
MPll::mb_pole
double mb_pole
Definition: MPll.h:295
MPll::initializeMPllParameters
std::vector< std::string > initializeMPllParameters()
A method for initializing the parameters necessary for MPll.
Definition: MPll.cpp:76
C_3
Definition: Doxygen/examples-src/myModel/src/myObservables.h:59
MPll::DeltaC9_KD
gslpp::complex DeltaC9_KD(double q2)
Definition: MPll.cpp:1090
Meson::getGegenalpha
const double & getGegenalpha(int i) const
A get method to get the Gegenbaur coefficient.
Definition: Meson.h:394
MPll::myF_2
std::unique_ptr< F_2 > myF_2
Definition: MPll.h:281
QCD::DOWN
Definition: QCD.h:325
FULLNLO
Definition: OrderScheme.h:37
StandardModel::getAle
double getAle() const
A get method to retrieve the fine-structure constant .
Definition: StandardModel.h:745
Flavour::getFlagUseDispersionRelation
bool getFlagUseDispersionRelation() const
Definition: Flavour.h:242
Flavour::getFlagFixedWCbtos
bool getFlagFixedWCbtos() const
Definition: Flavour.h:252
StandardModel::getCKM
CKM getCKM() const
A get method to retrieve the member object of type CKM.
Definition: StandardModel.h:876
QCD::lepton
lepton
An enum type for leptons.
Definition: QCD.h:310
StandardModel::getLeptons
Particle getLeptons(const QCD::lepton p) const
A get method to retrieve the member object of a lepton.
Definition: StandardModel.h:709