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