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