a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
MVll.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 "MVll.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_zeta.h>
15 #include <boost/bind.hpp>
16 #include <limits>
17 #include <TFitResult.h>
18 #include <gsl/gsl_sf_gegenbauer.h>
19 #include <gsl/gsl_sf_expint.h>
20 
21 MVll::MVll(const StandardModel& SM_i, QCD::meson meson_i, QCD::meson vector_i, QCD::lepton lep_i)
22 : mySM(SM_i), myF_1(new F_1()), myF_2(new F_2()),
23 N_cache(3, 0.),
24 V_cache(3, 0.),
25 A0_cache(3, 0.),
26 A1_cache(3, 0.),
27 T1_cache(3, 0.),
28 T2_cache(3, 0.),
29 k2_cache(2, 0.),
30 VL0_cache(3, 0.),
31 TL0_cache(3, 0.),
32 SL_cache(2, 0.),
33 Ycache(2, 0.),
34 H_V0cache(2, 0.),
35 H_V1cache(2, 0.),
36 H_V2cache(2, 0.),
37 H_Scache(2, 0.),
38 H_Pcache(4, 0.),
39 T_cache(5, 0.)
40 {
41  lep = lep_i;
42  meson = meson_i;
43  vectorM = vector_i;
44  dispersion = false;
45  FixedWCbtos = false;
46  mJ2 = 3.096 * 3.096;
47 
48  I0_updated = 0;
49  I1_updated = 0;
50  I2_updated = 0;
51  I3_updated = 0;
52  I4_updated = 0;
53  I5_updated = 0;
54  I6_updated = 0;
55  I7_updated = 0;
56  I8_updated = 0;
57  I9_updated = 0;
58  I10_updated = 0;
59  I11_updated = 0;
60 
61  VL1_updated = 0;
62  VL2_updated = 0;
63  TL1_updated = 0;
64  TL2_updated = 0;
65  VR1_updated = 0;
66  VR2_updated = 0;
67  TR1_updated = 0;
68  TR2_updated = 0;
69  VL0_updated = 0;
70  TL0_updated = 0;
71  VR0_updated = 0;
72  TR0_updated = 0;
73  SL_updated = 0;
74  SR_updated = 0;
75 
76  deltaTparpupdated = 0;
77  deltaTparmupdated = 0;
78  deltaTperpupdated = 0;
79 
80  w_sigma = gsl_integration_cquad_workspace_alloc(100);
81  // w_DTPPR = gsl_integration_cquad_workspace_alloc (100);
82  w_delta = gsl_integration_cquad_workspace_alloc(100);
83 
84  acc_Re_T_perp = gsl_interp_accel_alloc();
85  acc_Im_T_perp = gsl_interp_accel_alloc();
86  acc_Re_T_para = gsl_interp_accel_alloc();
87  acc_Im_T_para = gsl_interp_accel_alloc();
88 
89  spline_Re_T_perp = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
90  spline_Im_T_perp = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
91  spline_Re_T_para = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
92  spline_Im_T_para = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
93 
94 #if COMPUTECP
95  acc_Re_T_perp_conj = gsl_interp_accel_alloc();
96  acc_Im_T_perp_conj = gsl_interp_accel_alloc();
97  acc_Re_T_para_conj = gsl_interp_accel_alloc();
98  acc_Im_T_para_conj = gsl_interp_accel_alloc();
99 
100  spline_Re_T_perp_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
101  spline_Im_T_perp_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
102  spline_Re_T_para_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
103  spline_Im_T_para_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
104 #endif
105 
106  acc_Re_deltaC7_QCDF = gsl_interp_accel_alloc();
107  acc_Im_deltaC7_QCDF = gsl_interp_accel_alloc();
108  acc_Re_deltaC9_QCDF = gsl_interp_accel_alloc();
109  acc_Im_deltaC9_QCDF = gsl_interp_accel_alloc();
110 
111  spline_Re_deltaC7_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
112  spline_Im_deltaC7_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
113  spline_Re_deltaC9_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
114  spline_Im_deltaC9_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
115 
116 #if COMPUTECP
117  acc_Re_deltaC7_QCDF_conj = gsl_interp_accel_alloc();
118  acc_Im_deltaC7_QCDF_conj = gsl_interp_accel_alloc();
119  acc_Re_deltaC9_QCDF_conj = gsl_interp_accel_alloc();
120  acc_Im_deltaC9_QCDF_conj = gsl_interp_accel_alloc();
121 
122  spline_Re_deltaC7_QCDF_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
123  spline_Im_deltaC7_QCDF_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
124  spline_Re_deltaC9_QCDF_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
125  spline_Im_deltaC9_QCDF_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
126 #endif
127 
128  h_pole = false;
129 
130  M_PI2 = M_PI*M_PI;
131 
132  F87_1 = (4. / 3. * M_PI2 - 40. / 3.);
133  F87_2 = (32. / 9. * M_PI2 - 316. / 9.);
134  F87_3 = (200. / 27. * M_PI2 - 658. / 9.);
135 
136  F89_0 = (104. / 9. - 32. / 27. * M_PI2);
137  F89_1 = (1184. / 27. - 40. / 9. * M_PI2);
138  F89_2 = (-32. / 3. * M_PI2 + 14212. / 135.);
139  F89_3 = (-560. / 27. * M_PI2 + 193444. / 945.);
140 
141  CF = 4. / 3.;
142 
143 }
144 
145 MVll::~MVll()
146 {
147 }
148 
149 std::vector<std::string> MVll::initializeMVllParameters()
150 {
151  dispersion = mySM.getFlavour().getFlagUseDispersionRelation();
152  FixedWCbtos = mySM.getFlavour().getFlagFixedWCbtos();
153 
154 #if NFPOLARBASIS_MVLL
155  if (vectorM == StandardModel::PHI) mvllParameters = make_vector<std::string>()
156  << "a_0Vphi" << "a_1Vphi" << "a_2Vphi" << "MRV" << "a_0A0phi" << "a_1A0phi" << "a_2A0phi" << "MRA0"
157  << "a_0A1phi" << "a_1A1phi" << "a_2A1phi" << "MRA1" << "a_1A12phi" << "a_2A12phi" << "MRA12" /*a_0A12 and a_0T2 are not independent*/
158  << "a_0T1phi" << "a_1T1phi" << "a_2T1phi" << "MRT1" << "a_1T2phi" << "a_2T2phi" << "MRT2"
159  << "a_0T23phi" << "a_1T23phi" << "a_2T23phi" << "MRT23"
160  << "absh_0" << "absh_p" << "absh_m" << "argh_0" << "argh_p" << "argh_m"
161  << "absh_0_1" << "absh_p_1" << "absh_m_1" << "argh_0_1" << "argh_p_1" << "argh_m_1"
162  << "absh_p_2" << "absh_m_2" << "argh_p_2" << "argh_m_2" << "xs_phi" << "SU3_breaking_abs" << "SU3_breaking_arg";
163  else if (vectorM == StandardModel::K_star || vectorM == StandardModel::K_star_P) mvllParameters = make_vector<std::string>()
164  << "a_0V" << "a_1V" << "a_2V" << "MRV" << "a_0A0" << "a_1A0" << "a_2A0" << "MRA0"
165  << "a_0A1" << "a_1A1" << "a_2A1" << "MRA1" << "a_1A12" << "a_2A12" << "MRA12" /*a_0A12 and a_0T2 are not independent*/
166  << "a_0T1" << "a_1T1" << "a_2T1" << "MRT1" << "a_1T2" << "a_2T2" << "MRT2"
167  << "a_0T23" << "a_1T23" << "a_2T23" << "MRT23"
168  << "absh_0" << "absh_p" << "absh_m" << "argh_0" << "argh_p" << "argh_m"
169  << "absh_0_1" << "absh_p_1" << "absh_m_1" << "argh_0_1" << "argh_p_1" << "argh_m_1"
170  << "absh_p_2" << "absh_m_2" << "argh_p_2" << "argh_m_2";
171 #else
172  if (vectorM == StandardModel::PHI) mvllParameters = make_vector<std::string>()
173  << "a_0Vphi" << "a_1Vphi" << "a_2Vphi" << "MRV" << "a_0A0phi" << "a_1A0phi" << "a_2A0phi" << "MRA0"
174  << "a_0A1phi" << "a_1A1phi" << "a_2A1phi" << "MRA1" << "a_1A12phi" << "a_2A12phi" << "MRA12" /*a_0A12 and a_0T2 are not independent*/
175  << "a_0T1phi" << "a_1T1phi" << "a_2T1phi" << "MRT1" << "a_1T2phi" << "a_2T2phi" << "MRT2"
176  << "a_0T23phi" << "a_1T23phi" << "a_2T23phi" << "MRT23"
177  << "reh_0" << "reh_p" << "reh_m" << "imh_0" << "imh_p" << "imh_m"
178  << "reh_0_1" << "reh_p_1" << "reh_m_1" << "imh_0_1" << "imh_p_1" << "imh_m_1"
179  << "reh_p_2" << "reh_m_2" << "imh_p_2" << "imh_m_2" << "xs_phi" << "SU3_breaking_abs" << "SU3_breaking_arg";
180  else if (vectorM == StandardModel::K_star || vectorM == StandardModel::K_star_P) mvllParameters = make_vector<std::string>()
181  << "a_0V" << "a_1V" << "a_2V" << "MRV" << "a_0A0" << "a_1A0" << "a_2A0" << "MRA0"
182  << "a_0A1" << "a_1A1" << "a_2A1" << "MRA1" << "a_1A12" << "a_2A12" << "MRA12" /*a_0A12 and a_0T2 are not independent*/
183  << "a_0T1" << "a_1T1" << "a_2T1" << "MRT1" << "a_1T2" << "a_2T2" << "MRT2"
184  << "a_0T23" << "a_1T23" << "a_2T23" << "MRT23"
185  << "reh_0" << "reh_p" << "reh_m" << "imh_0" << "imh_p" << "imh_m"
186  << "reh_0_1" << "reh_p_1" << "reh_m_1" << "imh_0_1" << "imh_p_1" << "imh_m_1"
187  << "reh_p_2" << "reh_m_2" << "imh_p_2" << "imh_m_2";
188 #endif
189  else {
190  std::stringstream out;
191  out << vectorM;
192  throw std::runtime_error("MVll: vector " + out.str() + " not implemented");
193  }
194 
195  if (dispersion) {
196  mvllParameters.clear();
197  if (vectorM == StandardModel::PHI) mvllParameters = make_vector<std::string>()
198  << "a_0Vphi" << "a_1Vphi" << "a_2Vphi" << "MRV" << "a_0A0phi" << "a_1A0phi" << "a_2A0phi" << "MRA0"
199  << "a_0A1phi" << "a_1A1phi" << "a_2A1phi" << "MRA1" << "a_1A12phi" << "a_2A12phi" << "MRA12" /*a_0A12 and a_0T2 are not independent*/
200  << "a_0T1phi" << "a_1T1phi" << "a_2T1phi" << "MRT1" << "a_1T2phi" << "a_2T2phi" << "MRT2"
201  << "a_0T23phi" << "a_1T23phi" << "a_2T23phi" << "MRT23"
202  << "r1_1" << "r2_1" << "deltaC9_1" << "phDC9_1"
203  << "r1_2" << "r2_2" << "deltaC9_2" << "phDC9_2"
204  << "r1_3" << "r2_3" << "deltaC9_3" << "phDC9_3" << "xs_phi" << "SU3_breaking_abs" << "SU3_breaking_arg";
205  else if (vectorM == StandardModel::K_star || vectorM == StandardModel::K_star_P) mvllParameters = make_vector<std::string>()
206  << "a_0V" << "a_1V" << "a_2V" << "MRV" << "a_0A0" << "a_1A0" << "a_2A0" << "MRA0"
207  << "a_0A1" << "a_1A1" << "a_2A1" << "MRA1" << "a_1A12" << "a_2A12" << "MRA12" /*a_0A12 and a_0T2 are not independent*/
208  << "a_0T1" << "a_1T1" << "a_2T1" << "MRT1" << "a_1T2" << "a_2T2" << "MRT2"
209  << "a_0T23" << "a_1T23" << "a_2T23" << "MRT23"
210  << "r1_1" << "r2_1" << "deltaC9_1" << "phDC9_1"
211  << "r1_2" << "r2_2" << "deltaC9_2" << "phDC9_2"
212  << "r1_3" << "r2_3" << "deltaC9_3" << "phDC9_3";
213  }
214 
215  if (FixedWCbtos) mvllParameters.insert(mvllParameters.end(), { "C7_SM", "C9_SM", "C10_SM" });
216 
217  mySM.initializeMeson(meson);
218  mySM.initializeMeson(vectorM);
219  return mvllParameters;
220 }
221 
222 void MVll::updateParameters()
223 {
224  if (!mySM.getFlavour().getUpdateFlag(meson, vectorM, lep)) return;
225 
226 
227  GF = mySM.getGF();
228  ale = mySM.getAle();
229  Mlep = mySM.getLeptons(lep).getMass();
230  MM = mySM.getMesons(meson).getMass();
231  MV = mySM.getMesons(vectorM).getMass();
232  mu_b = mySM.getMub();
233  mu_h = sqrt(mu_b * .5); // From Beneke Neubert
234  Mb = mySM.getQuarks(QCD::BOTTOM).getMass(); // add the PS b mass
235  Mc = mySM.getQuarks(QCD::CHARM).getMass();
236  mb_pole = mySM.Mbar2Mp(Mb); /* Conversion to pole mass*/
237  mc_pole = mySM.Mbar2Mp(Mc, FULLNLO); /* Conversion to pole mass*/
238  Ms = mySM.getQuarks(QCD::STRANGE).getMass();
239  MW = mySM.Mw();
240  lambda_t = mySM.getCKM().computelamt_s();
241  lambda_u = mySM.getCKM().computelamu_s();
242  width = mySM.getMesons(meson).computeWidth();
243  alpha_s_mub = mySM.Als(mu_b);
244  fB = mySM.getMesons(meson).getDecayconst();
245  fpara = mySM.getMesons(vectorM).getDecayconst();
246  fperp = mySM.getMesons(vectorM).getDecayconst_p();
247 
248  switch (vectorM) {
249  case StandardModel::K_star:
250  case StandardModel::K_star_P:
251  a_0V = mySM.getOptionalParameter("a_0V");
252  a_1V = mySM.getOptionalParameter("a_1V");
253  a_2V = mySM.getOptionalParameter("a_2V");
254  MRV_2 = mySM.getOptionalParameter("MRV") * mySM.getOptionalParameter("MRV");
255 
256  a_0A0 = mySM.getOptionalParameter("a_0A0");
257  a_1A0 = mySM.getOptionalParameter("a_1A0");
258  a_2A0 = mySM.getOptionalParameter("a_2A0");
259  MRA0_2 = mySM.getOptionalParameter("MRA0") * mySM.getOptionalParameter("MRA0");
260 
261  a_0A1 = mySM.getOptionalParameter("a_0A1");
262  a_1A1 = mySM.getOptionalParameter("a_1A1");
263  a_2A1 = mySM.getOptionalParameter("a_2A1");
264  MRA1_2 = mySM.getOptionalParameter("MRA1") * mySM.getOptionalParameter("MRA1");
265 
266  a_0A12 = a_0A0 * (MM * MM - MV * MV) / (8. * MM * MV);
267  a_1A12 = mySM.getOptionalParameter("a_1A12");
268  a_2A12 = mySM.getOptionalParameter("a_2A12");
269  MRA12_2 = mySM.getOptionalParameter("MRA12") * mySM.getOptionalParameter("MRA12");
270 
271  a_0T1 = mySM.getOptionalParameter("a_0T1");
272  a_1T1 = mySM.getOptionalParameter("a_1T1");
273  a_2T1 = mySM.getOptionalParameter("a_2T1");
274  MRT1_2 = mySM.getOptionalParameter("MRT1") * mySM.getOptionalParameter("MRT1");
275 
276  a_0T2 = a_0T1;
277  a_1T2 = mySM.getOptionalParameter("a_1T2");
278  a_2T2 = mySM.getOptionalParameter("a_2T2");
279  MRT2_2 = mySM.getOptionalParameter("MRT2") * mySM.getOptionalParameter("MRT2");
280 
281  a_0T23 = mySM.getOptionalParameter("a_0T23");
282  a_1T23 = mySM.getOptionalParameter("a_1T23");
283  a_2T23 = mySM.getOptionalParameter("a_2T23");
284  MRT23_2 = mySM.getOptionalParameter("MRT23") * mySM.getOptionalParameter("MRT23");
285 
286 
287  if (vectorM == StandardModel::K_star) spectator_charge = mySM.getQuarks(QCD::DOWN).getCharge();
288  else spectator_charge = mySM.getQuarks(QCD::UP).getCharge();
289 
290  etaV = -1;
291  angmomV = 1.;
292 
293  b = 1;
294 
295  SU3_breaking = 1.;
296 
297  break;
298  case StandardModel::PHI:
299  a_0V = mySM.getOptionalParameter("a_0Vphi");
300  a_1V = mySM.getOptionalParameter("a_1Vphi");
301  a_2V = mySM.getOptionalParameter("a_2Vphi");
302  MRV_2 = mySM.getOptionalParameter("MRV") * mySM.getOptionalParameter("MRV");
303 
304  a_0A0 = mySM.getOptionalParameter("a_0A0phi");
305  a_1A0 = mySM.getOptionalParameter("a_1A0phi");
306  a_2A0 = mySM.getOptionalParameter("a_2A0phi");
307  MRA0_2 = mySM.getOptionalParameter("MRA0") * mySM.getOptionalParameter("MRA0");
308 
309  a_0A1 = mySM.getOptionalParameter("a_0A1phi");
310  a_1A1 = mySM.getOptionalParameter("a_1A1phi");
311  a_2A1 = mySM.getOptionalParameter("a_2A1phi");
312  MRA1_2 = mySM.getOptionalParameter("MRA1") * mySM.getOptionalParameter("MRA1");
313 
314  a_0A12 = a_0A0 * (MM * MM - MV * MV) / (8. * MM * MV);
315  a_1A12 = mySM.getOptionalParameter("a_1A12phi");
316  a_2A12 = mySM.getOptionalParameter("a_2A12phi");
317  MRA12_2 = mySM.getOptionalParameter("MRA12") * mySM.getOptionalParameter("MRA12");
318 
319  a_0T1 = mySM.getOptionalParameter("a_0T1phi");
320  a_1T1 = mySM.getOptionalParameter("a_1T1phi");
321  a_2T1 = mySM.getOptionalParameter("a_2T1phi");
322  MRT1_2 = mySM.getOptionalParameter("MRT1") * mySM.getOptionalParameter("MRT1");
323 
324  a_0T2 = a_0T1;
325  a_1T2 = mySM.getOptionalParameter("a_1T2phi");
326  a_2T2 = mySM.getOptionalParameter("a_2T2phi");
327  MRT2_2 = mySM.getOptionalParameter("MRT2") * mySM.getOptionalParameter("MRT2");
328 
329  a_0T23 = mySM.getOptionalParameter("a_0T23phi");
330  a_1T23 = mySM.getOptionalParameter("a_1T23phi");
331  a_2T23 = mySM.getOptionalParameter("a_2T23phi");
332  MRT23_2 = mySM.getOptionalParameter("MRT23") * mySM.getOptionalParameter("MRT23");
333 
334  spectator_charge = mySM.getQuarks(QCD::STRANGE).getCharge();
335 
336  ys = mySM.getMesons(QCD::B_S).getDgamma_gamma() / 2.;
337  xs = mySM.getOptionalParameter("xs_phi");
338 
339  etaV = -1;
340  angmomV = 1.;
341 
342  b = 1.; //0.489;
343 
344  SU3_breaking = 1. + gslpp::complex(mySM.getOptionalParameter("SU3_breaking_abs"),
345  mySM.getOptionalParameter("SU3_breaking_arg"), true);
346 
347  break;
348  default:
349  std::stringstream out;
350  out << vectorM;
351  throw std::runtime_error("MVll: vector " + out.str() + " not implemented");
352  }
353 
354  if (!dispersion) {
355 #if NFPOLARBASIS_MVLL
356  h_0[0] = gslpp::complex(mySM.getOptionalParameter("absh_0"), mySM.getOptionalParameter("argh_0"), true);
357  h_0[1] = gslpp::complex(mySM.getOptionalParameter("absh_p"), mySM.getOptionalParameter("argh_p"), true);
358  h_0[2] = gslpp::complex(mySM.getOptionalParameter("absh_m"), mySM.getOptionalParameter("argh_m"), true);
359 
360  h_1[0] = gslpp::complex(mySM.getOptionalParameter("absh_0_1"), mySM.getOptionalParameter("argh_0_1"), true);
361  h_1[1] = gslpp::complex(mySM.getOptionalParameter("absh_p_1"), mySM.getOptionalParameter("argh_p_1"), true);
362  h_1[2] = gslpp::complex(mySM.getOptionalParameter("absh_m_1"), mySM.getOptionalParameter("argh_m_1"), true);
363 
364  h_2[0] = 0.;
365  h_2[1] = gslpp::complex(mySM.getOptionalParameter("absh_p_2"), mySM.getOptionalParameter("argh_p_2"), true);
366  h_2[2] = gslpp::complex(mySM.getOptionalParameter("absh_m_2"), mySM.getOptionalParameter("argh_m_2"), true);
367 #else
368  h_0[0] = gslpp::complex(mySM.getOptionalParameter("reh_0"), mySM.getOptionalParameter("imh_0"), false);
369  h_0[1] = gslpp::complex(mySM.getOptionalParameter("reh_p"), mySM.getOptionalParameter("imh_p"), false);
370  h_0[2] = gslpp::complex(mySM.getOptionalParameter("reh_m"), mySM.getOptionalParameter("imh_m"), false);
371 
372  h_1[0] = gslpp::complex(mySM.getOptionalParameter("reh_0_1"), mySM.getOptionalParameter("imh_0_1"), false);
373  h_1[1] = gslpp::complex(mySM.getOptionalParameter("reh_p_1"), mySM.getOptionalParameter("imh_p_1"), false);
374  h_1[2] = gslpp::complex(mySM.getOptionalParameter("reh_m_1"), mySM.getOptionalParameter("imh_m_1"), false);
375 
376  h_2[0] = 0.;
377  h_2[1] = gslpp::complex(mySM.getOptionalParameter("reh_p_2"), mySM.getOptionalParameter("imh_p_2"), false);
378  h_2[2] = gslpp::complex(mySM.getOptionalParameter("reh_m_2"), mySM.getOptionalParameter("imh_m_2"), false);
379 #endif
380  } else {
381  h_0[0] = gslpp::complex(mySM.getOptionalParameter("r1_1"));
382  h_0[1] = gslpp::complex(mySM.getOptionalParameter("r1_2"));
383  h_0[2] = gslpp::complex(mySM.getOptionalParameter("r1_3"));
384 
385  h_1[0] = gslpp::complex(mySM.getOptionalParameter("r2_1"));
386  h_1[1] = gslpp::complex(mySM.getOptionalParameter("r2_2"));
387  h_1[2] = gslpp::complex(mySM.getOptionalParameter("r2_3"));
388 
389  h_2[0] = gslpp::complex(mySM.getOptionalParameter("deltaC9_1"));
390  h_2[1] = gslpp::complex(mySM.getOptionalParameter("deltaC9_2"));
391  h_2[2] = gslpp::complex(mySM.getOptionalParameter("deltaC9_3"));
392  exp_Phase[0] = exp(gslpp::complex::i() * mySM.getOptionalParameter("phDC9_1"));
393  exp_Phase[1] = exp(gslpp::complex::i() * mySM.getOptionalParameter("phDC9_2"));
394  exp_Phase[2] = exp(gslpp::complex::i() * mySM.getOptionalParameter("phDC9_3"));
395  }
396 
397  allcoeff = mySM.getFlavour().ComputeCoeffBMll(mu_b, lep); //check the mass scale, scheme fixed to NDR
398  allcoeffprime = mySM.getFlavour().ComputeCoeffprimeBMll(mu_b, lep); //check the mass scale, scheme fixed to NDR
399 
400  C_1 = ((*(allcoeff[LO]))(0) + (*(allcoeff[NLO]))(0));
401  C_1L_bar = (*(allcoeff[LO]))(0) / 2.;
402  C_2 = ((*(allcoeff[LO]))(1) + (*(allcoeff[NLO]))(1));
403  C_2L_bar = (*(allcoeff[LO]))(1) - (*(allcoeff[LO]))(0) / 6.;
404  C_3 = ((*(allcoeff[LO]))(2) + (*(allcoeff[NLO]))(2));
405  C_4 = ((*(allcoeff[LO]))(3) + (*(allcoeff[NLO]))(3));
406  C_5 = ((*(allcoeff[LO]))(4) + (*(allcoeff[NLO]))(4));
407  C_6 = ((*(allcoeff[LO]))(5) + (*(allcoeff[NLO]))(5));
408  C_8 = ((*(allcoeff[LO]))(7) + (*(allcoeff[NLO]))(7));
409  C_8L = (*(allcoeff[LO]))(7);
410  C_S = MW / Mb * (((*(allcoeff[LO]))(10) + (*(allcoeff[NLO]))(10)));
411  C_P = MW / Mb * (((*(allcoeff[LO]))(11) + (*(allcoeff[NLO]))(11)));
412  C_9p = (*(allcoeffprime[LO]))(8) + (*(allcoeffprime[NLO]))(8);
413  C_10p = (*(allcoeffprime[LO]))(9) + (*(allcoeffprime[NLO]))(9);
414  C_Sp = MW / Mb * ((*(allcoeffprime[LO]))(10) + (*(allcoeffprime[NLO]))(10));
415  C_Pp = MW / Mb * ((*(allcoeffprime[LO]))(11) + (*(allcoeffprime[NLO]))(11));
416 
417  if (FixedWCbtos) {
418  allcoeff_noSM = mySM.getFlavour().ComputeCoeffBMll(mu_b, lep, true); //check the mass scale, scheme fixed to NDR
419  C_7 = mySM.getOptionalParameter("C7_SM") + ((*(allcoeff_noSM[LO]))(6) + (*(allcoeff_noSM[NLO]))(6));
420  C_9 = mySM.getOptionalParameter("C9_SM") + ((*(allcoeff_noSM[LO]))(8) + (*(allcoeff_noSM[NLO]))(8));
421  C_10 = mySM.getOptionalParameter("C10_SM") + ((*(allcoeff_noSM[LO]))(9) + (*(allcoeff_noSM[NLO]))(9));
422  } else {
423  C_7 = ((*(allcoeff[LO]))(6) + (*(allcoeff[NLO]))(6));
424  C_9 = ((*(allcoeff[LO]))(8) + (*(allcoeff[NLO]))(8));
425  C_10 = ((*(allcoeff[LO]))(9) + (*(allcoeff[NLO]))(9));
426  }
427  C_7p = MsoMb * ((*(allcoeffprime[LO]))(6) + (*(allcoeffprime[NLO]))(6));
428  C_7p -= MsoMb * (C_7 + 1. / 3. * C_3 + 4 / 9 * C_4 + 20. / 3. * C_5 + 80. / 9. * C_6);
429 
430  allcoeffh = mySM.getFlavour().ComputeCoeffBMll(mu_h, lep); //check the mass scale, scheme fixed to NDR
431 
432  C_1Lh_bar = (*(allcoeffh[LO]))(0) / 2.;
433  C_2Lh_bar = (*(allcoeffh[LO]))(1) - (*(allcoeff[LO]))(0) / 6.;
434  C_8Lh = (*(allcoeffh[LO]))(7);
435 
436  checkCache();
437 
438  t_p = pow(MM + MV, 2.);
439  t_m = pow(MM - MV, 2.);
440  t_0 = t_p * (1. - sqrt(1. - t_m / t_p)); /*Modify it for Lattice*/
441  z_0 = (sqrt(t_p) - sqrt(t_p - t_0)) / (sqrt(t_p) + sqrt(t_p - t_0));
442  MMpMV = MM + MV;
443  MMpMV2 = MMpMV * MMpMV;
444  MMmMV = MM - MV;
445  MMmMV2 = MMmMV * MMmMV;
446  MM2 = MM*MM;
447  MM4 = MM2*MM2;
448  MV2 = MV*MV;
449  MV4 = MV2*MV2;
450  MMMV = MM*MV;
451  MM2mMV2 = MM2 - MV2;
452  MM2pMV2 = MM2 + MV2;
453  fourMV = 4. * MV;
454  twoMM2 = 2. * MM2;
455  twoMV2 = 2. * MV2;
456  onepMMoMV = (1. + MV / MM);
457  MM_MMpMV = MM * MMpMV;
458  twoMM_mbpms = 2. * MM * (Mb + Ms);
459  fourMM2 = 4. * MM2;
460  Mlep2 = Mlep*Mlep;
461  twoMlepMb = 2. * Mlep*Mb;
462  MboMW = Mb / MW;
463  MsoMb = Ms / Mb;
464  ninetysixM_PI3MM3 = 96. * M_PI * M_PI * M_PI * MM * MM*MM;
465  sixteenM_PI2 = 16. * M_PI2;
466  sixteenM_PI2MM2 = sixteenM_PI2 * MM*MM;
467  twoMboMM = 2 * Mb / MM;
468  H_0_pre = 8. / 27. + 4. / 9. * gslpp::complex::i() * M_PI;
469  H_0_WC = (C_3 + 4. / 3. * C_4 + 16. * C_5 + 64. / 3. * C_6);
470  H_c_WC = (4. / 3. * C_1 + C_2 + 6. * C_3 + 60. * C_5);
471  H_b_WC = (7. * C_3 + 4. / 3. * C_4 + 76. * C_5 + 64. / 3. * C_6);
472  mu_b2 = mu_b*mu_b;
473  Mc2 = Mc*Mc;
474  Mb2 = Mb*Mb;
475  fourMc2 = 4. * Mc2;
476  fourMb2 = 4. * Mb2;
477  logMc = log(Mc2 / mu_b2);
478  logMb = log(Mb2 / mu_b2);
479  fournineth = 4. / 9.;
480  half = 1. / 2.;
481  twothird = 2. / 3.;
482  ihalfMPI = gslpp::complex::i() * M_PI / 2.;
483  twoMM3 = 2. * MM2 * MM;
484  C2_inv = 1. / (2. * C_2.real());
485  gtilde_1_pre = -16. * pow(MM, 3.)*(MM + MV) * pow(M_PI, 2.);
486  gtilde_2_pre = -16. * pow(MM, 3.) * pow(M_PI, 2.) / MMpMV;
487  gtilde_3_pre = 64. * pow(MM, 3.) * pow(M_PI, 2.) * MV*MMpMV;
488  S_L_pre = (-2. * MM * (Mb + Ms));
489 
490  M_PI2osix = M_PI2 / 6.;
491  twoMM = 2. * MM;
492 
493  N_QCDF = M_PI2 / 3. * fB * fperp / MM;
494 
495  deltaT_0 = alpha_s_mub * CF / 4. / M_PI;
496  deltaT_1par = mySM.Als(mu_h) * CF / 4. * M_PI / 3. * mySM.getMesons(meson).getDecayconst() *
497  mySM.getMesons(vectorM).getDecayconst() / MM;
498  deltaT_1perp = mySM.Als(mu_h) * CF / 4. * M_PI / 3. * mySM.getMesons(meson).getDecayconst() *
499  mySM.getMesons(vectorM).getDecayconst_p() / MM;
500 
501  F87_0 = -32. / 9. * log(mu_b / Mb) + 8. / 27. * M_PI2 - 44. / 9. - 8. / 9. * gslpp::complex::i() * M_PI;
502 
503  NN = -(4. * GF * MM * ale * lambda_t) / (sqrt(2.)*4. * M_PI);
504  NN_conjugate = -(4. * GF * MM * ale * lambda_t.conjugate()) / (sqrt(2.)*4. * M_PI);
505 
506  std::map<std::pair<double, double>, unsigned int >::iterator it;
507 
508  if (I0_updated == 0) for (it = sigma0Cached.begin(); it != sigma0Cached.end(); ++it) it->second = 0;
509  if (I1_updated == 0) for (it = sigma1Cached.begin(); it != sigma1Cached.end(); ++it) it->second = 0;
510  if (I2_updated == 0) for (it = sigma2Cached.begin(); it != sigma2Cached.end(); ++it) it->second = 0;
511  if (I3_updated == 0) for (it = sigma3Cached.begin(); it != sigma3Cached.end(); ++it) it->second = 0;
512  if (I4_updated == 0) for (it = sigma4Cached.begin(); it != sigma4Cached.end(); ++it) it->second = 0;
513  if (I5_updated == 0) for (it = sigma5Cached.begin(); it != sigma5Cached.end(); ++it) it->second = 0;
514  if (I6_updated == 0) for (it = sigma6Cached.begin(); it != sigma6Cached.end(); ++it) it->second = 0;
515  if (I7_updated == 0) for (it = sigma7Cached.begin(); it != sigma7Cached.end(); ++it) it->second = 0;
516  if (I9_updated == 0) for (it = sigma9Cached.begin(); it != sigma9Cached.end(); ++it) it->second = 0;
517  if (I10_updated == 0) for (it = sigma10Cached.begin(); it != sigma10Cached.end(); ++it) it->second = 0;
518  if (I11_updated == 0) for (it = sigma11Cached.begin(); it != sigma11Cached.end(); ++it) it->second = 0;
519 
520  if (I0_updated == 0) for (it = delta0Cached.begin(); it != delta0Cached.end(); ++it) it->second = 0;
521  if (I1_updated == 0) for (it = delta1Cached.begin(); it != delta1Cached.end(); ++it) it->second = 0;
522  if (I2_updated == 0) for (it = delta2Cached.begin(); it != delta2Cached.end(); ++it) it->second = 0;
523  if (I3_updated == 0) for (it = delta3Cached.begin(); it != delta3Cached.end(); ++it) it->second = 0;
524  if (I11_updated == 0) for (it = delta11Cached.begin(); it != delta11Cached.end(); ++it) it->second = 0;
525 
526  std::map<double, unsigned int >::iterator iti;
527  if (deltaTparpupdated == 0) for (iti = deltaTparpCached.begin(); iti != deltaTparpCached.end(); ++iti) iti->second = 0;
528  if (deltaTparmupdated == 0) for (iti = deltaTparmCached.begin(); iti != deltaTparmCached.end(); ++iti) iti->second = 0;
529  if (deltaTperpupdated == 0) for (iti = deltaTparpCached.begin(); iti != deltaTparpCached.end(); ++iti) iti->second = 0;
530 
531  if (deltaTparpupdated * deltaTparmupdated == 0) for (it = I1Cached.begin(); it != I1Cached.end(); ++it) it->second = 0;
532 
533 #if SPLINE
534  if (mySM.getFlavour().getUpdateFlag(meson, vectorM, lep)) spline_QCDF_func();
535 #else
536  if (mySM.getFlavour().getUpdateFlag(meson, vectorM, lep)) fit_QCDF_func();
537 #endif
538 
539  mySM.getFlavour().setUpdateFlag(meson, vectorM, lep, false);
540 
541  return;
542 }
543 
544 void MVll::checkCache()
545 {
546 
547  if (MM == k2_cache(0) && MV == k2_cache(1)) {
548  k2_updated = 1;
549  z_updated = 1;
550  } else {
551  k2_updated = 0;
552  z_updated = 0;
553  k2_cache(0) = MM;
554  k2_cache(1) = MV;
555  }
556 
557  if (Mlep == beta_cache) {
558  beta_updated = 1;
559  } else {
560  beta_updated = 0;
561  beta_cache = Mlep;
562  }
563 
564  lambda_updated = k2_updated;
565  F_updated = lambda_updated * beta_updated;
566 
567  if (GF == N_cache(0) && ale == N_cache(1) && MM == N_cache(2) && lambda_t == Nc_cache) {
568  N_updated = 1;
569  } else {
570  N_updated = 0;
571  N_cache(0) = GF;
572  N_cache(1) = ale;
573  N_cache(2) = MM;
574  Nc_cache = lambda_t;
575  }
576 
577  if (a_0V == V_cache(0) && a_1V == V_cache(1) && a_2V == V_cache(2)) {
578  V_updated = V_updated * z_updated;
579  } else {
580  V_updated = 0;
581  V_cache(0) = a_0V;
582  V_cache(1) = a_1V;
583  V_cache(2) = a_2V;
584  }
585 
586  if (a_0A0 == A0_cache(0) && a_1A0 == A0_cache(1) && a_2A0 == A0_cache(2)) {
587  A0_updated = A0_updated * z_updated;
588  } else {
589  A0_updated = 0;
590  A0_cache(0) = a_0A0;
591  A0_cache(1) = a_1A0;
592  A0_cache(2) = a_2A0;
593  }
594 
595  if (a_0A1 == A1_cache(0) && a_1A1 == A1_cache(1) && a_2A1 == A1_cache(2)) {
596  A1_updated = A1_updated * z_updated;
597  } else {
598  A1_updated = 0;
599  A1_cache(0) = a_0A1;
600  A1_cache(1) = a_1A1;
601  A1_cache(2) = a_2A1;
602  }
603 
604  if (a_0T1 == T1_cache(0) && a_1T1 == T1_cache(1) && a_2T1 == T1_cache(2)) {
605  T1_updated = T1_updated * z_updated;
606  } else {
607  T1_updated = 0;
608  T1_cache(0) = a_0T1;
609  T1_cache(1) = a_1T1;
610  T1_cache(2) = a_2T1;
611  }
612 
613  if (a_0T2 == T2_cache(0) && a_1T2 == T2_cache(1) && a_2T2 == T2_cache(2)) {
614  T2_updated = T2_updated * z_updated;
615  } else {
616  T2_updated = 0;
617  T2_cache(0) = a_0T2;
618  T2_cache(1) = a_1T2;
619  T2_cache(2) = a_2T2;
620  }
621 
622  VL1_updated = k2_updated * lambda_updated * A1_updated * V_updated;
623  VL2_updated = VL1_updated;
624 
625  TL1_updated = k2_updated * lambda_updated * T1_updated * T2_updated;
626  TL2_updated = TL1_updated;
627 
628  VR1_updated = VL2_updated;
629  VR2_updated = VL1_updated;
630 
631  TR1_updated = TL2_updated;
632  TR2_updated = TL1_updated;
633 
634  if (Mb == SL_cache(0) && Ms == SL_cache(1)) {
635  Mb_Ms_updated = 1;
636  SL_updated = lambda_updated * A0_updated;
637  SR_updated = SL_updated;
638  } else {
639  Mb_Ms_updated = 0;
640  SL_updated = 0;
641  SR_updated = SL_updated;
642  SL_cache(0) = Mb;
643  SL_cache(1) = Ms;
644  }
645 
646  if (a_0A12 == VL0_cache(0) && a_1A12 == VL0_cache(1) && a_2A12 == VL0_cache(2)) {
647  VL0_updated = VL0_updated * z_updated;
648  VR0_updated = VL0_updated;
649  } else {
650  VL0_updated = 0;
651  VR0_updated = VL0_updated;
652  VL0_cache(0) = a_0A12;
653  VL0_cache(1) = a_1A12;
654  VL0_cache(2) = a_2A12;
655  }
656 
657  if (a_0T23 == TL0_cache(0) && a_1T23 == TL0_cache(1) && a_2T23 == TL0_cache(2)) {
658  TL0_updated = TL0_updated * z_updated;
659  TR0_updated = TL0_updated;
660  } else {
661  TL0_updated = 0;
662  TR0_updated = TL0_updated;
663  TL0_cache(0) = a_0T23;
664  TL0_cache(1) = a_1T23;
665  TL0_cache(2) = a_2T23;
666  }
667 
668 
669  if (C_1 == C_1_cache) {
670  C_1_updated = 1;
671  } else {
672  C_1_updated = 0;
673  C_1_cache = C_1;
674  }
675 
676  if (C_2 == C_2_cache) {
677  C_2_updated = 1;
678  } else {
679  C_2_updated = 0;
680  C_2_cache = C_2;
681  }
682 
683  if (C_3 == C_3_cache) {
684  C_3_updated = 1;
685  } else {
686  C_3_updated = 0;
687  C_3_cache = C_3;
688  }
689 
690  if (C_4 == C_4_cache) {
691  C_4_updated = 1;
692  } else {
693  C_4_updated = 0;
694  C_4_cache = C_4;
695  }
696 
697  if (C_5 == C_5_cache) {
698  C_5_updated = 1;
699  } else {
700  C_5_updated = 0;
701  C_5_cache = C_5;
702  }
703 
704  if (C_6 == C_6_cache) {
705  C_6_updated = 1;
706  } else {
707  C_6_updated = 0;
708  C_6_cache = C_6;
709  }
710 
711  if (C_7 == C_7_cache) {
712  C_7_updated = 1;
713  } else {
714  C_7_updated = 0;
715  C_7_cache = C_7;
716  }
717 
718  if (C_9 == C_9_cache) {
719  C_9_updated = 1;
720  } else {
721  C_9_updated = 0;
722  C_9_cache = C_9;
723  }
724 
725  if (C_10 == C_10_cache) {
726  C_10_updated = 1;
727  } else {
728  C_10_updated = 0;
729  C_10_cache = C_10;
730  }
731 
732  if (C_S == C_S_cache) {
733  C_S_updated = 1;
734  } else {
735  C_S_updated = 0;
736  C_S_cache = C_S;
737  }
738 
739  if (C_P == C_P_cache) {
740  C_P_updated = 1;
741  } else {
742  C_P_updated = 0;
743  C_P_cache = C_P;
744  }
745 
746  if (C_7p == C_7p_cache) {
747  C_7p_updated = 1;
748  } else {
749  C_7p_updated = 0;
750  C_7p_cache = C_7p;
751  }
752 
753  if (C_9p == C_9p_cache) {
754  C_9p_updated = 1;
755  } else {
756  C_9p_updated = 0;
757  C_9p_cache = C_9p;
758  }
759 
760  if (C_10p == C_10p_cache) {
761  C_10p_updated = 1;
762  } else {
763  C_10p_updated = 0;
764  C_10p_cache = C_10p;
765  }
766 
767  if (C_Sp == C_Sp_cache) {
768  C_Sp_updated = 1;
769  } else {
770  C_Sp_updated = 0;
771  C_Sp_cache = C_Sp;
772  }
773 
774  if (C_Pp == C_Pp_cache) {
775  C_Pp_updated = 1;
776  } else {
777  C_Pp_updated = 0;
778  C_Pp_cache = C_Pp;
779  }
780 
781  if (C_2Lh_bar == C_2Lh_cache) {
782  C_2Lh_updated = 1;
783  } else {
784  C_2Lh_updated = 0;
785  C_2Lh_cache = C_2Lh_bar;
786  }
787 
788  if (C_8Lh == C_8Lh_cache) {
789  C_8Lh_updated = 1;
790  } else {
791  C_8Lh_updated = 0;
792  C_8Lh_cache = C_8Lh;
793  }
794 
795  if (Mb == Ycache(0) && Mc == Ycache(1)) {
796  Yupdated = C_1_updated * C_2_updated * C_3_updated * C_4_updated * C_5_updated * C_6_updated;
797  } else {
798  Yupdated = 0;
799  Ycache(0) = Mb;
800  Ycache(1) = Mc;
801  }
802 
803  if (h_0[0] == h0Ccache[0] && h_1[0] == h0Ccache[1] && h_2[0] == h0Ccache[2] && SU3_breaking == h0Ccache[3]) {
804  h0_updated = 1;
805  } else {
806  h0_updated = 0;
807  h0Ccache[0] = h_0[0];
808  h0Ccache[1] = h_1[0];
809  h0Ccache[2] = h_2[0];
810  h0Ccache[3] = SU3_breaking;
811  }
812 
813  if (h_0[1] == h1Ccache[0] && h_1[1] == h1Ccache[1] && h_2[1] == h1Ccache[2] && SU3_breaking == h1Ccache[3]) {
814  h1_updated = 1;
815  } else {
816  h1_updated = 0;
817  h1Ccache[0] = h_0[1];
818  h1Ccache[1] = h_1[1];
819  h1Ccache[2] = h_2[1];
820  h1Ccache[3] = SU3_breaking;
821  }
822 
823  if (h_0[2] == h2Ccache[0] && h_1[2] == h2Ccache[1] && h_2[2] == h2Ccache[2] && SU3_breaking == h2Ccache[3]) {
824  h2_updated = 1;
825  } else {
826  h2_updated = 0;
827  h2Ccache[0] = h_0[2];
828  h2Ccache[1] = h_1[2];
829  h2Ccache[2] = h_2[2];
830  h2Ccache[3] = SU3_breaking;
831  }
832 
833  if (MM == H_V0cache(0) && Mb == H_V0cache(1)) {
834  H_V0updated = N_updated * C_9_updated * Yupdated * VL0_updated * C_9p_updated * VR0_updated * C_7_updated * TL0_updated * C_7p_updated * TR0_updated * h0_updated;
835  } else {
836  H_V0updated = 0;
837  H_V0cache(0) = MM;
838  H_V0cache(1) = Mb;
839  }
840 
841  if (MM == H_V1cache(0) && Mb == H_V1cache(1)) {
842  H_V1updated = N_updated * C_9_updated * Yupdated * VL1_updated * C_9p_updated * VR1_updated * C_7_updated * TL1_updated * C_7p_updated * TR1_updated * h1_updated;
843  } else {
844  H_V1updated = 0;
845  H_V1cache(0) = MM;
846  H_V1cache(1) = Mb;
847  }
848 
849  if (MM == H_V2cache(0) && Mb == H_V2cache(1)) {
850  H_V2updated = N_updated * C_9_updated * Yupdated * VL2_updated * C_9p_updated * VR2_updated * C_7_updated * TL2_updated * C_7p_updated * TR2_updated * h2_updated;
851  } else {
852  H_V2updated = 0;
853  H_V2cache(0) = MM;
854  H_V2cache(1) = Mb;
855  }
856 
857  H_A0updated = N_updated * C_10_updated * VL0_updated * C_10p_updated * VR0_updated;
858  H_A1updated = N_updated * C_10_updated * VL1_updated * C_10p_updated * VR1_updated;
859  H_A2updated = N_updated * C_10_updated * VL2_updated * C_10p_updated * VR2_updated;
860 
861  if (Mb == H_Scache(0) && MW == H_Scache(1)) {
862  H_Supdated = N_updated * C_S_updated * SL_updated * C_Sp_updated * SR_updated;
863  } else {
864  H_Supdated = 0;
865  H_Scache(0) = Mb;
866  H_Scache(1) = MW;
867  }
868 
869  if (Mb == H_Pcache(0) && MW == H_Pcache(1) && Mlep == H_Pcache(2) && Ms == H_Pcache(3)) {
870  H_Pupdated = N_updated * C_P_updated * SL_updated * C_Pp_updated * SR_updated * C_10_updated * C_10p_updated;
871  } else {
872  H_Pupdated = 0;
873  H_Pcache(0) = Mb;
874  H_Pcache(1) = MW;
875  H_Pcache(2) = Mlep;
876  H_Pcache(3) = Ms;
877 
878  }
879 
880  if (MM == T_cache(0) && Mb == T_cache(1) && Mc == T_cache(2) &&
881  mySM.getMesons(vectorM).getGegenalpha(0) == T_cache(3) && mySM.getMesons(vectorM).getGegenalpha(1) == T_cache(4)) {
882  T_updated = 1;
883  } else {
884  T_updated = 0;
885  T_cache(0) = MM;
886  T_cache(1) = Mb;
887  T_cache(2) = Mc;
888  T_cache(3) = mySM.getMesons(vectorM).getGegenalpha(0);
889  T_cache(4) = mySM.getMesons(vectorM).getGegenalpha(1);
890  }
891 
892  deltaTparpupdated = C_2Lh_updated * T_updated;
893  deltaTparmupdated = C_2Lh_updated * C_8Lh_updated * T_updated;
894  deltaTperpupdated = deltaTparpupdated;
895 
896  I0_updated = F_updated * H_V0updated * H_A0updated * H_Pupdated * beta_updated * H_Supdated * deltaTparmupdated;
897  I1_updated = F_updated * beta_updated * H_V1updated * H_V2updated * H_A1updated * H_A2updated * deltaTparmupdated;
898  I2_updated = F_updated * beta_updated * H_V0updated * H_A0updated * deltaTparmupdated;
899  I3_updated = F_updated * H_V1updated * H_V2updated * H_A1updated * H_A2updated * beta_updated * deltaTparmupdated;
900  I4_updated = F_updated * H_V1updated * H_V2updated * H_A1updated * H_A2updated * deltaTparmupdated;
901  I5_updated = F_updated * H_V0updated * H_V1updated * H_V2updated * H_A0updated * H_A1updated * H_A2updated * beta_updated * deltaTparmupdated;
902  I6_updated = F_updated * H_V1updated * H_V2updated * H_A0updated * H_A1updated * H_A2updated * H_V0updated * beta_updated * H_Supdated * deltaTparmupdated;
903  I7_updated = I4_updated * beta_updated;
904  I8_updated = F_updated * beta_updated * H_Supdated * H_V0updated * deltaTparmupdated;
905  I9_updated = I6_updated;
906  I10_updated = I5_updated;
907  I11_updated = I7_updated;
908 
909 }
910 
911 /*******************************************************************************
912  * Transverse Form Factors *
913  * ****************************************************************************/
914 
915 double MVll::FF_fit(double q2, double a_0, double a_1, double a_2, double MR_2)
916 {
917  return 1. / (1. - q2 / MR_2) * (a_0 + a_1 * (z(q2) - z_0) + a_2 * (z(q2) - z_0) * (z(q2) - z_0));
918 }
919 
920 double MVll::z(double q2)
921 {
922  return ( sqrt(t_p - q2) - sqrt(t_p - t_0)) / (sqrt(t_p - q2) + sqrt(t_p - t_0));
923 }
924 
925 double MVll::V(double q2)
926 {
927  return FF_fit(q2, a_0V, a_1V, a_2V, MRV_2);
928 }
929 
930 double MVll::A_0(double q2)
931 {
932  return FF_fit(q2, a_0A0, a_1A0, a_2A0, MRA0_2);
933 }
934 
935 double MVll::A_1(double q2)
936 {
937  return FF_fit(q2, a_0A1, a_1A1, a_2A1, MRA1_2);
938 }
939 
940 double MVll::A_2(double q2)
941 {
942  return (MMpMV2 * (MM2mMV2 - q2) * A_1(q2) - 16. * MM * MV2 * MMpMV * FF_fit(q2, a_0A12, a_1A12, a_2A12, MRA12_2)) / lambda(q2);
943 }
944 
945 double MVll::T_1(double q2)
946 {
947  return FF_fit(q2, a_0T1, a_1T1, a_2T1, MRT1_2);
948 }
949 
950 double MVll::T_2(double q2)
951 {
952  return FF_fit(q2, a_0T2, a_1T2, a_2T2, MRT2_2);
953 }
954 
955 double MVll::V_0t(double q2)
956 {
957  return fourMV / sqrt(q2) * FF_fit(q2, a_0A12, a_1A12, a_2A12, MRA12_2);
958 }
959 
960 double MVll::V_p(double q2)
961 {
962  return half * (onepMMoMV * A_1(q2) - sqrt(lambda(q2)) / (MM_MMpMV) * V(q2));
963 }
964 
965 double MVll::V_m(double q2)
966 {
967  return half * (onepMMoMV * A_1(q2) + sqrt(lambda(q2)) / (MM_MMpMV) * V(q2));
968 }
969 
970 double MVll::T_0t(double q2)
971 {
972  return 2 * sqrt(q2) * MV / MM_MMpMV * FF_fit(q2, a_0T23, a_1T23, a_2T23, MRT23_2);
973 }
974 
975 double MVll::T_p(double q2)
976 {
977  return (MM2mMV2 * T_2(q2) - sqrt(lambda(q2)) * T_1(q2)) / twoMM2;
978 }
979 
980 double MVll::T_m(double q2)
981 {
982  return (MM2mMV2 * T_2(q2) + sqrt(lambda(q2)) * T_1(q2)) / twoMM2;
983 }
984 
985 double MVll::S_L(double q2)
986 {
987  return -sqrt(lambda(q2)) / twoMM_mbpms * A_0(q2);
988 }
989 
990 /*******************************************************************************
991  * QCDF NLO *
992  * ****************************************************************************/
993 
994 gslpp::complex MVll::A_Seidel(double q2, double mb2)
995 {
996  double sh = q2 / mb2;
997  double z = (4. * mb2) / q2;
998  double lsh = log(sh);
999  gslpp::complex acsq = arccot((gslpp::complex)sqrt(z - 1.));
1000  double sh2 = sh*sh;
1001  double osh2 = (1. - sh)*(1. - sh);
1002  return (-(104.) / (243.) * log((mb2) / (mu_b2)) + (4. * sh) / (27. * (1. - sh)) * (dilog((gslpp::complex)sh) + lsh * log(1. - sh))
1003  + (1.) / (729. * osh2) * (6. * sh * (29. - 47. * sh) * lsh + 785. - 1600. * sh + 833. * sh * sh + 6. * M_PI * gslpp::complex::i() * (20. - 49. * sh + 47. * sh2))
1004  - (2.) / (243. * osh2 * (1. - sh)) * (2. * sqrt(z - 1.) * (-4. + 9. * sh - 15. * sh2 + 4. * sh2 * sh) * acsq + 9. * sh2 * sh * lsh * lsh + 18. * M_PI * gslpp::complex::i() * sh * (1. - 2. * sh) * lsh)
1005  + (2. * sh) / (243. * osh2 * osh2) * (36. * acsq * acsq + M_PI2 * (-4. + 9. * sh - 9. * sh2 + 3. * sh2 * sh)));
1006 }
1007 
1008 gslpp::complex MVll::B_Seidel(double q2, double mb2)
1009 {
1010  double sh = q2 / mb2;
1011  double z = (4. * mb2) / q2;
1012  double sqrt_z_m_1 = sqrt(z - 1.);
1013  gslpp::complex x1 = 0.5 + gslpp::complex::i() / 2. * sqrt_z_m_1;
1014  gslpp::complex x2 = 0.5 - gslpp::complex::i() / 2. * sqrt_z_m_1;
1015  gslpp::complex x3 = 0.5 + gslpp::complex::i() / (2. * sqrt_z_m_1);
1016  gslpp::complex x4 = 0.5 - gslpp::complex::i() / (2. * sqrt_z_m_1);
1017  gslpp::complex lx1 = log(x1);
1018  gslpp::complex lx2 = log(x2);
1019  gslpp::complex lx3 = log(x3);
1020  gslpp::complex lx4 = log(x4);
1021  gslpp::complex lx2_x1 = lx2 - lx1;
1022  gslpp::complex lzm1 = log(z - 1.);
1023  gslpp::complex acsq = arccot((gslpp::complex)sqrt_z_m_1);
1024  double sh2 = sh*sh;
1025  double lsh = log(sh);
1026  double osh2 = (1. - sh)*(1. - sh);
1027  double lmb_mu = log(mb2 / mu_b2);
1028  return (8. / (243. * sh) * ((4. - 34. * sh - 17. * M_PI * gslpp::complex::i() * sh) * lmb_mu + 8. * sh * lmb_mu * lmb_mu + 17. * sh * lsh * lmb_mu)
1029  + ((2. + sh) * sqrt_z_m_1) / (729. * sh) * (-48. * lmb_mu * acsq - 18. * M_PI * log(z - 1.) + 3. * gslpp::complex::i() * lzm1 * lzm1
1030  - 24. * gslpp::complex::i() * dilog(-x2 / x1) - 5. * M_PI2 * gslpp::complex::i()
1031  + 6. * gslpp::complex::i() * (-9. * lx1 * lx1 + lx2 * lx2 - 2. * lx4 * lx4 + 6. * lx1 * lx2 - 4. * lx1 * lx3 + 8. * lx1 * lx4)
1032  - 12. * M_PI * (2. * lx1 + lx3 + lx4)) - 2. / (243. * sh * (1 - sh)) * (4. * sh * (-8. + 17. * sh) * (dilog((gslpp::complex)sh) + lsh * log(1. - sh))
1033  + 3. * (2. + sh) * (3. - sh) * lx2_x1 * lx2_x1 + 12. * M_PI * (-6. - sh + sh2) * acsq) + 2. / (2187. * sh * osh2) * (-18. * sh * (120. - 211. * sh + 73. * sh2) * lsh
1034  - 288. - 8. * sh + 934. * sh2 - 692. * sh2 * sh + 18. * M_PI * gslpp::complex::i() * sh * (82. - 173. * sh + 73. * sh2))
1035  - 4. / (243. * sh * osh2 * (1 - sh)) * (-2. * sqrt_z_m_1 * (4. - 3. * sh - 18. * sh2 + 16. * sh2 * sh - 5. * sh2 * sh2) * acsq - 9. * sh * sh2 * lsh * lsh
1036  + 2. * M_PI * gslpp::complex::i() * sh * (8. - 33. * sh + 51. * sh2 - 17. * sh * sh2) * lsh) + 2. / (729. * sh * osh2 * osh2) * (72. * (3. - 8. * sh + 2. * sh2) * acsq * acsq
1037  - M_PI2 * (54. - 53. * sh - 286. * sh2 + 612. * sh * sh2 - 446. * sh2 * sh2 + 113. * sh2 * sh2 * sh)));
1038 }
1039 
1040 gslpp::complex MVll::C_Seidel(double q2)
1041 {
1042  return -(16.) / (81.) * log((q2) / (mu_b2)) + (428.) / (243.) - (64.) / (27.) * gsl_sf_zeta_int(3) + (16.) / (81.) * M_PI * gslpp::complex::i();
1043  /* gsl_sf_zeta_int returns a double */
1044 }
1045 
1046 gslpp::complex MVll::deltaC7_QCDF(double q2, bool conjugate, bool spline)
1047 {
1048 #if COMPUTECP && SPLINE
1049  if (spline && !conjugate) return gsl_spline_eval(spline_Re_deltaC7_QCDF, q2, acc_Re_deltaC7_QCDF);
1050  else if (spline && conjugate) return gsl_spline_eval(spline_Re_deltaC7_QCDF_conj, q2, acc_Re_deltaC7_QCDF_conj);
1051 #elif SPLINE
1052  if (spline) return gsl_spline_eval(spline_Re_deltaC7_QCDF, q2, acc_Re_deltaC7_QCDF);
1053 #endif
1054 
1055  double muh = mu_b / mb_pole;
1056  double z = mc_pole * mc_pole / mb_pole / mb_pole;
1057  double sh = q2 / mb_pole / mb_pole;
1058  double sh2 = sh*sh;
1059 
1060 #if FULLNLOQCDF_MVLL
1061  gslpp::complex A_Sdl = A_Seidel(q2, mb_pole*mb_pole); /* hep-ph/0403185v2.*/
1062  gslpp::complex Fu_17 = -A_Sdl; /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1063  gslpp::complex Fu_27 = 6. * A_Sdl; /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1064 #endif
1065  gslpp::complex F_17 = myF_1->F_17re(muh, z, sh, 20) + gslpp::complex::i() * myF_1->F_17im(muh, z, sh, 20); /*arXiv:0810.4077*/
1066  gslpp::complex F_27 = myF_2->F_27re(muh, z, sh, 20) + gslpp::complex::i() * myF_2->F_27im(muh, z, sh, 20); /*arXiv:0810.4077*/
1067  gslpp::complex F_87 = F87_0 + F87_1 * sh + F87_2 * sh2 + F87_3 * sh * sh2 - 8. / 9. * log(sh) * (sh + sh2 + sh * sh2);
1068 
1069  if (!conjugate) {
1070  gslpp::complex delta = C_1 * F_17 + C_2 * F_27;
1071  gslpp::complex delta_t = C_8 * F_87 + delta;
1072 #if FULLNLOQCDF_MVLL
1073  gslpp::complex delta_u = delta + C_1 * Fu_17 + C_2 * Fu_27;
1074  return -alpha_s_mub / (4. * M_PI) * (delta_t - lambda_u / lambda_t * delta_u);
1075 #else
1076  return -alpha_s_mub / (4. * M_PI) * delta_t;
1077 #endif
1078  } else {
1079  gslpp::complex delta = C_1.conjugate() * F_17 + C_2.conjugate() * F_27;
1080  gslpp::complex delta_t = C_8.conjugate() * F_87 + delta;
1081 #if FULLNLOQCDF_MVLL
1082  gslpp::complex delta_u = delta + C_1.conjugate() * Fu_17 + C_2.conjugate() * Fu_27;
1083  return -alpha_s_mub / (4. * M_PI) * (delta_t - (lambda_u / lambda_t).conjugate() * delta_u);
1084 #else
1085  return -alpha_s_mub / (4. * M_PI) * delta_t;
1086 #endif
1087  }
1088 }
1089 
1090 gslpp::complex MVll::deltaC9_QCDF(double q2, bool conjugate, bool spline)
1091 {
1092 #if COMPUTECP && SPLINE
1093  if (spline && !conjugate) return gsl_spline_eval(spline_Re_deltaC9_QCDF, q2, acc_Re_deltaC9_QCDF);
1094  else if (spline && conjugate) return gsl_spline_eval(spline_Re_deltaC9_QCDF_conj, q2, acc_Re_deltaC9_QCDF_conj);
1095 #elif SPLINE
1096  if (spline) return gsl_spline_eval(spline_Re_deltaC9_QCDF, q2, acc_Re_deltaC9_QCDF);
1097 #endif
1098 
1099  double muh = mu_b / mb_pole;
1100  double z = mc_pole * mc_pole / mb_pole / mb_pole;
1101  double sh = q2 / mb_pole / mb_pole;
1102  double sh2 = sh*sh;
1103 
1104 #if FULLNLOQCDF_MVLL
1105  gslpp::complex B_Sdl = B_Seidel(q2, mb_pole*mb_pole); /* hep-ph/0403185v2.*/
1106  gslpp::complex C_Sdl = C_Seidel(q2); /* hep-ph/0403185v2.*/
1107  gslpp::complex Fu_19 = -(B_Sdl + 4. * C_Sdl); /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1108  gslpp::complex Fu_29 = -(-6. * B_Sdl + 3. * C_Sdl); /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1109 #endif
1110  gslpp::complex F_19 = myF_1->F_19re(muh, z, sh, 20) + gslpp::complex::i() * myF_1->F_19im(muh, z, sh, 20); /*arXiv:0810.4077*/
1111  gslpp::complex F_29 = myF_2->F_29re(muh, z, sh, 20) + gslpp::complex::i() * myF_2->F_29im(muh, z, sh, 20); /*arXiv:0810.4077*/
1112  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));
1113 
1114  if (!conjugate) {
1115  gslpp::complex delta = C_1 * F_19 + C_2 * F_29;
1116  gslpp::complex delta_t = C_8 * F_89 + delta;
1117 #if FULLNLOQCDF_MVLL
1118  gslpp::complex delta_u = delta + C_1 * Fu_19 + C_2 * Fu_29;
1119  return -alpha_s_mub / (4. * M_PI) * (delta_t - lambda_u / lambda_t * delta_u);
1120 #else
1121  return -alpha_s_mub / (4. * M_PI) * delta_t;
1122 #endif
1123  } else {
1124  gslpp::complex delta = C_1.conjugate() * F_19 + C_2.conjugate() * F_29;
1125  gslpp::complex delta_t = C_8.conjugate() * F_89 + delta;
1126 #if FULLNLOQCDF_MVLL
1127  gslpp::complex delta_u = delta + C_1.conjugate() * Fu_19 + C_2.conjugate() * Fu_29;
1128  return -alpha_s_mub / (4. * M_PI) * (delta_t - (lambda_u / lambda_t).conjugate() * delta_u);
1129 #else
1130  return -alpha_s_mub / (4. * M_PI) * delta_t;
1131 #endif
1132  }
1133 }
1134 
1135 gslpp::complex MVll::Cq34(bool conjugate)
1136 {
1137  gslpp::complex T_t = C_3 + 4. / 3. * (C_4 + 12. * C_5 + 16. * C_6);
1138  gslpp::complex T_u = 0.; /* 0 for K*0, phi*/
1139  if (meson == QCD::B_P) T_u = -3. * C_2;
1140  else if (vectorM == QCD::PHI) T_t = T_t + 6. * (C_3 + 10. * C_5);
1141  if (!conjugate) return T_t + lambda_u / lambda_t * T_u;
1142  else return T_t + (lambda_u / lambda_t).conjugate() * T_u;
1143 }
1144 
1145 gslpp::complex MVll::T_para_minus_WA(bool conjugate)
1146 {
1147  return -spectator_charge * 4. * MM / mb_pole * Cq34(conjugate);
1148 }
1149 
1150 gslpp::complex MVll::T_perp_WA_1()
1151 {
1152  return -spectator_charge * 4. / mb_pole * (C_3 + 4. / 3. * (C_4 + 3. * C_5 + 4. * C_6));
1153 }
1154 
1155 gslpp::complex MVll::T_perp_WA_2(bool conjugate)
1156 {
1157  return spectator_charge * 2. / mb_pole * Cq34(conjugate);
1158 }
1159 
1160 gslpp::complex MVll::T_perp_plus_O8(double q2, double u)
1161 {
1162  double ubar = 1. - u;
1163  double ed = -1. / 3.;
1164 
1165  return -(alpha_s_mub / (3. * M_PI))*4. * ed * C_8 / (u + ubar * q2 / MM2);
1166 }
1167 
1168 gslpp::complex MVll::T_para_minus_O8(double q2, double u)
1169 {
1170  double ubar = 1. - u;
1171 
1172  return (alpha_s_mub / (3. * M_PI))*spectator_charge * 8. * C_8 / (ubar + u * q2 / MM2);
1173 }
1174 
1175 gslpp::complex MVll::t_perp(double q2, double u, double m2)
1176 {
1177  double EV = (MM2 - q2 + MV2) / (2. * MM);
1178  double ubar = 1. - u;
1179 
1180  return (2. * MM) / (ubar * EV) * I1(q2, u, m2) + q2 / (ubar * ubar * EV * EV) * B0diff(q2, u, m2);
1181 
1182 }
1183 
1184 gslpp::complex MVll::t_para(double q2, double u, double m2)
1185 {
1186  double EV = (MM2pMV2 - q2) / (2. * MM);
1187  double ubar = 1. - u;
1188  return (2. * MM) / (ubar * EV) * I1(q2, u, m2) + (ubar * MM2 + u * q2) / (ubar * ubar * EV * EV) * B0diff(q2, u, m2);
1189 }
1190 
1191 gslpp::complex MVll::I1(double q2, double u, double m2)
1192 {
1193  if (m2 == 0.) return 1.;
1194 
1195  ubar = 1. - u;
1196  xp = 0.5 + sqrt(0.25 - ((gslpp::complex) m2) / (ubar * MM2 + u * q2));
1197  xm = 0.5 - sqrt(0.25 - ((gslpp::complex) m2) / (ubar * MM2 + u * q2));
1198  yp = 0.5 + sqrt(0.25 - ((gslpp::complex) m2) / q2);
1199  ym = 0.5 - sqrt(0.25 - ((gslpp::complex) m2) / q2);
1200  L1xp = log(1. - 1. / xp) * log(1. - xp) - M_PI2osix + dilog(xp / (xp - 1.));
1201  L1xm = log(1. - 1. / xm) * log(1. - xm) - M_PI2osix + dilog(xm / (xm - 1.));
1202  L1yp = log(1. - 1. / yp) * log(1. - yp) - M_PI2osix + dilog(yp / (yp - 1.));
1203  L1ym = log(1. - 1. / ym) * log(1. - ym) - M_PI2osix + dilog(ym / (ym - 1.));
1204 
1205  return 1. + 2. * m2 / ubar / (MM2 - q2)*(L1xp + L1xm - L1yp - L1ym);
1206 }
1207 
1208 gslpp::complex MVll::B0diff(double q2, double u, double m2)
1209 {
1210  double ubar = 1. - u;
1211 
1212  if (m2 == 0.) return -log((gslpp::complex)(-(2. / q2))) + log((gslpp::complex)(-(2. / (q2 * u + MM2 * ubar))));
1213  else return B0(ubar * MM2 + u * q2, m2) - B0(q2, m2);
1214 }
1215 
1216 gslpp::complex MVll::B0(double s, double m2)
1217 {
1218  if (4. * m2 / s == 1.) return gslpp::complex(0.);
1219  else return -2. * sqrt(4. * (m2 - gslpp::complex::i()*1.e-10) / s - 1.) * arctan(1. / sqrt(4. * (m2 - gslpp::complex::i()*1.e-10) / s - 1.));
1220 }
1221 
1222 gslpp::complex MVll::h_func(double s, double m2)
1223 {
1224  if (m2 == 0.) return 8. / 27. + 4. * gslpp::complex::i() * M_PI / 9. + 8. * log(mu_b) / 9. - 4. * log(s) / 9.;
1225  if (s == 0.) return -4. / 9. * (1. + log(m2 / mu_b / mu_b));
1226 
1227  double z = 4 * m2 / s;
1228  gslpp::complex term;
1229  if (z > 1) term = atan(1. / sqrt(z - 1.));
1230  else term = log((1. + sqrt(1. - z)) / sqrt(z)) - ihalfMPI;
1231 
1232  return -4. / 9. * log(m2 / mu_b / mu_b) + 8. / 27. + 4. / 9. * z - 4. / 9. * (2. + z) * sqrt(std::abs(z - 1.)) * term;
1233 
1234 }
1235 
1236 gslpp::complex MVll::T_perp_plus_QSS(double q2, double u, bool conjugate)
1237 {
1238  gslpp::complex t_perp_mc = t_perp(q2, u, mc_pole * mc_pole);
1239  double eu = 0.666666667;
1240 #if FULLNLOQCDF_MVLL
1241  gslpp::complex t_perp_mb = t_perp(q2, u, mb_pole*mb_pole);
1242  gslpp::complex t_perp_0 = t_perp(q2, u, 0.);
1243  double ed = -0.333333333;
1244 
1245  gslpp::complex T_t = (eu * t_perp_mc * (-C_1 / 6. + C_2 + 6. * C_6)
1246  + ed * t_perp_mb * (C_3 - C_4/6. + 16. * C_5 + 10. * C_6/3. + 4. * mb_pole / MM * (-C_3 + C_4/6. - 4. * C_5 + 2. * C_6/3.))
1247  + ed * t_perp_0 * (C_3 - C_4/6. + 16. * C_5 - 8. * C_6/3.));
1248 
1249  gslpp::complex T_u = eu * (t_perp_mc - t_perp_0)*(C_2 - C_1 / 6.);
1250 
1251  if (!conjugate) return alpha_s_mub / (3. * M_PI) * MM / (2. * mb_pole)*(T_t + lambda_u / lambda_t * T_u);
1252  else return alpha_s_mub / (3. * M_PI) * MM / (2. * mb_pole)*(T_t + (lambda_u / lambda_t).conjugate() * T_u);
1253 #else
1254  return alpha_s_mub / (3. * M_PI) * MM / (2. * mb_pole)*(eu * t_perp_mc * (-C_1 / 6. + C_2 + 6. * C_6));
1255 #endif
1256 }
1257 
1258 gslpp::complex MVll::T_para_plus_QSS(double q2, double u, bool conjugate)
1259 {
1260  gslpp::complex t_para_mc = t_para(q2, u, mc_pole * mc_pole);
1261  double eu = 0.666666667;
1262 #if FULLNLOQCDF_MVLL
1263  gslpp::complex t_para_mb = t_para(q2, u, mb_pole*mb_pole);
1264  gslpp::complex t_para_0 = t_para(q2, u, 0.);
1265  double ed = -0.333333333;
1266 
1267  gslpp::complex T_t = (eu * t_para_mc * (-C_1 / 6. + C_2 + 6. * C_6)
1268  + ed * t_para_mb * (C_3 - C_4/6. + 16.*C_5 + 10.*C_6/3.)
1269  + ed * t_para_0 * (C_3 - C_4/6. + 16.*C_5 - 8.*C_6/3.));
1270 
1271  gslpp::complex T_u = eu * (t_para_mc - t_para_0) * (C_2 - C_1/6.);
1272 
1273  if (!conjugate) return alpha_s_mub / (3. * M_PI) * MM / mb_pole * (T_t + lambda_u / lambda_t * T_u);
1274  else return alpha_s_mub / (3. * M_PI) * MM / mb_pole * (T_t + (lambda_u / lambda_t).conjugate() * T_u);
1275 #else
1276  return alpha_s_mub / (3. * M_PI) * MM / mb_pole * (eu * t_para_mc * (-C_1 / 6. + C_2 + 6. * C_6));
1277 #endif
1278 }
1279 
1280 gslpp::complex MVll::T_para_minus_QSS(double q2, double u, bool conjugate)
1281 {
1282  double ubar = 1. - u;
1283  gslpp::complex h_mc = h_func(ubar * MM2 + u*q2, mc_pole * mc_pole);
1284 #if FULLNLOQCDF_MVLL
1285  gslpp::complex h_mb = h_func(ubar*MM2 + u*q2, mb_pole*mb_pole);
1286  gslpp::complex h_0 = h_func(ubar*MM2 + u*q2, 0);
1287 
1288  gslpp::complex T_t = (h_mc * (-C_1 / 6. + C_2 + C_4 + 10. * C_6)
1289  + h_mb * (C_3 + 5.*C_4/6. + 16.*C_5 + 22.*C_6/3.)
1290  + h_0 * (C_3 + 17.*C_4/6. + 16.*C_5 + 82.*C_6/3.)
1291  - 8./27. * (-15.*C_4/2. + 12.*C_5 - 32.*C_6));
1292 
1293  gslpp::complex T_u = (h_mc - h_0)*(C_2 - C_1/6.);
1294 
1295  if (!conjugate) return alpha_s_mub / (3. * M_PI) * spectator_charge * 6. * MM / mb_pole * (T_t + lambda_u / lambda_t * T_u);
1296  else return alpha_s_mub / (3. * M_PI) * spectator_charge * 6. * MM / mb_pole * (T_t + (lambda_u / lambda_t).conjugate() * T_u);
1297 #else
1298  return alpha_s_mub / (3. * M_PI) * spectator_charge * 6. * MM / mb_pole * (h_mc * (-C_1 / 6. + C_2 + C_4 + 10. * C_6));
1299 #endif
1300 }
1301 
1302 double MVll::phi_V(double u)
1303 {
1304  return 6. * u * (1. - u) * (1. + mySM.getMesons(vectorM).getGegenalpha(0) * gsl_sf_gegenpoly_1(3. / 2., (2. * u - 1.)) + mySM.getMesons(vectorM).getGegenalpha(1) * gsl_sf_gegenpoly_2(3. / 2., (2. * u - 1.)));
1305 }
1306 
1307 gslpp::complex MVll::lambda_B_minus(double q2)
1308 {
1309  double w0 = mySM.getMesons(meson).getLambdaM();
1310  return 1. / (exp(-q2 / MM / w0) / w0 * (-gsl_sf_expint_Ei(q2 / MM / w0) + gslpp::complex::i() * M_PI));
1311 }
1312 
1313 double MVll::T_perp_real(double q2, double u, bool conjugate)
1314 {
1315  gslpp::complex T_amp = N_QCDF / mySM.getMesons(meson).getLambdaM() * phi_V(u) * (T_perp_plus_O8(q2, u) + T_perp_plus_QSS(q2, u, conjugate));
1316 #if FULLNLOQCDF_MVLL
1317  double ubar = 1. - u;
1318 
1319  T_amp += N_QCDF/(ubar + u*q2/MM2) * phi_V(u) * T_perp_WA_1()
1320  + N_QCDF/mySM.getMesons(meson).getLambdaM() * fpara/fperp * MV/(1. - q2/MM2) * T_perp_WA_2(conjugate);
1321  /*last term proportional to T_perp_WA_2 is a constant but is included in the integral because u is integrated over the range [0,1]*/
1322 #endif
1323  return T_amp.real();
1324 }
1325 
1326 double MVll::T_perp_imag(double q2, double u, bool conjugate)
1327 {
1328  gslpp::complex T_amp = N_QCDF / mySM.getMesons(meson).getLambdaM() * phi_V(u) * (T_perp_plus_O8(q2, u) + T_perp_plus_QSS(q2, u, conjugate));
1329 #if FULLNLOQCDF_MVLL
1330  double ubar = 1. - u;
1331 
1332  T_amp += N_QCDF/(ubar + u*q2/MM2) * phi_V(u) * T_perp_WA_1()
1333  + N_QCDF/mySM.getMesons(meson).getLambdaM() * fpara/fperp * MV/(1. - q2/MM2) * T_perp_WA_2(conjugate);
1334  /*last term proportional to T_perp_WA_2 is a constant but is included in the integral because u is integrated over the range [0,1]*/
1335 #endif
1336  return T_amp.imag();
1337 }
1338 
1339 double MVll::T_para_real(double q2, double u, bool conjugate)
1340 {
1341  double N = N_QCDF * (MV / ((MM2pMV2 - q2) / (2. * MM)));
1342 
1343  gslpp::complex T_amp = (N / lambda_B_minus(q2) * (T_para_minus_O8(q2, u) + T_para_minus_QSS(q2, u, conjugate))
1344  + N / mySM.getMesons(meson).getLambdaM() * T_para_plus_QSS(q2, u, conjugate)) * phi_V(u);
1345 #if FULLNLOQCDF_MVLL
1346  T_amp += N / lambda_B_minus(q2) * T_para_minus_WA(conjugate)* phi_V(u);
1347 #endif
1348  return sqrt(q2) * T_amp.real();
1349 }
1350 
1351 double MVll::T_para_imag(double q2, double u, bool conjugate)
1352 {
1353  double N = N_QCDF * (MV / ((MM2pMV2 - q2) / (2. * MM)));
1354 
1355  gslpp::complex T_amp = (N / lambda_B_minus(q2) * (/* + */T_para_minus_O8(q2, u) + T_para_minus_QSS(q2, u, conjugate))
1356  + N / mySM.getMesons(meson).getLambdaM() * T_para_plus_QSS(q2, u, conjugate)) * phi_V(u);
1357 #if FULLNLOQCDF_MVLL
1358  T_amp += N / lambda_B_minus(q2) * T_para_minus_WA(conjugate) * phi_V(u);
1359 #endif
1360  return sqrt(q2) * T_amp.imag();
1361 }
1362 
1363 double MVll::T_perp_real(double q2, bool conjugate)
1364 {
1365  FS = convertToGslFunction(boost::bind(&MVll::T_perp_real, &(*this), q2, _1, conjugate));
1366  gsl_integration_cquad(&FS, 0., 1., 1.e-2, 1.e-1, w_sigma, &avaSigma, &errSigma, NULL);
1367 
1368  return avaSigma;
1369 }
1370 
1371 double MVll::T_perp_imag(double q2, bool conjugate)
1372 {
1373  FS = convertToGslFunction(boost::bind(&MVll::T_perp_imag, &(*this), q2, _1, conjugate));
1374  gsl_integration_cquad(&FS, 0., 1., 1.e-2, 1.e-1, w_sigma, &avaSigma, &errSigma, NULL);
1375 
1376  return avaSigma;
1377 }
1378 
1379 double MVll::T_para_real(double q2, bool conjugate)
1380 {
1381  FS = convertToGslFunction(boost::bind(&MVll::T_para_real, &(*this), q2, _1, conjugate));
1382  gsl_integration_cquad(&FS, 0., 1., 1.e-2, 1.e-1, w_sigma, &avaSigma, &errSigma, NULL);
1383 
1384  return avaSigma;
1385 }
1386 
1387 double MVll::T_para_imag(double q2, bool conjugate)
1388 {
1389  FS = convertToGslFunction(boost::bind(&MVll::T_para_imag, &(*this), q2, _1, conjugate));
1390  gsl_integration_cquad(&FS, 0., 1., 1.e-2, 1.e-1, w_sigma, &avaSigma, &errSigma, NULL);
1391 
1392  return avaSigma;
1393 }
1394 
1395 double MVll::QCDF_fit_func(double* x, double* p)
1396 {
1397  return p[0] + p[1] * x[0] + p[2] * x[0] * x[0] + p[3] * x[0] * x[0] * x[0] + p[4] * x[0] * x[0] * x[0] * x[0] + p[5] * x[0] * x[0] * x[0] * x[0] * x[0] + p[6] * x[0] * x[0] * x[0] * x[0] * x[0] * x[0];
1398 }
1399 
1400 void MVll::fit_QCDF_func()
1401 {
1402  int dim = 0;
1403  for (double i = 0.001; i < 8.6; i += 0.5) {
1404  myq2.push_back(i);
1405  Re_T_perp.push_back(T_perp_real(i, false));
1406  Im_T_perp.push_back(T_perp_imag(i, false));
1407  Re_T_para.push_back(T_para_real(i, false));
1408  Im_T_para.push_back(T_para_imag(i, false));
1409 
1410 #if COMPUTECP
1411  Re_T_perp_conj.push_back(T_perp_real(i, true));
1412  Im_T_perp_conj.push_back(T_perp_imag(i, true));
1413  Re_T_para_conj.push_back(T_para_real(i, true));
1414  Im_T_para_conj.push_back(T_para_imag(i, true));
1415 #endif
1416  dim++;
1417  }
1418 
1419  gr1 = TGraph(dim, myq2.data(), Re_T_perp.data());
1420  QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Re_T_perp");
1421  Re_T_perp_res = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
1422  Re_T_perp.clear();
1423 
1424  gr1 = TGraph(dim, myq2.data(), Im_T_perp.data());
1425  QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Im_T_perp");
1426  Im_T_perp_res = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
1427  Im_T_perp.clear();
1428 
1429  gr1 = TGraph(dim, myq2.data(), Re_T_para.data());
1430  QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Re_T_para");
1431  Re_T_para_res = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
1432  Re_T_para.clear();
1433 
1434  gr1 = TGraph(dim, myq2.data(), Im_T_para.data());
1435  QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Im_T_para");
1436  Im_T_para_res = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
1437  Im_T_para.clear();
1438 
1439 #if COMPUTECP
1440  gr1 = TGraph(dim, myq2.data(), Re_T_perp_conj.data());
1441  QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Re_T_perp_conj");
1442  Re_T_perp_res_conj = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
1443  Re_T_perp_conj.clear();
1444 
1445  gr1 = TGraph(dim, myq2.data(), Im_T_perp_conj.data());
1446  QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Im_T_perp_conj");
1447  Im_T_perp_res_conj = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
1448  Im_T_perp_conj.clear();
1449 
1450  gr1 = TGraph(dim, myq2.data(), Re_T_para_conj.data());
1451  QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Re_T_para_conj");
1452  Re_T_para_res_conj = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
1453  Re_T_para_conj.clear();
1454 
1455  gr1 = TGraph(dim, myq2.data(), Im_T_para_conj.data());
1456  QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Im_T_para_conj");
1457  Im_T_para_res_conj = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
1458  Im_T_para_conj.clear();
1459 #endif
1460 
1461  myq2.clear();
1462 }
1463 
1464 void MVll::spline_QCDF_func()
1465 {
1466  int dim = GSL_INTERP_DIM;
1467  int dim_DC = GSL_INTERP_DIM_DC;
1468  double min = 0.001;
1469  double interval = (9.9 - min) / ((double) dim);
1470  double interval_DC = (9.9 - min) / ((double) dim_DC);
1471  double q2_spline[dim];
1472  double fq2_Re_T_perp[dim], fq2_Im_T_perp[dim], fq2_Re_T_para[dim], fq2_Im_T_para[dim];
1473  double q2_spline_DC[dim_DC];
1474  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];
1475 #if COMPUTECP
1476  double fq2_Re_T_perp_conj[dim], fq2_Im_T_perp_conj[dim], fq2_Re_T_para_conj[dim], fq2_Im_T_para_conj[dim];
1477  double fq2_Re_deltaC7_QCDF_conj[dim_DC], fq2_Im_deltaC7_QCDF_conj[dim_DC], fq2_Re_deltaC9_QCDF_conj[dim_DC], fq2_Im_deltaC9_QCDF_conj[dim_DC];
1478 #endif
1479 
1480  for (int i = 0; i < dim; i++) {
1481  q2_spline[i] = min + (double) i*interval;
1482  fq2_Re_T_perp[i] = T_perp_real(q2_spline[i], false);
1483  fq2_Im_T_perp[i] = T_perp_imag(q2_spline[i], false);
1484  fq2_Re_T_para[i] = T_para_real(q2_spline[i], false);
1485  fq2_Im_T_para[i] = T_para_imag(q2_spline[i], false);
1486 
1487 #if COMPUTECP
1488  fq2_Re_T_perp_conj[i] = T_perp_real(q2_spline[i], true);
1489  fq2_Im_T_perp_conj[i] = T_perp_imag(q2_spline[i], true);
1490  fq2_Re_T_para_conj[i] = T_para_real(q2_spline[i], true);
1491  fq2_Im_T_para_conj[i] = T_para_imag(q2_spline[i], true);
1492 #endif
1493  }
1494  for (int i = 0; i < dim_DC; i++) {
1495  q2_spline_DC[i] = min + (double) i*interval_DC;
1496  fq2_Re_deltaC7_QCDF[i] = deltaC7_QCDF(q2_spline_DC[i], false).real();
1497  fq2_Im_deltaC7_QCDF[i] = deltaC7_QCDF(q2_spline_DC[i], false).imag();
1498  fq2_Re_deltaC9_QCDF[i] = deltaC9_QCDF(q2_spline_DC[i], false).real();
1499  fq2_Im_deltaC9_QCDF[i] = deltaC9_QCDF(q2_spline_DC[i], false).imag();
1500 
1501 #if COMPUTECP
1502  fq2_Re_deltaC7_QCDF_conj[i] = deltaC7_QCDF(q2_spline_DC[i], true).real();
1503  fq2_Im_deltaC7_QCDF_conj[i] = deltaC7_QCDF(q2_spline_DC[i], true).imag();
1504  fq2_Re_deltaC9_QCDF_conj[i] = deltaC9_QCDF(q2_spline_DC[i], true).real();
1505  fq2_Im_deltaC9_QCDF_conj[i] = deltaC9_QCDF(q2_spline_DC[i], true).imag();
1506 #endif
1507  }
1508 
1509  gsl_spline_init(spline_Re_T_perp, q2_spline, fq2_Re_T_perp, dim);
1510  gsl_spline_init(spline_Im_T_perp, q2_spline, fq2_Im_T_perp, dim);
1511  gsl_spline_init(spline_Re_T_para, q2_spline, fq2_Re_T_para, dim);
1512  gsl_spline_init(spline_Im_T_para, q2_spline, fq2_Im_T_para, dim);
1513 
1514  gsl_spline_init(spline_Re_deltaC7_QCDF, q2_spline_DC, fq2_Re_deltaC7_QCDF, dim_DC);
1515  gsl_spline_init(spline_Im_deltaC7_QCDF, q2_spline_DC, fq2_Im_deltaC7_QCDF, dim_DC);
1516  gsl_spline_init(spline_Re_deltaC9_QCDF, q2_spline_DC, fq2_Re_deltaC9_QCDF, dim_DC);
1517  gsl_spline_init(spline_Im_deltaC9_QCDF, q2_spline_DC, fq2_Im_deltaC9_QCDF, dim_DC);
1518 
1519 #if COMPUTECP
1520  gsl_spline_init(spline_Re_T_perp_conj, q2_spline, fq2_Re_T_perp_conj, dim);
1521  gsl_spline_init(spline_Im_T_perp_conj, q2_spline, fq2_Im_T_perp_conj, dim);
1522  gsl_spline_init(spline_Re_T_para_conj, q2_spline, fq2_Re_T_para_conj, dim);
1523  gsl_spline_init(spline_Im_T_para_conj, q2_spline, fq2_Im_T_para_conj, dim);
1524 
1525  gsl_spline_init(spline_Re_deltaC7_QCDF_conj, q2_spline_DC, fq2_Re_deltaC7_QCDF_conj, dim_DC);
1526  gsl_spline_init(spline_Im_deltaC7_QCDF_conj, q2_spline_DC, fq2_Im_deltaC7_QCDF_conj, dim_DC);
1527  gsl_spline_init(spline_Re_deltaC9_QCDF_conj, q2_spline_DC, fq2_Re_deltaC9_QCDF_conj, dim_DC);
1528  gsl_spline_init(spline_Im_deltaC9_QCDF_conj, q2_spline_DC, fq2_Im_deltaC9_QCDF_conj, dim_DC);
1529 #endif
1530 
1531 }
1532 
1533 gslpp::complex MVll::T_minus(double q2, bool conjugate)
1534 {
1535 #if COMPUTECP && SPLINE
1536  if (!conjugate) return -2. * MM * mb_pole / q2 * (1. - q2 / MM2) * (gsl_spline_eval(spline_Re_T_perp, q2, acc_Re_T_perp) + gslpp::complex::i() * gsl_spline_eval(spline_Im_T_perp, q2, acc_Im_T_perp));
1537  else return -2. * MM * mb_pole / q2 * (1. - q2 / MM2) * (gsl_spline_eval(spline_Re_T_perp_conj, q2, acc_Re_T_perp_conj) + gslpp::complex::i() * gsl_spline_eval(spline_Im_T_perp_conj, q2, acc_Im_T_perp_conj));
1538 #elif SPLINE
1539  return -2. * MM * mb_pole / q2 * (1. - q2 / MM2) * (gsl_spline_eval(spline_Re_T_perp, q2, acc_Re_T_perp) + gslpp::complex::i() * gsl_spline_eval(spline_Im_T_perp, q2, acc_Im_T_perp));
1540 #endif
1541 
1542 #if COMPUTECP && !SPLINE
1543  if (!conjugate) return -2. * MM * mb_pole / q2 * (1. - q2 / MM2) * (QCDF_fit_func(&q2, const_cast<double *> (Re_T_perp_res->GetParams())) + gslpp::complex::i() * QCDF_fit_func(&q2, const_cast<double *> (Im_T_perp_res->GetParams())));
1544  else return -2. * MM * mb_pole / q2 * (1. - q2 / MM2) * (QCDF_fit_func(&q2, const_cast<double *> (Re_T_perp_res_conj->GetParams())) + gslpp::complex::i() * QCDF_fit_func(&q2, const_cast<double *> (Im_T_perp_res_conj->GetParams())));
1545 #elif !SPLINE
1546  return -2. * MM * mb_pole / q2 * (1. - q2 / MM2) * (QCDF_fit_func(&q2, const_cast<double *> (Re_T_perp_res->GetParams())) + gslpp::complex::i() * QCDF_fit_func(&q2, const_cast<double *> (Im_T_perp_res->GetParams())));
1547 #endif
1548 }
1549 
1550 gslpp::complex MVll::T_0(double q2, bool conjugate)
1551 {
1552 #if COMPUTECP && SPLINE
1553  if (!conjugate) return -(1. - q2 / MM2)* (1. - q2 / MM2) * MM * mb_pole / sqrt(q2) * (gsl_spline_eval(spline_Re_T_para, q2, acc_Re_T_para) + gslpp::complex::i() * gsl_spline_eval(spline_Im_T_para, q2, acc_Im_T_para));
1554  else return -(1. - q2 / MM2)* (1. - q2 / MM2) * MM * mb_pole / sqrt(q2) * (gsl_spline_eval(spline_Re_T_para_conj, q2, acc_Re_T_para_conj) + gslpp::complex::i() * gsl_spline_eval(spline_Im_T_para_conj, q2, acc_Im_T_para_conj));
1555 #elif SPLINE
1556  return -(1. - q2 / MM2)* (1. - q2 / MM2) * MM * mb_pole / sqrt(q2) * (gsl_spline_eval(spline_Re_T_para, q2, acc_Re_T_para) + gslpp::complex::i() * gsl_spline_eval(spline_Im_T_para, q2, acc_Im_T_para));
1557 #endif
1558 
1559 #if COMPUTECP && !SPLINE
1560  if (!conjugate) return -(1. - q2 / MM2)* (1. - q2 / MM2) * MM * mb_pole / sqrt(q2) * (QCDF_fit_func(&q2, const_cast<double *> (Re_T_para_res->GetParams())) + gslpp::complex::i() * QCDF_fit_func(&q2, const_cast<double *> (Im_T_para_res->GetParams())));
1561  else return -(1. - q2 / MM2)* (1. - q2 / MM2) * MM * mb_pole / sqrt(q2) * (QCDF_fit_func(&q2, const_cast<double *> (Re_T_para_res_conj->GetParams())) + gslpp::complex::i() * QCDF_fit_func(&q2, const_cast<double *> (Im_T_para_res_conj->GetParams())));
1562 #elif !SPLINE
1563  return -(1. - q2 / MM2)* (1. - q2 / MM2) * MM * mb_pole / sqrt(q2) * (QCDF_fit_func(&q2, const_cast<double *> (Re_T_para_res->GetParams())) + gslpp::complex::i() * QCDF_fit_func(&q2, const_cast<double *> (Im_T_para_res->GetParams())));
1564 #endif
1565 }
1566 
1567 /*******************************************************************************
1568  * Helicity amplitudes *
1569  * ****************************************************************************/
1570 gslpp::complex MVll::H(double q2, double m2, double mu2)
1571 {
1572  double x = 4. * m2 / q2;
1573  gslpp::complex par;
1574 
1575  if (x > 1.) par = sqrt(x - 1.) * atan(1. / sqrt(x - 1.));
1576  else par = sqrt(1. - x) * (log((1. + sqrt(1. - x)) / sqrt(x)) - ihalfMPI);
1577 
1578  return -fournineth * (log(m2 / mu2) - twothird - x) - fournineth * (2. + x) * par;
1579 }
1580 
1581 gslpp::complex MVll::H_0(double q2)
1582 {
1583  return (H_0_pre - fournineth * log(q2 / mu_b2));
1584 }
1585 
1586 gslpp::complex MVll::Y(double q2)
1587 {
1588  return -half * H_0(q2) * H_0_WC + H(q2, mc_pole*mc_pole, mu_b2) * H_c_WC - half * H(q2, mb_pole*mb_pole, mu_b2) * H_b_WC;
1589 }
1590 
1591 gslpp::complex MVll::funct_g(double q2)
1592 {
1593  if (q2 < 4. * Mc * Mc)
1594  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.)));
1595  else
1596  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));
1597 }
1598 
1599 gslpp::complex MVll::DeltaC9_KD(double q2, int com)
1600 {
1601  return ((h_0[com] * (1. - 1. / q2) + h_2[com] / q2) / (1. + h_1[com] * (1. - q2) / mJ2) - (3. * (-0.267) + 1.117) * funct_g(q2))*exp_Phase[com];
1602  /* C_1 = -0.267 and C_2 = 1.117 in KMPW */
1603 }
1604 
1605 gslpp::complex MVll::h_lambda(int hel, double q2)
1606 {
1607  if (!dispersion) {
1608  if (h_pole == true) return SU3_breaking * (h_0[hel]+(1. - h_2[hel]) * q2 * (h_1[hel] - h_0[hel]) / (q2 - h_2[hel]));
1609  else if (hel == 1) return SU3_breaking * (h_0[1] + h_1[1] * q2 + h_2[1] * q2 * q2 + (twoMboMM * h_0[2] * T_p(q2) + h_1[2] * q2 / MM2 * V_p(q2)) / sixteenM_PI2);
1610  else if (hel == 2) return SU3_breaking * (twoMboMM * h_0[2] * T_m(q2) + h_1[2] * q2 / MM2 * V_m(q2)) / sixteenM_PI2 + h_2[2] * q2 * q2;
1611  else return SU3_breaking * ((h_0[hel] + h_1[hel] * q2) * sqrt(q2) + (twoMboMM * h_0[2] * T_0t(q2) + h_1[2] * q2 * V_0t(q2) / MM2) / sixteenM_PI2);
1612  } else {
1613  if (hel == 0) return SU3_breaking * (-sqrt(q2) / (MM2 * 16. * M_PI * M_PI) * ((MMpMV2 * (MM2mMV2 - q2) * A_1(q2) * DeltaC9_KD(q2, 1) - lambda(q2) * A_2(q2) * DeltaC9_KD(q2, 2)) / (4. * MV * MM * MMpMV)));
1614  else if (hel == 1) {
1615  if (q2 == 0.) return SU3_breaking * (-1. / (MM2 * 16. * M_PI * M_PI) * (
1616  (MMpMV * A_1(0.)) / (2. * MM) * ((-h_0[1] + h_2[1]) / (1. + h_1[1] / mJ2)) * exp_Phase[1]
1617  - sqrt(lambda(0.)) / (2. * MM * MMpMV) * V(0.) * ((-h_0[0] + h_2[0]) / (1. + h_1[0] / mJ2)) * exp_Phase[1]));
1618  else return SU3_breaking * (-q2 / (MM2 * 16. * M_PI * M_PI) * ((MMpMV * A_1(q2)) / (2. * MM) * DeltaC9_KD(q2, 1) - sqrt(lambda(q2)) / (2. * MM * MMpMV) * V(q2) * DeltaC9_KD(q2, 0)));
1619  } else {
1620  if (q2 == 0.) return SU3_breaking * (-1. / (MM2 * 16. * M_PI * M_PI) *
1621  ((MMpMV * A_1(0.)) / (2. * MM) * ((-h_0[1] + h_2[1]) / (1. + h_1[1] / mJ2)) * exp_Phase[1]
1622  + sqrt(lambda(0.)) / (2. * MM * MMpMV) * V(0.) * ((-h_0[0] + h_2[0]) / (1. + h_1[0] / mJ2)) * exp_Phase[1]));
1623  else return SU3_breaking * (-q2 / (MM2 * 16. * M_PI * M_PI) * ((MMpMV * A_1(q2)) / (2. * MM) * DeltaC9_KD(q2, 1) + sqrt(lambda(q2)) / (2. * MM * MMpMV) * V(q2) * DeltaC9_KD(q2, 0)));
1624  }
1625  }
1626 }
1627 
1628 gslpp::complex MVll::H_V_0(double q2, bool bar)
1629 {
1630  if (!bar) return -gslpp::complex::i() * NN * (((C_9 + deltaC9_QCDF(q2, !bar, SPLINE) /*+ fDeltaC9_0(q2)*/ + Y(q2)) - etaV * pow(-1, angmomV) * C_9p) * V_0t(q2) + T_0(q2, !bar) + MM2 / q2 * (twoMboMM * (C_7 + deltaC7_QCDF(q2, !bar, SPLINE) - etaV * pow(-1, angmomV) * C_7p) * T_0t(q2) - sixteenM_PI2 * h_lambda(0, q2)));
1631  return -gslpp::complex::i() * NN_conjugate * (((C_9.conjugate() + deltaC9_QCDF(q2, bar, SPLINE) /*+ fDeltaC9_0(q2)*/ + Y(q2)) - etaV * pow(-1, angmomV) * C_9p.conjugate()) * V_0t(q2) + T_0(q2, bar) + MM2 / q2 * (twoMboMM * (C_7.conjugate() + deltaC7_QCDF(q2, bar, SPLINE) - etaV * pow(-1, angmomV) * C_7p.conjugate()) * T_0t(q2) - sixteenM_PI2 * h_lambda(0, q2)));
1632 
1633 }
1634 
1635 gslpp::complex MVll::H_V_p(double q2, bool bar)
1636 {
1637  if (!bar) return -gslpp::complex::i() * NN * (((C_9 + deltaC9_QCDF(q2, !bar, SPLINE) /*+ fDeltaC9_p(q2)*/ + Y(q2)) * V_p(q2) - etaV * pow(-1, angmomV) * C_9p * V_m(q2)) + MM2 / q2 * (twoMboMM * ((C_7 + deltaC7_QCDF(q2, !bar, SPLINE)) * T_p(q2) - etaV * pow(-1, angmomV) * C_7p * T_m(q2)) - sixteenM_PI2 * h_lambda(1, q2)));
1638  return -gslpp::complex::i() * NN_conjugate * (((C_9.conjugate() + deltaC9_QCDF(q2, bar, SPLINE) /*+ fDeltaC9_p(q2)*/ + Y(q2)) * V_p(q2) - etaV * pow(-1, angmomV) * C_9p.conjugate() * V_m(q2)) + MM2 / q2 * (twoMboMM * ((C_7.conjugate() + deltaC7_QCDF(q2, bar, SPLINE)) * T_p(q2) - etaV * pow(-1, angmomV) * C_7p.conjugate() * T_m(q2)) - sixteenM_PI2 * h_lambda(1, q2)));
1639 }
1640 
1641 gslpp::complex MVll::H_V_m(double q2, bool bar)
1642 {
1643  if (!bar) return -gslpp::complex::i() * NN * (((C_9 + deltaC9_QCDF(q2, !bar, SPLINE) /*+ fDeltaC9_m(q2)*/ + Y(q2)) * V_m(q2) - etaV * pow(-1, angmomV) * C_9p * V_p(q2)) + T_minus(q2, !bar) + MM2 / q2 * (twoMboMM * ((C_7 + deltaC7_QCDF(q2, !bar, SPLINE)) * T_m(q2) - etaV * pow(-1, angmomV) * C_7p * T_p(q2)) - sixteenM_PI2 * h_lambda(2, q2)));
1644  return -gslpp::complex::i() * NN_conjugate * (((C_9.conjugate() + deltaC9_QCDF(q2, bar, SPLINE) /*+ fDeltaC9_m(q2)*/ + Y(q2)) * V_m(q2) - etaV * pow(-1, angmomV) * C_9p.conjugate() * V_p(q2)) + T_minus(q2, bar) + MM2 / q2 * (twoMboMM * ((C_7.conjugate() + deltaC7_QCDF(q2, bar, SPLINE)) * T_m(q2) - etaV * pow(-1, angmomV) * C_7p.conjugate() * T_p(q2)) - sixteenM_PI2 * h_lambda(2, q2)));
1645 }
1646 
1647 gslpp::complex MVll::H_A_0(double q2, bool bar)
1648 {
1649  if (!bar) return gslpp::complex::i() * NN * (-C_10 + etaV * pow(-1, angmomV) * C_10p) * V_0t(q2);
1650  return gslpp::complex::i() * NN_conjugate * (-C_10.conjugate() + etaV * pow(-1, angmomV) * C_10p.conjugate()) * V_0t(q2);
1651 }
1652 
1653 gslpp::complex MVll::H_A_p(double q2, bool bar)
1654 {
1655  if (!bar) return gslpp::complex::i() * NN * (-C_10 * V_p(q2) + etaV * pow(-1, angmomV) * C_10p * V_m(q2));
1656  return gslpp::complex::i() * NN_conjugate * (-C_10.conjugate() * V_p(q2) + etaV * pow(-1, angmomV) * C_10p.conjugate() * V_m(q2));
1657 }
1658 
1659 gslpp::complex MVll::H_A_m(double q2, bool bar)
1660 {
1661  if (!bar) return gslpp::complex::i() * NN * (-C_10 * V_m(q2) + etaV * pow(-1, angmomV) * C_10p * V_p(q2));
1662  return gslpp::complex::i() * NN_conjugate * (-C_10.conjugate() * V_m(q2) + etaV * pow(-1, angmomV) * C_10p.conjugate() * V_p(q2));
1663 }
1664 
1665 gslpp::complex MVll::H_S(double q2, bool bar)
1666 {
1667  if (!bar) return gslpp::complex::i() * NN * MboMW * (C_S - etaV * pow(-1, angmomV) * C_Sp) * S_L(q2);
1668  return gslpp::complex::i() * NN_conjugate * MboMW * (C_S.conjugate() - etaV * pow(-1, angmomV) * C_Sp.conjugate()) * S_L(q2);
1669 }
1670 
1671 gslpp::complex MVll::H_P(double q2, bool bar)
1672 {
1673  if (!bar) return gslpp::complex::i() * NN * (MboMW * (C_P - etaV * pow(-1, angmomV) * C_Pp) + twoMlepMb / q2 * (C_10 * (1. + etaV * pow(-1, angmomV) * MsoMb) - C_10p * (etaV * pow(-1, angmomV) + MsoMb))) * S_L(q2);
1674  return gslpp::complex::i() * NN_conjugate * (MboMW * (C_P.conjugate() - etaV * pow(-1, angmomV) * C_Pp.conjugate()) + twoMlepMb / q2 * (C_10.conjugate()*(1. + etaV * pow(-1, angmomV) * MsoMb) - C_10p.conjugate()*(etaV * pow(-1, angmomV) + MsoMb))) * S_L(q2);
1675 }
1676 
1677 /*******************************************************************************
1678  * Angular coefficients *
1679  * ****************************************************************************/
1680 double MVll::k2(double q2)
1681 {
1682  return (MM4 + q2 * q2 + MV4 - twoMV2 * q2 - twoMM2 * (q2 + MV2)) / fourMM2;
1683 }
1684 
1685 double MVll::beta(double q2)
1686 {
1687  return sqrt(1. - 4. * Mlep2 / q2);
1688 }
1689 
1690 double MVll::beta2(double q2)
1691 {
1692  return 1. - 4. * Mlep2 / q2;
1693 }
1694 
1695 double MVll::lambda(double q2)
1696 {
1697  return (MM4 + q2 * q2 + MV4 - twoMV2 * q2 - twoMM2 * (q2 + MV2));
1698 }
1699 
1700 double MVll::F(double q2, double b_i)
1701 {
1702  return sqrt(lambda(q2)) * beta(q2) * q2 * b_i / (ninetysixM_PI3MM3);
1703 }
1704 
1705 double MVll::I_1c(double q2, bool bar)
1706 {
1707  return F(q2, b)*((H_V_0(q2, bar).abs2() + H_A_0(q2, bar).abs2()) / 2. + H_P(q2, bar).abs2() + 2. * Mlep2 / q2 * (H_V_0(q2, bar).abs2()
1708  - H_A_0(q2, bar).abs2()) + beta2(q2) * H_S(q2, bar).abs2());
1709 }
1710 
1711 double MVll::I_1s(double q2, bool bar)
1712 {
1713  return F(q2, b)*((beta2(q2) + 2.) / 8. * (H_V_p(q2, bar).abs2() + H_V_m(q2, bar).abs2() + H_A_p(q2, bar).abs2() + H_A_m(q2, bar).abs2()) +
1714  Mlep2 / q2 * (H_V_p(q2, bar).abs2() + H_V_m(q2, bar).abs2() - H_A_p(q2, bar).abs2() - H_A_m(q2, bar).abs2()));
1715 }
1716 
1717 double MVll::I_2c(double q2, bool bar)
1718 {
1719  return -F(q2, b) * beta2(q2) / 2. * (H_V_0(q2, bar).abs2() + H_A_0(q2, bar).abs2());
1720 }
1721 
1722 double MVll::I_2s(double q2, bool bar)
1723 {
1724  return F(q2, b) * beta2(q2) / 8. * (H_V_p(q2, bar).abs2() + H_V_m(q2, bar).abs2() + H_A_p(q2, bar).abs2() + H_A_m(q2, bar).abs2());
1725 }
1726 
1727 double MVll::I_3(double q2, bool bar)
1728 {
1729  return -F(q2, b) * beta2(q2) / 2. * ((H_V_p(q2, bar) * H_V_m(q2, bar).conjugate()).real() + (H_A_p(q2, bar) * H_A_m(q2, bar).conjugate()).real());
1730 }
1731 
1732 double MVll::I_4(double q2, bool bar)
1733 {
1734  return F(q2, b) * beta2(q2) / 4. * (((H_V_m(q2, bar) + H_V_p(q2, bar)) * H_V_0(q2, bar).conjugate()).real() + ((H_A_m(q2, bar) + H_A_p(q2, bar)) * H_A_0(q2, bar).conjugate()).real());
1735 }
1736 
1737 double MVll::I_5(double q2, bool bar)
1738 {
1739  return F(q2, b)*(beta(q2) / 2. * (((H_V_m(q2, bar) - H_V_p(q2, bar)) * H_A_0(q2, bar).conjugate()).real() + ((H_A_m(q2, bar) - H_A_p(q2, bar)) * H_V_0(q2, bar).conjugate()).real()) -
1740  beta(q2) * Mlep / sqrt(q2)*(H_S(q2, bar).conjugate()*(H_V_p(q2, bar) + H_V_m(q2, bar))).real());
1741 }
1742 
1743 double MVll::I_6s(double q2, bool bar)
1744 {
1745  return F(q2, b) * beta(q2)*(H_V_m(q2, bar)*(H_A_m(q2, bar).conjugate()) - H_V_p(q2, bar)*(H_A_p(q2, bar).conjugate())).real();
1746 }
1747 
1748 double MVll::I_6c(double q2, bool bar)
1749 {
1750  return 4. * F(q2, b) * beta(q2) * Mlep / sqrt(q2)*(H_S(q2, bar).conjugate() * H_V_0(q2, bar)).real();
1751 }
1752 
1753 double MVll::I_7(double q2, bool bar)
1754 {
1755  return F(q2, b)*(beta(q2) / 2. * (((H_V_m(q2, bar) + H_V_p(q2, bar)) * H_A_0(q2, bar).conjugate()).imag() + ((H_A_m(q2, bar) + H_A_p(q2, bar)) * H_V_0(q2, bar).conjugate()).imag()) -
1756  beta(q2) * Mlep / sqrt(q2)*(H_S(q2, bar).conjugate()*(H_V_m(q2, bar) - H_V_p(q2, bar))).imag());
1757 }
1758 
1759 double MVll::I_8(double q2, bool bar)
1760 {
1761  return F(q2, b) * beta2(q2) / 4. * (((H_V_m(q2, bar) - H_V_p(q2, bar)) * H_V_0(q2, bar).conjugate()).imag() + ((H_A_m(q2, bar) - H_A_p(q2, bar)) * H_A_0(q2, bar).conjugate()).imag());
1762 }
1763 
1764 double MVll::I_9(double q2, bool bar)
1765 {
1766  return F(q2, b) * beta2(q2) / 2. * ((H_V_p(q2, bar) * H_V_m(q2, bar).conjugate()).imag() + (H_A_p(q2, bar) * H_A_m(q2, bar).conjugate()).imag());
1767 }
1768 
1769 double MVll::h_1c(double q2, bool bar)
1770 {
1771  return F(q2, b)*((H_V_0(q2, bar).abs2() + H_A_0(q2, bar).abs2()) + 2. * H_P(q2, bar).abs2() + 4. * Mlep2 / q2 * (H_V_0(q2, bar).abs2()
1772  - H_A_0(q2, bar).abs2()) - 2. * beta2(q2) * H_S(q2, bar).abs2());
1773 }
1774 
1775 double MVll::h_1s(double q2, bool bar)
1776 {
1777  return F(q2, b)*((beta2(q2) + 2.) / 2. * ((H_V_p(q2, bar) * H_V_m(q2, bar).conjugate()).real()
1778  + (H_A_p(q2, bar) * H_A_m(q2, bar).conjugate()).real()) +
1779  4. * Mlep2 / q2 * ((H_V_p(q2, bar) * H_V_m(q2, bar).conjugate()).real()
1780  - (H_A_p(q2, bar) * H_A_m(q2, bar).conjugate()).real()));
1781 }
1782 
1783 double MVll::h_2c(double q2, bool bar)
1784 {
1785  return -F(q2, b) * beta2(q2) * (H_V_0(q2, bar).abs2() + H_A_0(q2, bar).abs2());
1786 }
1787 
1788 double MVll::h_2s(double q2, bool bar)
1789 {
1790  return F(q2, b) * beta2(q2) / 2. * ((H_V_p(q2, bar) * H_V_m(q2, bar).conjugate()).real()
1791  + (H_A_p(q2, bar) * H_A_m(q2, bar).conjugate()).real());
1792 }
1793 
1794 double MVll::h_3(double q2, bool bar)
1795 {
1796  return -F(q2, b) * beta2(q2) / 2. * (H_V_p(q2, bar).abs2() + H_V_m(q2, bar).abs2()
1797  + H_A_p(q2, bar).abs2() + H_A_m(q2, bar).abs2());
1798 }
1799 
1800 double MVll::h_4(double q2, bool bar)
1801 {
1802  return F(q2, b) * beta2(q2) / 2. * (((H_V_m(q2, bar) + H_V_p(q2, bar)) * H_V_0(q2, bar).conjugate()).real()
1803  + ((H_A_m(q2, bar) + H_A_p(q2, bar)) * H_A_0(q2, bar).conjugate()).real());
1804 }
1805 
1806 double MVll::h_7(double q2, bool bar)
1807 {
1808  return F(q2, b)*(beta(q2) * (((H_V_m(q2, bar) + H_V_p(q2, bar)) * H_A_0(q2, bar).conjugate()).imag()
1809  + ((H_A_m(q2, bar) + H_A_p(q2, bar)) * H_V_0(q2, bar).conjugate()).imag()) -
1810  beta(q2) * 2. * Mlep / sqrt(q2)*(H_S(q2, bar).conjugate()*(H_V_m(q2, bar) - H_V_p(q2, bar))).imag());
1811 }
1812 
1813 double MVll::s_5(double q2, bool bar)
1814 {
1815  return beta(q2) * (2. * Mlep * ((H_V_m(q2, bar) + H_V_p(q2, bar)) * F(q2, b) * H_S(q2, bar).conjugate()).imag() / sqrt(q2)
1816  - F(q2, b)*((H_A_m(q2, bar) - H_A_p(q2, bar)) * H_V_0(q2, bar).conjugate()
1817  + (H_V_m(q2, bar) - H_V_p(q2, bar)) * H_A_0(q2, bar).conjugate()).imag());
1818 }
1819 
1820 double MVll::s_6s(double q2, bool bar)
1821 {
1822  return 2. * beta(q2) * F(q2, b) * (H_A_p(q2, bar) * H_V_m(q2, bar).conjugate() + H_V_p(q2, bar) * H_A_m(q2, bar).conjugate()).imag();
1823 }
1824 
1825 double MVll::s_6c(double q2, bool bar)
1826 {
1827  return -8. * beta(q2) * Mlep * (H_V_0(q2, bar) * F(q2, b) * H_S(q2, bar).conjugate()).imag() / sqrt(q2);
1828 }
1829 
1830 double MVll::s_8(double q2, bool bar)
1831 {
1832  return beta2(q2) * F(q2, b) * ((H_A_m(q2, bar) - H_A_p(q2, bar)) * H_A_0(q2, bar).conjugate()
1833  + (H_V_m(q2, bar) - H_V_p(q2, bar)) * H_V_0(q2, bar).conjugate()).real() / 2.;
1834 }
1835 
1836 double MVll::s_9(double q2, bool bar)
1837 {
1838  return beta2(q2) * F(q2, b) * (H_A_p(q2, bar).abs2() - H_A_m(q2, bar).abs2()
1839  + H_V_p(q2, bar).abs2() - H_V_m(q2, bar).abs2()) / 2.;
1840 }
1841 
1842 double MVll::integrateSigma(int i, double q_min, double q_max)
1843 {
1844  updateParameters();
1845 
1846  std::pair<double, double > qbin = std::make_pair(q_min, q_max);
1847 
1848  old_handler = gsl_set_error_handler_off();
1849 
1850  switch (i) {
1851  case 0:
1852  if (sigma0Cached[qbin] == 0) {
1853  FS = convertToGslFunction(boost::bind(&MVll::getSigma1c, &(*this), _1));
1854  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();
1855  cacheSigma0[qbin] = avaSigma;
1856  sigma0Cached[qbin] = 1;
1857  }
1858  return cacheSigma0[qbin];
1859  break;
1860  case 1:
1861  if (sigma1Cached[qbin] == 0) {
1862  FS = convertToGslFunction(boost::bind(&MVll::getSigma1s, &(*this), _1));
1863  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();
1864  cacheSigma1[qbin] = avaSigma;
1865  sigma1Cached[qbin] = 1;
1866  }
1867  return cacheSigma1[qbin];
1868  break;
1869  case 2:
1870  if (sigma2Cached[qbin] == 0) {
1871  FS = convertToGslFunction(boost::bind(&MVll::getSigma2c, &(*this), _1));
1872  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();
1873  cacheSigma2[qbin] = avaSigma;
1874  sigma2Cached[qbin] = 1;
1875  }
1876  return cacheSigma2[qbin];
1877  break;
1878  case 3:
1879  if (sigma3Cached[qbin] == 0) {
1880  FS = convertToGslFunction(boost::bind(&MVll::getSigma2s, &(*this), _1));
1881  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();
1882  cacheSigma3[qbin] = avaSigma;
1883  sigma3Cached[qbin] = 1;
1884  }
1885  return cacheSigma3[qbin];
1886  break;
1887  case 4:
1888  if (sigma4Cached[qbin] == 0) {
1889  FS = convertToGslFunction(boost::bind(&MVll::getSigma3, &(*this), _1));
1890  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();
1891  cacheSigma4[qbin] = avaSigma;
1892  sigma4Cached[qbin] = 1;
1893  }
1894  return cacheSigma4[qbin];
1895  break;
1896  case 5:
1897  if (sigma5Cached[qbin] == 0) {
1898  FS = convertToGslFunction(boost::bind(&MVll::getSigma4, &(*this), _1));
1899  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();
1900  cacheSigma5[qbin] = avaSigma;
1901  sigma5Cached[qbin] = 1;
1902  }
1903  return cacheSigma5[qbin];
1904  break;
1905  case 6:
1906  if (sigma6Cached[qbin] == 0) {
1907  FS = convertToGslFunction(boost::bind(&MVll::getSigma5, &(*this), _1));
1908  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();
1909  cacheSigma6[qbin] = avaSigma;
1910  sigma6Cached[qbin] = 1;
1911  }
1912  return cacheSigma6[qbin];
1913  break;
1914  case 7:
1915  if (sigma7Cached[qbin] == 0) {
1916  FS = convertToGslFunction(boost::bind(&MVll::getSigma6s, &(*this), _1));
1917  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();
1918  cacheSigma7[qbin] = avaSigma;
1919  sigma7Cached[qbin] = 1;
1920  }
1921  return cacheSigma7[qbin];
1922  break;
1923  case 9:
1924  if (sigma9Cached[qbin] == 0) {
1925  FS = convertToGslFunction(boost::bind(&MVll::getSigma7, &(*this), _1));
1926  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();
1927  cacheSigma9[qbin] = avaSigma;
1928  sigma9Cached[qbin] = 1;
1929  }
1930  return cacheSigma9[qbin];
1931  break;
1932  case 10:
1933  if (sigma10Cached[qbin] == 0) {
1934  FS = convertToGslFunction(boost::bind(&MVll::getSigma8, &(*this), _1));
1935  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();
1936  cacheSigma10[qbin] = avaSigma;
1937  sigma10Cached[qbin] = 1;
1938  }
1939  return cacheSigma10[qbin];
1940  break;
1941  case 11:
1942  if (sigma11Cached[qbin] == 0) {
1943  FS = convertToGslFunction(boost::bind(&MVll::getSigma9, &(*this), _1));
1944  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();
1945  cacheSigma11[qbin] = avaSigma;
1946  sigma11Cached[qbin] = 1;
1947  }
1948  return cacheSigma11[qbin];
1949  break;
1950  default:
1951  std::stringstream out;
1952  out << i;
1953  throw std::runtime_error("MVll::integrateSigma: index " + out.str() + " not implemented");
1954  }
1955 
1956  gsl_set_error_handler(old_handler);
1957 
1958 }
1959 
1960 double MVll::getSigma(int i, double q_2)
1961 {
1962  updateParameters();
1963 
1964  switch (i) {
1965  case 0:
1966  return getSigma1c(q_2);
1967  break;
1968  case 1:
1969  return getSigma1s(q_2);
1970  break;
1971  case 2:
1972  return getSigma2c(q_2);
1973  break;
1974  case 3:
1975  return getSigma2s(q_2);
1976  break;
1977  case 4:
1978  return getSigma3(q_2);
1979  break;
1980  case 5:
1981  return getSigma4(q_2);
1982  break;
1983  case 6:
1984  return getSigma5(q_2);
1985  break;
1986  case 7:
1987  return getSigma6s(q_2);
1988  break;
1989  case 9:
1990  return getSigma7(q_2);
1991  break;
1992  case 10:
1993  return getSigma8(q_2);
1994  break;
1995  case 11:
1996  return getSigma9(q_2);
1997  break;
1998  default:
1999  std::stringstream out;
2000  out << i;
2001  throw std::runtime_error("MVll::getSigma: index " + out.str() + " not implemented");
2002  }
2003 }
2004 
2005 double MVll::integrateDelta(int i, double q_min, double q_max)
2006 {
2007  updateParameters();
2008 
2009  std::pair<double, double > qbin = std::make_pair(q_min, q_max);
2010 
2011  old_handler = gsl_set_error_handler_off();
2012 
2013  switch (i) {
2014  case 0:
2015  if (delta0Cached[qbin] == 0) {
2016  FD = convertToGslFunction(boost::bind(&MVll::getDelta1c, &(*this), _1));
2017  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();
2018  cacheDelta0[qbin] = avaDelta;
2019  delta0Cached[qbin] = 1;
2020  }
2021  return cacheDelta0[qbin];
2022  break;
2023  case 1:
2024  if (delta1Cached[qbin] == 0) {
2025  FD = convertToGslFunction(boost::bind(&MVll::getDelta1s, &(*this), _1));
2026  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();
2027  cacheDelta1[qbin] = avaDelta;
2028  delta1Cached[qbin] = 1;
2029  }
2030  return cacheDelta1[qbin];
2031  break;
2032  case 2:
2033  if (delta2Cached[qbin] == 0) {
2034  FD = convertToGslFunction(boost::bind(&MVll::getDelta2c, &(*this), _1));
2035  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();
2036  cacheDelta2[qbin] = avaDelta;
2037  delta2Cached[qbin] = 1;
2038  }
2039  return cacheDelta2[qbin];
2040  break;
2041  case 3:
2042  if (delta3Cached[qbin] == 0) {
2043  FD = convertToGslFunction(boost::bind(&MVll::getDelta2s, &(*this), _1));
2044  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();
2045  cacheDelta3[qbin] = avaDelta;
2046  delta3Cached[qbin] = 1;
2047  }
2048  return cacheDelta3[qbin];
2049  break;
2050  case 6:
2051  if (delta6Cached[qbin] == 0) {
2052  FD = convertToGslFunction(boost::bind(&MVll::getDelta5, &(*this), _1));
2053  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();
2054  cacheDelta6[qbin] = avaDelta;
2055  delta6Cached[qbin] = 1;
2056  }
2057  return cacheDelta6[qbin];
2058  break;
2059  case 7:
2060  if (delta7Cached[qbin] == 0) {
2061  FD = convertToGslFunction(boost::bind(&MVll::getDelta6s, &(*this), _1));
2062  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();
2063  cacheDelta7[qbin] = avaDelta;
2064  delta7Cached[qbin] = 1;
2065  }
2066  return cacheDelta7[qbin];
2067  break;
2068  case 8:
2069  if (delta8Cached[qbin] == 0) {
2070  FD = convertToGslFunction(boost::bind(&MVll::getDelta6c, &(*this), _1));
2071  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();
2072  cacheDelta8[qbin] = avaDelta;
2073  delta8Cached[qbin] = 1;
2074  }
2075  return cacheDelta8[qbin];
2076  break;
2077  case 10:
2078  if (delta10Cached[qbin] == 0) {
2079  FD = convertToGslFunction(boost::bind(&MVll::getDelta8, &(*this), _1));
2080  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();
2081  cacheDelta10[qbin] = avaDelta;
2082  delta10Cached[qbin] = 1;
2083  }
2084  return cacheDelta10[qbin];
2085  break;
2086  case 11:
2087  if (delta11Cached[qbin] == 0) {
2088  FD = convertToGslFunction(boost::bind(&MVll::getDelta9, &(*this), _1));
2089  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();
2090  cacheDelta11[qbin] = avaDelta;
2091  delta11Cached[qbin] = 1;
2092  }
2093  return cacheDelta11[qbin];
2094  break;
2095  default:
2096  std::stringstream out;
2097  out << i;
2098  throw std::runtime_error("MVll::integrateDelta: index " + out.str() + " not implemented");
2099  }
2100 
2101  gsl_set_error_handler(old_handler);
2102 
2103 }
MVll::H_A_0
gslpp::complex H_A_0(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:1644
MVll::H_A_m
gslpp::complex H_A_m(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:1656
std_make_vector.h
QCD::BOTTOM
Definition: QCD.h:329
gslpp::arccot
complex arccot(const complex &z)
Definition: gslpp_complex.cpp:519
MVll::GF
double GF
Definition: MVll.h:738
gslpp::arctan
complex arctan(const complex &z)
Definition: gslpp_complex.cpp:492
MVll::vectorM
QCD::meson vectorM
Definition: MVll.h:729
MVll::initializeMVllParameters
std::vector< std::string > initializeMVllParameters()
A method for initializing the parameters necessary for MVll.
Definition: MVll.cpp:149
gslpp::dilog
complex dilog(const complex &z)
Definition: gslpp_complex.cpp:370
MVll::MVll
MVll(const StandardModel &SM_i, QCD::meson meson_i, QCD::meson vector_i, QCD::lepton lep_i)
Constructor.
Definition: MVll.cpp:21
MVll::MV
double MV
Definition: MVll.h:742
QCD::B_S
Definition: QCD.h:345
MVll::H_P
gslpp::complex H_P(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:1668
MVll::h_lambda
gslpp::complex h_lambda(int hel, double q2)
The non-pertubative ccbar contributions to the helicity amplitudes.
Definition: MVll.cpp:1602
F_2.h
CKM::computelamt_s
gslpp::complex computelamt_s() const
The product of the CKM elements .
Definition: CKM.cpp:136
MVll::ale
double ale
Definition: MVll.h:739
MVll::H_V_p
gslpp::complex H_V_p(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:1632
QCD::initializeMeson
void initializeMeson(QCD::meson meson_i) const
A method to initialize a meson.
Definition: QCD.cpp:236
MVll::exp_Phase
gslpp::complex exp_Phase[3]
Definition: MVll.h:736
make_vector
Definition: std_make_vector.h:15
LO
Definition: OrderScheme.h:33
QCD::UP
Definition: QCD.h:324
StandardModel.h
F_1.h
MVll::xs
double xs
Definition: MVll.h:753
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
MVll::getSigma
double getSigma(int i, double q_2)
The value of from to .
Definition: MVll.cpp:1956
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
gslpp::complex::abs2
double abs2() const
Definition: gslpp_complex.cpp:86
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:474
MVll::mc_pole
double mc_pole
Definition: MVll.h:748
MVll::mb_pole
double mb_pole
Definition: MVll.h:747
MVll::H_V_0
gslpp::complex H_V_0(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:1625
MVll::lep
QCD::lepton lep
Definition: MVll.h:727
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
MVll::mu_h
double mu_h
Definition: MVll.h:745
MVll::mySM
const StandardModel & mySM
Definition: MVll.h:726
MVll::myF_1
std::unique_ptr< F_1 > myF_1
Definition: MVll.h:731
MVll.h
CKM::computelamu_s
gslpp::complex computelamu_s() const
The product of the CKM elements .
Definition: CKM.cpp:146
MVll::H_A_p
gslpp::complex H_A_p(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:1650
MVll::etaV
int etaV
Definition: MVll.h:755
gslpp::complex::imag
const double & imag() const
Definition: gslpp_complex.cpp:59
MVll::H_V_m
gslpp::complex H_V_m(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:1638
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
MVll::ys
double ys
Definition: MVll.h:752
QCD::K_star
Definition: QCD.h:348
QCD::meson
meson
An enum type for mesons.
Definition: QCD.h:336
MVll::Mlep
double Mlep
Definition: MVll.h:740
QCD::B_P
Definition: QCD.h:344
F_2
Definition: F_2.h:15
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
MVll::mJ2
double mJ2
Definition: MVll.h:735
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
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
MVll::meson
QCD::meson meson
Definition: MVll.h:728
Meson::getDecayconst
const double & getDecayconst() const
A get method for the decay constant of the meson.
Definition: Meson.h:360
MVll::beta
double beta(double q2)
The factor used in the angular coefficients .
Definition: MVll.cpp:1681
MVll::~MVll
virtual ~MVll()
Destructor.
Definition: MVll.cpp:145
MVll::H_S
gslpp::complex H_S(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:1662
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
Meson::getDecayconst_p
const double & getDecayconst_p() const
A get method for the perpendicular decay constant of a vector meson.
Definition: Meson.h:378
StandardModel::getFlavour
const Flavour & getFlavour() const
Definition: StandardModel.h:1006
MVll::integrateSigma
double integrateSigma(int i, double q_min, double q_max)
The integral of from to .
Definition: MVll.cpp:1838
C_4
Definition: Doxygen/examples-src/myModel/src/myObservables.h:76
F_1
Definition: F_1.h:15
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
MVll::FixedWCbtos
bool FixedWCbtos
Definition: MVll.h:734
QCD::K_star_P
Definition: QCD.h:349
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::STRANGE
Definition: QCD.h:327
gslpp::complex::real
const double & real() const
Definition: gslpp_complex.cpp:53
MVll::width
double width
Definition: MVll.h:751
QCD::PHI
Definition: QCD.h:347
MVll::Mb
double Mb
Definition: MVll.h:743
MVll::myF_2
std::unique_ptr< F_2 > myF_2
Definition: MVll.h:732
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
MVll::angmomV
double angmomV
Definition: MVll.h:754
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
NLO
Definition: OrderScheme.h:34
Meson::getDgamma_gamma
const double & getDgamma_gamma() const
Definition: Meson.h:411
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
MVll::mu_b
double mu_b
Definition: MVll.h:744
MVll::spectator_charge
double spectator_charge
Definition: MVll.h:750
MVll::integrateDelta
double integrateDelta(int i, double q_min, double q_max)
The integral of from to .
Definition: MVll.cpp:2001
C_3
Definition: Doxygen/examples-src/myModel/src/myObservables.h:59
Meson::getGegenalpha
const double & getGegenalpha(int i) const
A get method to get the Gegenbaur coefficient.
Definition: Meson.h:394
MVll::MM
double MM
Definition: MVll.h:741
MVll::dispersion
bool dispersion
Definition: MVll.h:733
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
MVll::Ms
double Ms
Definition: MVll.h:749
Flavour::getFlagUseDispersionRelation
bool getFlagUseDispersionRelation() const
Definition: Flavour.h:242
Flavour::getFlagFixedWCbtos
bool getFlagFixedWCbtos() const
Definition: Flavour.h:252
MVll::mvllParameters
std::vector< std::string > mvllParameters
Definition: MVll.h:730
StandardModel::getCKM
CKM getCKM() const
A get method to retrieve the member object of type CKM.
Definition: StandardModel.h:876
MVll::Mc
double Mc
Definition: MVll.h:746
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