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