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