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/bind.hpp>
16#include <limits>
17#include <TFitResult.h>
18#include <gsl/gsl_sf_gegenbauer.h>
19#include <gsl/gsl_sf_expint.h>
20using namespace boost::placeholders;
21
22MVll::MVll(const StandardModel& SM_i, QCD::meson meson_i, QCD::meson vector_i, QCD::lepton lep_i)
23: mySM(SM_i), myF_1(new F_1()), myF_2(new F_2()),
24N_cache(3, 0.),
25V_cache(3, 0.),
26A0_cache(3, 0.),
27A1_cache(3, 0.),
28T1_cache(3, 0.),
29T2_cache(3, 0.),
30k2_cache(2, 0.),
31VL0_cache(3, 0.),
32TL0_cache(3, 0.),
33SL_cache(2, 0.),
34Ycache(2, 0.),
35H_V0cache(2, 0.),
36H_V1cache(2, 0.),
37H_V2cache(2, 0.),
38H_Scache(2, 0.),
39H_Pcache(4, 0.),
40Itree_cache(3, 0.),
41T_cache(5, 0.)
42{
43 lep = lep_i;
44 meson = meson_i;
45 vectorM = vector_i;
46 dispersion = false;
47 zExpansion = false;
48 FixedWCbtos = false;
49 NeutrinoTree_flag = false;
50 MVll_DM_flag = false;
51 mJpsi = 3.0969;
52 mJ2 = mJpsi * mJpsi;
53 mPsi2S = 3.6861;
54 mPsi2S2 = mPsi2S * mPsi2S;
55 mD2 = 1.8648 * 1.8648;
56
57 I0_updated = 0;
58 I1_updated = 0;
59 I2_updated = 0;
60 I3_updated = 0;
61 I4_updated = 0;
62 I5_updated = 0;
63 I6_updated = 0;
64 I7_updated = 0;
65 I8_updated = 0;
66 I9_updated = 0;
67 I10_updated = 0;
68 I11_updated = 0;
69 Itree_updated = 0;
70
71 VL1_updated = 0;
72 VL2_updated = 0;
73 TL1_updated = 0;
74 TL2_updated = 0;
75 VR1_updated = 0;
76 VR2_updated = 0;
77 TR1_updated = 0;
78 TR2_updated = 0;
79 VL0_updated = 0;
80 TL0_updated = 0;
81 VR0_updated = 0;
82 TR0_updated = 0;
83 SL_updated = 0;
84 SR_updated = 0;
85
86 deltaTparpupdated = 0;
87 deltaTparmupdated = 0;
88 deltaTperpupdated = 0;
89
90 w_sigma = gsl_integration_cquad_workspace_alloc(100);
91 // w_DTPPR = gsl_integration_cquad_workspace_alloc (100);
92 w_sigmaTree = gsl_integration_cquad_workspace_alloc(100);
93 w_delta = gsl_integration_cquad_workspace_alloc(100);
94
95 acc_Re_T_perp = gsl_interp_accel_alloc();
96 acc_Im_T_perp = gsl_interp_accel_alloc();
97 acc_Re_T_para = gsl_interp_accel_alloc();
98 acc_Im_T_para = gsl_interp_accel_alloc();
99
100 spline_Re_T_perp = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
101 spline_Im_T_perp = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
102 spline_Re_T_para = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
103 spline_Im_T_para = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
104
105#if COMPUTECP
106 acc_Re_T_perp_conj = gsl_interp_accel_alloc();
107 acc_Im_T_perp_conj = gsl_interp_accel_alloc();
108 acc_Re_T_para_conj = gsl_interp_accel_alloc();
109 acc_Im_T_para_conj = gsl_interp_accel_alloc();
110
111 spline_Re_T_perp_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
112 spline_Im_T_perp_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
113 spline_Re_T_para_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
114 spline_Im_T_para_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM);
115#endif
116
117 acc_Re_deltaC7_QCDF = gsl_interp_accel_alloc();
118 acc_Im_deltaC7_QCDF = gsl_interp_accel_alloc();
119 acc_Re_deltaC9_QCDF = gsl_interp_accel_alloc();
120 acc_Im_deltaC9_QCDF = gsl_interp_accel_alloc();
121
122 spline_Re_deltaC7_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
123 spline_Im_deltaC7_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
124 spline_Re_deltaC9_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
125 spline_Im_deltaC9_QCDF = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
126
127#if COMPUTECP
128 acc_Re_deltaC7_QCDF_conj = gsl_interp_accel_alloc();
129 acc_Im_deltaC7_QCDF_conj = gsl_interp_accel_alloc();
130 acc_Re_deltaC9_QCDF_conj = gsl_interp_accel_alloc();
131 acc_Im_deltaC9_QCDF_conj = gsl_interp_accel_alloc();
132
133 spline_Re_deltaC7_QCDF_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
134 spline_Im_deltaC7_QCDF_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
135 spline_Re_deltaC9_QCDF_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
136 spline_Im_deltaC9_QCDF_conj = gsl_spline_alloc(gsl_interp_cspline, GSL_INTERP_DIM_DC);
137#endif
138
139 h_pole = false;
140
141 M_PI2 = M_PI*M_PI;
142
143 F87_1 = (4. / 3. * M_PI2 - 40. / 3.);
144 F87_2 = (32. / 9. * M_PI2 - 316. / 9.);
145 F87_3 = (200. / 27. * M_PI2 - 658. / 9.);
146
147 F89_0 = (104. / 9. - 32. / 27. * M_PI2);
148 F89_1 = (1184. / 27. - 40. / 9. * M_PI2);
149 F89_2 = (-32. / 3. * M_PI2 + 14212. / 135.);
150 F89_3 = (-560. / 27. * M_PI2 + 193444. / 945.);
151
152 CF = 4. / 3.;
153
154}
155
156MVll::~MVll()
157{
158}
159
160std::vector<std::string> MVll::initializeMVllParameters()
161{
162 dispersion = mySM.getFlavour().getFlagUseDispersionRelation();
163 zExpansion = mySM.getFlavour().getFlagUsezExpansion();
164 FixedWCbtos = mySM.getFlavour().getFlagFixedWCbtos();
165 NeutrinoTree_flag = mySM.getFlavour().getFlagNeutrinoTree();
166 MVll_DM_flag = mySM.getFlavour().getFlagMVll_DM();
167
168#if NFPOLARBASIS_MVLL
169 if (vectorM == StandardModel::PHI)
170 if (MVll_DM_flag) mvllParameters = make_vector<std::string>()
171 << "a_0fphi" << "a_1fphi" << "a_2fphi" << "MRf" << "a_0gphi" << "a_1gphi" << "a_2gphi" << "MRg"
172 << "a_1F1phi" << "a_2F1phi" << "MRF1" << "a_1F2phi" << "a_2F2phi" << "MRF2" /*a_0F1 and a_0F2 are not independent*/
173 << "a_0T1phi" << "a_1T1phi" << "a_2T1phi" << "MRT1" << "a_1T2phi" << "a_2T2phi" << "MRT2"
174 << "a_1T0phi" << "a_2T0phi" << "MRT0" /*a_0T0 and a_0T2 are not independent*/
175 << "Chi1minus" << "Chi1plus" << "Chi0plus" << "Chi0minus" << "ChiTT" << "ChiBB"
176 << "absh_0" << "absh_p" << "absh_m" << "argh_0" << "argh_p" << "argh_m"
177 << "absh_0_1" << "absh_p_1" << "absh_m_1" << "argh_0_1" << "argh_p_1" << "argh_m_1"
178 << "absh_p_2" << "absh_m_2" << "argh_p_2" << "argh_m_2" << "xs_phi" << "SU3_breaking_abs" << "SU3_breaking_arg";
179 else mvllParameters = make_vector<std::string>()
180 << "a_0Vphi" << "a_1Vphi" << "a_2Vphi" << "MRV" << "a_0A0phi" << "a_1A0phi" << "a_2A0phi" << "MRA0"
181 << "a_0A1phi" << "a_1A1phi" << "a_2A1phi" << "MRA1" << "a_1A12phi" << "a_2A12phi" << "MRA12" /*a_0A12 and a_0T2 are not independent*/
182 << "a_0T1phi" << "a_1T1phi" << "a_2T1phi" << "MRT1" << "a_1T2phi" << "a_2T2phi" << "MRT2"
183 << "a_0T23phi" << "a_1T23phi" << "a_2T23phi" << "MRT23"
184 << "absh_0" << "absh_p" << "absh_m" << "argh_0" << "argh_p" << "argh_m"
185 << "absh_0_1" << "absh_p_1" << "absh_m_1" << "argh_0_1" << "argh_p_1" << "argh_m_1"
186 << "absh_p_2" << "absh_m_2" << "argh_p_2" << "argh_m_2" << "xs_phi" << "SU3_breaking_abs" << "SU3_breaking_arg";
187 else if (vectorM == StandardModel::K_star || vectorM == StandardModel::K_star_P)
188 if (MVll_DM_flag) mvllParameters = make_vector<std::string>()
189 << "a_0f" << "a_1f" << "a_2f" << "MRf" << "a_0g" << "a_1g" << "a_2g" << "MRg"
190 << "a_1F1" << "a_2F1" << "MRF1" << "a_1F2" << "a_2F2" << "MRF2" /*a_0F1 and a_0F2 are not independent*/
191 << "a_0T1" << "a_1T1" << "a_2T1" << "MRT1" << "a_1T2" << "a_2T2" << "MRT2"
192 << "a_1T0" << "a_2T0" << "MRT0" /*a_0T0 and a_0T2 are not independent*/
193 << "Chi1minus" << "Chi1plus" << "Chi0plus" << "Chi0minus" << "ChiTT" << "ChiBB"
194 << "absh_0" << "absh_p" << "absh_m" << "argh_0" << "argh_p" << "argh_m"
195 << "absh_0_1" << "absh_p_1" << "absh_m_1" << "argh_0_1" << "argh_p_1" << "argh_m_1"
196 << "absh_p_2" << "absh_m_2" << "argh_p_2" << "argh_m_2";
197 else mvllParameters = make_vector<std::string>()
198 << "a_0V" << "a_1V" << "a_2V" << "MRV" << "a_0A0" << "a_1A0" << "a_2A0" << "MRA0"
199 << "a_0A1" << "a_1A1" << "a_2A1" << "MRA1" << "a_1A12" << "a_2A12" << "MRA12" /*a_0A12 and a_0T2 are not independent*/
200 << "a_0T1" << "a_1T1" << "a_2T1" << "MRT1" << "a_1T2" << "a_2T2" << "MRT2"
201 << "a_0T23" << "a_1T23" << "a_2T23" << "MRT23"
202 << "absh_0" << "absh_p" << "absh_m" << "argh_0" << "argh_p" << "argh_m"
203 << "absh_0_1" << "absh_p_1" << "absh_m_1" << "argh_0_1" << "argh_p_1" << "argh_m_1"
204 << "absh_p_2" << "absh_m_2" << "argh_p_2" << "argh_m_2";
205#else
206 if (vectorM == StandardModel::PHI)
207 if (MVll_DM_flag) mvllParameters = make_vector<std::string>()
208 << "a_0fphi" << "a_1fphi" << "a_2fphi" << "MRf" << "a_0gphi" << "a_1gphi" << "a_2gphi" << "MRg"
209 << "a_1F1phi" << "a_2F1phi" << "MRF1" << "a_1F2phi" << "a_2F2phi" << "MRF2" /*a_0F1 and a_0F2 are not independent*/
210 << "a_0T1phi" << "a_1T1phi" << "a_2T1phi" << "MRT1" << "a_1T2phi" << "a_2T2phi" << "MRT2"
211 << "a_1T0phi" << "a_2T0phi" << "MRT0" /*a_0T0 and a_0T2 are not independent*/
212 << "Chi1minus" << "Chi1plus" << "Chi0plus" << "Chi0minus" << "ChiTT" << "ChiBB"
213 << "reh_0" << "reh_p" << "reh_m" << "imh_0" << "imh_p" << "imh_m"
214 << "reh_0_1" << "reh_p_1" << "reh_m_1" << "imh_0_1" << "imh_p_1" << "imh_m_1"
215 << "reh_p_2" << "reh_m_2" << "imh_p_2" << "imh_m_2" << "xs_phi" << "SU3_breaking_abs" << "SU3_breaking_arg";
216 else mvllParameters = make_vector<std::string>()
217 << "a_0Vphi" << "a_1Vphi" << "a_2Vphi" << "MRV" << "a_0A0phi" << "a_1A0phi" << "a_2A0phi" << "MRA0"
218 << "a_0A1phi" << "a_1A1phi" << "a_2A1phi" << "MRA1" << "a_1A12phi" << "a_2A12phi" << "MRA12" /*a_0A12 and a_0T2 are not independent*/
219 << "a_0T1phi" << "a_1T1phi" << "a_2T1phi" << "MRT1" << "a_1T2phi" << "a_2T2phi" << "MRT2"
220 << "a_0T23phi" << "a_1T23phi" << "a_2T23phi" << "MRT23"
221 << "reh_0" << "reh_p" << "reh_m" << "imh_0" << "imh_p" << "imh_m"
222 << "reh_0_1" << "reh_p_1" << "reh_m_1" << "imh_0_1" << "imh_p_1" << "imh_m_1"
223 << "reh_p_2" << "reh_m_2" << "imh_p_2" << "imh_m_2" << "xs_phi" << "SU3_breaking_abs" << "SU3_breaking_arg";
224 else if (vectorM == StandardModel::K_star || vectorM == StandardModel::K_star_P)
225 if (MVll_DM_flag) mvllParameters = make_vector<std::string>()
226 << "a_0f" << "a_1f" << "a_2f" << "MRf" << "a_0g" << "a_1g" << "a_2g" << "MRg"
227 << "a_1F1" << "a_2F1" << "MRF1" << "a_1F2" << "a_2F2" << "MRF2" /*a_0F1 and a_0F2 are not independent*/
228 << "a_0T1" << "a_1T1" << "a_2T1" << "MRT1" << "a_1T2" << "a_2T2" << "MRT2"
229 << "a_1T0" << "a_2T0" << "MRT0" /*a_0T0 and a_0T2 are not independent*/
230 << "Chi1minus" << "Chi1plus" << "Chi0plus" << "Chi0minus" << "ChiTT" << "ChiBB"
231 << "reh_0" << "reh_p" << "reh_m" << "imh_0" << "imh_p" << "imh_m"
232 << "reh_0_1" << "reh_p_1" << "reh_m_1" << "imh_0_1" << "imh_p_1" << "imh_m_1"
233 << "reh_p_2" << "reh_m_2" << "imh_p_2" << "imh_m_2";
234 else mvllParameters = make_vector<std::string>()
235 << "a_0V" << "a_1V" << "a_2V" << "MRV" << "a_0A0" << "a_1A0" << "a_2A0" << "MRA0"
236 << "a_0A1" << "a_1A1" << "a_2A1" << "MRA1" << "a_1A12" << "a_2A12" << "MRA12" /*a_0A12 and a_0T2 are not independent*/
237 << "a_0T1" << "a_1T1" << "a_2T1" << "MRT1" << "a_1T2" << "a_2T2" << "MRT2"
238 << "a_0T23" << "a_1T23" << "a_2T23" << "MRT23"
239 << "reh_0" << "reh_p" << "reh_m" << "imh_0" << "imh_p" << "imh_m"
240 << "reh_0_1" << "reh_p_1" << "reh_m_1" << "imh_0_1" << "imh_p_1" << "imh_m_1"
241 << "reh_p_2" << "reh_m_2" << "imh_p_2" << "imh_m_2";
242#endif
243 else {
244 std::stringstream out;
245 out << vectorM;
246 throw std::runtime_error("MVll: vector " + out.str() + " not implemented");
247 }
248
249 if (dispersion) {
250 mvllParameters.clear();
251 if (vectorM == StandardModel::PHI)
252 if (MVll_DM_flag) mvllParameters = make_vector<std::string>()
253 << "a_0fphi" << "a_1fphi" << "a_2fphi" << "MRf" << "a_0gphi" << "a_1gphi" << "a_2gphi" << "MRg"
254 << "a_1F1phi" << "a_2F1phi" << "MRF1" << "a_1F2phi" << "a_2F2phi" << "MRF2" /*a_0F1 and a_0F2 are not independent*/
255 << "a_0T1phi" << "a_1T1phi" << "a_2T1phi" << "MRT1" << "a_1T2phi" << "a_2T2phi" << "MRT2"
256 << "a_1T0phi" << "a_2T0phi" << "MRT0" /*a_0T0 and a_0T2 are not independent*/
257 << "Chi1minus" << "Chi1plus" << "Chi0plus" << "Chi0minus" << "ChiTT" << "ChiBB"
258 << "r1_1" << "r2_1" << "deltaC9_1" << "phDC9_1"
259 << "r1_2" << "r2_2" << "deltaC9_2" << "phDC9_2"
260 << "r1_3" << "r2_3" << "deltaC9_3" << "phDC9_3" << "xs_phi" << "SU3_breaking_abs" << "SU3_breaking_arg";
261 else mvllParameters = make_vector<std::string>()
262 << "a_0Vphi" << "a_1Vphi" << "a_2Vphi" << "MRV" << "a_0A0phi" << "a_1A0phi" << "a_2A0phi" << "MRA0"
263 << "a_0A1phi" << "a_1A1phi" << "a_2A1phi" << "MRA1" << "a_1A12phi" << "a_2A12phi" << "MRA12" /*a_0A12 and a_0T2 are not independent*/
264 << "a_0T1phi" << "a_1T1phi" << "a_2T1phi" << "MRT1" << "a_1T2phi" << "a_2T2phi" << "MRT2"
265 << "a_0T23phi" << "a_1T23phi" << "a_2T23phi" << "MRT23"
266 << "r1_1" << "r2_1" << "deltaC9_1" << "phDC9_1"
267 << "r1_2" << "r2_2" << "deltaC9_2" << "phDC9_2"
268 << "r1_3" << "r2_3" << "deltaC9_3" << "phDC9_3" << "xs_phi" << "SU3_breaking_abs" << "SU3_breaking_arg";
269 else if (vectorM == StandardModel::K_star || vectorM == StandardModel::K_star_P)
270 if (MVll_DM_flag) mvllParameters = make_vector<std::string>()
271 << "a_0f" << "a_1f" << "a_2f" << "MRf" << "a_0g" << "a_1g" << "a_2g" << "MRg"
272 << "a_1F1" << "a_2F1" << "MRF1" << "a_1F2" << "a_2F2" << "MRF2" /*a_0F1 and a_0F2 are not independent*/
273 << "a_0T1" << "a_1T1" << "a_2T1" << "MRT1" << "a_1T2" << "a_2T2" << "MRT2"
274 << "a_1T0" << "a_2T0" << "MRT0" /*a_0T0 and a_0T2 are not independent*/
275 << "Chi1minus" << "Chi1plus" << "Chi0plus" << "Chi0minus" << "ChiTT" << "ChiBB"
276 << "r1_1" << "r2_1" << "deltaC9_1" << "phDC9_1"
277 << "r1_2" << "r2_2" << "deltaC9_2" << "phDC9_2"
278 << "r1_3" << "r2_3" << "deltaC9_3" << "phDC9_3";
279 else mvllParameters = make_vector<std::string>()
280 << "a_0V" << "a_1V" << "a_2V" << "MRV" << "a_0A0" << "a_1A0" << "a_2A0" << "MRA0"
281 << "a_0A1" << "a_1A1" << "a_2A1" << "MRA1" << "a_1A12" << "a_2A12" << "MRA12" /*a_0A12 and a_0T2 are not independent*/
282 << "a_0T1" << "a_1T1" << "a_2T1" << "MRT1" << "a_1T2" << "a_2T2" << "MRT2"
283 << "a_0T23" << "a_1T23" << "a_2T23" << "MRT23"
284 << "r1_1" << "r2_1" << "deltaC9_1" << "phDC9_1"
285 << "r1_2" << "r2_2" << "deltaC9_2" << "phDC9_2"
286 << "r1_3" << "r2_3" << "deltaC9_3" << "phDC9_3";
287 }
288
289 if (zExpansion) {
290 mvllParameters.clear();
291 if (vectorM == StandardModel::PHI)
292 if (MVll_DM_flag) mvllParameters = make_vector<std::string>()
293 << "a_0fphi" << "a_1fphi" << "a_2fphi" << "MRf" << "a_0gphi" << "a_1gphi" << "a_2gphi" << "MRg"
294 << "a_1F1phi" << "a_2F1phi" << "MRF1" << "a_1F2phi" << "a_2F2phi" << "MRF2" /*a_0F1 and a_0F2 are not independent*/
295 << "a_0T1phi" << "a_1T1phi" << "a_2T1phi" << "MRT1" << "a_1T2phi" << "a_2T2phi" << "MRT2"
296 << "a_1T0phi" << "a_2T0phi" << "MRT0" /*a_0T0 and a_0T2 are not independent*/
297 << "Chi1minus" << "Chi1plus" << "Chi0plus" << "Chi0minus" << "ChiTT" << "ChiBB"
298 << "DeltaC9" << "DeltaC10"
299 << "re_beta_0_0" << "re_beta_0_1" << "re_beta_0_2" << "re_beta_0_3" << "re_beta_0_4" << "re_beta_0_5" << "re_beta_0_6"
300 << "im_beta_0_0" << "im_beta_0_1" << "im_beta_0_2" << "im_beta_0_3" << "im_beta_0_4" << "im_beta_0_5" << "im_beta_0_6"
301 << "re_beta_1_0" << "re_beta_1_1" << "re_beta_1_2" << "re_beta_1_3" << "re_beta_1_4" << "re_beta_1_5" << "re_beta_1_6"
302 << "im_beta_1_0" << "im_beta_1_1" << "im_beta_1_2" << "im_beta_1_3" << "im_beta_1_4" << "im_beta_1_5" << "im_beta_1_6"
303 << "re_beta_2_0" << "re_beta_2_1" << "re_beta_2_2" << "re_beta_2_3" << "re_beta_2_4" << "re_beta_2_5" << "re_beta_2_6"
304 << "im_beta_2_0" << "im_beta_2_1" << "im_beta_2_2" << "im_beta_2_3" << "im_beta_2_4" << "im_beta_2_5" << "im_beta_2_6"
305 << "xs_phi" << "SU3_breaking_abs" << "SU3_breaking_arg";
306 else mvllParameters = make_vector<std::string>()
307 << "a_0Vphi" << "a_1Vphi" << "a_2Vphi" << "MRV" << "a_0A0phi" << "a_1A0phi" << "a_2A0phi" << "MRA0"
308 << "a_0A1phi" << "a_1A1phi" << "a_2A1phi" << "MRA1" << "a_1A12phi" << "a_2A12phi" << "MRA12" /*a_0A12 and a_0T2 are not independent*/
309 << "a_0T1phi" << "a_1T1phi" << "a_2T1phi" << "MRT1" << "a_1T2phi" << "a_2T2phi" << "MRT2"
310 << "a_0T23phi" << "a_1T23phi" << "a_2T23phi" << "MRT23" << "DeltaC9" << "DeltaC10"
311 << "re_beta_0_0" << "re_beta_0_1" << "re_beta_0_2" << "re_beta_0_3" << "re_beta_0_4" << "re_beta_0_5" << "re_beta_0_6"
312 << "im_beta_0_0" << "im_beta_0_1" << "im_beta_0_2" << "im_beta_0_3" << "im_beta_0_4" << "im_beta_0_5" << "im_beta_0_6"
313 << "re_beta_1_0" << "re_beta_1_1" << "re_beta_1_2" << "re_beta_1_3" << "re_beta_1_4" << "re_beta_1_5" << "re_beta_1_6"
314 << "im_beta_1_0" << "im_beta_1_1" << "im_beta_1_2" << "im_beta_1_3" << "im_beta_1_4" << "im_beta_1_5" << "im_beta_1_6"
315 << "re_beta_2_0" << "re_beta_2_1" << "re_beta_2_2" << "re_beta_2_3" << "re_beta_2_4" << "re_beta_2_5" << "re_beta_2_6"
316 << "im_beta_2_0" << "im_beta_2_1" << "im_beta_2_2" << "im_beta_2_3" << "im_beta_2_4" << "im_beta_2_5" << "im_beta_2_6"
317 << "xs_phi" << "SU3_breaking_abs" << "SU3_breaking_arg";
318 else if (vectorM == StandardModel::K_star || vectorM == StandardModel::K_star_P)
319 if (MVll_DM_flag) mvllParameters = make_vector<std::string>()
320 << "a_0f" << "a_1f" << "a_2f" << "MRf" << "a_0g" << "a_1g" << "a_2g" << "MRg"
321 << "a_1F1" << "a_2F1" << "MRF1" << "a_1F2" << "a_2F2" << "MRF2" /*a_0F1 and a_0F2 are not independent*/
322 << "a_0T1" << "a_1T1" << "a_2T1" << "MRT1" << "a_1T2" << "a_2T2" << "MRT2"
323 << "a_1T0" << "a_2T0" << "MRT0" /*a_0T0 and a_0T2 are not independent*/
324 << "Chi1minus" << "Chi1plus" << "Chi0plus" << "Chi0minus" << "ChiTT" << "ChiBB"
325 << "DeltaC9" << "DeltaC10"
326 << "re_beta_0_0" << "re_beta_0_1" << "re_beta_0_2" << "re_beta_0_3" << "re_beta_0_4" << "re_beta_0_5" << "re_beta_0_6"
327 << "im_beta_0_0" << "im_beta_0_1" << "im_beta_0_2" << "im_beta_0_3" << "im_beta_0_4" << "im_beta_0_5" << "im_beta_0_6"
328 << "re_beta_1_0" << "re_beta_1_1" << "re_beta_1_2" << "re_beta_1_3" << "re_beta_1_4" << "re_beta_1_5" << "re_beta_1_6"
329 << "im_beta_1_0" << "im_beta_1_1" << "im_beta_1_2" << "im_beta_1_3" << "im_beta_1_4" << "im_beta_1_5" << "im_beta_1_6"
330 << "re_beta_2_0" << "re_beta_2_1" << "re_beta_2_2" << "re_beta_2_3" << "re_beta_2_4" << "re_beta_2_5" << "re_beta_2_6"
331 << "im_beta_2_0" << "im_beta_2_1" << "im_beta_2_2" << "im_beta_2_3" << "im_beta_2_4" << "im_beta_2_5" << "im_beta_2_6";
332 else mvllParameters = make_vector<std::string>()
333 << "a_0V" << "a_1V" << "a_2V" << "MRV" << "a_0A0" << "a_1A0" << "a_2A0" << "MRA0"
334 << "a_0A1" << "a_1A1" << "a_2A1" << "MRA1" << "a_1A12" << "a_2A12" << "MRA12" /*a_0A12 and a_0T2 are not independent*/
335 << "a_0T1" << "a_1T1" << "a_2T1" << "MRT1" << "a_1T2" << "a_2T2" << "MRT2"
336 << "a_0T23" << "a_1T23" << "a_2T23" << "MRT23" << "DeltaC9" << "DeltaC10"
337 << "re_beta_0_0" << "re_beta_0_1" << "re_beta_0_2" << "re_beta_0_3" << "re_beta_0_4" << "re_beta_0_5" << "re_beta_0_6"
338 << "im_beta_0_0" << "im_beta_0_1" << "im_beta_0_2" << "im_beta_0_3" << "im_beta_0_4" << "im_beta_0_5" << "im_beta_0_6"
339 << "re_beta_1_0" << "re_beta_1_1" << "re_beta_1_2" << "re_beta_1_3" << "re_beta_1_4" << "re_beta_1_5" << "re_beta_1_6"
340 << "im_beta_1_0" << "im_beta_1_1" << "im_beta_1_2" << "im_beta_1_3" << "im_beta_1_4" << "im_beta_1_5" << "im_beta_1_6"
341 << "re_beta_2_0" << "re_beta_2_1" << "re_beta_2_2" << "re_beta_2_3" << "re_beta_2_4" << "re_beta_2_5" << "re_beta_2_6"
342 << "im_beta_2_0" << "im_beta_2_1" << "im_beta_2_2" << "im_beta_2_3" << "im_beta_2_4" << "im_beta_2_5" << "im_beta_2_6";
343 }
344
345 if (FixedWCbtos)
346 if (lep != QCD::NEUTRINO_1) mvllParameters.insert(mvllParameters.end(), { "C7_SM", "C9_SM", "C10_SM" });
347 else mvllParameters.insert(mvllParameters.end(), { "CLnunu_SM" });
348
349 mySM.initializeMeson(meson);
350 mySM.initializeMeson(vectorM);
351 return mvllParameters;
352}
353
354void MVll::updateParameters()
355{
356 if (!mySM.getFlavour().getUpdateFlag(meson, vectorM, lep)) return;
357
358
359 GF = mySM.getGF();
360 ale = mySM.getAle();
361 if (lep == QCD::NEUTRINO_1){
362 Mlep = 0.;
363 }
364 else{
365 Mlep = mySM.getLeptons(lep).getMass();
366 }
367
368 MM = mySM.getMesons(meson).getMass();
369 MV = mySM.getMesons(vectorM).getMass();
370 mu_b = mySM.getMub();
371 mu_h = sqrt(mu_b * .5); // From Beneke Neubert
372 Mb = mySM.getQuarks(QCD::BOTTOM).getMass(); // add the PS b mass
373 Mc = mySM.getQuarks(QCD::CHARM).getMass();
374 mb_pole = mySM.Mbar2Mp(Mb, QCD::BOTTOM); /* Conversion to pole mass*/
375 mc_pole = mySM.Mbar2Mp(Mc, QCD::CHARM); /* Conversion to pole mass*/
376 Ms = mySM.getQuarks(QCD::STRANGE).getMass();
377 MW = mySM.Mw();
378 lambda_t = mySM.getCKM().computelamt_s();
379 lambda_u = mySM.getCKM().computelamu_s();
380 width = mySM.getMesons(meson).computeWidth();
381 alpha_s_mub = mySM.Als(mu_b);
382 fB = mySM.getMesons(meson).getDecayconst();
383 fpara = mySM.getMesons(vectorM).getDecayconst();
384 fperp = mySM.getMesons(vectorM).getDecayconst_p();
385
386 t_p = pow(MM + MV, 2.);
387 t_m = pow(MM - MV, 2.);
388 rV = MV/MM;
389 MM2 = MM*MM;
390 MM3 = MM2*MM;
391
392 switch (vectorM) {
393 case StandardModel::K_star:
394 case StandardModel::K_star_P:
395 if (MVll_DM_flag) {
396 Chi1minus = mySM.getOptionalParameter("Chi1minus"); //0.000623174575;
397 Chi1plus = mySM.getOptionalParameter("Chi1plus"); //0.000543940610;
398 Chi0plus = mySM.getOptionalParameter("Chi0plus"); //0.0142;
399 Chi0minus = mySM.getOptionalParameter("Chi0minus"); //0.0138586514;
400 ChiTT = mySM.getOptionalParameter("ChiTT"); //0.0454644444;
401 ChiBB = mySM.getOptionalParameter("ChiBB"); //0.0423069792;
402
403 a_0f = mySM.getOptionalParameter("a_0f");
404 a_1f = mySM.getOptionalParameter("a_1f");
405 a_2f = mySM.getOptionalParameter("a_2f");
406 MRf_2 = mySM.getOptionalParameter("MRf") * mySM.getOptionalParameter("MRf");
407
408 a_0g = mySM.getOptionalParameter("a_0g");
409 a_1g = mySM.getOptionalParameter("a_1g");
410 a_2g = mySM.getOptionalParameter("a_2g");
411 MRg_2 = mySM.getOptionalParameter("MRg") * mySM.getOptionalParameter("MRg");
412
413 a_1F1 = mySM.getOptionalParameter("a_1F1");
414 a_2F1 = mySM.getOptionalParameter("a_2F1");
415 MRF1_2 = mySM.getOptionalParameter("MRF1") * mySM.getOptionalParameter("MRF1");
416 a_0F1 = f_DM(t_m,a_0f,a_1f,a_2f,MRf_2)*MM*(1. - rV)*phi_F1(t_m, MRF1_2) - a_1F1*z_DM(t_m) - a_2F1*z_DM(t_m)*z_DM(t_m);
417
418 a_1F2 = mySM.getOptionalParameter("a_1F2");
419 a_2F2 = mySM.getOptionalParameter("a_2F2");
420 MRF2_2 = mySM.getOptionalParameter("MRF2") * mySM.getOptionalParameter("MRF2");
421 a_0F2 = F1_DM(0.,a_0F1,a_1F1,a_2F1,MRF1_2)*2./MM2/(1. - rV*rV)*phi_F2(0., MRF2_2) - a_1F2*z_DM(0.) - a_2F2*z_DM(0.)*z_DM(0.);
422
423 a_0T1 = mySM.getOptionalParameter("a_0T1");
424 a_1T1 = mySM.getOptionalParameter("a_1T1");
425 a_2T1 = mySM.getOptionalParameter("a_2T1");
426 MRT1_2 = mySM.getOptionalParameter("MRT1") * mySM.getOptionalParameter("MRT1");
427
428 a_1T2 = mySM.getOptionalParameter("a_1T2");
429 a_2T2 = mySM.getOptionalParameter("a_2T2");
430 MRT2_2 = mySM.getOptionalParameter("MRT2") * mySM.getOptionalParameter("MRT2");
431 a_0T2 = T1_DM(0.,a_0T1,a_1T1,a_2T1,MRT1_2)*phi_T2(0., MRT2_2) - a_1T2*z_DM(0.) - a_2T2*z_DM(0.)*z_DM(0.);
432
433 a_1T0 = mySM.getOptionalParameter("a_1T0");
434 a_2T0 = mySM.getOptionalParameter("a_2T0");
435 MRT0_2 = mySM.getOptionalParameter("MRT0") * mySM.getOptionalParameter("MRT0");
436 a_0T0 = T2_DM(t_m,a_0T2,a_1T2,a_2T2,MRT2_2)*phi_T0(t_m, MRT0_2) - a_1T0*z_DM(t_m) - a_2T0*z_DM(t_m)*z_DM(t_m);
437
438 unitarity_bound_f_F1 = pow(a_0f,2) + pow(a_1f,2) + pow(a_2f,2) + pow(a_0F1,2) + pow(a_1F1,2) + pow(a_2F1,2);
439 unitarity_bound_g = pow(a_0g,2) + pow(a_1g,2) + pow(a_2g,2);
440 unitarity_bound_F2 = pow(a_0F2,2) + pow(a_1F2,2) + pow(a_2F2,2);
441 unitarity_bound_T1 = pow(a_0T1,2) + pow(a_1T1,2) + pow(a_2T1,2);
442 unitarity_bound_T2_T0 = pow(a_0T2,2) + pow(a_1T2,2) + pow(a_2T2,2) + pow(a_0T0,2) + pow(a_1T0,2) + pow(a_2T0,2);
443 } else {
444 a_0V = mySM.getOptionalParameter("a_0V");
445 a_1V = mySM.getOptionalParameter("a_1V");
446 a_2V = mySM.getOptionalParameter("a_2V");
447 MRV_2 = mySM.getOptionalParameter("MRV") * mySM.getOptionalParameter("MRV");
448
449 a_0A0 = mySM.getOptionalParameter("a_0A0");
450 a_1A0 = mySM.getOptionalParameter("a_1A0");
451 a_2A0 = mySM.getOptionalParameter("a_2A0");
452 MRA0_2 = mySM.getOptionalParameter("MRA0") * mySM.getOptionalParameter("MRA0");
453
454 a_0A1 = mySM.getOptionalParameter("a_0A1");
455 a_1A1 = mySM.getOptionalParameter("a_1A1");
456 a_2A1 = mySM.getOptionalParameter("a_2A1");
457 MRA1_2 = mySM.getOptionalParameter("MRA1") * mySM.getOptionalParameter("MRA1");
458
459 a_0A12 = a_0A0 * (MM * MM - MV * MV) / (8. * MM * MV);
460 a_1A12 = mySM.getOptionalParameter("a_1A12");
461 a_2A12 = mySM.getOptionalParameter("a_2A12");
462 MRA12_2 = mySM.getOptionalParameter("MRA12") * mySM.getOptionalParameter("MRA12");
463
464 a_0T1 = mySM.getOptionalParameter("a_0T1");
465 a_1T1 = mySM.getOptionalParameter("a_1T1");
466 a_2T1 = mySM.getOptionalParameter("a_2T1");
467 MRT1_2 = mySM.getOptionalParameter("MRT1") * mySM.getOptionalParameter("MRT1");
468
469 a_0T2 = a_0T1;
470 a_1T2 = mySM.getOptionalParameter("a_1T2");
471 a_2T2 = mySM.getOptionalParameter("a_2T2");
472 MRT2_2 = mySM.getOptionalParameter("MRT2") * mySM.getOptionalParameter("MRT2");
473
474 a_0T23 = mySM.getOptionalParameter("a_0T23");
475 a_1T23 = mySM.getOptionalParameter("a_1T23");
476 a_2T23 = mySM.getOptionalParameter("a_2T23");
477 MRT23_2 = mySM.getOptionalParameter("MRT23") * mySM.getOptionalParameter("MRT23");
478 }
479
480 if (vectorM == StandardModel::K_star) spectator_charge = mySM.getQuarks(QCD::DOWN).getCharge();
481 else spectator_charge = mySM.getQuarks(QCD::UP).getCharge();
482
483 etaV = -1;
484 angmomV = 1.;
485
486 b = 1.;
487
488 SU3_breaking = 1.;
489
490 break;
491 case StandardModel::PHI:
492 if (MVll_DM_flag) {
493 Chi1minus = mySM.getOptionalParameter("Chi1minus"); //0.000623174575;
494 Chi1plus = mySM.getOptionalParameter("Chi1plus"); //0.000543940610;
495 Chi0plus = mySM.getOptionalParameter("Chi0plus"); //0.0142;
496 Chi0minus = mySM.getOptionalParameter("Chi0minus"); //0.0138586514;
497 ChiTT = mySM.getOptionalParameter("ChiTT"); //0.0454644444;
498 ChiBB = mySM.getOptionalParameter("ChiBB"); //0.0423069792;
499
500 a_0f = mySM.getOptionalParameter("a_0fphi");
501 a_1f = mySM.getOptionalParameter("a_1fphi");
502 a_2f = mySM.getOptionalParameter("a_2fphi");
503 MRf_2 = mySM.getOptionalParameter("MRf") * mySM.getOptionalParameter("MRf");
504
505 a_0g = mySM.getOptionalParameter("a_0gphi");
506 a_1g = mySM.getOptionalParameter("a_1gphi");
507 a_2g = mySM.getOptionalParameter("a_2gphi");
508 MRg_2 = mySM.getOptionalParameter("MRg") * mySM.getOptionalParameter("MRg");
509
510 a_1F1 = mySM.getOptionalParameter("a_1F1phi");
511 a_2F1 = mySM.getOptionalParameter("a_2F1phi");
512 MRF1_2 = mySM.getOptionalParameter("MRF1") * mySM.getOptionalParameter("MRF1");
513 a_0F1 = f_DM(t_m,a_0f,a_1f,a_2f,MRf_2)*MM*(1. - rV)*phi_F1(t_m, MRF1_2) - a_1F1*z_DM(t_m) - a_2F1*z_DM(t_m)*z_DM(t_m);
514
515 a_1F2 = mySM.getOptionalParameter("a_1F2phi");
516 a_2F2 = mySM.getOptionalParameter("a_2F2phi");
517 MRF2_2 = mySM.getOptionalParameter("MRF2") * mySM.getOptionalParameter("MRF2");
518 a_0F2 = F1_DM(0.,a_0F1,a_1F1,a_2F1,MRF1_2)*2./MM2/(1. - rV*rV)*phi_F2(0., MRF2_2) - a_1F2*z_DM(0.) - a_2F2*z_DM(0.)*z_DM(0.);
519
520 a_0T1 = mySM.getOptionalParameter("a_0T1phi");
521 a_1T1 = mySM.getOptionalParameter("a_1T1phi");
522 a_2T1 = mySM.getOptionalParameter("a_2T1phi");
523 MRT1_2 = mySM.getOptionalParameter("MRT1") * mySM.getOptionalParameter("MRT1");
524
525 a_1T2 = mySM.getOptionalParameter("a_1T2phi");
526 a_2T2 = mySM.getOptionalParameter("a_2T2phi");
527 MRT2_2 = mySM.getOptionalParameter("MRT2") * mySM.getOptionalParameter("MRT2");
528 a_0T2 = T1_DM(0.,a_0T1,a_1T1,a_2T1,MRT1_2)*phi_T2(0., MRT2_2) - a_1T2*z_DM(0.) - a_2T2*z_DM(0.)*z_DM(0.);
529
530 a_1T0 = mySM.getOptionalParameter("a_1T0phi");
531 a_2T0 = mySM.getOptionalParameter("a_2T0phi");
532 MRT0_2 = mySM.getOptionalParameter("MRT0") * mySM.getOptionalParameter("MRT0");
533 a_0T0 = T2_DM(t_m,a_0T2,a_1T2,a_2T2,MRT2_2)*phi_T0(t_m, MRT0_2) - a_1T0*z_DM(t_m) - a_2T0*z_DM(t_m)*z_DM(t_m);
534
535 unitarity_bound_f_F1 = pow(a_0f,2) + pow(a_1f,2) + pow(a_2f,2) + pow(a_0F1,2) + pow(a_1F1,2) + pow(a_2F1,2);
536 unitarity_bound_g = pow(a_0g,2) + pow(a_1g,2) + pow(a_2g,2);
537 unitarity_bound_F2 = pow(a_0F2,2) + pow(a_1F2,2) + pow(a_2F2,2);
538 unitarity_bound_T1 = pow(a_0T1,2) + pow(a_1T1,2) + pow(a_2T1,2);
539 unitarity_bound_T2_T0 = pow(a_0T2,2) + pow(a_1T2,2) + pow(a_2T2,2) + pow(a_0T0,2) + pow(a_1T0,2) + pow(a_2T0,2);
540 } else {
541 a_0V = mySM.getOptionalParameter("a_0Vphi");
542 a_1V = mySM.getOptionalParameter("a_1Vphi");
543 a_2V = mySM.getOptionalParameter("a_2Vphi");
544 MRV_2 = mySM.getOptionalParameter("MRV") * mySM.getOptionalParameter("MRV");
545
546 a_0A0 = mySM.getOptionalParameter("a_0A0phi");
547 a_1A0 = mySM.getOptionalParameter("a_1A0phi");
548 a_2A0 = mySM.getOptionalParameter("a_2A0phi");
549 MRA0_2 = mySM.getOptionalParameter("MRA0") * mySM.getOptionalParameter("MRA0");
550
551 a_0A1 = mySM.getOptionalParameter("a_0A1phi");
552 a_1A1 = mySM.getOptionalParameter("a_1A1phi");
553 a_2A1 = mySM.getOptionalParameter("a_2A1phi");
554 MRA1_2 = mySM.getOptionalParameter("MRA1") * mySM.getOptionalParameter("MRA1");
555
556 a_0A12 = a_0A0 * (MM * MM - MV * MV) / (8. * MM * MV);
557 a_1A12 = mySM.getOptionalParameter("a_1A12phi");
558 a_2A12 = mySM.getOptionalParameter("a_2A12phi");
559 MRA12_2 = mySM.getOptionalParameter("MRA12") * mySM.getOptionalParameter("MRA12");
560
561 a_0T1 = mySM.getOptionalParameter("a_0T1phi");
562 a_1T1 = mySM.getOptionalParameter("a_1T1phi");
563 a_2T1 = mySM.getOptionalParameter("a_2T1phi");
564 MRT1_2 = mySM.getOptionalParameter("MRT1") * mySM.getOptionalParameter("MRT1");
565
566 a_0T2 = a_0T1;
567 a_1T2 = mySM.getOptionalParameter("a_1T2phi");
568 a_2T2 = mySM.getOptionalParameter("a_2T2phi");
569 MRT2_2 = mySM.getOptionalParameter("MRT2") * mySM.getOptionalParameter("MRT2");
570
571 a_0T23 = mySM.getOptionalParameter("a_0T23phi");
572 a_1T23 = mySM.getOptionalParameter("a_1T23phi");
573 a_2T23 = mySM.getOptionalParameter("a_2T23phi");
574 MRT23_2 = mySM.getOptionalParameter("MRT23") * mySM.getOptionalParameter("MRT23");
575 }
576
577 spectator_charge = mySM.getQuarks(QCD::STRANGE).getCharge();
578
579 ys = mySM.getMesons(QCD::B_S).getDgamma_gamma() / 2.;
580 xs = mySM.getOptionalParameter("xs_phi");
581
582 etaV = -1;
583 angmomV = 1.;
584
585 b = 1.; //0.489;
586
587 SU3_breaking = 1. + gslpp::complex(mySM.getOptionalParameter("SU3_breaking_abs"),
588 mySM.getOptionalParameter("SU3_breaking_arg"), true);
589
590 break;
591 default:
592 std::stringstream out;
593 out << vectorM;
594 throw std::runtime_error("MVll: vector " + out.str() + " not implemented");
595 }
596
597 if (zExpansion) {
598 beta_0[0] = gslpp::complex(mySM.getOptionalParameter("re_beta_0_0"), mySM.getOptionalParameter("im_beta_0_0"), false);
599 beta_0[1] = gslpp::complex(mySM.getOptionalParameter("re_beta_0_1"), mySM.getOptionalParameter("im_beta_0_1"), false);
600 beta_0[2] = gslpp::complex(mySM.getOptionalParameter("re_beta_0_2"), mySM.getOptionalParameter("im_beta_0_2"), false);
601 beta_0[3] = gslpp::complex(mySM.getOptionalParameter("re_beta_0_3"), mySM.getOptionalParameter("im_beta_0_3"), false);
602 beta_0[4] = gslpp::complex(mySM.getOptionalParameter("re_beta_0_4"), mySM.getOptionalParameter("im_beta_0_4"), false);
603 beta_0[5] = gslpp::complex(mySM.getOptionalParameter("re_beta_0_5"), mySM.getOptionalParameter("im_beta_0_5"), false);
604 beta_0[6] = gslpp::complex(mySM.getOptionalParameter("re_beta_0_6"), mySM.getOptionalParameter("im_beta_0_6"), false);
605
606 beta_1[0] = gslpp::complex(mySM.getOptionalParameter("re_beta_1_0"), mySM.getOptionalParameter("im_beta_1_0"), false);
607 beta_1[1] = gslpp::complex(mySM.getOptionalParameter("re_beta_1_1"), mySM.getOptionalParameter("im_beta_1_1"), false);
608 beta_1[2] = gslpp::complex(mySM.getOptionalParameter("re_beta_1_2"), mySM.getOptionalParameter("im_beta_1_2"), false);
609 beta_1[3] = gslpp::complex(mySM.getOptionalParameter("re_beta_1_3"), mySM.getOptionalParameter("im_beta_1_3"), false);
610 beta_1[4] = gslpp::complex(mySM.getOptionalParameter("re_beta_1_4"), mySM.getOptionalParameter("im_beta_1_4"), false);
611 beta_1[5] = gslpp::complex(mySM.getOptionalParameter("re_beta_1_5"), mySM.getOptionalParameter("im_beta_1_5"), false);
612 beta_1[6] = gslpp::complex(mySM.getOptionalParameter("re_beta_1_6"), mySM.getOptionalParameter("im_beta_1_6"), false);
613
614 beta_2[0] = gslpp::complex(mySM.getOptionalParameter("re_beta_2_0"), mySM.getOptionalParameter("im_beta_2_0"), false);
615 beta_2[1] = gslpp::complex(mySM.getOptionalParameter("re_beta_2_1"), mySM.getOptionalParameter("im_beta_2_1"), false);
616 beta_2[2] = gslpp::complex(mySM.getOptionalParameter("re_beta_2_2"), mySM.getOptionalParameter("im_beta_2_2"), false);
617 beta_2[3] = gslpp::complex(mySM.getOptionalParameter("re_beta_2_3"), mySM.getOptionalParameter("im_beta_2_3"), false);
618 beta_2[4] = gslpp::complex(mySM.getOptionalParameter("re_beta_2_4"), mySM.getOptionalParameter("im_beta_2_4"), false);
619 beta_2[5] = gslpp::complex(mySM.getOptionalParameter("re_beta_2_5"), mySM.getOptionalParameter("im_beta_2_5"), false);
620 beta_2[6] = gslpp::complex(mySM.getOptionalParameter("re_beta_2_6"), mySM.getOptionalParameter("im_beta_2_6"), false);
621
622 DeltaC9 = mySM.getOptionalParameter("DeltaC9");
623 DeltaC10 = mySM.getOptionalParameter("DeltaC10");
624 } else if (dispersion) {
625 h_0[0] = gslpp::complex(mySM.getOptionalParameter("r1_1"));
626 h_0[1] = gslpp::complex(mySM.getOptionalParameter("r1_2"));
627 h_0[2] = gslpp::complex(mySM.getOptionalParameter("r1_3"));
628
629 h_1[0] = gslpp::complex(mySM.getOptionalParameter("r2_1"));
630 h_1[1] = gslpp::complex(mySM.getOptionalParameter("r2_2"));
631 h_1[2] = gslpp::complex(mySM.getOptionalParameter("r2_3"));
632
633 h_2[0] = gslpp::complex(mySM.getOptionalParameter("deltaC9_1"));
634 h_2[1] = gslpp::complex(mySM.getOptionalParameter("deltaC9_2"));
635 h_2[2] = gslpp::complex(mySM.getOptionalParameter("deltaC9_3"));
636 exp_Phase[0] = exp(gslpp::complex::i() * mySM.getOptionalParameter("phDC9_1"));
637 exp_Phase[1] = exp(gslpp::complex::i() * mySM.getOptionalParameter("phDC9_2"));
638 exp_Phase[2] = exp(gslpp::complex::i() * mySM.getOptionalParameter("phDC9_3"));
639 } else {
640#if NFPOLARBASIS_MVLL
641 h_0[0] = gslpp::complex(mySM.getOptionalParameter("absh_0"), mySM.getOptionalParameter("argh_0"), true);
642 h_0[1] = gslpp::complex(mySM.getOptionalParameter("absh_p"), mySM.getOptionalParameter("argh_p"), true);
643 h_0[2] = gslpp::complex(mySM.getOptionalParameter("absh_m"), mySM.getOptionalParameter("argh_m"), true);
644
645 h_1[0] = gslpp::complex(mySM.getOptionalParameter("absh_0_1"), mySM.getOptionalParameter("argh_0_1"), true);
646 h_1[1] = gslpp::complex(mySM.getOptionalParameter("absh_p_1"), mySM.getOptionalParameter("argh_p_1"), true);
647 h_1[2] = gslpp::complex(mySM.getOptionalParameter("absh_m_1"), mySM.getOptionalParameter("argh_m_1"), true);
648
649 h_2[0] = 0.;
650 h_2[1] = gslpp::complex(mySM.getOptionalParameter("absh_p_2"), mySM.getOptionalParameter("argh_p_2"), true);
651 h_2[2] = gslpp::complex(mySM.getOptionalParameter("absh_m_2"), mySM.getOptionalParameter("argh_m_2"), true);
652#else
653 h_0[0] = gslpp::complex(mySM.getOptionalParameter("reh_0"), mySM.getOptionalParameter("imh_0"), false);
654 h_0[1] = gslpp::complex(mySM.getOptionalParameter("reh_p"), mySM.getOptionalParameter("imh_p"), false);
655 h_0[2] = gslpp::complex(mySM.getOptionalParameter("reh_m"), mySM.getOptionalParameter("imh_m"), false);
656
657 h_1[0] = gslpp::complex(mySM.getOptionalParameter("reh_0_1"), mySM.getOptionalParameter("imh_0_1"), false);
658 h_1[1] = gslpp::complex(mySM.getOptionalParameter("reh_p_1"), mySM.getOptionalParameter("imh_p_1"), false);
659 h_1[2] = gslpp::complex(mySM.getOptionalParameter("reh_m_1"), mySM.getOptionalParameter("imh_m_1"), false);
660
661 h_2[0] = 0.;
662 h_2[1] = gslpp::complex(mySM.getOptionalParameter("reh_p_2"), mySM.getOptionalParameter("imh_p_2"), false);
663 h_2[2] = gslpp::complex(mySM.getOptionalParameter("reh_m_2"), mySM.getOptionalParameter("imh_m_2"), false);
664#endif
665 }
666 sqrt3 = sqrt(3.);
667
668 if (lep == QCD::NEUTRINO_1){
669 VusVub_abs2 = (mySM.getCKM().computelamu_s() * mySM.getCKM().computelamu_s().conjugate()).abs();
670 GF4 = GF * GF * GF * GF;
671 fM2 = mySM.getMesons(meson).getDecayconst() * mySM.getMesons(meson).getDecayconst();
672 fV2 = mySM.getMesons(vectorM).getDecayconst() * mySM.getMesons(vectorM).getDecayconst();
673 mtau = mySM.getLeptons(QCD::TAU).getMass();
674 mtau2 = mtau * mtau;
675 //from PDG 2024 tau lifetime: need SM prediction
676 Gammatau = HCUT / 0.2903;
677
678 allcoeff_nu = mySM.getFlavour().ComputeCoeffsnunu(QCD::NEUTRINO_1);
679 C_R_nunu_e = ((*(allcoeff_nu[LO]))(1) + (*(allcoeff_nu[NLO]))(1) + (*(allcoeff_nu[NLO_QED11]))(1));
680 if (FixedWCbtos) {
681 allcoeff_noSM_nu = mySM.getFlavour().ComputeCoeffsnunu(QCD::NEUTRINO_1,true); //check the mass scale, scheme fixed to NDR
682 C_L_nunu_e = mySM.getOptionalParameter("CLnunu_SM") + ((*(allcoeff_noSM_nu[LO]))(0) + (*(allcoeff_noSM_nu[NLO]))(0) + (*(allcoeff_noSM_nu[NLO_QED11]))(0));
683 } else
684 C_L_nunu_e = ((*(allcoeff_nu[LO]))(0) + (*(allcoeff_nu[NLO]))(0) + (*(allcoeff_nu[NLO_QED11]))(0));
685
686 allcoeff_nu = mySM.getFlavour().ComputeCoeffsnunu(QCD::NEUTRINO_2);
687 C_R_nunu_mu = ((*(allcoeff_nu[LO]))(1) + (*(allcoeff_nu[NLO]))(1) + (*(allcoeff_nu[NLO_QED11]))(1));
688 if (FixedWCbtos) {
689 allcoeff_noSM_nu = mySM.getFlavour().ComputeCoeffsnunu(QCD::NEUTRINO_2,true); //check the mass scale, scheme fixed to NDR
690 C_L_nunu_mu = mySM.getOptionalParameter("CLnunu_SM") + ((*(allcoeff_noSM_nu[LO]))(0) + (*(allcoeff_noSM_nu[NLO]))(0) + (*(allcoeff_noSM_nu[NLO_QED11]))(0));
691 } else
692 C_L_nunu_mu = ((*(allcoeff_nu[LO]))(0) + (*(allcoeff_nu[NLO]))(0) + (*(allcoeff_nu[NLO_QED11]))(0));
693
694 allcoeff_nu = mySM.getFlavour().ComputeCoeffsnunu(QCD::NEUTRINO_3);
695 C_R_nunu_tau = ((*(allcoeff_nu[LO]))(1) + (*(allcoeff_nu[NLO]))(1) + (*(allcoeff_nu[NLO_QED11]))(1));
696 if (FixedWCbtos) {
697 allcoeff_noSM_nu = mySM.getFlavour().ComputeCoeffsnunu(QCD::NEUTRINO_3,true); //check the mass scale, scheme fixed to NDR
698 C_L_nunu_tau = mySM.getOptionalParameter("CLnunu_SM") + ((*(allcoeff_noSM_nu[LO]))(0) + (*(allcoeff_noSM_nu[NLO]))(0) + (*(allcoeff_noSM_nu[NLO_QED11]))(0));
699 } else
700 C_L_nunu_tau = ((*(allcoeff_nu[LO]))(0) + (*(allcoeff_nu[NLO]))(0) + (*(allcoeff_nu[NLO_QED11]))(0));
701
702 C_L_nunu = sqrt(C_L_nunu_e * C_L_nunu_e + C_L_nunu_mu * C_L_nunu_mu + C_L_nunu_tau * C_L_nunu_tau);
703 C_R_nunu = sqrt(C_R_nunu_e * C_R_nunu_e + C_R_nunu_mu * C_R_nunu_mu + C_R_nunu_tau * C_R_nunu_tau);
704 }
705 else{
706 allcoeff = mySM.getFlavour().ComputeCoeffBMll(mu_b, lep); //check the mass scale, scheme fixed to NDR
707 allcoeffprime = mySM.getFlavour().ComputeCoeffprimeBMll(mu_b, lep); //check the mass scale, scheme fixed to NDR
708
709 C_1 = ((*(allcoeff[LO]))(0) + (*(allcoeff[NLO]))(0));
710 C_1L_bar = (*(allcoeff[LO]))(0) / 2.;
711 C_2 = ((*(allcoeff[LO]))(1) + (*(allcoeff[NLO]))(1));
712 C_2L_bar = (*(allcoeff[LO]))(1) - (*(allcoeff[LO]))(0) / 6.;
713 C_3 = ((*(allcoeff[LO]))(2) + (*(allcoeff[NLO]))(2));
714 C_4 = ((*(allcoeff[LO]))(3) + (*(allcoeff[NLO]))(3));
715 C_5 = ((*(allcoeff[LO]))(4) + (*(allcoeff[NLO]))(4));
716 C_6 = ((*(allcoeff[LO]))(5) + (*(allcoeff[NLO]))(5));
717 C_8 = ((*(allcoeff[LO]))(7) + (*(allcoeff[NLO]))(7));
718 C_8L = (*(allcoeff[LO]))(7);
719 C_S = MW / Mb * (((*(allcoeff[LO]))(10) + (*(allcoeff[NLO]))(10)));
720 C_P = MW / Mb * (((*(allcoeff[LO]))(11) + (*(allcoeff[NLO]))(11)));
721 C_9p = (*(allcoeffprime[LO]))(8) + (*(allcoeffprime[NLO]))(8);
722 C_10p = (*(allcoeffprime[LO]))(9) + (*(allcoeffprime[NLO]))(9);
723 C_Sp = MW / Mb * ((*(allcoeffprime[LO]))(10) + (*(allcoeffprime[NLO]))(10));
724 C_Pp = MW / Mb * ((*(allcoeffprime[LO]))(11) + (*(allcoeffprime[NLO]))(11));
725
726 if (FixedWCbtos) {
727 allcoeff_noSM = mySM.getFlavour().ComputeCoeffBMll(mu_b, lep, true); //check the mass scale, scheme fixed to NDR
728 C_7 = mySM.getOptionalParameter("C7_SM") + ((*(allcoeff_noSM[LO]))(6) + (*(allcoeff_noSM[NLO]))(6));
729 C_9 = mySM.getOptionalParameter("C9_SM") + ((*(allcoeff_noSM[LO]))(8) + (*(allcoeff_noSM[NLO]))(8));
730 C_10 = mySM.getOptionalParameter("C10_SM") + ((*(allcoeff_noSM[LO]))(9) + (*(allcoeff_noSM[NLO]))(9));
731 } else {
732 C_7 = ((*(allcoeff[LO]))(6) + (*(allcoeff[NLO]))(6));
733 C_9 = ((*(allcoeff[LO]))(8) + (*(allcoeff[NLO]))(8));
734 C_10 = ((*(allcoeff[LO]))(9) + (*(allcoeff[NLO]))(9));
735 }
736 C_7p = MsoMb * ((*(allcoeffprime[LO]))(6) + (*(allcoeffprime[NLO]))(6));
737 C_7p -= MsoMb * (C_7 + 1. / 3. * C_3 + 4 / 9 * C_4 + 20. / 3. * C_5 + 80. / 9. * C_6);
738
739 allcoeffh = mySM.getFlavour().ComputeCoeffBMll(mu_h, lep); //check the mass scale, scheme fixed to NDR
740
741 C_1Lh_bar = (*(allcoeffh[LO]))(0) / 2.;
742 C_2Lh_bar = (*(allcoeffh[LO]))(1) - (*(allcoeff[LO]))(0) / 6.;
743 C_8Lh = (*(allcoeffh[LO]))(7);
744
745 if (zExpansion) {
746 C_9 += DeltaC9;
747 C_10 += DeltaC10;
748 }
749 }
750
751 checkCache();
752
753 t_0 = t_p * (1. - sqrt(1. - t_m / t_p)); /*Modify it for Lattice*/
754 z_0 = (sqrt(t_p) - sqrt(t_p - t_0)) / (sqrt(t_p) + sqrt(t_p - t_0));
755 s_p = 4. * mD2;
756 // s_0 = 4.;
757 s_0 = s_p - sqrt(s_p * (s_p - mPsi2S2));
758 Q2 = - Mb*Mb;
759 chiOPE = 0.000181;
760 twoalphaBtoKst = 2.276;
761 rho_0 = 0.7977;
762 rho_1 = -0.8298;
763 rho_2 = 0.8372;
764 rho_3 = -0.8396;
765 rho_4 = 0.8406;
766 rho_5 = -0.8412;
767 onemrho_0_2 = 1. - rho_0*rho_0;
768 onemrho_1_2 = 1. - rho_1*rho_1;
769 onemrho_2_2 = 1. - rho_2*rho_2;
770 onemrho_3_2 = 1. - rho_3*rho_3;
771 onemrho_4_2 = 1. - rho_4*rho_4;
772 onemrho_5_2 = 1. - rho_5*rho_5;
773 MMpMV = MM + MV;
774 MMpMV2 = MMpMV * MMpMV;
775 MMmMV = MM - MV;
776 MMmMV2 = MMmMV * MMmMV;
777 MM4 = MM2*MM2;
778 MV2 = MV*MV;
779 MV4 = MV2*MV2;
780 MMMV = MM*MV;
781 MM2mMV2 = MM2 - MV2;
782 MM2pMV2 = MM2 + MV2;
783 fourMV = 4. * MV;
784 twoMM2 = 2. * MM2;
785 twoMV2 = 2. * MV2;
786 onepMMoMV = (1. + MV / MM);
787 MM_MMpMV = MM * MMpMV;
788 twoMM_mbpms = 2. * MM * (Mb + Ms);
789 fourMM2 = 4. * MM2;
790 Mlep2 = Mlep*Mlep;
791 twoMlepMb = 2. * Mlep*Mb;
792 MboMW = Mb / MW;
793 MsoMb = Ms / Mb;
794 ninetysixM_PI3MM3 = 96. * M_PI * M_PI * M_PI * MM * MM*MM;
795 sixteenM_PI2 = 16. * M_PI2;
796 sixteenM_PI2MM2 = sixteenM_PI2 * MM*MM;
797 twoMboMM = 2 * Mb / MM;
798 H_0_pre = 8. / 27. + 4. / 9. * gslpp::complex::i() * M_PI;
799 H_0_WC = (C_3 + 4. / 3. * C_4 + 16. * C_5 + 64. / 3. * C_6);
800 H_c_WC = (4. / 3. * C_1 + C_2 + 6. * C_3 + 60. * C_5);
801 H_b_WC = (7. * C_3 + 4. / 3. * C_4 + 76. * C_5 + 64. / 3. * C_6);
802 mu_b2 = mu_b*mu_b;
803 Mc2 = Mc*Mc;
804 Mb2 = Mb*Mb;
805 fourMc2 = 4. * Mc2;
806 fourMb2 = 4. * Mb2;
807 logMc = log(Mc2 / mu_b2);
808 logMb = log(Mb2 / mu_b2);
809 fournineth = 4. / 9.;
810 half = 1. / 2.;
811 twothird = 2. / 3.;
812 ihalfMPI = gslpp::complex::i() * M_PI / 2.;
813 twoMM3 = 2. * MM2 * MM;
814 C2_inv = 1. / (2. * C_2.real());
815 gtilde_1_pre = -16. * pow(MM, 3.)*(MM + MV) * pow(M_PI, 2.);
816 gtilde_2_pre = -16. * pow(MM, 3.) * pow(M_PI, 2.) / MMpMV;
817 gtilde_3_pre = 64. * pow(MM, 3.) * pow(M_PI, 2.) * MV*MMpMV;
818 S_L_pre = (-2. * MM * (Mb + Ms));
819
820 M_PI2osix = M_PI2 / 6.;
821 twoMM = 2. * MM;
822
823 N_QCDF = M_PI2 / 3. * fB * fperp / MM;
824
825 deltaT_0 = alpha_s_mub * CF / 4. / M_PI;
826 deltaT_1par = mySM.Als(mu_h) * CF / 4. * M_PI / 3. * mySM.getMesons(meson).getDecayconst() *
827 mySM.getMesons(vectorM).getDecayconst() / MM;
828 deltaT_1perp = mySM.Als(mu_h) * CF / 4. * M_PI / 3. * mySM.getMesons(meson).getDecayconst() *
829 mySM.getMesons(vectorM).getDecayconst_p() / MM;
830
831 F87_0 = -32. / 9. * log(mu_b / Mb) + 8. / 27. * M_PI2 - 44. / 9. - 8. / 9. * gslpp::complex::i() * M_PI;
832
833 NN = -(4. * GF * MM * ale * lambda_t) / (sqrt(2.)*4. * M_PI);
834 NN_conjugate = -(4. * GF * MM * ale * lambda_t.conjugate()) / (sqrt(2.)*4. * M_PI);
835
836 std::map<std::pair<double, double>, unsigned int >::iterator it;
837
838 if (I0_updated == 0) for (it = sigma0Cached.begin(); it != sigma0Cached.end(); ++it) it->second = 0;
839 if (I1_updated == 0) for (it = sigma1Cached.begin(); it != sigma1Cached.end(); ++it) it->second = 0;
840 if (I2_updated == 0) for (it = sigma2Cached.begin(); it != sigma2Cached.end(); ++it) it->second = 0;
841 if (I3_updated == 0) for (it = sigma3Cached.begin(); it != sigma3Cached.end(); ++it) it->second = 0;
842 if (I4_updated == 0) for (it = sigma4Cached.begin(); it != sigma4Cached.end(); ++it) it->second = 0;
843 if (I5_updated == 0) for (it = sigma5Cached.begin(); it != sigma5Cached.end(); ++it) it->second = 0;
844 if (I6_updated == 0) for (it = sigma6Cached.begin(); it != sigma6Cached.end(); ++it) it->second = 0;
845 if (I7_updated == 0) for (it = sigma7Cached.begin(); it != sigma7Cached.end(); ++it) it->second = 0;
846 if (I8_updated == 0) for (it = sigma8Cached.begin(); it != sigma8Cached.end(); ++it) it->second = 0;
847 if (I9_updated == 0) for (it = sigma9Cached.begin(); it != sigma9Cached.end(); ++it) it->second = 0;
848 if (I10_updated == 0) for (it = sigma10Cached.begin(); it != sigma10Cached.end(); ++it) it->second = 0;
849 if (I11_updated == 0) for (it = sigma11Cached.begin(); it != sigma11Cached.end(); ++it) it->second = 0;
850
851 if (I0_updated == 0) for (it = delta0Cached.begin(); it != delta0Cached.end(); ++it) it->second = 0;
852 if (I1_updated == 0) for (it = delta1Cached.begin(); it != delta1Cached.end(); ++it) it->second = 0;
853 if (I2_updated == 0) for (it = delta2Cached.begin(); it != delta2Cached.end(); ++it) it->second = 0;
854 if (I3_updated == 0) for (it = delta3Cached.begin(); it != delta3Cached.end(); ++it) it->second = 0;
855 if (I11_updated == 0) for (it = delta11Cached.begin(); it != delta11Cached.end(); ++it) it->second = 0;
856
857 if (Itree_updated) for (it = sigmaTreeCached.begin(); it != sigmaTreeCached.end(); ++it) it->second = 0;
858
859 std::map<double, unsigned int >::iterator iti;
860 if (deltaTparpupdated == 0) for (iti = deltaTparpCached.begin(); iti != deltaTparpCached.end(); ++iti) iti->second = 0;
861 if (deltaTparmupdated == 0) for (iti = deltaTparmCached.begin(); iti != deltaTparmCached.end(); ++iti) iti->second = 0;
862 if (deltaTperpupdated == 0) for (iti = deltaTparpCached.begin(); iti != deltaTparpCached.end(); ++iti) iti->second = 0;
863
864 if (deltaTparpupdated * deltaTparmupdated == 0) for (it = I1Cached.begin(); it != I1Cached.end(); ++it) it->second = 0;
865
866#if SPLINE
867 if (mySM.getFlavour().getUpdateFlag(meson, vectorM, lep)) spline_QCDF_func();
868#else
869 if (mySM.getFlavour().getUpdateFlag(meson, vectorM, lep)) fit_QCDF_func();
870#endif
871
872 mySM.getFlavour().setUpdateFlag(meson, vectorM, lep, false);
873
874 /*
875 std::cout << "MVll: meson type: " << vectorM << std::endl;
876 std::cout << "MM: " << MM << std::endl;
877 std::cout << "MV: " << MV << std::endl;
878
879 std::cout << "a_0F1: " << a_0F1 << std::endl;
880 std::cout << "a_0F2: " << a_0F2 << std::endl;
881 std::cout << "a_0T0: " << a_0T0 << std::endl;
882 std::cout << "a_0T2: " << a_0T2 << std::endl;
883
884 std::cout << "f_DM(4.): " << f_DM(4., a_0f, a_1f, a_2f, MRf_2) << std::endl;
885 std::cout << "g_DM(4.): " << g_DM(4., a_0g, a_1g, a_2g, MRg_2) << std::endl;
886 std::cout << "F1_DM(4.): " << F1_DM(4., a_0F1, a_1F1, a_2F1, MRF1_2) << std::endl;
887 std::cout << "F2_DM(4.): " << F2_DM(4., a_0F2, a_1F2, a_2F2, MRF2_2) << std::endl;
888 std::cout << "T0_DM(4.): " << T0_DM(4., a_0T0, a_1T0, a_2T0, MRT0_2) << std::endl;
889 std::cout << "T1_DM(4.): " << T1_DM(4., a_0T1, a_1T1, a_2T1, MRT1_2) << std::endl;
890 std::cout << "T2_DM(4.): " << T2_DM(4., a_0T2, a_1T2, a_2T2, MRT2_2) << std::endl << std::endl;
891
892 std::cout << "V(1.): " << V(1.) << std::endl;
893 std::cout << "A_0(1.): " << A_0(1.) << std::endl;
894 std::cout << "A_1(1.): " << A_1(1.) << std::endl;
895 std::cout << "A_2(1.): " << A_2(1.) << std::endl;
896 std::cout << "T_1(1.): " << T_1(1.) << std::endl;
897 std::cout << "T_2(1.): " << T_2(1.) << std::endl;
898 std::cout << "V_p(1.): " << V_p(1.) << std::endl;
899 std::cout << "V_m(1.): " << V_m(1.) << std::endl;
900 std::cout << "V_0t(1.): " << V_0t(1.) << std::endl;
901 std::cout << "T_p(1.): " << T_p(1.) << std::endl;
902 std::cout << "T_m(1.): " << T_m(1.) << std::endl;
903 std::cout << "T_0t(1.): " << T_0t(1.) << std::endl << std::endl;
904
905 std::cout << "V(4.): " << V(4.) << std::endl;
906 std::cout << "A_0(4.): " << A_0(4.) << std::endl;
907 std::cout << "A_1(4.): " << A_1(4.) << std::endl;
908 std::cout << "A_2(4.): " << A_2(4.) << std::endl;
909 std::cout << "T_1(4.): " << T_1(4.) << std::endl;
910 std::cout << "T_2(4.): " << T_2(4.) << std::endl;
911 std::cout << "V_p(4.): " << V_p(4.) << std::endl;
912 std::cout << "V_m(4.): " << V_m(4.) << std::endl;
913 std::cout << "V_0t(4.): " << V_0t(4.) << std::endl;
914 std::cout << "T_p(4.): " << T_p(4.) << std::endl;
915 std::cout << "T_m(4.): " << T_m(4.) << std::endl;
916 std::cout << "T_0t(4.): " << T_0t(4.) << std::endl << std::endl;
917
918 std::cout << "V(8.): " << V(8.) << std::endl;
919 std::cout << "A_0(8.): " << A_0(8.) << std::endl;
920 std::cout << "A_1(8.): " << A_1(8.) << std::endl;
921 std::cout << "A_2(8.): " << A_2(8.) << std::endl;
922 std::cout << "T_1(8.): " << T_1(8.) << std::endl;
923 std::cout << "T_2(8.): " << T_2(8.) << std::endl;
924 std::cout << "V_p(8.): " << V_p(8.) << std::endl;
925 std::cout << "V_m(8.): " << V_m(8.) << std::endl;
926 std::cout << "V_0t(8.): " << V_0t(8.) << std::endl;
927 std::cout << "T_p(8.): " << T_p(8.) << std::endl;
928 std::cout << "T_m(8.): " << T_m(8.) << std::endl;
929 std::cout << "T_0t(8.): " << T_0t(8.) << std::endl << std::endl;
930 */
931
932 return;
933}
934
935void MVll::checkCache()
936{
937
938 if (MM == k2_cache(0) && MV == k2_cache(1)) {
939 k2_updated = 1;
940 z_updated = 1;
941 } else {
942 k2_updated = 0;
943 z_updated = 0;
944 k2_cache(0) = MM;
945 k2_cache(1) = MV;
946 }
947
948 if (Mlep == beta_cache) {
949 beta_updated = 1;
950 } else {
951 beta_updated = 0;
952 beta_cache = Mlep;
953 }
954
955 lambda_updated = k2_updated;
956 F_updated = lambda_updated * beta_updated;
957
958 if (GF == N_cache(0) && ale == N_cache(1) && MM == N_cache(2) && lambda_t == Nc_cache) {
959 N_updated = 1;
960 } else {
961 N_updated = 0;
962 N_cache(0) = GF;
963 N_cache(1) = ale;
964 N_cache(2) = MM;
965 Nc_cache = lambda_t;
966 }
967 if (MVll_DM_flag) {
968 if (a_0g == V_cache(0) && a_1g == V_cache(1) && a_2g == V_cache(2)) {
969 V_updated = V_updated * z_updated;
970 } else {
971 V_updated = 0;
972 V_cache(0) = a_0g;
973 V_cache(1) = a_1g;
974 V_cache(2) = a_2g;
975 }
976
977 if (a_0F2 == A0_cache(0) && a_1F2 == A0_cache(1) && a_2F2 == A0_cache(2)) {
978 A0_updated = A0_updated * z_updated;
979 } else {
980 A0_updated = 0;
981 A0_cache(0) = a_0F2;
982 A0_cache(1) = a_1F2;
983 A0_cache(2) = a_2F2;
984 }
985
986 if (a_0f == A1_cache(0) && a_1f == A1_cache(1) && a_2f == A1_cache(2)) {
987 A1_updated = A1_updated * z_updated;
988 } else {
989 A1_updated = 0;
990 A1_cache(0) = a_0f;
991 A1_cache(1) = a_1f;
992 A1_cache(2) = a_2f;
993 }
994
995 if (a_0T1 == T1_cache(0) && a_1T1 == T1_cache(1) && a_2T1 == T1_cache(2)) {
996 T1_updated = T1_updated * z_updated;
997 } else {
998 T1_updated = 0;
999 T1_cache(0) = a_0T1;
1000 T1_cache(1) = a_1T1;
1001 T1_cache(2) = a_2T1;
1002 }
1003
1004 if (a_0T2 == T2_cache(0) && a_1T2 == T2_cache(1) && a_2T2 == T2_cache(2)) {
1005 T2_updated = T2_updated * z_updated;
1006 } else {
1007 T2_updated = 0;
1008 T2_cache(0) = a_0T2;
1009 T2_cache(1) = a_1T2;
1010 T2_cache(2) = a_2T2;
1011 }
1012 } else {
1013 if (a_0V == V_cache(0) && a_1V == V_cache(1) && a_2V == V_cache(2)) {
1014 V_updated = V_updated * z_updated;
1015 } else {
1016 V_updated = 0;
1017 V_cache(0) = a_0V;
1018 V_cache(1) = a_1V;
1019 V_cache(2) = a_2V;
1020 }
1021
1022 if (a_0A0 == A0_cache(0) && a_1A0 == A0_cache(1) && a_2A0 == A0_cache(2)) {
1023 A0_updated = A0_updated * z_updated;
1024 } else {
1025 A0_updated = 0;
1026 A0_cache(0) = a_0A0;
1027 A0_cache(1) = a_1A0;
1028 A0_cache(2) = a_2A0;
1029 }
1030
1031 if (a_0A1 == A1_cache(0) && a_1A1 == A1_cache(1) && a_2A1 == A1_cache(2)) {
1032 A1_updated = A1_updated * z_updated;
1033 } else {
1034 A1_updated = 0;
1035 A1_cache(0) = a_0A1;
1036 A1_cache(1) = a_1A1;
1037 A1_cache(2) = a_2A1;
1038 }
1039
1040 if (a_0T1 == T1_cache(0) && a_1T1 == T1_cache(1) && a_2T1 == T1_cache(2)) {
1041 T1_updated = T1_updated * z_updated;
1042 } else {
1043 T1_updated = 0;
1044 T1_cache(0) = a_0T1;
1045 T1_cache(1) = a_1T1;
1046 T1_cache(2) = a_2T1;
1047 }
1048
1049 if (a_0T2 == T2_cache(0) && a_1T2 == T2_cache(1) && a_2T2 == T2_cache(2)) {
1050 T2_updated = T2_updated * z_updated;
1051 } else {
1052 T2_updated = 0;
1053 T2_cache(0) = a_0T2;
1054 T2_cache(1) = a_1T2;
1055 T2_cache(2) = a_2T2;
1056 }
1057 }
1058
1059 VL1_updated = k2_updated * lambda_updated * A1_updated * V_updated;
1060 VL2_updated = VL1_updated;
1061
1062 TL1_updated = k2_updated * lambda_updated * T1_updated * T2_updated;
1063 TL2_updated = TL1_updated;
1064
1065 VR1_updated = VL2_updated;
1066 VR2_updated = VL1_updated;
1067
1068 TR1_updated = TL2_updated;
1069 TR2_updated = TL1_updated;
1070
1071 if (Mb == SL_cache(0) && Ms == SL_cache(1)) {
1072 Mb_Ms_updated = 1;
1073 SL_updated = lambda_updated * A0_updated;
1074 SR_updated = SL_updated;
1075 } else {
1076 Mb_Ms_updated = 0;
1077 SL_updated = 0;
1078 SR_updated = SL_updated;
1079 SL_cache(0) = Mb;
1080 SL_cache(1) = Ms;
1081 }
1082
1083 if (MVll_DM_flag) {
1084 if (a_0F1 == VL0_cache(0) && a_1F1 == VL0_cache(1) && a_2F1 == VL0_cache(2)) {
1085 VL0_updated = VL0_updated * z_updated;
1086 VR0_updated = VL0_updated;
1087 } else {
1088 VL0_updated = 0;
1089 VR0_updated = VL0_updated;
1090 VL0_cache(0) = a_0F1;
1091 VL0_cache(1) = a_1F1;
1092 VL0_cache(2) = a_2F1;
1093 }
1094
1095 if (a_0T0 == TL0_cache(0) && a_1T0 == TL0_cache(1) && a_2T0 == TL0_cache(2)) {
1096 TL0_updated = TL0_updated * z_updated;
1097 TR0_updated = TL0_updated;
1098 } else {
1099 TL0_updated = 0;
1100 TR0_updated = TL0_updated;
1101 TL0_cache(0) = a_0T0;
1102 TL0_cache(1) = a_1T0;
1103 TL0_cache(2) = a_2T0;
1104 }
1105 } else {
1106 if (a_0A12 == VL0_cache(0) && a_1A12 == VL0_cache(1) && a_2A12 == VL0_cache(2)) {
1107 VL0_updated = VL0_updated * z_updated;
1108 VR0_updated = VL0_updated;
1109 } else {
1110 VL0_updated = 0;
1111 VR0_updated = VL0_updated;
1112 VL0_cache(0) = a_0A12;
1113 VL0_cache(1) = a_1A12;
1114 VL0_cache(2) = a_2A12;
1115 }
1116
1117 if (a_0T23 == TL0_cache(0) && a_1T23 == TL0_cache(1) && a_2T23 == TL0_cache(2)) {
1118 TL0_updated = TL0_updated * z_updated;
1119 TR0_updated = TL0_updated;
1120 } else {
1121 TL0_updated = 0;
1122 TR0_updated = TL0_updated;
1123 TL0_cache(0) = a_0T23;
1124 TL0_cache(1) = a_1T23;
1125 TL0_cache(2) = a_2T23;
1126 }
1127 }
1128
1129
1130 if (C_1 == C_1_cache) {
1131 C_1_updated = 1;
1132 } else {
1133 C_1_updated = 0;
1134 C_1_cache = C_1;
1135 }
1136
1137 if (C_2 == C_2_cache) {
1138 C_2_updated = 1;
1139 } else {
1140 C_2_updated = 0;
1141 C_2_cache = C_2;
1142 }
1143
1144 if (C_3 == C_3_cache) {
1145 C_3_updated = 1;
1146 } else {
1147 C_3_updated = 0;
1148 C_3_cache = C_3;
1149 }
1150
1151 if (C_4 == C_4_cache) {
1152 C_4_updated = 1;
1153 } else {
1154 C_4_updated = 0;
1155 C_4_cache = C_4;
1156 }
1157
1158 if (C_5 == C_5_cache) {
1159 C_5_updated = 1;
1160 } else {
1161 C_5_updated = 0;
1162 C_5_cache = C_5;
1163 }
1164
1165 if (C_6 == C_6_cache) {
1166 C_6_updated = 1;
1167 } else {
1168 C_6_updated = 0;
1169 C_6_cache = C_6;
1170 }
1171
1172 if (C_7 == C_7_cache) {
1173 C_7_updated = 1;
1174 } else {
1175 C_7_updated = 0;
1176 C_7_cache = C_7;
1177 }
1178
1179 if (C_9 == C_9_cache) {
1180 C_9_updated = 1;
1181 } else {
1182 C_9_updated = 0;
1183 C_9_cache = C_9;
1184 }
1185
1186 if (C_10 == C_10_cache) {
1187 C_10_updated = 1;
1188 } else {
1189 C_10_updated = 0;
1190 C_10_cache = C_10;
1191 }
1192
1193 if (C_S == C_S_cache) {
1194 C_S_updated = 1;
1195 } else {
1196 C_S_updated = 0;
1197 C_S_cache = C_S;
1198 }
1199
1200 if (C_P == C_P_cache) {
1201 C_P_updated = 1;
1202 } else {
1203 C_P_updated = 0;
1204 C_P_cache = C_P;
1205 }
1206
1207 if (C_7p == C_7p_cache) {
1208 C_7p_updated = 1;
1209 } else {
1210 C_7p_updated = 0;
1211 C_7p_cache = C_7p;
1212 }
1213
1214 if (C_9p == C_9p_cache) {
1215 C_9p_updated = 1;
1216 } else {
1217 C_9p_updated = 0;
1218 C_9p_cache = C_9p;
1219 }
1220
1221 if (C_10p == C_10p_cache) {
1222 C_10p_updated = 1;
1223 } else {
1224 C_10p_updated = 0;
1225 C_10p_cache = C_10p;
1226 }
1227
1228 if (C_Sp == C_Sp_cache) {
1229 C_Sp_updated = 1;
1230 } else {
1231 C_Sp_updated = 0;
1232 C_Sp_cache = C_Sp;
1233 }
1234
1235 if (C_Pp == C_Pp_cache) {
1236 C_Pp_updated = 1;
1237 } else {
1238 C_Pp_updated = 0;
1239 C_Pp_cache = C_Pp;
1240 }
1241
1242 if (C_2Lh_bar == C_2Lh_cache) {
1243 C_2Lh_updated = 1;
1244 } else {
1245 C_2Lh_updated = 0;
1246 C_2Lh_cache = C_2Lh_bar;
1247 }
1248
1249 if (C_8Lh == C_8Lh_cache) {
1250 C_8Lh_updated = 1;
1251 } else {
1252 C_8Lh_updated = 0;
1253 C_8Lh_cache = C_8Lh;
1254 }
1255
1256 if (C_L_nunu == C_L_nunu_cache) {
1257 C_L_nunu_updated = 1;
1258 } else {
1259 C_L_nunu_updated = 0;
1260 C_L_nunu_cache = C_L_nunu;
1261 }
1262
1263 if (C_R_nunu == C_R_nunu_cache) {
1264 C_R_nunu_updated = 1;
1265 } else {
1266 C_R_nunu_updated = 0;
1267 C_R_nunu_cache = C_R_nunu;
1268 }
1269
1270 if (Mb == Ycache(0) && Mc == Ycache(1)) {
1271 Yupdated = C_1_updated * C_2_updated * C_3_updated * C_4_updated * C_5_updated * C_6_updated;
1272 } else {
1273 Yupdated = 0;
1274 Ycache(0) = Mb;
1275 Ycache(1) = Mc;
1276 }
1277
1278 if (zExpansion) {
1279 if (beta_0[0] == beta0Ccache[0] && beta_0[1] == beta0Ccache[1] && beta_0[2] == beta0Ccache[2] && beta_0[3] == beta0Ccache[3]
1280 && beta_0[4] == beta0Ccache[4] && beta_0[5] == beta0Ccache[5] && beta_0[6] == beta0Ccache[6] && SU3_breaking == beta0Ccache[7]) {
1281 h0_updated = 1;
1282 } else {
1283 h0_updated = 0;
1284 beta0Ccache[0] = beta_0[0];
1285 beta0Ccache[1] = beta_0[1];
1286 beta0Ccache[2] = beta_0[2];
1287 beta0Ccache[3] = beta_0[3];
1288 beta0Ccache[4] = beta_0[4];
1289 beta0Ccache[5] = beta_0[5];
1290 beta0Ccache[6] = beta_0[6];
1291 beta0Ccache[7] = SU3_breaking;
1292 }
1293
1294 if (beta_1[0] == beta1Ccache[0] && beta_1[1] == beta1Ccache[1] && beta_1[2] == beta1Ccache[2] && beta_1[3] == beta1Ccache[3]
1295 && beta_1[4] == beta1Ccache[4] && beta_1[5] == beta1Ccache[5] && beta_1[6] == beta1Ccache[6] && SU3_breaking == beta1Ccache[7]) {
1296 h1_updated = 1;
1297 } else {
1298 h1_updated = 0;
1299 beta1Ccache[0] = beta_1[0];
1300 beta1Ccache[1] = beta_1[1];
1301 beta1Ccache[2] = beta_1[2];
1302 beta1Ccache[3] = beta_1[3];
1303 beta1Ccache[4] = beta_1[4];
1304 beta1Ccache[5] = beta_1[5];
1305 beta1Ccache[6] = beta_1[6];
1306 beta1Ccache[7] = SU3_breaking;
1307 }
1308
1309 if (beta_2[0] == beta2Ccache[0] && beta_2[1] == beta2Ccache[1] && beta_2[2] == beta2Ccache[2] && beta_2[3] == beta2Ccache[3]
1310 && beta_2[4] == beta2Ccache[4] && beta_2[5] == beta2Ccache[5] && beta_2[6] == beta2Ccache[6] && SU3_breaking == beta2Ccache[7]) {
1311 h2_updated = 1;
1312 } else {
1313 h2_updated = 0;
1314 beta2Ccache[0] = beta_2[0];
1315 beta2Ccache[1] = beta_2[1];
1316 beta2Ccache[2] = beta_2[2];
1317 beta2Ccache[3] = beta_2[3];
1318 beta2Ccache[4] = beta_2[4];
1319 beta2Ccache[5] = beta_2[5];
1320 beta2Ccache[6] = beta_2[6];
1321 beta2Ccache[7] = SU3_breaking;
1322 }
1323 } else {
1324 if (h_0[0] == h0Ccache[0] && h_1[0] == h0Ccache[1] && h_2[0] == h0Ccache[2] && SU3_breaking == h0Ccache[3]) {
1325 h0_updated = 1;
1326 } else {
1327 h0_updated = 0;
1328 h0Ccache[0] = h_0[0];
1329 h0Ccache[1] = h_1[0];
1330 h0Ccache[2] = h_2[0];
1331 h0Ccache[3] = SU3_breaking;
1332 }
1333
1334 if (h_0[1] == h1Ccache[0] && h_1[1] == h1Ccache[1] && h_2[1] == h1Ccache[2] && SU3_breaking == h1Ccache[3]) {
1335 h1_updated = 1;
1336 } else {
1337 h1_updated = 0;
1338 h1Ccache[0] = h_0[1];
1339 h1Ccache[1] = h_1[1];
1340 h1Ccache[2] = h_2[1];
1341 h1Ccache[3] = SU3_breaking;
1342 }
1343
1344 if (h_0[2] == h2Ccache[0] && h_1[2] == h2Ccache[1] && h_2[2] == h2Ccache[2] && SU3_breaking == h2Ccache[3]) {
1345 h2_updated = 1;
1346 } else {
1347 h2_updated = 0;
1348 h2Ccache[0] = h_0[2];
1349 h2Ccache[1] = h_1[2];
1350 h2Ccache[2] = h_2[2];
1351 h2Ccache[3] = SU3_breaking;
1352 }
1353 }
1354
1355 if (lep == QCD::NEUTRINO_1){
1356 H_V0updated = N_updated * VL0_updated * C_L_nunu_updated * C_R_nunu_updated * VR0_updated;
1357 H_V1updated = N_updated * VL1_updated * C_L_nunu_updated * C_R_nunu_updated * VR1_updated;
1358 H_V2updated = N_updated * VL2_updated * C_L_nunu_updated * C_R_nunu_updated * VR2_updated;
1359 H_A0updated = N_updated * VL0_updated * C_L_nunu_updated * C_R_nunu_updated * VR0_updated;
1360 H_A1updated = N_updated * VL1_updated * C_L_nunu_updated * C_R_nunu_updated * VR1_updated;
1361 H_A2updated = N_updated * VL2_updated * C_L_nunu_updated * C_R_nunu_updated * VR2_updated;
1362 } else {
1363 if (MM == H_V0cache(0) && Mb == H_V0cache(1)) {
1364 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;
1365 } else {
1366 H_V0updated = 0;
1367 H_V0cache(0) = MM;
1368 H_V0cache(1) = Mb;
1369 }
1370
1371 if (MM == H_V1cache(0) && Mb == H_V1cache(1)) {
1372 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;
1373 } else {
1374 H_V1updated = 0;
1375 H_V1cache(0) = MM;
1376 H_V1cache(1) = Mb;
1377 }
1378
1379 if (MM == H_V2cache(0) && Mb == H_V2cache(1)) {
1380 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;
1381 } else {
1382 H_V2updated = 0;
1383 H_V2cache(0) = MM;
1384 H_V2cache(1) = Mb;
1385 }
1386
1387 H_A0updated = N_updated * C_10_updated * VL0_updated * C_10p_updated * VR0_updated;
1388 H_A1updated = N_updated * C_10_updated * VL1_updated * C_10p_updated * VR1_updated;
1389 H_A2updated = N_updated * C_10_updated * VL2_updated * C_10p_updated * VR2_updated;
1390 }
1391
1392 if (Mb == H_Scache(0) && MW == H_Scache(1)) {
1393 H_Supdated = N_updated * C_S_updated * SL_updated * C_Sp_updated * SR_updated;
1394 } else {
1395 H_Supdated = 0;
1396 H_Scache(0) = Mb;
1397 H_Scache(1) = MW;
1398 }
1399
1400 if (Mb == H_Pcache(0) && MW == H_Pcache(1) && Mlep == H_Pcache(2) && Ms == H_Pcache(3)) {
1401 H_Pupdated = N_updated * C_P_updated * SL_updated * C_Pp_updated * SR_updated * C_10_updated * C_10p_updated;
1402 } else {
1403 H_Pupdated = 0;
1404 H_Pcache(0) = Mb;
1405 H_Pcache(1) = MW;
1406 H_Pcache(2) = Mlep;
1407 H_Pcache(3) = Ms;
1408
1409 }
1410
1411 if (MM == T_cache(0) && Mb == T_cache(1) && Mc == T_cache(2) &&
1412 mySM.getMesons(vectorM).getGegenalpha(0) == T_cache(3) && mySM.getMesons(vectorM).getGegenalpha(1) == T_cache(4)) {
1413 T_updated = 1;
1414 } else {
1415 T_updated = 0;
1416 T_cache(0) = MM;
1417 T_cache(1) = Mb;
1418 T_cache(2) = Mc;
1419 T_cache(3) = mySM.getMesons(vectorM).getGegenalpha(0);
1420 T_cache(4) = mySM.getMesons(vectorM).getGegenalpha(1);
1421 }
1422
1423 deltaTparpupdated = C_2Lh_updated * T_updated;
1424 deltaTparmupdated = C_2Lh_updated * C_8Lh_updated * T_updated;
1425 deltaTperpupdated = deltaTparpupdated;
1426
1427 I0_updated = F_updated * H_V0updated * H_A0updated * H_Pupdated * beta_updated * H_Supdated * deltaTparmupdated;
1428 I1_updated = F_updated * beta_updated * H_V1updated * H_V2updated * H_A1updated * H_A2updated * deltaTparmupdated;
1429 I2_updated = F_updated * beta_updated * H_V0updated * H_A0updated * deltaTparmupdated;
1430 I3_updated = F_updated * H_V1updated * H_V2updated * H_A1updated * H_A2updated * beta_updated * deltaTparmupdated;
1431 I4_updated = F_updated * H_V1updated * H_V2updated * H_A1updated * H_A2updated * deltaTparmupdated;
1432 I5_updated = F_updated * H_V0updated * H_V1updated * H_V2updated * H_A0updated * H_A1updated * H_A2updated * beta_updated * deltaTparmupdated;
1433 I6_updated = F_updated * H_V1updated * H_V2updated * H_A0updated * H_A1updated * H_A2updated * H_V0updated * beta_updated * H_Supdated * deltaTparmupdated;
1434 I7_updated = I4_updated * beta_updated;
1435 I8_updated = F_updated * beta_updated * H_Supdated * H_V0updated * deltaTparmupdated;
1436 I9_updated = I6_updated;
1437 I10_updated = I5_updated;
1438 I11_updated = I7_updated;
1439
1440 if (MM2 == Itree_cache(0) && mtau2 == Itree_cache(1) && MV2 == Itree_cache(2)) {
1441 Itree_updated = 1;
1442 } else {
1443 Itree_updated = 0;
1444 Itree_cache(0) = MM2;
1445 Itree_cache(1) = mtau2;
1446 Itree_cache(2) = MV2;
1447 }
1448
1449}
1450
1451/*******************************************************************************
1452 * Transverse Form Factors *
1453 * ****************************************************************************/
1454
1455double MVll::FF_fit(double q2, double a_0, double a_1, double a_2, double MR_2)
1456{
1457 return 1. / (1. - q2 / MR_2) * (a_0 + a_1 * (z(q2) - z_0) + a_2 * (z(q2) - z_0) * (z(q2) - z_0));
1458}
1459
1460double MVll::z(double q2)
1461{
1462 return ( sqrt(t_p - q2) - sqrt(t_p - t_0)) / (sqrt(t_p - q2) + sqrt(t_p - t_0));
1463}
1464
1465double MVll::z_DM(double q2)
1466{
1467 return (sqrt(t_p - q2) - sqrt(t_p - t_m)) / (sqrt(t_p - q2) + sqrt(t_p - t_m));
1468}
1469
1470double MVll::phi_f(double q2, double MRf_2)
1471{
1472 double z = z_DM(q2);
1473 double z_M = z_DM(MRf_2);
1474
1475 return 4.*rV/MM2*sqrt(2./3./Chi1plus/M_PI) * (1. + z)*pow(1. - z,1.5)/pow((1. + rV)*(1. - z)+2.*sqrt(rV)*(1. + z),4) * (z - z_M)/(1. - z_M*z);
1476}
1477
1478double MVll::phi_g(double q2, double MRg_2)
1479{
1480 double z = z_DM(q2);
1481 double z_M = z_DM(MRg_2);
1482
1483 return 16.*rV*rV*sqrt(2./3./Chi1minus/M_PI) * (1. + z)*(1. + z)*pow(1. - z,-0.5)/pow((1. + rV)*(1. - z)+2.*sqrt(rV)*(1. + z),4) * (z - z_M)/(1. - z_M*z);
1484}
1485
1486double MVll::phi_F1(double q2, double MRF1_2)
1487{
1488 double z = z_DM(q2);
1489 double z_M = z_DM(MRF1_2);
1490
1491 return 2.*rV/MM3*sqrt(4./3./Chi1plus/M_PI) * (1. + z)*pow(1. - z,2.5)/pow((1. + rV)*(1. - z)+2.*sqrt(rV)*(1. + z),5) * (z - z_M)/(1. - z_M*z);
1492}
1493
1494double MVll::phi_F2(double q2, double MRF2_2)
1495{
1496 double z = z_DM(q2);
1497 double z_M = z_DM(MRF2_2);
1498
1499 return 8.*rV*rV*sqrt(4./Chi0minus/M_PI) * (1. + z)*(1. + z)*pow(1. - z,-0.5)/pow((1. + rV)*(1. - z)+2.*sqrt(rV)*(1. + z),4) * (z - z_M)/(1. - z_M*z);
1500}
1501
1502double MVll::phi_T0(double q2, double MRT0_2)
1503{
1504 double z = z_DM(q2);
1505 double z_M = z_DM(MRT0_2);
1506
1507 return 2.*rV*(1. + rV)/MM*sqrt(4./3./ChiBB/M_PI) * (1. + z)*pow(1. - z,1.5)/pow((1. + rV)*(1. - z)+2.*sqrt(rV)*(1. + z),4) * (z - z_M)/(1. - z_M*z);
1508}
1509
1510double MVll::phi_T1(double q2, double MRT1_2)
1511{
1512 double z = z_DM(q2);
1513 double z_M = z_DM(MRT1_2);
1514
1515 return 32.*rV*rV/MM*sqrt(2./3./ChiTT/M_PI) * (1. + z)*(1. + z)*pow(1. - z,0.5)/pow((1. + rV)*(1. - z)+2.*sqrt(rV)*(1. + z),5) * (z - z_M)/(1. - z_M*z);
1516}
1517
1518double MVll::phi_T2(double q2, double MRT2_2)
1519{
1520 double z = z_DM(q2);
1521 double z_M = z_DM(MRT2_2);
1522
1523 return 4.*rV*(1. - rV*rV)/MM*sqrt(2./3./ChiBB/M_PI) * (1. + z)*pow(1. - z,2.5)/pow((1. + rV)*(1. - z)+2.*sqrt(rV)*(1. + z),5) * (z - z_M)/(1. - z_M*z);
1524}
1525
1526double MVll::f_DM(double q2, double a_0f, double a_1f, double a_2f, double MRf_2)
1527{
1528 double z = z_DM(q2);
1529 return (a_0f + a_1f*z + a_2f*z*z) / phi_f(q2, MRf_2);
1530}
1531
1532double MVll::g_DM(double q2, double a_0g, double a_1g, double a_2g, double MRg_2)
1533{
1534 double z = z_DM(q2);
1535 return (a_0g + a_1g*z + a_2g*z*z) / phi_g(q2, MRg_2);
1536}
1537
1538double MVll::F1_DM(double q2, double a_0F1, double a_1F1, double a_2F1, double MRF1_2)
1539{
1540 double z = z_DM(q2);
1541 return (a_0F1 + a_1F1*z + a_2F1*z*z) / phi_F1(q2, MRF1_2);
1542}
1543
1544double MVll::F2_DM(double q2, double a_0F2, double a_1F2, double a_2F2, double MRF2_2)
1545{
1546 double z = z_DM(q2);
1547 return (a_0F2 + a_1F2*z + a_2F2*z*z) / phi_F2(q2, MRF2_2);
1548}
1549
1550double MVll::T0_DM(double q2, double a_0T0, double a_1T0, double a_2T0, double MRT0_2)
1551{
1552 double z = z_DM(q2);
1553 return (a_0T0 + a_1T0*z + a_2T0*z*z) / phi_T0(q2, MRT0_2);
1554}
1555
1556double MVll::T1_DM(double q2, double a_0T1, double a_1T1, double a_2T1, double MRT1_2)
1557{
1558 double z = z_DM(q2);
1559 return (a_0T1 + a_1T1*z + a_2T1*z*z) / phi_T1(q2, MRT1_2);
1560}
1561
1562double MVll::T2_DM(double q2, double a_0T2, double a_1T2, double a_2T2, double MRT2_2)
1563{
1564 double z = z_DM(q2);
1565 return (a_0T2 + a_1T2*z + a_2T2*z*z) / phi_T2(q2, MRT2_2);
1566}
1567
1568double MVll::V(double q2)
1569{
1570 if (MVll_DM_flag) {
1571 return g_DM(q2, a_0g, a_1g, a_2g, MRg_2)*MMpMV/2.;
1572 } else {
1573 return FF_fit(q2, a_0V, a_1V, a_2V, MRV_2);
1574 }
1575}
1576
1577double MVll::A_0(double q2)
1578{
1579 if (MVll_DM_flag) {
1580 return F2_DM(q2, a_0F2, a_1F2, a_2F2, MRF2_2)/2.;
1581 } else {
1582 return FF_fit(q2, a_0A0, a_1A0, a_2A0, MRA0_2);
1583 }
1584}
1585
1586double MVll::A_1(double q2)
1587{
1588 if (MVll_DM_flag) {
1589 return f_DM(q2, a_0f, a_1f, a_2f, MRf_2)/MMpMV;
1590 } else {
1591 return FF_fit(q2, a_0A1, a_1A1, a_2A1, MRA1_2);
1592 }
1593}
1594
1595double MVll::A_2(double q2)
1596{
1597 double A12 = 0.;
1598 if (MVll_DM_flag) {
1599 A12 = F1_DM(q2, a_0F1, a_1F1, a_2F1, MRF1_2)/MMMV/8.;
1600 } else {
1601 A12 = FF_fit(q2, a_0A12, a_1A12, a_2A12, MRA12_2);
1602 }
1603
1604 return (MMpMV2 * (MM2mMV2 - q2) * A_1(q2) - 16. * MM * MV2 * MMpMV * A12) / lambda(q2);
1605}
1606
1607double MVll::T_1(double q2)
1608{
1609 if (MVll_DM_flag) {
1610 return T1_DM(q2, a_0T1, a_1T1, a_2T1, MRT1_2);
1611 } else {
1612 return FF_fit(q2, a_0T1, a_1T1, a_2T1, MRT1_2);
1613 }
1614}
1615
1616double MVll::T_2(double q2)
1617{
1618 if (MVll_DM_flag) {
1619 return T2_DM(q2, a_0T2, a_1T2, a_2T2, MRT2_2);
1620 } else {
1621 return FF_fit(q2, a_0T2, a_1T2, a_2T2, MRT2_2);
1622 }
1623}
1624
1625double MVll::V_0t(double q2)
1626{
1627 double A12 = 0.;
1628 if (MVll_DM_flag) {
1629 A12 = F1_DM(q2, a_0F1, a_1F1, a_2F1, MRF1_2)/MMMV/8.;
1630 } else {
1631 A12 = FF_fit(q2, a_0A12, a_1A12, a_2A12, MRA12_2);
1632 }
1633
1634 return fourMV / sqrt(q2) * A12;
1635}
1636
1637double MVll::V_p(double q2)
1638{
1639 return half * (onepMMoMV * A_1(q2) - sqrt(lambda(q2)) / (MM_MMpMV) * V(q2));
1640}
1641
1642double MVll::V_m(double q2)
1643{
1644 return half * (onepMMoMV * A_1(q2) + sqrt(lambda(q2)) / (MM_MMpMV) * V(q2));
1645}
1646
1647double MVll::T_0t(double q2)
1648{
1649 double T23 = 0.;
1650 if (MVll_DM_flag) {
1651 T23 = T0_DM(q2, a_0T0, a_1T0, a_2T0, MRT0_2)*MMpMV*MMpMV/4./MM;
1652 } else {
1653 T23 = FF_fit(q2, a_0T23, a_1T23, a_2T23, MRT23_2);
1654 }
1655
1656 return 2 * sqrt(q2) * MV / MM_MMpMV * T23;
1657}
1658
1659double MVll::T_p(double q2)
1660{
1661 return (MM2mMV2 * T_2(q2) - sqrt(lambda(q2)) * T_1(q2)) / twoMM2;
1662}
1663
1664double MVll::T_m(double q2)
1665{
1666 return (MM2mMV2 * T_2(q2) + sqrt(lambda(q2)) * T_1(q2)) / twoMM2;
1667}
1668
1669double MVll::S_L(double q2)
1670{
1671 return -sqrt(lambda(q2)) / twoMM_mbpms * A_0(q2);
1672}
1673
1674/*******************************************************************************
1675 * QCDF NLO *
1676 * ****************************************************************************/
1677
1678gslpp::complex MVll::A_Seidel(double q2, double mb2)
1679{
1680 double sh = q2 / mb2;
1681 double z = (4. * mb2) / q2;
1682 double lsh = log(sh);
1683 gslpp::complex acsq = arccot((gslpp::complex)sqrt(z - 1.));
1684 double sh2 = sh*sh;
1685 double osh2 = (1. - sh)*(1. - sh);
1686 return (-(104.) / (243.) * log((mb2) / (mu_b2)) + (4. * sh) / (27. * (1. - sh)) * (dilog((gslpp::complex)sh) + lsh * log(1. - sh))
1687 + (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))
1688 - (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)
1689 + (2. * sh) / (243. * osh2 * osh2) * (36. * acsq * acsq + M_PI2 * (-4. + 9. * sh - 9. * sh2 + 3. * sh2 * sh)));
1690}
1691
1692gslpp::complex MVll::B_Seidel(double q2, double mb2)
1693{
1694 double sh = q2 / mb2;
1695 double z = (4. * mb2) / q2;
1696 double sqrt_z_m_1 = sqrt(z - 1.);
1697 gslpp::complex x1 = 0.5 + gslpp::complex::i() / 2. * sqrt_z_m_1;
1698 gslpp::complex x2 = 0.5 - gslpp::complex::i() / 2. * sqrt_z_m_1;
1699 gslpp::complex x3 = 0.5 + gslpp::complex::i() / (2. * sqrt_z_m_1);
1700 gslpp::complex x4 = 0.5 - gslpp::complex::i() / (2. * sqrt_z_m_1);
1701 gslpp::complex lx1 = log(x1);
1702 gslpp::complex lx2 = log(x2);
1703 gslpp::complex lx3 = log(x3);
1704 gslpp::complex lx4 = log(x4);
1705 gslpp::complex lx2_x1 = lx2 - lx1;
1706 gslpp::complex lzm1 = log(z - 1.);
1707 gslpp::complex acsq = arccot((gslpp::complex)sqrt_z_m_1);
1708 double sh2 = sh*sh;
1709 double lsh = log(sh);
1710 double osh2 = (1. - sh)*(1. - sh);
1711 double lmb_mu = log(mb2 / mu_b2);
1712 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)
1713 + ((2. + sh) * sqrt_z_m_1) / (729. * sh) * (-48. * lmb_mu * acsq - 18. * M_PI * log(z - 1.) + 3. * gslpp::complex::i() * lzm1 * lzm1
1714 - 24. * gslpp::complex::i() * dilog(-x2 / x1) - 5. * M_PI2 * gslpp::complex::i()
1715 + 6. * gslpp::complex::i() * (-9. * lx1 * lx1 + lx2 * lx2 - 2. * lx4 * lx4 + 6. * lx1 * lx2 - 4. * lx1 * lx3 + 8. * lx1 * lx4)
1716 - 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))
1717 + 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
1718 - 288. - 8. * sh + 934. * sh2 - 692. * sh2 * sh + 18. * M_PI * gslpp::complex::i() * sh * (82. - 173. * sh + 73. * sh2))
1719 - 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
1720 + 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
1721 - M_PI2 * (54. - 53. * sh - 286. * sh2 + 612. * sh * sh2 - 446. * sh2 * sh2 + 113. * sh2 * sh2 * sh)));
1722}
1723
1724gslpp::complex MVll::C_Seidel(double q2)
1725{
1726 return -(16.) / (81.) * log((q2) / (mu_b2)) + (428.) / (243.) - (64.) / (27.) * gsl_sf_zeta_int(3) + (16.) / (81.) * M_PI * gslpp::complex::i();
1727 /* gsl_sf_zeta_int returns a double */
1728}
1729
1730gslpp::complex MVll::deltaC7_QCDF(double q2, bool conjugate, bool spline)
1731{
1732 if (zExpansion)
1733 return 0.;
1734 else {
1735 #if COMPUTECP && SPLINE
1736 if (spline && !conjugate) return gsl_spline_eval(spline_Re_deltaC7_QCDF, q2, acc_Re_deltaC7_QCDF);
1737 else if (spline && conjugate) return gsl_spline_eval(spline_Re_deltaC7_QCDF_conj, q2, acc_Re_deltaC7_QCDF_conj);
1738 #elif SPLINE
1739 if (spline) return gsl_spline_eval(spline_Re_deltaC7_QCDF, q2, acc_Re_deltaC7_QCDF);
1740 #endif
1741
1742 double muh = mu_b / mb_pole;
1743 double z = mc_pole * mc_pole / mb_pole / mb_pole;
1744 double sh = q2 / mb_pole / mb_pole;
1745 double sh2 = sh*sh;
1746
1747 #if FULLNLOQCDF_MVLL
1748 gslpp::complex A_Sdl = A_Seidel(q2, mb_pole*mb_pole); /* hep-ph/0403185v2.*/
1749 gslpp::complex Fu_17 = -A_Sdl; /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1750 gslpp::complex Fu_27 = 6. * A_Sdl; /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1751 #endif
1752 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*/
1753 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*/
1754 gslpp::complex F_87 = F87_0 + F87_1 * sh + F87_2 * sh2 + F87_3 * sh * sh2 - 8. / 9. * log(sh) * (sh + sh2 + sh * sh2);
1755
1756 if (!conjugate) {
1757 gslpp::complex delta = C_1 * F_17 + C_2 * F_27;
1758 gslpp::complex delta_t = C_8 * F_87 + delta;
1759 #if FULLNLOQCDF_MVLL
1760 gslpp::complex delta_u = delta + C_1 * Fu_17 + C_2 * Fu_27;
1761 return -alpha_s_mub / (4. * M_PI) * (delta_t - lambda_u / lambda_t * delta_u);
1762 #else
1763 return -alpha_s_mub / (4. * M_PI) * delta_t;
1764 #endif
1765 } else {
1766 gslpp::complex delta = C_1.conjugate() * F_17 + C_2.conjugate() * F_27;
1767 gslpp::complex delta_t = C_8.conjugate() * F_87 + delta;
1768 #if FULLNLOQCDF_MVLL
1769 gslpp::complex delta_u = delta + C_1.conjugate() * Fu_17 + C_2.conjugate() * Fu_27;
1770 return -alpha_s_mub / (4. * M_PI) * (delta_t - (lambda_u / lambda_t).conjugate() * delta_u);
1771 #else
1772 return -alpha_s_mub / (4. * M_PI) * delta_t;
1773 #endif
1774 }
1775 }
1776}
1777
1778gslpp::complex MVll::deltaC9_QCDF(double q2, bool conjugate, bool spline)
1779{
1780 if (zExpansion)
1781 return 0.;
1782 else {
1783 #if COMPUTECP && SPLINE
1784 if (spline && !conjugate) return gsl_spline_eval(spline_Re_deltaC9_QCDF, q2, acc_Re_deltaC9_QCDF);
1785 else if (spline && conjugate) return gsl_spline_eval(spline_Re_deltaC9_QCDF_conj, q2, acc_Re_deltaC9_QCDF_conj);
1786 #elif SPLINE
1787 if (spline) return gsl_spline_eval(spline_Re_deltaC9_QCDF, q2, acc_Re_deltaC9_QCDF);
1788 #endif
1789
1790 double muh = mu_b / mb_pole;
1791 double z = mc_pole * mc_pole / mb_pole / mb_pole;
1792 double sh = q2 / mb_pole / mb_pole;
1793 double sh2 = sh*sh;
1794
1795 #if FULLNLOQCDF_MVLL
1796 gslpp::complex B_Sdl = B_Seidel(q2, mb_pole*mb_pole); /* hep-ph/0403185v2.*/
1797 gslpp::complex C_Sdl = C_Seidel(q2); /* hep-ph/0403185v2.*/
1798 gslpp::complex Fu_19 = -(B_Sdl + 4. * C_Sdl); /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1799 gslpp::complex Fu_29 = -(-6. * B_Sdl + 3. * C_Sdl); /* sign different from hep-ph/0403185v2 but consistent with hep-ph/0412400 */
1800 #endif
1801 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*/
1802 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*/
1803 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));
1804
1805 if (!conjugate) {
1806 gslpp::complex delta = C_1 * F_19 + C_2 * F_29;
1807 gslpp::complex delta_t = C_8 * F_89 + delta;
1808 #if FULLNLOQCDF_MVLL
1809 gslpp::complex delta_u = delta + C_1 * Fu_19 + C_2 * Fu_29;
1810 return -alpha_s_mub / (4. * M_PI) * (delta_t - lambda_u / lambda_t * delta_u);
1811 #else
1812 return -alpha_s_mub / (4. * M_PI) * delta_t;
1813 #endif
1814 } else {
1815 gslpp::complex delta = C_1.conjugate() * F_19 + C_2.conjugate() * F_29;
1816 gslpp::complex delta_t = C_8.conjugate() * F_89 + delta;
1817 #if FULLNLOQCDF_MVLL
1818 gslpp::complex delta_u = delta + C_1.conjugate() * Fu_19 + C_2.conjugate() * Fu_29;
1819 return -alpha_s_mub / (4. * M_PI) * (delta_t - (lambda_u / lambda_t).conjugate() * delta_u);
1820 #else
1821 return -alpha_s_mub / (4. * M_PI) * delta_t;
1822 #endif
1823 }
1824 }
1825}
1826
1827gslpp::complex MVll::Cq34(bool conjugate)
1828{
1829 gslpp::complex T_t = C_3 + 4. / 3. * (C_4 + 12. * C_5 + 16. * C_6);
1830 gslpp::complex T_u = 0.; /* 0 for K*0, phi*/
1831 if (meson == QCD::B_P) T_u = -3. * C_2;
1832 else if (vectorM == QCD::PHI) T_t = T_t + 6. * (C_3 + 10. * C_5);
1833 if (!conjugate) return T_t + lambda_u / lambda_t * T_u;
1834 else return T_t + (lambda_u / lambda_t).conjugate() * T_u;
1835}
1836
1837gslpp::complex MVll::T_para_minus_WA(bool conjugate)
1838{
1839 return -spectator_charge * 4. * MM / mb_pole * Cq34(conjugate);
1840}
1841
1842gslpp::complex MVll::T_perp_WA_1()
1843{
1844 return -spectator_charge * 4. / mb_pole * (C_3 + 4. / 3. * (C_4 + 3. * C_5 + 4. * C_6));
1845}
1846
1847gslpp::complex MVll::T_perp_WA_2(bool conjugate)
1848{
1849 return spectator_charge * 2. / mb_pole * Cq34(conjugate);
1850}
1851
1852gslpp::complex MVll::T_perp_plus_O8(double q2, double u)
1853{
1854 double ubar = 1. - u;
1855 double ed = -1. / 3.;
1856
1857 return -(alpha_s_mub / (3. * M_PI))*4. * ed * C_8 / (u + ubar * q2 / MM2);
1858}
1859
1860gslpp::complex MVll::T_para_minus_O8(double q2, double u)
1861{
1862 double ubar = 1. - u;
1863
1864 return (alpha_s_mub / (3. * M_PI))*spectator_charge * 8. * C_8 / (ubar + u * q2 / MM2);
1865}
1866
1867gslpp::complex MVll::t_perp(double q2, double u, double m2)
1868{
1869 double EV = (MM2 - q2 + MV2) / (2. * MM);
1870 double ubar = 1. - u;
1871
1872 return (2. * MM) / (ubar * EV) * I1(q2, u, m2) + q2 / (ubar * ubar * EV * EV) * B0diff(q2, u, m2);
1873
1874}
1875
1876gslpp::complex MVll::t_para(double q2, double u, double m2)
1877{
1878 double EV = (MM2pMV2 - q2) / (2. * MM);
1879 double ubar = 1. - u;
1880 return (2. * MM) / (ubar * EV) * I1(q2, u, m2) + (ubar * MM2 + u * q2) / (ubar * ubar * EV * EV) * B0diff(q2, u, m2);
1881}
1882
1883gslpp::complex MVll::I1(double q2, double u, double m2)
1884{
1885 if (m2 == 0.) return 1.;
1886
1887 ubar = 1. - u;
1888 xp = 0.5 + sqrt(0.25 - ((gslpp::complex) m2) / (ubar * MM2 + u * q2));
1889 xm = 0.5 - sqrt(0.25 - ((gslpp::complex) m2) / (ubar * MM2 + u * q2));
1890 yp = 0.5 + sqrt(0.25 - ((gslpp::complex) m2) / q2);
1891 ym = 0.5 - sqrt(0.25 - ((gslpp::complex) m2) / q2);
1892 L1xp = log(1. - 1. / xp) * log(1. - xp) - M_PI2osix + dilog(xp / (xp - 1.));
1893 L1xm = log(1. - 1. / xm) * log(1. - xm) - M_PI2osix + dilog(xm / (xm - 1.));
1894 L1yp = log(1. - 1. / yp) * log(1. - yp) - M_PI2osix + dilog(yp / (yp - 1.));
1895 L1ym = log(1. - 1. / ym) * log(1. - ym) - M_PI2osix + dilog(ym / (ym - 1.));
1896
1897 return 1. + 2. * m2 / ubar / (MM2 - q2)*(L1xp + L1xm - L1yp - L1ym);
1898}
1899
1900gslpp::complex MVll::B0diff(double q2, double u, double m2)
1901{
1902 double ubar = 1. - u;
1903
1904 if (m2 == 0.) return -log((gslpp::complex)(-(2. / q2))) + log((gslpp::complex)(-(2. / (q2 * u + MM2 * ubar))));
1905 else return B0(ubar * MM2 + u * q2, m2) - B0(q2, m2);
1906}
1907
1908gslpp::complex MVll::B0(double s, double m2)
1909{
1910 if (4. * m2 / s == 1.) return gslpp::complex(0.);
1911 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.));
1912}
1913
1914gslpp::complex MVll::h_func(double s, double m2)
1915{
1916 if (m2 == 0.) return 8. / 27. + 4. * gslpp::complex::i() * M_PI / 9. + 8. * log(mu_b) / 9. - 4. * log(s) / 9.;
1917 if (s == 0.) return -4. / 9. * (1. + log(m2 / mu_b / mu_b));
1918
1919 double z = 4 * m2 / s;
1920 gslpp::complex term;
1921 if (z > 1) term = atan(1. / sqrt(z - 1.));
1922 else term = log((1. + sqrt(1. - z)) / sqrt(z)) - ihalfMPI;
1923
1924 return -4. / 9. * log(m2 / mu_b / mu_b) + 8. / 27. + 4. / 9. * z - 4. / 9. * (2. + z) * sqrt(std::abs(z - 1.)) * term;
1925
1926}
1927
1928gslpp::complex MVll::T_perp_plus_QSS(double q2, double u, bool conjugate)
1929{
1930 gslpp::complex t_perp_mc = t_perp(q2, u, mc_pole * mc_pole);
1931 double eu = 0.666666667;
1932#if FULLNLOQCDF_MVLL
1933 gslpp::complex t_perp_mb = t_perp(q2, u, mb_pole*mb_pole);
1934 gslpp::complex t_perp_0 = t_perp(q2, u, 0.);
1935 double ed = -0.333333333;
1936
1937 gslpp::complex T_t = (eu * t_perp_mc * (-C_1 / 6. + C_2 + 6. * C_6)
1938 + 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.))
1939 + ed * t_perp_0 * (C_3 - C_4/6. + 16. * C_5 - 8. * C_6/3.));
1940
1941 gslpp::complex T_u = eu * (t_perp_mc - t_perp_0)*(C_2 - C_1 / 6.);
1942
1943 if (!conjugate) return alpha_s_mub / (3. * M_PI) * MM / (2. * mb_pole)*(T_t + lambda_u / lambda_t * T_u);
1944 else return alpha_s_mub / (3. * M_PI) * MM / (2. * mb_pole)*(T_t + (lambda_u / lambda_t).conjugate() * T_u);
1945#else
1946 return alpha_s_mub / (3. * M_PI) * MM / (2. * mb_pole)*(eu * t_perp_mc * (-C_1 / 6. + C_2 + 6. * C_6));
1947#endif
1948}
1949
1950gslpp::complex MVll::T_para_plus_QSS(double q2, double u, bool conjugate)
1951{
1952 gslpp::complex t_para_mc = t_para(q2, u, mc_pole * mc_pole);
1953 double eu = 0.666666667;
1954#if FULLNLOQCDF_MVLL
1955 gslpp::complex t_para_mb = t_para(q2, u, mb_pole*mb_pole);
1956 gslpp::complex t_para_0 = t_para(q2, u, 0.);
1957 double ed = -0.333333333;
1958
1959 gslpp::complex T_t = (eu * t_para_mc * (-C_1 / 6. + C_2 + 6. * C_6)
1960 + ed * t_para_mb * (C_3 - C_4/6. + 16.*C_5 + 10.*C_6/3.)
1961 + ed * t_para_0 * (C_3 - C_4/6. + 16.*C_5 - 8.*C_6/3.));
1962
1963 gslpp::complex T_u = eu * (t_para_mc - t_para_0) * (C_2 - C_1/6.);
1964
1965 if (!conjugate) return alpha_s_mub / (3. * M_PI) * MM / mb_pole * (T_t + lambda_u / lambda_t * T_u);
1966 else return alpha_s_mub / (3. * M_PI) * MM / mb_pole * (T_t + (lambda_u / lambda_t).conjugate() * T_u);
1967#else
1968 return alpha_s_mub / (3. * M_PI) * MM / mb_pole * (eu * t_para_mc * (-C_1 / 6. + C_2 + 6. * C_6));
1969#endif
1970}
1971
1972gslpp::complex MVll::T_para_minus_QSS(double q2, double u, bool conjugate)
1973{
1974 double ubar = 1. - u;
1975 gslpp::complex h_mc = h_func(ubar * MM2 + u*q2, mc_pole * mc_pole);
1976#if FULLNLOQCDF_MVLL
1977 gslpp::complex h_mb = h_func(ubar*MM2 + u*q2, mb_pole*mb_pole);
1978 gslpp::complex h_0 = h_func(ubar*MM2 + u*q2, 0);
1979
1980 gslpp::complex T_t = (h_mc * (-C_1 / 6. + C_2 + C_4 + 10. * C_6)
1981 + h_mb * (C_3 + 5.*C_4/6. + 16.*C_5 + 22.*C_6/3.)
1982 + h_0 * (C_3 + 17.*C_4/6. + 16.*C_5 + 82.*C_6/3.)
1983 - 8./27. * (-15.*C_4/2. + 12.*C_5 - 32.*C_6));
1984
1985 gslpp::complex T_u = (h_mc - h_0)*(C_2 - C_1/6.);
1986
1987 if (!conjugate) return alpha_s_mub / (3. * M_PI) * spectator_charge * 6. * MM / mb_pole * (T_t + lambda_u / lambda_t * T_u);
1988 else return alpha_s_mub / (3. * M_PI) * spectator_charge * 6. * MM / mb_pole * (T_t + (lambda_u / lambda_t).conjugate() * T_u);
1989#else
1990 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));
1991#endif
1992}
1993
1994double MVll::phi_V(double u)
1995{
1996 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.)));
1997}
1998
1999gslpp::complex MVll::lambda_B_minus(double q2)
2000{
2001 double w0 = mySM.getMesons(meson).getLambdaM();
2002 return 1. / (exp(-q2 / MM / w0) / w0 * (-gsl_sf_expint_Ei(q2 / MM / w0) + gslpp::complex::i() * M_PI));
2003}
2004
2005double MVll::T_perp_real(double q2, double u, bool conjugate)
2006{
2007 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));
2008#if FULLNLOQCDF_MVLL
2009 double ubar = 1. - u;
2010
2011 T_amp += N_QCDF/(ubar + u*q2/MM2) * phi_V(u) * T_perp_WA_1()
2012 + N_QCDF/mySM.getMesons(meson).getLambdaM() * fpara/fperp * MV/(1. - q2/MM2) * T_perp_WA_2(conjugate);
2013 /*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]*/
2014#endif
2015 return T_amp.real();
2016}
2017
2018double MVll::T_perp_imag(double q2, double u, bool conjugate)
2019{
2020 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));
2021#if FULLNLOQCDF_MVLL
2022 double ubar = 1. - u;
2023
2024 T_amp += N_QCDF/(ubar + u*q2/MM2) * phi_V(u) * T_perp_WA_1()
2025 + N_QCDF/mySM.getMesons(meson).getLambdaM() * fpara/fperp * MV/(1. - q2/MM2) * T_perp_WA_2(conjugate);
2026 /*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]*/
2027#endif
2028 return T_amp.imag();
2029}
2030
2031double MVll::T_para_real(double q2, double u, bool conjugate)
2032{
2033 double N = N_QCDF * (MV / ((MM2pMV2 - q2) / (2. * MM)));
2034
2035 gslpp::complex T_amp = (N / lambda_B_minus(q2) * (T_para_minus_O8(q2, u) + T_para_minus_QSS(q2, u, conjugate))
2036 + N / mySM.getMesons(meson).getLambdaM() * T_para_plus_QSS(q2, u, conjugate)) * phi_V(u);
2037#if FULLNLOQCDF_MVLL
2038 T_amp += N / lambda_B_minus(q2) * T_para_minus_WA(conjugate)* phi_V(u);
2039#endif
2040 return sqrt(q2) * T_amp.real();
2041}
2042
2043double MVll::T_para_imag(double q2, double u, bool conjugate)
2044{
2045 double N = N_QCDF * (MV / ((MM2pMV2 - q2) / (2. * MM)));
2046
2047 gslpp::complex T_amp = (N / lambda_B_minus(q2) * (/* + */T_para_minus_O8(q2, u) + T_para_minus_QSS(q2, u, conjugate))
2048 + N / mySM.getMesons(meson).getLambdaM() * T_para_plus_QSS(q2, u, conjugate)) * phi_V(u);
2049#if FULLNLOQCDF_MVLL
2050 T_amp += N / lambda_B_minus(q2) * T_para_minus_WA(conjugate) * phi_V(u);
2051#endif
2052 return sqrt(q2) * T_amp.imag();
2053}
2054
2055double MVll::T_perp_real(double q2, bool conjugate)
2056{
2057 FS = convertToGslFunction(bind(&MVll::T_perp_real, &(*this), q2, _1, conjugate));
2058 gsl_integration_cquad(&FS, 0., 1., 1.e-2, 1.e-1, w_sigma, &avaSigma, &errSigma, NULL);
2059
2060 return avaSigma;
2061}
2062
2063double MVll::T_perp_imag(double q2, bool conjugate)
2064{
2065 FS = convertToGslFunction(bind(&MVll::T_perp_imag, &(*this), q2, _1, conjugate));
2066 gsl_integration_cquad(&FS, 0., 1., 1.e-2, 1.e-1, w_sigma, &avaSigma, &errSigma, NULL);
2067
2068 return avaSigma;
2069}
2070
2071double MVll::T_para_real(double q2, bool conjugate)
2072{
2073 FS = convertToGslFunction(bind(&MVll::T_para_real, &(*this), q2, _1, conjugate));
2074 gsl_integration_cquad(&FS, 0., 1., 1.e-2, 1.e-1, w_sigma, &avaSigma, &errSigma, NULL);
2075
2076 return avaSigma;
2077}
2078
2079double MVll::T_para_imag(double q2, bool conjugate)
2080{
2081 FS = convertToGslFunction(bind(&MVll::T_para_imag, &(*this), q2, _1, conjugate));
2082 gsl_integration_cquad(&FS, 0., 1., 1.e-2, 1.e-1, w_sigma, &avaSigma, &errSigma, NULL);
2083
2084 return avaSigma;
2085}
2086
2087double MVll::QCDF_fit_func(double* x, double* p)
2088{
2089 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];
2090}
2091
2092void MVll::fit_QCDF_func()
2093{
2094 int dim = 0;
2095 for (double i = 0.001; i < 8.6; i += 0.5) {
2096 myq2.push_back(i);
2097 Re_T_perp.push_back(T_perp_real(i, false));
2098 Im_T_perp.push_back(T_perp_imag(i, false));
2099 Re_T_para.push_back(T_para_real(i, false));
2100 Im_T_para.push_back(T_para_imag(i, false));
2101
2102#if COMPUTECP
2103 Re_T_perp_conj.push_back(T_perp_real(i, true));
2104 Im_T_perp_conj.push_back(T_perp_imag(i, true));
2105 Re_T_para_conj.push_back(T_para_real(i, true));
2106 Im_T_para_conj.push_back(T_para_imag(i, true));
2107#endif
2108 dim++;
2109 }
2110
2111 gr1 = TGraph(dim, myq2.data(), Re_T_perp.data());
2112 QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Re_T_perp");
2113 Re_T_perp_res = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
2114 Re_T_perp.clear();
2115
2116 gr1 = TGraph(dim, myq2.data(), Im_T_perp.data());
2117 QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Im_T_perp");
2118 Im_T_perp_res = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
2119 Im_T_perp.clear();
2120
2121 gr1 = TGraph(dim, myq2.data(), Re_T_para.data());
2122 QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Re_T_para");
2123 Re_T_para_res = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
2124 Re_T_para.clear();
2125
2126 gr1 = TGraph(dim, myq2.data(), Im_T_para.data());
2127 QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Im_T_para");
2128 Im_T_para_res = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
2129 Im_T_para.clear();
2130
2131#if COMPUTECP
2132 gr1 = TGraph(dim, myq2.data(), Re_T_perp_conj.data());
2133 QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Re_T_perp_conj");
2134 Re_T_perp_res_conj = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
2135 Re_T_perp_conj.clear();
2136
2137 gr1 = TGraph(dim, myq2.data(), Im_T_perp_conj.data());
2138 QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Im_T_perp_conj");
2139 Im_T_perp_res_conj = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
2140 Im_T_perp_conj.clear();
2141
2142 gr1 = TGraph(dim, myq2.data(), Re_T_para_conj.data());
2143 QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Re_T_para_conj");
2144 Re_T_para_res_conj = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
2145 Re_T_para_conj.clear();
2146
2147 gr1 = TGraph(dim, myq2.data(), Im_T_para_conj.data());
2148 QCDFfit = TF1("QCDFfit", this, &MVll::QCDF_fit_func, 0.001, 8.51, 7, "MVll", "Im_T_para_conj");
2149 Im_T_para_res_conj = gr1.Fit(&QCDFfit, "SQN0+rob=0.99");
2150 Im_T_para_conj.clear();
2151#endif
2152
2153 myq2.clear();
2154}
2155
2156void MVll::spline_QCDF_func()
2157{
2158 int dim = GSL_INTERP_DIM;
2159 int dim_DC = GSL_INTERP_DIM_DC;
2160 double min = 0.001;
2161 double interval = (9.9 - min) / ((double) dim);
2162 double interval_DC = (9.9 - min) / ((double) dim_DC);
2163 double q2_spline[dim];
2164 double fq2_Re_T_perp[dim], fq2_Im_T_perp[dim], fq2_Re_T_para[dim], fq2_Im_T_para[dim];
2165 double q2_spline_DC[dim_DC];
2166 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];
2167#if COMPUTECP
2168 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];
2169 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];
2170#endif
2171
2172 for (int i = 0; i < dim; i++) {
2173 q2_spline[i] = min + (double) i*interval;
2174 fq2_Re_T_perp[i] = T_perp_real(q2_spline[i], false);
2175 fq2_Im_T_perp[i] = T_perp_imag(q2_spline[i], false);
2176 fq2_Re_T_para[i] = T_para_real(q2_spline[i], false);
2177 fq2_Im_T_para[i] = T_para_imag(q2_spline[i], false);
2178
2179#if COMPUTECP
2180 fq2_Re_T_perp_conj[i] = T_perp_real(q2_spline[i], true);
2181 fq2_Im_T_perp_conj[i] = T_perp_imag(q2_spline[i], true);
2182 fq2_Re_T_para_conj[i] = T_para_real(q2_spline[i], true);
2183 fq2_Im_T_para_conj[i] = T_para_imag(q2_spline[i], true);
2184#endif
2185 }
2186 for (int i = 0; i < dim_DC; i++) {
2187 q2_spline_DC[i] = min + (double) i*interval_DC;
2188 fq2_Re_deltaC7_QCDF[i] = deltaC7_QCDF(q2_spline_DC[i], false).real();
2189 fq2_Im_deltaC7_QCDF[i] = deltaC7_QCDF(q2_spline_DC[i], false).imag();
2190 fq2_Re_deltaC9_QCDF[i] = deltaC9_QCDF(q2_spline_DC[i], false).real();
2191 fq2_Im_deltaC9_QCDF[i] = deltaC9_QCDF(q2_spline_DC[i], false).imag();
2192
2193#if COMPUTECP
2194 fq2_Re_deltaC7_QCDF_conj[i] = deltaC7_QCDF(q2_spline_DC[i], true).real();
2195 fq2_Im_deltaC7_QCDF_conj[i] = deltaC7_QCDF(q2_spline_DC[i], true).imag();
2196 fq2_Re_deltaC9_QCDF_conj[i] = deltaC9_QCDF(q2_spline_DC[i], true).real();
2197 fq2_Im_deltaC9_QCDF_conj[i] = deltaC9_QCDF(q2_spline_DC[i], true).imag();
2198#endif
2199 }
2200
2201 gsl_spline_init(spline_Re_T_perp, q2_spline, fq2_Re_T_perp, dim);
2202 gsl_spline_init(spline_Im_T_perp, q2_spline, fq2_Im_T_perp, dim);
2203 gsl_spline_init(spline_Re_T_para, q2_spline, fq2_Re_T_para, dim);
2204 gsl_spline_init(spline_Im_T_para, q2_spline, fq2_Im_T_para, dim);
2205
2206 gsl_spline_init(spline_Re_deltaC7_QCDF, q2_spline_DC, fq2_Re_deltaC7_QCDF, dim_DC);
2207 gsl_spline_init(spline_Im_deltaC7_QCDF, q2_spline_DC, fq2_Im_deltaC7_QCDF, dim_DC);
2208 gsl_spline_init(spline_Re_deltaC9_QCDF, q2_spline_DC, fq2_Re_deltaC9_QCDF, dim_DC);
2209 gsl_spline_init(spline_Im_deltaC9_QCDF, q2_spline_DC, fq2_Im_deltaC9_QCDF, dim_DC);
2210
2211#if COMPUTECP
2212 gsl_spline_init(spline_Re_T_perp_conj, q2_spline, fq2_Re_T_perp_conj, dim);
2213 gsl_spline_init(spline_Im_T_perp_conj, q2_spline, fq2_Im_T_perp_conj, dim);
2214 gsl_spline_init(spline_Re_T_para_conj, q2_spline, fq2_Re_T_para_conj, dim);
2215 gsl_spline_init(spline_Im_T_para_conj, q2_spline, fq2_Im_T_para_conj, dim);
2216
2217 gsl_spline_init(spline_Re_deltaC7_QCDF_conj, q2_spline_DC, fq2_Re_deltaC7_QCDF_conj, dim_DC);
2218 gsl_spline_init(spline_Im_deltaC7_QCDF_conj, q2_spline_DC, fq2_Im_deltaC7_QCDF_conj, dim_DC);
2219 gsl_spline_init(spline_Re_deltaC9_QCDF_conj, q2_spline_DC, fq2_Re_deltaC9_QCDF_conj, dim_DC);
2220 gsl_spline_init(spline_Im_deltaC9_QCDF_conj, q2_spline_DC, fq2_Im_deltaC9_QCDF_conj, dim_DC);
2221#endif
2222
2223}
2224
2225gslpp::complex MVll::T_minus(double q2, bool conjugate)
2226{
2227 if (zExpansion)
2228 return 0.;
2229 else {
2230 #if COMPUTECP && SPLINE
2231 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));
2232 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));
2233 #elif SPLINE
2234 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));
2235 #endif
2236
2237 #if COMPUTECP && !SPLINE
2238 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())));
2239 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())));
2240 #elif !SPLINE
2241 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())));
2242 #endif
2243 }
2244
2245}
2246
2247gslpp::complex MVll::T_0(double q2, bool conjugate)
2248{
2249 if (zExpansion)
2250 return 0.;
2251 else {
2252 #if COMPUTECP && SPLINE
2253 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));
2254 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));
2255 #elif SPLINE
2256 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));
2257 #endif
2258
2259 #if COMPUTECP && !SPLINE
2260 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())));
2261 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())));
2262 #elif !SPLINE
2263 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())));
2264 #endif
2265 }
2266}
2267
2268/*******************************************************************************
2269 * Helicity amplitudes *
2270 * ****************************************************************************/
2271gslpp::complex MVll::H(double q2, double m2, double mu2)
2272{
2273 double x = 4. * m2 / q2;
2274 gslpp::complex par;
2275
2276 if (x > 1.) par = sqrt(x - 1.) * atan(1. / sqrt(x - 1.));
2277 else par = sqrt(1. - x) * (log((1. + sqrt(1. - x)) / sqrt(x)) - ihalfMPI);
2278
2279 return -fournineth * (log(m2 / mu2) - twothird - x) - fournineth * (2. + x) * par;
2280}
2281
2282gslpp::complex MVll::H_0(double q2)
2283{
2284 return (H_0_pre - fournineth * log(q2 / mu_b2));
2285}
2286
2287gslpp::complex MVll::Y(double q2)
2288{
2289 if (zExpansion)
2290 return 0.;
2291 else
2292 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;
2293}
2294
2295gslpp::complex MVll::funct_g(double q2)
2296{
2297 if (q2 < 4. * Mc * Mc)
2298 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.)));
2299 else
2300 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));
2301}
2302
2303gslpp::complex MVll::DeltaC9_KD(double q2, int com)
2304{
2305 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];
2306 /* C_1 = -0.267 and C_2 = 1.117 in KMPW */
2307}
2308
2309gslpp::complex MVll::zh(double q2)
2310{
2311 return ( sqrt(s_p - q2) - sqrt(s_p - s_0)) / (sqrt(s_p - q2) + sqrt(s_p - s_0));
2312}
2313
2314gslpp::complex MVll::P(double q2)
2315{
2316 gslpp::complex facmj2 = ( zh(q2) - zh(mJ2) ) / ( 1. - zh(q2)*zh(mJ2).conjugate() );
2317 if(fabs(q2 - mJ2)< 1.e-5) facmj2 = 1/(4.*(mJ2 - s_p));
2318 gslpp::complex facmPsi2S2 = ( zh(q2) - zh(mPsi2S2) ) / ( 1. - zh(q2)*zh(mPsi2S2).conjugate() );
2319 if(fabs(q2 - mPsi2S2)< 1.e-5) facmPsi2S2 = 1/(4.*(mPsi2S2 - s_p));
2320 // at the pole it returns directly the residue, i.e. Lim_{q2->mres2} P(q2)/(q2-mres2)
2321 return facmj2*facmPsi2S2;
2322}
2323
2324gslpp::complex MVll::Phi_1(double q2)
2325{
2326 return zh(q2)*1. - rho_0*1;
2327}
2328
2329gslpp::complex MVll::Phi_1_st(double q2)
2330{
2331 return 1. - rho_0*zh(q2)*1.;
2332}
2333
2334gslpp::complex MVll::Phi_2(double q2)
2335{
2336 return zh(q2)*Phi_1(q2) - rho_1*Phi_1_st(q2);
2337}
2338
2339gslpp::complex MVll::Phi_2_st(double q2)
2340{
2341 return Phi_1_st(q2) - rho_1*zh(q2)*Phi_1(q2);
2342}
2343
2344gslpp::complex MVll::Phi_3(double q2)
2345{
2346 return zh(q2)*Phi_2(q2) - rho_2*Phi_2_st(q2);
2347}
2348
2349gslpp::complex MVll::Phi_3_st(double q2)
2350{
2351 return Phi_2_st(q2) - rho_2*zh(q2)*Phi_2(q2);
2352}
2353
2354gslpp::complex MVll::Phi_4(double q2)
2355{
2356 return zh(q2)*Phi_3(q2) - rho_3*Phi_3_st(q2);
2357}
2358
2359gslpp::complex MVll::Phi_4_st(double q2)
2360{
2361 return Phi_3_st(q2) - rho_3*zh(q2)*Phi_3(q2);
2362}
2363
2364gslpp::complex MVll::Phi_5(double q2)
2365{
2366 return zh(q2)*Phi_4(q2) - rho_4*Phi_4_st(q2);
2367}
2368
2369gslpp::complex MVll::Phi_5_st(double q2)
2370{
2371 return Phi_4_st(q2) - rho_4*zh(q2)*Phi_4(q2);
2372}
2373
2374gslpp::complex MVll::Phi_6(double q2)
2375{
2376 return zh(q2)*Phi_5(q2) - rho_5*Phi_5_st(q2);
2377}
2378
2379gslpp::complex MVll::Phi_6_st(double q2)
2380{
2381 return Phi_5_st(q2) - rho_5*zh(q2)*Phi_5(q2);
2382}
2383
2384gslpp::complex MVll::p0()
2385{
2386 return 1. / sqrt(twoalphaBtoKst);
2387}
2388
2389gslpp::complex MVll::p1(double q2)
2390{
2391 return Phi_1(q2) / sqrt(twoalphaBtoKst * onemrho_0_2);
2392}
2393
2394gslpp::complex MVll::p2(double q2)
2395{
2396 return Phi_2(q2) / sqrt(twoalphaBtoKst * onemrho_0_2 * onemrho_1_2);
2397}
2398
2399gslpp::complex MVll::p3(double q2)
2400{
2401 return Phi_3(q2) / sqrt(twoalphaBtoKst * onemrho_0_2 * onemrho_1_2 * onemrho_2_2);
2402}
2403
2404gslpp::complex MVll::p4(double q2)
2405{
2406 return Phi_4(q2) / sqrt(twoalphaBtoKst * onemrho_0_2 * onemrho_1_2 * onemrho_2_2 * onemrho_3_2);
2407}
2408
2409gslpp::complex MVll::p5(double q2)
2410{
2411 return Phi_5(q2) / sqrt(twoalphaBtoKst * onemrho_0_2 * onemrho_1_2 * onemrho_2_2 * onemrho_3_2 * onemrho_4_2);
2412}
2413
2414gslpp::complex MVll::p6(double q2)
2415{
2416 return Phi_6(q2) / sqrt(twoalphaBtoKst * onemrho_0_2 * onemrho_1_2 * onemrho_2_2 * onemrho_3_2 * onemrho_4_2 * onemrho_5_2);
2417}
2418
2419gslpp::complex MVll::phi_1(double q2)
2420{
2421 return - sqrt( 2.*sqrt((4.*mD2-Q2)*(4.*mD2-s_0)) + 8.*mD2 - Q2 - s_0 ) / ( 2.*sqrt((4.*mD2-Q2)*(4.*mD2-s_0)) + 8.*mD2 + Q2*(zh(q2)-1.) - s_0*(zh(q2)+1.) ) ;
2422}
2423
2424gslpp::complex MVll::phi_2(double q2)
2425{
2426 gslpp::complex zhm1_2 = (zh(q2)-1.)*(zh(q2)-1.);
2427 gslpp::complex zhp1_2 = (zh(q2)+1.)*(zh(q2)+1.);
2428
2429 return sqrt( MM4*zhm1_2*zhm1_2 - 2.*MM2*zhm1_2*(-16.*mD2*zh(q2) + MV2*zhm1_2 + s_0*zhp1_2) + (16.*mD2*zh(q2) + MV2*zhm1_2 - s_0*zhp1_2)*(16.*mD2*zh(q2) + MV2*zhm1_2 - s_0*zhp1_2) );
2430}
2431
2432gslpp::complex MVll::phi_3(double q2)
2433{
2434 return sqrt( 8.*mD2 + 4.*sqrt(4.*mD2*mD2 - mD2*s_0) - s_0 ) / ( -8.*mD2 - 4.*sqrt(4.*mD2*mD2 - mD2*s_0) + s_0*(zh(q2)+1.) ) ;
2435}
2436
2437gslpp::complex MVll::phi_4(double q2)
2438{
2439 return 1. / sqrt( s_0*(zh(q2)+1.)*(zh(q2)+1.) - 16.*mD2*zh(q2) ) ;
2440}
2441
2442gslpp::complex MVll::DeltaC9_zExpansion(double q2, int tran)
2443{
2444 gslpp::complex z = zh(q2);
2445
2446 gslpp::complex invpref = 4.*M_PI*sqrt(2.*(4.*mD2-s_0)/3./chiOPE)*sqrt(1+zh(q2)) * P(q2);
2447
2448 if (tran == 0) {
2449 invpref *= MM4 * pow(1.-zh(q2),4.5) * phi_1(q2)*phi_1(q2)*phi_1(q2) * sqrt(phi_2(q2)) * phi_3(q2)*phi_3(q2) * phi_4(q2)*phi_4(q2);
2450
2451 return 1./invpref * (beta_0[0] + beta_0[1]*z + beta_0[2]*z*z + beta_0[3]*z*z*z + beta_0[4]*z*z*z*z + beta_0[5]*z*z*z*z*z + beta_0[6]*z*z*z*z*z*z);
2452 } else if (tran == 1) { // parallel
2453 invpref *= MM2*MM * pow(1.-zh(q2),3.5) * phi_1(q2)*phi_1(q2)*phi_1(q2) * sqrt(phi_2(q2)) * phi_3(q2)*phi_3(q2)*phi_3(q2);
2454
2455 return 1./invpref * (beta_1[0] + beta_1[1]*z + beta_1[2]*z*z + beta_1[3]*z*z*z + beta_1[4]*z*z*z*z + beta_1[5]*z*z*z*z*z + beta_1[6]*z*z*z*z*z*z);
2456 } else { // perpendicular
2457 invpref *= MM2*MM * pow(1.-zh(q2),3.5) * phi_1(q2)*phi_1(q2)*phi_1(q2) * sqrt(phi_2(q2)) * phi_3(q2)*phi_3(q2)*phi_3(q2);
2458
2459 return 1./invpref * (beta_2[0] + beta_2[1]*z + beta_2[2]*z*z + beta_2[3]*z*z*z + beta_2[4]*z*z*z*z + beta_2[5]*z*z*z*z*z + beta_2[6]*z*z*z*z*z*z);
2460 }
2461}
2462
2463gslpp::complex MVll::h_lambda(int hel, double q2)
2464{
2465 if (zExpansion) {
2466 if (hel == 0)
2467 return DeltaC9_zExpansion(q2, 0) * MM / sqrt(q2);
2468 else if (hel == 1)
2469 return (DeltaC9_zExpansion(q2, 1) - DeltaC9_zExpansion(q2, 2)) / sqrt(2.);
2470 else
2471 return (DeltaC9_zExpansion(q2, 1) + DeltaC9_zExpansion(q2, 2)) / sqrt(2.);
2472 } else if (dispersion) {
2473 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)));
2474 else if (hel == 1) {
2475 if (q2 == 0.) return SU3_breaking * (-1. / (MM2 * 16. * M_PI * M_PI) * (
2476 (MMpMV * A_1(0.)) / (2. * MM) * ((-h_0[1] + h_2[1]) / (1. + h_1[1] / mJ2)) * exp_Phase[1]
2477 - sqrt(lambda(0.)) / (2. * MM * MMpMV) * V(0.) * ((-h_0[0] + h_2[0]) / (1. + h_1[0] / mJ2)) * exp_Phase[0]));
2478 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)));
2479 } else {
2480 if (q2 == 0.) return SU3_breaking * (-1. / (MM2 * 16. * M_PI * M_PI) *
2481 ((MMpMV * A_1(0.)) / (2. * MM) * ((-h_0[1] + h_2[1]) / (1. + h_1[1] / mJ2)) * exp_Phase[1]
2482 + sqrt(lambda(0.)) / (2. * MM * MMpMV) * V(0.) * ((-h_0[0] + h_2[0]) / (1. + h_1[0] / mJ2)) * exp_Phase[0]));
2483 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)));
2484 }
2485 } else {
2486 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]));
2487 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);
2488 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;
2489 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);
2490 }
2491}
2492
2493double MVll::Delta_C9_zExp(int hel)
2494{
2495 if (hel == 0)
2496 return beta_0[3].real()*(-26.55265491727846*a_0A12)/a_0A12/a_0A12 +
2497 beta_0[2].real()*(-60.539167428104925*a_0A12)/a_0A12/a_0A12 +
2498 beta_0[1].real()*(-138.02728217972742*a_0A12)/a_0A12/a_0A12 +
2499 beta_0[0].real()*(-314.6975988486678*a_0A12)/a_0A12/a_0A12;
2500 else if (hel == 1)
2501 return (63.24357991272575*a_0A1 - 293.67248647811704*a_0V + 66.1650673421469*a_1A1 - 46.966706577539846*a_1V)*(beta_1[0].real() - beta_2[0].real())/
2502 (1.*a_0A1 - 0.7098414384537659*a_0V)/(1.*a_0A1 - 0.7098414384537659*a_0V) +
2503 (-119.89709952961475*a_0A1 - 24.007514603707598*a_0V + 29.020190982985117*a_1A1 - 20.59973411156516*a_1V)*(beta_1[1].real() - beta_2[1].real())/
2504 (1.*a_0A1 - 0.7098414384537659*a_0V)/(1.*a_0A1 - 0.7098414384537659*a_0V) +
2505 (-117.34075946812884*a_0A1 + 35.43498229234759*a_0V + 12.728340172828181*a_1A1 - 9.035103297409211*a_1V)*(beta_1[2].real() - beta_2[2].real())/
2506 (1.*a_0A1 - 0.7098414384537659*a_0V)/(1.*a_0A1 - 0.7098414384537659*a_0V) +
2507 (-79.86709064819027*a_0A1 + 35.702158475408076*a_0V + 5.582687021261181*a_1A1 - 3.962822585609206*a_1V)*(beta_1[3].real() - beta_2[3].real())/
2508 (1.*a_0A1 - 0.7098414384537659*a_0V)/(1.*a_0A1 - 0.7098414384537659*a_0V);
2509 else
2510 return (63.24357991272575*a_0A1 + 293.67248647811704*a_0V + 66.1650673421469*a_1A1 + 46.966706577539846*a_1V)*(beta_1[0].real() + beta_2[0].real())/
2511 (1.*a_0A1 + 0.7098414384537659*a_0V)/(1.*a_0A1 + 0.7098414384537659*a_0V) +
2512 (-119.89709952961475*a_0A1 + 24.007514603707598*a_0V + 29.020190982985117*a_1A1 + 20.59973411156516*a_1V)*(beta_1[1].real() + beta_2[1].real())/
2513 (1.*a_0A1 + 0.7098414384537659*a_0V)/(1.*a_0A1 + 0.7098414384537659*a_0V) +
2514 (-117.34075946812884*a_0A1 - 35.43498229234759*a_0V + 12.728340172828181*a_1A1 + 9.035103297409211*a_1V)*(beta_1[2].real() + beta_2[2].real())/
2515 (1.*a_0A1 + 0.7098414384537659*a_0V)/(1.*a_0A1 + 0.7098414384537659*a_0V) +
2516 (-79.86709064819027*a_0A1 - 35.702158475408076*a_0V + 5.582687021261181*a_1A1 + 3.962822585609206*a_1V)*(beta_1[3].real() - beta_2[3].real())/
2517 (1.*a_0A1 + 0.7098414384537659*a_0V)/(1.*a_0A1 + 0.7098414384537659*a_0V);
2518}
2519
2520gslpp::complex MVll::H_V_0(double q2, bool bar)
2521{
2522 if (lep == QCD::NEUTRINO_1) {
2523 if (!bar) return -gslpp::complex::i() * NN * (C_L_nunu - etaV * pow(-1, angmomV) * C_R_nunu) * V_0t(q2);
2524 return -gslpp::complex::i() * NN_conjugate * (C_L_nunu.conjugate() - etaV * pow(-1, angmomV) * C_R_nunu.conjugate()) * V_0t(q2);
2525 }
2526 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)));
2527 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)));
2528
2529}
2530
2531gslpp::complex MVll::H_V_p(double q2, bool bar)
2532{
2533 if (lep == QCD::NEUTRINO_1) {
2534 if (!bar) return -gslpp::complex::i() * NN * (C_L_nunu * V_p(q2) - etaV * pow(-1, angmomV) * C_R_nunu * V_m(q2));
2535 return -gslpp::complex::i() * NN_conjugate * (C_L_nunu.conjugate() * V_p(q2) - etaV * pow(-1, angmomV) * C_R_nunu.conjugate() * V_m(q2));
2536 }
2537 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)));
2538 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)));
2539}
2540
2541gslpp::complex MVll::H_V_m(double q2, bool bar)
2542{
2543 if (lep == QCD::NEUTRINO_1) {
2544 if (!bar) return -gslpp::complex::i() * NN * (C_L_nunu * V_m(q2) - etaV * pow(-1, angmomV) * C_R_nunu * V_p(q2));
2545 return -gslpp::complex::i() * NN_conjugate * (C_L_nunu.conjugate() * V_m(q2) - etaV * pow(-1, angmomV) * C_R_nunu.conjugate() * V_p(q2));
2546 }
2547 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)));
2548 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)));
2549}
2550
2551gslpp::complex MVll::H_A_0(double q2, bool bar)
2552{
2553 if (lep == QCD::NEUTRINO_1) {
2554 if (!bar) return -gslpp::complex::i() * NN * (C_L_nunu - etaV * pow(-1, angmomV) * C_R_nunu) * V_0t(q2);
2555 return -gslpp::complex::i() * NN_conjugate * (C_L_nunu.conjugate() - etaV * pow(-1, angmomV) * C_R_nunu.conjugate()) * V_0t(q2);
2556 }
2557 if (!bar) return gslpp::complex::i() * NN * (-C_10 + etaV * pow(-1, angmomV) * C_10p) * V_0t(q2);
2558 return gslpp::complex::i() * NN_conjugate * (-C_10.conjugate() + etaV * pow(-1, angmomV) * C_10p.conjugate()) * V_0t(q2);
2559}
2560
2561gslpp::complex MVll::H_A_p(double q2, bool bar)
2562{
2563 if (lep == QCD::NEUTRINO_1) {
2564 if (!bar) return -gslpp::complex::i() * NN * (C_L_nunu * V_p(q2) - etaV * pow(-1, angmomV) * C_R_nunu * V_m(q2));
2565 return -gslpp::complex::i() * NN_conjugate * (C_L_nunu.conjugate() * V_p(q2) - etaV * pow(-1, angmomV) * C_R_nunu.conjugate() * V_m(q2));
2566 }
2567 if (!bar) return gslpp::complex::i() * NN * (-C_10 * V_p(q2) + etaV * pow(-1, angmomV) * C_10p * V_m(q2));
2568 return gslpp::complex::i() * NN_conjugate * (-C_10.conjugate() * V_p(q2) + etaV * pow(-1, angmomV) * C_10p.conjugate() * V_m(q2));
2569}
2570
2571gslpp::complex MVll::H_A_m(double q2, bool bar)
2572{
2573 if (lep == QCD::NEUTRINO_1) {
2574 if (!bar) return -gslpp::complex::i() * NN * (C_L_nunu * V_m(q2) - etaV * pow(-1, angmomV) * C_R_nunu * V_p(q2));
2575 return -gslpp::complex::i() * NN_conjugate * (C_L_nunu.conjugate() * V_m(q2) - etaV * pow(-1, angmomV) * C_R_nunu.conjugate() * V_p(q2));
2576 }
2577 if (!bar) return gslpp::complex::i() * NN * (-C_10 * V_m(q2) + etaV * pow(-1, angmomV) * C_10p * V_p(q2));
2578 return gslpp::complex::i() * NN_conjugate * (-C_10.conjugate() * V_m(q2) + etaV * pow(-1, angmomV) * C_10p.conjugate() * V_p(q2));
2579}
2580
2581gslpp::complex MVll::H_S(double q2, bool bar)
2582{
2583 if (lep == QCD::NEUTRINO_1) return 0.;
2584
2585 if (!bar) return gslpp::complex::i() * NN * MboMW * (C_S - etaV * pow(-1, angmomV) * C_Sp) * S_L(q2);
2586 return gslpp::complex::i() * NN_conjugate * MboMW * (C_S.conjugate() - etaV * pow(-1, angmomV) * C_Sp.conjugate()) * S_L(q2);
2587}
2588
2589gslpp::complex MVll::H_P(double q2, bool bar)
2590{
2591 if (lep == QCD::NEUTRINO_1) return 0.;
2592
2593 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);
2594 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);
2595}
2596
2597gslpp::complex MVll::AmpMVpsi_zExpansion(double mpsi, int tran)
2598{
2599 updateParameters();
2600
2601 // amplitude at charmonium resonance, i.e. q2 = mJ2 or mPsi2S2
2602 double q2 = mpsi*mpsi;
2603 double fpsi = 0.;
2604 // decay constant of the charmonium state estimated from EXP decay width in e+ e-
2605 if(fabs(mpsi - mJpsi) <1.e-5){
2606 double Gammaepm = 5.971/100.*(92.6*1e-6);
2607 fpsi = sqrt(Gammaepm*(3.*sqrt(q2))/(4.*M_PI*ale*ale)/(4./9.));
2608 }
2609 else if(fabs(mpsi - mPsi2S)< 1.e-5){
2610 double Gammaepm = 7.93/1000.*(294.*1e-6);
2611 fpsi = sqrt(Gammaepm*(3.*sqrt(q2))/(4.*M_PI*ale*ale)/(4./9.));
2612 }
2613 else{
2614 return 0.;
2615 }
2616 gslpp::complex Norm = GF*lambda_t.conjugate()*sqrt(sqrt(lambda(q2))/(2.*M_PI*MM))*MM*MM/sqrt(q2)/fpsi;
2617 if(tran == 0) Norm *= MM/sqrt(q2);
2618 return Norm*DeltaC9_zExpansion(q2,tran);
2619}
2620
2621/*******************************************************************************
2622 * Angular coefficients *
2623 * ****************************************************************************/
2624double MVll::k2(double q2)
2625{
2626 return (MM4 + q2 * q2 + MV4 - twoMV2 * q2 - twoMM2 * (q2 + MV2)) / fourMM2;
2627}
2628
2629double MVll::beta(double q2)
2630{
2631 return sqrt(1. - 4. * Mlep2 / q2);
2632}
2633
2634double MVll::beta2(double q2)
2635{
2636 return 1. - 4. * Mlep2 / q2;
2637}
2638
2639double MVll::lambda(double q2)
2640{
2641 return (MM4 + q2 * q2 + MV4 - twoMV2 * q2 - twoMM2 * (q2 + MV2));
2642}
2643
2644double MVll::F(double q2, double b_i)
2645{
2646 return sqrt(lambda(q2)) * beta(q2) * q2 * b_i / (ninetysixM_PI3MM3);
2647}
2648
2649double MVll::I_1c(double q2, bool bar)
2650{
2651 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()
2652 - H_A_0(q2, bar).abs2()) + beta2(q2) * H_S(q2, bar).abs2());
2653}
2654
2655double MVll::I_1s(double q2, bool bar)
2656{
2657 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()) +
2658 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()));
2659}
2660
2661double MVll::I_2c(double q2, bool bar)
2662{
2663 return -F(q2, b) * beta2(q2) / 2. * (H_V_0(q2, bar).abs2() + H_A_0(q2, bar).abs2());
2664}
2665
2666double MVll::I_2s(double q2, bool bar)
2667{
2668 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());
2669}
2670
2671double MVll::I_3(double q2, bool bar)
2672{
2673 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());
2674}
2675
2676double MVll::I_4(double q2, bool bar)
2677{
2678 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());
2679}
2680
2681double MVll::I_5(double q2, bool bar)
2682{
2683 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()) -
2684 beta(q2) * Mlep / sqrt(q2)*(H_S(q2, bar).conjugate()*(H_V_p(q2, bar) + H_V_m(q2, bar))).real());
2685}
2686
2687double MVll::I_6s(double q2, bool bar)
2688{
2689 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();
2690}
2691
2692double MVll::I_6c(double q2, bool bar)
2693{
2694 return 4. * F(q2, b) * beta(q2) * Mlep / sqrt(q2)*(H_S(q2, bar).conjugate() * H_V_0(q2, bar)).real();
2695}
2696
2697double MVll::I_7(double q2, bool bar)
2698{
2699 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()) -
2700 beta(q2) * Mlep / sqrt(q2)*(H_S(q2, bar).conjugate()*(H_V_m(q2, bar) - H_V_p(q2, bar))).imag());
2701}
2702
2703double MVll::I_8(double q2, bool bar)
2704{
2705 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());
2706}
2707
2708double MVll::I_9(double q2, bool bar)
2709{
2710 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());
2711}
2712
2713double MVll::h_1c(double q2, bool bar)
2714{
2715 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()
2716 - H_A_0(q2, bar).abs2()) - 2. * beta2(q2) * H_S(q2, bar).abs2());
2717}
2718
2719double MVll::h_1s(double q2, bool bar)
2720{
2721 return F(q2, b)*((beta2(q2) + 2.) / 2. * ((H_V_p(q2, bar) * H_V_m(q2, bar).conjugate()).real()
2722 + (H_A_p(q2, bar) * H_A_m(q2, bar).conjugate()).real()) +
2723 4. * Mlep2 / q2 * ((H_V_p(q2, bar) * H_V_m(q2, bar).conjugate()).real()
2724 - (H_A_p(q2, bar) * H_A_m(q2, bar).conjugate()).real()));
2725}
2726
2727double MVll::h_2c(double q2, bool bar)
2728{
2729 return -F(q2, b) * beta2(q2) * (H_V_0(q2, bar).abs2() + H_A_0(q2, bar).abs2());
2730}
2731
2732double MVll::h_2s(double q2, bool bar)
2733{
2734 return F(q2, b) * beta2(q2) / 2. * ((H_V_p(q2, bar) * H_V_m(q2, bar).conjugate()).real()
2735 + (H_A_p(q2, bar) * H_A_m(q2, bar).conjugate()).real());
2736}
2737
2738double MVll::h_3(double q2, bool bar)
2739{
2740 return -F(q2, b) * beta2(q2) / 2. * (H_V_p(q2, bar).abs2() + H_V_m(q2, bar).abs2()
2741 + H_A_p(q2, bar).abs2() + H_A_m(q2, bar).abs2());
2742}
2743
2744double MVll::h_4(double q2, bool bar)
2745{
2746 return F(q2, b) * beta2(q2) / 2. * (((H_V_m(q2, bar) + H_V_p(q2, bar)) * H_V_0(q2, bar).conjugate()).real()
2747 + ((H_A_m(q2, bar) + H_A_p(q2, bar)) * H_A_0(q2, bar).conjugate()).real());
2748}
2749
2750double MVll::h_7(double q2, bool bar)
2751{
2752 return F(q2, b)*(beta(q2) * (((H_V_m(q2, bar) + H_V_p(q2, bar)) * H_A_0(q2, bar).conjugate()).imag()
2753 + ((H_A_m(q2, bar) + H_A_p(q2, bar)) * H_V_0(q2, bar).conjugate()).imag()) -
2754 beta(q2) * 2. * Mlep / sqrt(q2)*(H_S(q2, bar).conjugate()*(H_V_m(q2, bar) - H_V_p(q2, bar))).imag());
2755}
2756
2757double MVll::s_5(double q2, bool bar)
2758{
2759 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)
2760 - F(q2, b)*((H_A_m(q2, bar) - H_A_p(q2, bar)) * H_V_0(q2, bar).conjugate()
2761 + (H_V_m(q2, bar) - H_V_p(q2, bar)) * H_A_0(q2, bar).conjugate()).imag());
2762}
2763
2764double MVll::s_6s(double q2, bool bar)
2765{
2766 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();
2767}
2768
2769double MVll::s_6c(double q2, bool bar)
2770{
2771 return -8. * beta(q2) * Mlep * (H_V_0(q2, bar) * F(q2, b) * H_S(q2, bar).conjugate()).imag() / sqrt(q2);
2772}
2773
2774double MVll::s_8(double q2, bool bar)
2775{
2776 return beta2(q2) * F(q2, b) * ((H_A_m(q2, bar) - H_A_p(q2, bar)) * H_A_0(q2, bar).conjugate()
2777 + (H_V_m(q2, bar) - H_V_p(q2, bar)) * H_V_0(q2, bar).conjugate()).real() / 2.;
2778}
2779
2780double MVll::s_9(double q2, bool bar)
2781{
2782 return beta2(q2) * F(q2, b) * (H_A_p(q2, bar).abs2() - H_A_m(q2, bar).abs2()
2783 + H_V_p(q2, bar).abs2() - H_V_m(q2, bar).abs2()) / 2.;
2784}
2785
2786double MVll::integrateSigma(int i, double q_min, double q_max)
2787{
2788 updateParameters();
2789
2790 std::pair<double, double > qbin = std::make_pair(q_min, q_max);
2791
2792 old_handler = gsl_set_error_handler_off();
2793
2794 switch (i) {
2795 case 0:
2796 if (sigma0Cached[qbin] == 0) {
2797 FS = convertToGslFunction(bind(&MVll::getSigma1c, &(*this), _1));
2798 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();
2799 cacheSigma0[qbin] = avaSigma;
2800 sigma0Cached[qbin] = 1;
2801 }
2802 return cacheSigma0[qbin];
2803 break;
2804 case 1:
2805 if (sigma1Cached[qbin] == 0) {
2806 FS = convertToGslFunction(bind(&MVll::getSigma1s, &(*this), _1));
2807 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();
2808 cacheSigma1[qbin] = avaSigma;
2809 sigma1Cached[qbin] = 1;
2810 }
2811 return cacheSigma1[qbin];
2812 break;
2813 case 2:
2814 if (sigma2Cached[qbin] == 0) {
2815 FS = convertToGslFunction(bind(&MVll::getSigma2c, &(*this), _1));
2816 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();
2817 cacheSigma2[qbin] = avaSigma;
2818 sigma2Cached[qbin] = 1;
2819 }
2820 return cacheSigma2[qbin];
2821 break;
2822 case 3:
2823 if (sigma3Cached[qbin] == 0) {
2824 FS = convertToGslFunction(bind(&MVll::getSigma2s, &(*this), _1));
2825 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();
2826 cacheSigma3[qbin] = avaSigma;
2827 sigma3Cached[qbin] = 1;
2828 }
2829 return cacheSigma3[qbin];
2830 break;
2831 case 4:
2832 if (sigma4Cached[qbin] == 0) {
2833 FS = convertToGslFunction(bind(&MVll::getSigma3, &(*this), _1));
2834 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();
2835 cacheSigma4[qbin] = avaSigma;
2836 sigma4Cached[qbin] = 1;
2837 }
2838 return cacheSigma4[qbin];
2839 break;
2840 case 5:
2841 if (sigma5Cached[qbin] == 0) {
2842 FS = convertToGslFunction(bind(&MVll::getSigma4, &(*this), _1));
2843 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();
2844 cacheSigma5[qbin] = avaSigma;
2845 sigma5Cached[qbin] = 1;
2846 }
2847 return cacheSigma5[qbin];
2848 break;
2849 case 6:
2850 if (sigma6Cached[qbin] == 0) {
2851 FS = convertToGslFunction(bind(&MVll::getSigma5, &(*this), _1));
2852 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();
2853 cacheSigma6[qbin] = avaSigma;
2854 sigma6Cached[qbin] = 1;
2855 }
2856 return cacheSigma6[qbin];
2857 break;
2858 case 7:
2859 if (sigma7Cached[qbin] == 0) {
2860 FS = convertToGslFunction(bind(&MVll::getSigma6s, &(*this), _1));
2861 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();
2862 cacheSigma7[qbin] = avaSigma;
2863 sigma7Cached[qbin] = 1;
2864 }
2865 return cacheSigma7[qbin];
2866 break;
2867 case 8:
2868 if (sigma8Cached[qbin] == 0) {
2869 FS = convertToGslFunction(bind(&MVll::getSigma6c, &(*this), _1));
2870 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();
2871 cacheSigma8[qbin] = avaSigma;
2872 sigma8Cached[qbin] = 1;
2873 }
2874 return cacheSigma8[qbin];
2875 break;
2876 case 9:
2877 if (sigma9Cached[qbin] == 0) {
2878 FS = convertToGslFunction(bind(&MVll::getSigma7, &(*this), _1));
2879 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();
2880 cacheSigma9[qbin] = avaSigma;
2881 sigma9Cached[qbin] = 1;
2882 }
2883 return cacheSigma9[qbin];
2884 break;
2885 case 10:
2886 if (sigma10Cached[qbin] == 0) {
2887 FS = convertToGslFunction(bind(&MVll::getSigma8, &(*this), _1));
2888 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();
2889 cacheSigma10[qbin] = avaSigma;
2890 sigma10Cached[qbin] = 1;
2891 }
2892 return cacheSigma10[qbin];
2893 break;
2894 case 11:
2895 if (sigma11Cached[qbin] == 0) {
2896 FS = convertToGslFunction(bind(&MVll::getSigma9, &(*this), _1));
2897 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();
2898 cacheSigma11[qbin] = avaSigma;
2899 sigma11Cached[qbin] = 1;
2900 }
2901 return cacheSigma11[qbin];
2902 break;
2903 default:
2904 std::stringstream out;
2905 out << i;
2906 throw std::runtime_error("MVll::integrateSigma: index " + out.str() + " not implemented");
2907 }
2908
2909 gsl_set_error_handler(old_handler);
2910
2911}
2912
2913double MVll::getSigma(int i, double q_2)
2914{
2915 updateParameters();
2916
2917 switch (i) {
2918 case 0:
2919 return getSigma1c(q_2);
2920 break;
2921 case 1:
2922 return getSigma1s(q_2);
2923 break;
2924 case 2:
2925 return getSigma2c(q_2);
2926 break;
2927 case 3:
2928 return getSigma2s(q_2);
2929 break;
2930 case 4:
2931 return getSigma3(q_2);
2932 break;
2933 case 5:
2934 return getSigma4(q_2);
2935 break;
2936 case 6:
2937 return getSigma5(q_2);
2938 break;
2939 case 7:
2940 return getSigma6s(q_2);
2941 break;
2942 case 8:
2943 return getSigma6c(q_2);
2944 break;
2945 case 9:
2946 return getSigma7(q_2);
2947 break;
2948 case 10:
2949 return getSigma8(q_2);
2950 break;
2951 case 11:
2952 return getSigma9(q_2);
2953 break;
2954 default:
2955 std::stringstream out;
2956 out << i;
2957 throw std::runtime_error("MVll::getSigma: index " + out.str() + " not implemented");
2958 }
2959}
2960
2961double MVll::integrateDelta(int i, double q_min, double q_max)
2962{
2963 updateParameters();
2964
2965 std::pair<double, double > qbin = std::make_pair(q_min, q_max);
2966
2967 old_handler = gsl_set_error_handler_off();
2968
2969 switch (i) {
2970 case 0:
2971 if (delta0Cached[qbin] == 0) {
2972 FD = convertToGslFunction(bind(&MVll::getDelta1c, &(*this), _1));
2973 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();
2974 cacheDelta0[qbin] = avaDelta;
2975 delta0Cached[qbin] = 1;
2976 }
2977 return cacheDelta0[qbin];
2978 break;
2979 case 1:
2980 if (delta1Cached[qbin] == 0) {
2981 FD = convertToGslFunction(bind(&MVll::getDelta1s, &(*this), _1));
2982 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();
2983 cacheDelta1[qbin] = avaDelta;
2984 delta1Cached[qbin] = 1;
2985 }
2986 return cacheDelta1[qbin];
2987 break;
2988 case 2:
2989 if (delta2Cached[qbin] == 0) {
2990 FD = convertToGslFunction(bind(&MVll::getDelta2c, &(*this), _1));
2991 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();
2992 cacheDelta2[qbin] = avaDelta;
2993 delta2Cached[qbin] = 1;
2994 }
2995 return cacheDelta2[qbin];
2996 break;
2997 case 3:
2998 if (delta3Cached[qbin] == 0) {
2999 FD = convertToGslFunction(bind(&MVll::getDelta2s, &(*this), _1));
3000 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();
3001 cacheDelta3[qbin] = avaDelta;
3002 delta3Cached[qbin] = 1;
3003 }
3004 return cacheDelta3[qbin];
3005 break;
3006 case 6:
3007 if (delta6Cached[qbin] == 0) {
3008 FD = convertToGslFunction(bind(&MVll::getDelta5, &(*this), _1));
3009 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();
3010 cacheDelta6[qbin] = avaDelta;
3011 delta6Cached[qbin] = 1;
3012 }
3013 return cacheDelta6[qbin];
3014 break;
3015 case 7:
3016 if (delta7Cached[qbin] == 0) {
3017 FD = convertToGslFunction(bind(&MVll::getDelta6s, &(*this), _1));
3018 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();
3019 cacheDelta7[qbin] = avaDelta;
3020 delta7Cached[qbin] = 1;
3021 }
3022 return cacheDelta7[qbin];
3023 break;
3024 case 8:
3025 if (delta8Cached[qbin] == 0) {
3026 FD = convertToGslFunction(bind(&MVll::getDelta6c, &(*this), _1));
3027 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();
3028 cacheDelta8[qbin] = avaDelta;
3029 delta8Cached[qbin] = 1;
3030 }
3031 return cacheDelta8[qbin];
3032 break;
3033 case 10:
3034 if (delta10Cached[qbin] == 0) {
3035 FD = convertToGslFunction(bind(&MVll::getDelta8, &(*this), _1));
3036 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();
3037 cacheDelta10[qbin] = avaDelta;
3038 delta10Cached[qbin] = 1;
3039 }
3040 return cacheDelta10[qbin];
3041 break;
3042 case 11:
3043 if (delta11Cached[qbin] == 0) {
3044 FD = convertToGslFunction(bind(&MVll::getDelta9, &(*this), _1));
3045 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();
3046 cacheDelta11[qbin] = avaDelta;
3047 delta11Cached[qbin] = 1;
3048 }
3049 return cacheDelta11[qbin];
3050 break;
3051 default:
3052 std::stringstream out;
3053 out << i;
3054 throw std::runtime_error("MVll::integrateDelta: index " + out.str() + " not implemented");
3055 }
3056
3057 gsl_set_error_handler(old_handler);
3058
3059}
3060double MVll::integrateSigmaTree(double q_min, double q_max)
3061{
3062 if (lep != QCD::NEUTRINO_1 or meson != QCD::B_P or !NeutrinoTree_flag) return 0.;
3063
3064 updateParameters();
3065
3066 //phase space limit where tree-level contribution is relevant (0908.1174)
3067 double q_cut = (mtau2 - MV2) * (MM2 - mtau2) / mtau2;
3068 if (q_max >= q_cut) {
3069 if (q_min == 0.) return getintegratedSigmaTree();
3070 q_max = q_cut;
3071 }
3072
3073 double prefactor = mySM.getMesons(meson).getLifetime() / HCUT * GF4 * VusVub_abs2 * fV2 * fM2 / (64. * M_PI2 * MM3 * Gammatau) * mtau2 * mtau;
3074
3075 std::pair<double, double > qbin = std::make_pair(q_min, q_max);
3076
3077 old_handler = gsl_set_error_handler_off();
3078
3079 if (sigmaTreeCached[qbin] == 0) {
3080 FD = convertToGslFunction(bind(&MVll::SigmaTree, &(*this), _1));
3081 if (gsl_integration_cquad(&FD, q_min, q_max, 1.e-2, 1.e-1, w_sigmaTree, &avaSigmaTree, &errSigmaTree, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
3082 cacheSigmaTree[qbin] = avaSigmaTree;
3083 sigmaTreeCached[qbin] = 1;
3084 }
3085 return prefactor * cacheSigmaTree[qbin];
3086
3087 gsl_set_error_handler(old_handler);
3088}
3089
3090double MVll::SigmaTree(double q2)
3091{
3092 return (MM2 - mtau2) * (mtau2 - MV2) - q2 * (mtau2 - 2. * MV2);
3093}
3094
3095double MVll::getintegratedSigmaTree()
3096{
3097 return mySM.getMesons(meson).getLifetime() / HCUT * GF4 * VusVub_abs2 * fV2 * fM2 / (128. * M_PI2 * MM3 * Gammatau) * mtau * (mtau2 - MV2) * (mtau2 - MV2) * (MM2 - mtau2) * (MM2 - mtau2) * (1. + 2.* MV2 / mtau2);
3098}
F_1.h
F_2.h
LO
@ LO
Definition: OrderScheme.h:34
NLO
@ NLO
Definition: OrderScheme.h:35
NLO_QED11
@ NLO_QED11
Definition: OrderScheme.h:59
CKM::computelamt_s
const gslpp::complex computelamt_s() const
The product of the CKM elements .
Definition: CKM.cpp:174
CKM::computelamu_s
const gslpp::complex computelamu_s() const
The product of the CKM elements .
Definition: CKM.cpp:184
F_1
Definition: F_1.h:15
F_2
Definition: F_2.h:15
Flavour::getFlagMVll_DM
bool getFlagMVll_DM() const
Definition: Flavour.h:379
Flavour::getFlagUseDispersionRelation
bool getFlagUseDispersionRelation() const
Definition: Flavour.h:343
Flavour::ComputeCoeffsnunu
gslpp::vector< gslpp::complex > ** ComputeCoeffsnunu(QCD::lepton lepton=QCD::NOLEPTON, bool noSM=false) const
Definition: Flavour.cpp:170
Flavour::getFlagFixedWCbtos
bool getFlagFixedWCbtos() const
Definition: Flavour.h:363
Flavour::getFlagNeutrinoTree
bool getFlagNeutrinoTree() const
Definition: Flavour.h:383
Flavour::setUpdateFlag
void setUpdateFlag(QCD::meson meson_i, QCD::meson meson_j, QCD::lepton lep_i, bool updated_i) const
sets the update flag for the initial and final state dependent object for .
Definition: Flavour.cpp:315
Flavour::ComputeCoeffprimeBMll
gslpp::vector< gslpp::complex > ** ComputeCoeffprimeBMll(double mu, QCD::lepton lepton, schemes scheme=NDR) const
Computes the chirality flipped Wilson coefficient for the process .
Definition: Flavour.cpp:205
Flavour::getUpdateFlag
bool getUpdateFlag(QCD::meson meson_i, QCD::meson meson_j, QCD::lepton lep_i) const
gets the update flag for the initial and final state dependent object for .
Definition: Flavour.cpp:340
Flavour::ComputeCoeffBMll
gslpp::vector< gslpp::complex > ** ComputeCoeffBMll(double mu, QCD::lepton lepton, bool noSM=false, schemes scheme=NDR) const
Computes the Wilson coefficient for the process .
Definition: Flavour.cpp:195
Flavour::getFlagUsezExpansion
bool getFlagUsezExpansion() const
Definition: Flavour.h:347
MVll::T_para_minus_WA
gslpp::complex T_para_minus_WA(bool conjugate)
QCDF Correction from various BFS paper (hep-ph/0412400). Weak Annihilation.
Definition: MVll.cpp:1837
MVll::deltaC7_QCDF
gslpp::complex deltaC7_QCDF(double q2, bool conjugate, bool spline=false)
QCDF Correction from various BFS papers (hep-ph/0403185, hep-ph/0412400) and Greub et....
Definition: MVll.cpp:1730
MVll::FixedWCbtos
bool FixedWCbtos
Definition: MVll.h:838
MVll::mvllParameters
std::vector< std::string > mvllParameters
Definition: MVll.h:833
MVll::mySM
const StandardModel & mySM
Definition: MVll.h:829
MVll::xs
double xs
Definition: MVll.h:861
MVll::mu_h
double mu_h
Definition: MVll.h:853
MVll::zExpansion
bool zExpansion
Definition: MVll.h:837
MVll::phi_V
double phi_V(double u)
QCDF Correction from various BFS paper (hep-ph/0106067).Vector meson distribution amplitude.
Definition: MVll.cpp:1994
MVll::spline_QCDF_func
void spline_QCDF_func()
Definition: MVll.cpp:2156
MVll::t_para
gslpp::complex t_para(double q2, double u, double m2)
QCDF Correction from various BFS paper (hep-ph/0106067). Part of 4 quark operator contribution.
Definition: MVll.cpp:1876
MVll::B_Seidel
gslpp::complex B_Seidel(double q2, double mb2)
Definition: MVll.cpp:1692
MVll::MVll_DM_flag
bool MVll_DM_flag
Definition: MVll.h:840
MVll::H_A_p
gslpp::complex H_A_p(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:2561
MVll::T_perp_plus_QSS
gslpp::complex T_perp_plus_QSS(double q2, double u, bool conjugate)
QCDF Correction from various BFS paper (hep-ph/0106067). 4 quark operator contribution.
Definition: MVll.cpp:1928
MVll::ale
double ale
Definition: MVll.h:847
MVll::T_para_real
double T_para_real(double q2, double u, bool conjugate)
QCDF Correction from various BFS papers (hep-ph/0106067, hep-ph/0412400). Total.
Definition: MVll.cpp:2031
MVll::T_perp_WA_1
gslpp::complex T_perp_WA_1()
QCDF Correction from various BFS paper (hep-ph/0412400). Weak Annihilation.
Definition: MVll.cpp:1842
MVll::deltaC9_QCDF
gslpp::complex deltaC9_QCDF(double q2, bool conjugate, bool spline=false)
QCDF Correction from various BFS papers (hep-ph/0403185, hep-ph/0412400) and Greub et....
Definition: MVll.cpp:1778
MVll::Mb
double Mb
Definition: MVll.h:851
MVll::myF_2
std::unique_ptr< F_2 > myF_2
Definition: MVll.h:835
MVll::Cq34
gslpp::complex Cq34(bool conjugate)
QCDF Correction from various BFS paper (hep-ph/0412400). Part of Weak Annihilation.
Definition: MVll.cpp:1827
MVll::QCDF_fit_func
double QCDF_fit_func(double *x, double *p)
Definition: MVll.cpp:2087
MVll::mPsi2S2
double mPsi2S2
Definition: MVll.h:842
MVll::MM
double MM
Definition: MVll.h:849
MVll::T_perp_real
double T_perp_real(double q2, double u, bool conjugate)
QCDF Correction from various BFS papers (hep-ph/0106067, hep-ph/0412400). Total.
Definition: MVll.cpp:2005
MVll::T_para_plus_QSS
gslpp::complex T_para_plus_QSS(double q2, double u, bool conjugate)
QCDF Correction from various BFS paper (hep-ph/0106067). 4 quark operator contribution.
Definition: MVll.cpp:1950
MVll::T_para_minus_O8
gslpp::complex T_para_minus_O8(double q2, double u)
QCDF Correction from various BFS paper (hep-ph/0106067). Chromomagnetic dipole contribution contribut...
Definition: MVll.cpp:1860
MVll::C_Seidel
gslpp::complex C_Seidel(double q2)
Definition: MVll.cpp:1724
MVll::H_S
gslpp::complex H_S(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:2581
MVll::mD2
double mD2
Definition: MVll.h:843
MVll::initializeMVllParameters
std::vector< std::string > initializeMVllParameters()
A method for initializing the parameters necessary for MVll.
Definition: MVll.cpp:160
MVll::myF_1
std::unique_ptr< F_1 > myF_1
Definition: MVll.h:834
MVll::integrateDelta
double integrateDelta(int i, double q_min, double q_max)
The integral of from to .
Definition: MVll.cpp:2961
MVll::width
double width
Definition: MVll.h:859
MVll::alpha_s_mub
double alpha_s_mub
Definition: MVll.h:864
MVll::H_V_0
gslpp::complex H_V_0(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:2520
MVll::meson
QCD::meson meson
Definition: MVll.h:831
MVll::T_para_imag
double T_para_imag(double q2, double u, bool conjugate)
QCDF Correction from various BFS papers (hep-ph/0106067, hep-ph/0412400). Total.
Definition: MVll.cpp:2043
MVll::~MVll
virtual ~MVll()
Destructor.
Definition: MVll.cpp:156
MVll::fit_QCDF_func
void fit_QCDF_func()
Definition: MVll.cpp:2092
MVll::T_perp_imag
double T_perp_imag(double q2, double u, bool conjugate)
QCDF Correction from various BFS papers (hep-ph/0106067, hep-ph/0412400). Total.
Definition: MVll.cpp:2018
MVll::dispersion
bool dispersion
Definition: MVll.h:836
MVll::h_func
gslpp::complex h_func(double s, double m2)
Definition: MVll.cpp:1914
MVll::GF
double GF
Definition: MVll.h:846
MVll::T_minus
gslpp::complex T_minus(double q2, bool conjugate)
Definition: MVll.cpp:2225
MVll::getSigma
double getSigma(int i, double q_2)
The value of from to .
Definition: MVll.cpp:2913
MVll::etaV
int etaV
Definition: MVll.h:863
MVll::H_V_p
gslpp::complex H_V_p(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:2531
MVll::lambda_B_minus
gslpp::complex lambda_B_minus(double q2)
Definition: MVll.cpp:1999
MVll::T_0
gslpp::complex T_0(double q2, bool conjugate)
Definition: MVll.cpp:2247
MVll::Ms
double Ms
Definition: MVll.h:857
MVll::mPsi2S
double mPsi2S
Definition: MVll.h:842
MVll::h_lambda
gslpp::complex h_lambda(int hel, double q2)
The non-pertubative ccbar contributions to the helicity amplitudes.
Definition: MVll.cpp:2463
MVll::exp_Phase
gslpp::complex exp_Phase[3]
Definition: MVll.h:844
MVll::mJpsi
double mJpsi
Definition: MVll.h:841
MVll::MV
double MV
Definition: MVll.h:850
MVll::getintegratedSigmaTree
double getintegratedSigmaTree()
The integral of from 0 to .
Definition: MVll.cpp:3095
MVll::integrateSigmaTree
double integrateSigmaTree(double q_min, double q_max)
The integral of from to (arxiv/2301.06990)
Definition: MVll.cpp:3060
MVll::T_para_minus_QSS
gslpp::complex T_para_minus_QSS(double q2, double u, bool conjugate)
QCDF Correction from various BFS paper (hep-ph/0106067). 4 quark operator contribution.
Definition: MVll.cpp:1972
MVll::mc_pole
double mc_pole
Definition: MVll.h:856
MVll::angmomV
double angmomV
Definition: MVll.h:862
MVll::T_perp_WA_2
gslpp::complex T_perp_WA_2(bool conjugate)
QCDF Correction from various BFS paper (hep-ph/0412400). Weak Annihilation.
Definition: MVll.cpp:1847
MVll::FF_fit
double FF_fit(double q2, double a_0, double a_1, double a_2, double MR2)
The fit function from , .
Definition: MVll.cpp:1455
MVll::t_perp
gslpp::complex t_perp(double q2, double u, double m2)
QCDF Correction from various BFS paper (hep-ph/0106067). Part of 4 quark operator contribution.
Definition: MVll.cpp:1867
MVll::MVll
MVll(const StandardModel &SM_i, QCD::meson meson_i, QCD::meson vector_i, QCD::lepton lep_i)
Constructor.
Definition: MVll.cpp:22
MVll::vectorM
QCD::meson vectorM
Definition: MVll.h:832
MVll::H_V_m
gslpp::complex H_V_m(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:2541
MVll::A_Seidel
gslpp::complex A_Seidel(double q2, double mb2)
Definition: MVll.cpp:1678
MVll::AmpMVpsi_zExpansion
gslpp::complex AmpMVpsi_zExpansion(double mpsi, int tran)
Polarization amplitudes for M to V psi, Eq. B.16 of arXiv:2206.03797.
Definition: MVll.cpp:2597
MVll::spectator_charge
double spectator_charge
Definition: MVll.h:858
MVll::Mlep
double Mlep
Definition: MVll.h:848
MVll::B0diff
gslpp::complex B0diff(double q2, double u, double m2)
Definition: MVll.cpp:1900
MVll::Delta_C9_zExp
double Delta_C9_zExp(int hel)
The non-pertubative ccbar contributions to the helicity amplitudes.
Definition: MVll.cpp:2493
MVll::H_A_0
gslpp::complex H_A_0(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:2551
MVll::SigmaTree
double SigmaTree(double q2)
Definition: MVll.cpp:3090
MVll::lep
QCD::lepton lep
Definition: MVll.h:830
MVll::I1
gslpp::complex I1(double q2, double u, double m2)
Definition: MVll.cpp:1883
MVll::B0
gslpp::complex B0(double s, double m2)
Definition: MVll.cpp:1908
MVll::H_A_m
gslpp::complex H_A_m(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:2571
MVll::T_perp_plus_O8
gslpp::complex T_perp_plus_O8(double q2, double u)
QCDF Correction from various BFS paper (hep-ph/0106067). Chromomagnetic dipole contribution contribut...
Definition: MVll.cpp:1852
MVll::H_P
gslpp::complex H_P(double q2, bool bar)
The helicity amplitude .
Definition: MVll.cpp:2589
MVll::Mc
double Mc
Definition: MVll.h:854
MVll::integrateSigma
double integrateSigma(int i, double q_min, double q_max)
The integral of from to .
Definition: MVll.cpp:2786
MVll::NeutrinoTree_flag
bool NeutrinoTree_flag
Definition: MVll.h:839
MVll::mb_pole
double mb_pole
Definition: MVll.h:855
MVll::beta
double beta(double q2)
The factor used in the angular coefficients .
Definition: MVll.cpp:2629
MVll::mu_b
double mu_b
Definition: MVll.h:852
MVll::ys
double ys
Definition: MVll.h:860
MVll::mJ2
double mJ2
Definition: MVll.h:841
Meson::getLambdaM
const double & getLambdaM() const
Definition: Meson.h:402
Meson::getDecayconst_p
const double & getDecayconst_p() const
A get method for the perpendicular decay constant of a vector meson.
Definition: Meson.h:378
Meson::getDgamma_gamma
const double & getDgamma_gamma() const
Definition: Meson.h:411
Meson::computeWidth
double computeWidth() const
A method to compute the width of the meson from its lifetime.
Definition: Meson.cpp:521
Meson::getDecayconst
const double & getDecayconst() const
A get method for the decay constant of the meson.
Definition: Meson.h:360
Meson::getGegenalpha
const double & getGegenalpha(int i) const
A get method to get the Gegenbaur coefficient.
Definition: Meson.h:394
Meson::getLifetime
double getLifetime() const
A get method for the lifetime of the meson.
Definition: Meson.h:351
Particle::getMass
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
Particle::getCharge
double getCharge() const
A get method to access the particle charge.
Definition: Particle.h:97
QCD::meson
meson
An enum type for mesons.
Definition: QCD.h:336
QCD::PHI
@ PHI
Definition: QCD.h:348
QCD::K_star
@ K_star
Definition: QCD.h:349
QCD::B_P
@ B_P
Definition: QCD.h:345
QCD::K_star_P
@ K_star_P
Definition: QCD.h:350
QCD::B_S
@ B_S
Definition: QCD.h:346
QCD::getOptionalParameter
const double getOptionalParameter(std::string name) const
A method to get parameters that are specific to only one set of observables.
Definition: QCD.h:450
QCD::getMesons
const Meson & getMesons(const QCD::meson m) const
A get method to access a meson as an object of the type Meson.
Definition: QCD.h:526
QCD::UP
@ UP
Definition: QCD.h:324
QCD::BOTTOM
@ BOTTOM
Definition: QCD.h:329
QCD::DOWN
@ DOWN
Definition: QCD.h:325
QCD::STRANGE
@ STRANGE
Definition: QCD.h:327
QCD::CHARM
@ CHARM
Definition: QCD.h:326
QCD::lepton
lepton
An enum type for leptons.
Definition: QCD.h:310
QCD::NEUTRINO_2
@ NEUTRINO_2
Definition: QCD.h:313
QCD::NEUTRINO_1
@ NEUTRINO_1
Definition: QCD.h:311
QCD::NEUTRINO_3
@ NEUTRINO_3
Definition: QCD.h:315
QCD::TAU
@ TAU
Definition: QCD.h:316
QCD::getQuarks
const Particle & getQuarks(const QCD::quark q) const
A get method to access a quark as an object of the type Particle.
Definition: QCD.h:536
QCD::initializeMeson
void initializeMeson(QCD::meson meson_i) const
A method to initialize a meson.
Definition: QCD.cpp:280
QCD::getMub
const double getMub() const
A get method to access the threshold between five- and four-flavour theory in GeV.
Definition: QCD.h:582
QCD::Mbar2Mp
const double Mbar2Mp(const double mbar, const quark q, const orders order=FULLNNLO) const
Converts the mass to the pole mass.
Definition: QCD.cpp:1552
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:521
StandardModel::getLeptons
const Particle & getLeptons(const QCD::lepton p) const
A get method to retrieve the member object of a lepton.
Definition: StandardModel.h:756
StandardModel::getCKM
const CKM & getCKM() const
A get method to retrieve the member object of type CKM.
Definition: StandardModel.h:952
StandardModel::Mw
virtual const double Mw() const
The SM prediction for the -boson mass in the on-shell scheme, .
Definition: StandardModel/src/StandardModel.cpp:1025
StandardModel::getFlavour
const Flavour & getFlavour() const
Definition: StandardModel.h:1065
StandardModel::getGF
const double getGF() const
A get method to retrieve the Fermi constant .
Definition: StandardModel.h:792
StandardModel::Als
const double Als(const double mu, const orders order, const bool Nf_thr, const bool qed_flag) const
The running QCD coupling in the scheme including QED corrections.
Definition: StandardModel.h:1120
StandardModel::getAle
const double getAle() const
A get method to retrieve the fine-structure constant .
Definition: StandardModel.h:801
h_0
A class for the correction in .
Definition: MVllObservables.h:1952
s
Test Observable.
unitarity_bound_F2
A class for the unitarity constraints on form factors .
Definition: MVllObservables.h:1733
unitarity_bound_T1
A class for the unitarity constraints on form factors .
Definition: MVllObservables.h:1767
unitarity_bound_T2_T0
A class for the unitarity constraints on form factors and .
Definition: MVllObservables.h:1801
unitarity_bound_f_F1
A class for the unitarity constraints on form factors and .
Definition: MVllObservables.h:1665
unitarity_bound_g
A class for the unitarity constraints on form factors .
Definition: MVllObservables.h:1699