a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
MPlnu.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 "MPlnu.h"
10#include "std_make_vector.h"
11#include "gslpp_function_adapter.h"
12#include <gsl/gsl_sf_zeta.h>
13#include <boost/bind/bind.hpp>
14#include <limits>
15#include <TFitResult.h>
16#include <gsl/gsl_sf_gegenbauer.h>
17#include <gsl/gsl_sf_expint.h>
18#include <limits>
19using namespace boost::placeholders;
20
21MPlnu::MPlnu(const StandardModel& SM_i, QCD::meson meson_i, QCD::meson pseudoscalar_i, QCD::lepton lep_i)
22: mySM(SM_i), ig(ROOT::Math::Integration::kADAPTIVESINGULAR, ROOT::Math::Integration::kGAUSS51), wf()
23{
24 lep = lep_i;
25 meson = meson_i;
26 pseudoscalarM = pseudoscalar_i;
27 CLNflag = false;
28 BGLflag = false;
29 DMflag = false;
30 btocNPpmflag = false;
31 NPanalysis = false;
32
33 checkcache_int_tau = false;
34 checkcache_int_mu = false;
35 checkcache_int_el = false;
36
37 double max_double = std::numeric_limits<double>::max();
38
39 fplusz0_cache = max_double;
40 rho1to2_cache = max_double;
41 N_0_cache = max_double;
42 alpha_0_cache = max_double;
43 alpha_p_cache = max_double;
44 beta_0_cache = max_double;
45 beta_p_cache = max_double;
46 gamma_0_cache = max_double;
47 gamma_p_cache = max_double;
48 af0_1_cache = max_double;
49 af0_2_cache = max_double;
50 afplus_0_cache = max_double;
51 afplus_1_cache = max_double;
52 afplus_2_cache = max_double;
53
54#if NBGL == 3
55 af0_3_cache = max_double;
56 afplus_3_cache = max_double;
57#endif
58
59 afplus1_cache = max_double;
60 afplus2_cache = max_double;
61 afplus3_cache = max_double;
62 afplus4_cache = max_double;
63 afplus5_cache = max_double;
64 afplus6_cache = max_double;
65 afplus7_cache = max_double;
66 afplus8_cache = max_double;
67 afplus9_cache = max_double;
68 afplus10_cache = max_double;
69 afzero1_cache = max_double;
70 afzero2_cache = max_double;
71 afzero3_cache = max_double;
72 afzero4_cache = max_double;
73 afzero5_cache = max_double;
74 afzero6_cache = max_double;
75 afzero7_cache = max_double;
76 afzero8_cache = max_double;
77 afzero9_cache = max_double;
78 afzero10_cache = max_double;
79 bfplus1_cache = max_double;
80 bfplus2_cache = max_double;
81 bfplus3_cache = max_double;
82 bfplus4_cache = max_double;
83 bfplus5_cache = max_double;
84 bfplus6_cache = max_double;
85 bfplus7_cache = max_double;
86 bfplus8_cache = max_double;
87 bfplus9_cache = max_double;
88 bfplus10_cache = max_double;
89 bfzero1_cache = max_double;
90 bfzero2_cache = max_double;
91 bfzero3_cache = max_double;
92 bfzero4_cache = max_double;
93 bfzero5_cache = max_double;
94 bfzero6_cache = max_double;
95 bfzero7_cache = max_double;
96 bfzero8_cache = max_double;
97 bfzero9_cache = max_double;
98 bfzero10_cache = max_double;
99
100 CS_cache = max_double;
101 CSp_cache = max_double;
102 CP_cache = max_double;
103 CPp_cache = max_double;
104 CV_cache = max_double;
105 CVp_cache = max_double;
106 CA_cache = max_double;
107 CAp_cache = max_double;
108 CT_cache = max_double;
109 CTp_cache = max_double;
110
111 ig.SetRelTolerance(1.E-6); // set relative tolerance
112 ig.SetAbsTolerance(1.E-6); // set absoulte tolerance
113
114}
115
117{}
118
119std::vector<std::string> MPlnu::initializeMPlnuParameters()
120{
124 btocNPpmflag = (mySM.getModelName().compare("RealWeakEFTCCPM") == 0);
125 NPanalysis = (mySM.getModelName().compare("RealWeakEFTCCPM") == 0 || mySM.getModelName().compare("RealWeakEFTCC") == 0);
126 if (CLNflag + BGLflag + DMflag != true) throw std::runtime_error("MPlnu: Set only one among CLNflag, BGLflag, DMflag to true");
127 else mplnuParameters = make_vector<std::string>();
128 if (CLNflag) {
129 mplnuParameters.clear();
130 if (pseudoscalarM == StandardModel::D_P) mplnuParameters = make_vector<std::string>()
131 << "fplusz0" << "rho1to2"
132 << "N_0" << "alpha_0" << "alpha_p" << "beta_0" << "beta_p" << "gamma_0" << "gamma_p";
133 }
134 else if (BGLflag) {
135 mplnuParameters.clear();
136 if (pseudoscalarM == StandardModel::D_P) mplnuParameters = make_vector<std::string>()
137 << "af0_1" << "af0_2" << "afplus_0" << "afplus_1" << "afplus_2"
138#if NBGL == 3
139 << "af0_3" << "afplus_3"
140#endif
141 << "mBc1m_1" << "mBc1m_2" << "mBc1m_3" << "mBc1m_4"
142 << "mBc0p_1" << "mBc0p_2" << "chitildeT" << "chiL" << "n_I";
143 }
144 else if (DMflag) {
145 mplnuParameters.clear();
146 mplnuParameters = make_vector<std::string>()
147 << "afplus1" << "afplus2" << "afplus3" << "afplus4" << "afplus5"
148 << "afplus6" << "afplus7" << "afplus8" << "afplus9" << "afplus10"
149 << "afzero1" << "afzero2" << "afzero3" << "afzero4" << "afzero5"
150 << "afzero6" << "afzero7" << "afzero8" << "afzero9" << "afzero10"
151 << "bfplus1" << "bfplus2" << "bfplus3" << "bfplus4" << "bfplus5"
152 << "bfplus6" << "bfplus7" << "bfplus8" << "bfplus9" << "bfplus10"
153 << "bfzero1" << "bfzero2" << "bfzero3" << "bfzero4" << "bfzero5"
154 << "bfzero6" << "bfzero7" << "bfzero8" << "bfzero9" << "bfzero10";
155 }
156 else {
157 std::stringstream out;
158 out << pseudoscalarM;
159 throw std::runtime_error("MPlnu: vector " + out.str() + " not implemented");
160 }
161
164 return mplnuParameters;
165}
166
167void MPlnu::updateParameters()
168{
170
172 Mnu = 0.; // neutrinos assumed to be massless
176 w0 = (MM*MM+MP*MP)/(2.*MM*MP);
177 RV = 2.*sqrt(MM*MP)/(MM+MP);
178 mu_b = MM; // mySM.getMub();
179 Mb = mySM.getQuarks(QCD::BOTTOM).getMass(); // add the PS b mass
180 Mc = mySM.getQuarks(QCD::CHARM).getMass(); // add the PS b mass
181 ale_mub = mySM.Ale(mu_b,FULLNLO);
182 /* Amplitude propto 4*GF*Vij/sqrt(2) & kinematics requires 1/(2^9 pi^3 MB^3) */
183 amplsq_factor = 1./(64.*M_PI*M_PI*M_PI*MM*MM*MM);
184 q2min = Mlep*Mlep;
185 q2max = (MM-MP)*(MM-MP);
186
187 /* SM + NP Wilson coefficients */
188 gslpp::complex norm = 4./sqrt(2.);
189 gslpp::vector<gslpp::complex> ** allcoeff_bclnu = mySM.getFlavour().ComputeCoeffdiujlknu(2,1,0,mu_b);
190 CV = (*(allcoeff_bclnu[LO]))(0)/norm*(1.+ale_mub/M_PI*log(mySM.getMz()/mu_b))/2.;
191 CA = -CV;
192 CVp = (*(allcoeff_bclnu[LO]))(1)/norm/2.;
193 CAp = -CVp;
194 CS = (*(allcoeff_bclnu[LO]))(2)/norm/2.;
195 CSp = (*(allcoeff_bclnu[LO]))(3)/norm/2.;
196 CP = -CS;
197 CPp = -CSp;
198 C7 = 0.;
199 C7p = 0.;
200 CT = (*(allcoeff_bclnu[LO]))(4)/norm/2.;
201 CTp = 0.;
202
203 if (NPanalysis) {
204 if (lep == StandardModel::TAU) {
205 if (btocNPpmflag) {
206 CV += (mySM.getCCC3() - mySM.getCCC4()) / 2. / M_SQRT2;
207 CVp = (mySM.getCCC3() + mySM.getCCC4()) / 2. / M_SQRT2;
208 CA -= (mySM.getCCC3() - mySM.getCCC4()) / 2. / M_SQRT2;
209 CAp = -(mySM.getCCC3() + mySM.getCCC4()) / 2. / M_SQRT2;
210 CS = (mySM.getCCC1() - mySM.getCCC2()) / 2. / M_SQRT2;
211 CSp = (mySM.getCCC1() + mySM.getCCC2()) / 2. / M_SQRT2;
212 CP = -(mySM.getCCC1() - mySM.getCCC2()) / 2. / M_SQRT2;
213 CPp = -(mySM.getCCC1() + mySM.getCCC2()) / 2. / M_SQRT2;
214 CTp = mySM.getCCC5();
215 } else {
216 CV += mySM.getCCC3() / 2.;
217 CVp = mySM.getCCC4() / 2.;
218 CA -= mySM.getCCC3() / 2.;
219 CAp = -mySM.getCCC4() / 2.;
220 CS = mySM.getCCC1() / 2.;
221 CSp = mySM.getCCC2() / 2.;
222 CP = -mySM.getCCC1() / 2.;
223 CPp = -mySM.getCCC2() / 2.;
224 CTp = mySM.getCCC5();
225 }
226 }
227 }
228
229 switch (pseudoscalarM) {
231 if (CLNflag) {
232 fplusz0 = mySM.getOptionalParameter("fplusz0");
233 rho1to2 = mySM.getOptionalParameter("rho1to2");
234 N_0 = mySM.getOptionalParameter("N_0");
235 alpha_0 = mySM.getOptionalParameter("alpha_0");
236 alpha_p = mySM.getOptionalParameter("alpha_p");
237 beta_0 = mySM.getOptionalParameter("beta_0");
238 beta_p = mySM.getOptionalParameter("beta_p");
239 gamma_0 = mySM.getOptionalParameter("gamma_0");
240 gamma_p = mySM.getOptionalParameter("gamma_p");
241 }
242 else if (BGLflag) {
243 af0_1 = mySM.getOptionalParameter("af0_1");
244 af0_2 = mySM.getOptionalParameter("af0_2");
245 afplus_0 = mySM.getOptionalParameter("afplus_0");
246 afplus_1 = mySM.getOptionalParameter("afplus_1");
247 afplus_2 = mySM.getOptionalParameter("afplus_2");
248#if NBGL == 3
249 af0_3 = mySM.getOptionalParameter("af0_3");
250 afplus_3 = mySM.getOptionalParameter("afplus_3");
251#endif
252 mBc1m_1 = mySM.getOptionalParameter("mBc1m_1");
253 mBc1m_2 = mySM.getOptionalParameter("mBc1m_2");
254 mBc1m_3 = mySM.getOptionalParameter("mBc1m_3");
255 mBc1m_4 = mySM.getOptionalParameter("mBc1m_4");
256 mBc0p_1 = mySM.getOptionalParameter("mBc0p_1");
257 mBc0p_2 = mySM.getOptionalParameter("mBc0p_2");
258 chitildeT = mySM.getOptionalParameter("chitildeT");
259 chiL = mySM.getOptionalParameter("chiL");
260 nI = mySM.getOptionalParameter("n_I");
261
262 z1m_1 = sqrt((MM+MP)*(MM+MP)-mBc1m_1*mBc1m_1)-sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP));
263 z1m_1 /= (sqrt((MM+MP)*(MM+MP)-mBc1m_1*mBc1m_1)+sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP)));
264 z1m_2 = sqrt((MM+MP)*(MM+MP)-mBc1m_2*mBc1m_2)-sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP));
265 z1m_2 /= (sqrt((MM+MP)*(MM+MP)-mBc1m_2*mBc1m_2)+sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP)));
266 z1m_3 = sqrt((MM+MP)*(MM+MP)-mBc1m_3*mBc1m_3)-sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP));
267 z1m_3 /= (sqrt((MM+MP)*(MM+MP)-mBc1m_3*mBc1m_3)+sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP)));
268
269 z0p_1 = sqrt((MM+MP)*(MM+MP)-mBc0p_1*mBc0p_1)-sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP));
270 z0p_1 /= (sqrt((MM+MP)*(MM+MP)-mBc0p_1*mBc0p_1)+sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP)));
271 z0p_2 = sqrt((MM+MP)*(MM+MP)-mBc0p_2*mBc0p_2)-sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP));
272 z0p_2 /= (sqrt((MM+MP)*(MM+MP)-mBc0p_2*mBc0p_2)+sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP)));
273 }
274 else if (DMflag) {
275 afplus1 = mySM.getOptionalParameter("afplus1");
276 afplus2 = mySM.getOptionalParameter("afplus2");
277 afplus3 = mySM.getOptionalParameter("afplus3");
278 afplus4 = mySM.getOptionalParameter("afplus4");
279 afplus5 = mySM.getOptionalParameter("afplus5");
280 afplus6 = mySM.getOptionalParameter("afplus6");
281 afplus7 = mySM.getOptionalParameter("afplus7");
282 afplus8 = mySM.getOptionalParameter("afplus8");
283 afplus9 = mySM.getOptionalParameter("afplus9");
284 afplus10 = mySM.getOptionalParameter("afplus10");
285 afzero1 = mySM.getOptionalParameter("afzero1");
286 afzero2 = mySM.getOptionalParameter("afzero2");
287 afzero3 = mySM.getOptionalParameter("afzero3");
288 afzero4 = mySM.getOptionalParameter("afzero4");
289 afzero5 = mySM.getOptionalParameter("afzero5");
290 afzero6 = mySM.getOptionalParameter("afzero6");
291 afzero7 = mySM.getOptionalParameter("afzero7");
292 afzero8 = mySM.getOptionalParameter("afzero8");
293 afzero9 = mySM.getOptionalParameter("afzero9");
294 afzero10 = mySM.getOptionalParameter("afzero10");
295 bfplus1 = mySM.getOptionalParameter("bfplus1");
296 bfplus2 = mySM.getOptionalParameter("bfplus2");
297 bfplus3 = mySM.getOptionalParameter("bfplus3");
298 bfplus4 = mySM.getOptionalParameter("bfplus4");
299 bfplus5 = mySM.getOptionalParameter("bfplus5");
300 bfplus6 = mySM.getOptionalParameter("bfplus6");
301 bfplus7 = mySM.getOptionalParameter("bfplus7");
302 bfplus8 = mySM.getOptionalParameter("bfplus8");
303 bfplus9 = mySM.getOptionalParameter("bfplus9");
304 bfplus10 = mySM.getOptionalParameter("bfplus10");
305 bfzero1 = mySM.getOptionalParameter("bfzero1");
306 bfzero2 = mySM.getOptionalParameter("bfzero2");
307 bfzero3 = mySM.getOptionalParameter("bfzero3");
308 bfzero4 = mySM.getOptionalParameter("bfzero4");
309 bfzero5 = mySM.getOptionalParameter("bfzero5");
310 bfzero6 = mySM.getOptionalParameter("bfzero6");
311 bfzero7 = mySM.getOptionalParameter("bfzero7");
312 bfzero8 = mySM.getOptionalParameter("bfzero8");
313 bfzero9 = mySM.getOptionalParameter("bfzero9");
314 bfzero10 = mySM.getOptionalParameter("bfzero10");
315 }
316 else{ };
317 break;
318 default:
319 std::stringstream out;
320 out << pseudoscalarM;
321 throw std::runtime_error("MPlnu: vector " + out.str() + " not implemented");
322 }
323
324 if ((fplusz0 != fplusz0_cache) || (rho1to2 != rho1to2_cache)
325 || (N_0 != N_0_cache)
326 || (alpha_0 != alpha_0_cache) || (alpha_p != alpha_p_cache)
327 || (beta_0 != beta_0_cache) || (beta_p != beta_p_cache)
328 || (gamma_0 != gamma_0_cache) || (gamma_p != gamma_p_cache)
329 || (af0_1 != af0_1_cache) || (af0_2 != af0_2_cache)
330 || (afplus_0 != afplus_0_cache) || (afplus_1 != afplus_1_cache)
331 || (afplus_2 != afplus_2_cache)
332#if NBGL == 3
333 || (af0_3 != af0_3_cache) || (afplus_3 != afplus_3_cache)
334#endif
335 || (afplus1 != afplus1_cache) || (afplus2 != afplus2_cache)
336 || (afplus3 != afplus3_cache) || (afplus4 != afplus4_cache)
337 || (afplus5 != afplus5_cache) || (afplus6 != afplus6_cache)
338 || (afplus7 != afplus7_cache) || (afplus8 != afplus8_cache)
339 || (afplus9 != afplus9_cache) || (afplus10 != afplus10_cache)
340 || (afzero1 != afzero1_cache) || (afzero2 != afzero2_cache)
341 || (afzero3 != afzero3_cache) || (afzero4 != afzero4_cache)
342 || (afzero5 != afzero5_cache) || (afzero6 != afzero6_cache)
343 || (afzero7 != afzero7_cache) || (afzero8 != afzero8_cache)
344 || (afzero9 != afzero9_cache) || (afzero10 != afzero10_cache)
345 || (bfplus1 != bfplus1_cache) || (bfplus2 != bfplus2_cache)
346 || (bfplus3 != bfplus3_cache) || (bfplus4 != bfplus4_cache)
347 || (bfplus5 != bfplus5_cache) || (bfplus6 != bfplus6_cache)
348 || (bfplus7 != bfplus7_cache) || (bfplus8 != bfplus8_cache)
349 || (bfplus9 != bfplus9_cache) || (bfplus10 != bfplus10_cache)
350 || (bfzero1 != bfzero1_cache) || (bfzero2 != bfzero2_cache)
351 || (bfzero3 != bfzero3_cache) || (bfzero4 != bfzero4_cache)
352 || (bfzero5 != bfzero5_cache) || (bfzero6 != bfzero6_cache)
353 || (bfzero7 != bfzero7_cache) || (bfzero8 != bfzero8_cache)
354 || (bfzero9 != bfzero9_cache) || (bfzero10 != bfzero10_cache)
355 || (CS != CS_cache) || (CSp != CSp_cache)
356 || (CP != CP_cache) || (CPp != CPp_cache)
357 || (CV != CV_cache) || (CVp != CVp_cache)
358 || (CA != CA_cache) || (CAp != CAp_cache)
359 || (CT != CT_cache) || (CTp != CTp_cache)) {
360 checkcache_int_tau = false;
361 checkcache_int_mu = false;
362 checkcache_int_el = false;
363 }
364
365 if (!checkcache_int_tau || !checkcache_int_mu || !checkcache_int_el) {
366 if (lep == StandardModel::TAU) {
367 cached_intJ1_tau = integrateJ(1, q2min, q2max);
368 cached_intJ3_tau = integrateJ(3, q2min, q2max);
369 cached_intJ2_tau = 0.;
370 // cached_intJ2_tau = integrateJ(2,q2min,q2max);
371 checkcache_int_tau = true;
372 }
373 if (lep == StandardModel::MU) {
374 cached_intJ1_mu = integrateJ(1, q2min, q2max);
375 cached_intJ3_mu = integrateJ(3, q2min, q2max);
376 cached_intJ2_mu = 0.;
377 // cached_intJ2_mu = integrateJ(2,q2min,q2max);
378 checkcache_int_mu = true;
379 }
381 cached_intJ1_el = integrateJ(1, q2min, q2max);
382 cached_intJ3_el = integrateJ(3, q2min, q2max);
383 cached_intJ2_el = 0.;
384 // cached_intJ2_el = integrateJ(2,q2min,q2max);
385 checkcache_int_el = true;
386 }
387 }
388 if (CLNflag) {
389 fplusz0_cache = fplusz0;
390 rho1to2_cache = rho1to2;
391 N_0_cache = N_0;
392 alpha_0_cache = alpha_0;
393 alpha_p_cache = alpha_p;
394 beta_0_cache = beta_0;
395 beta_p_cache = beta_p;
396 gamma_0_cache = gamma_0;
397 gamma_p_cache = gamma_p;
398 af0_0 = 0.;
399 }
400 else if (BGLflag) {
401 af0_1_cache = af0_1;
402 af0_2_cache = af0_2;
403 afplus_0_cache = afplus_0;
404 afplus_1_cache = afplus_1;
405 afplus_2_cache = afplus_2;
406#if NBGL == 3
407 af0_3_cache = af0_3;
408 afplus_3_cache = afplus_3;
409#endif
410 /* f+(q2=0) = f0(q2=0) */
411 z0 = (sqrt(w0+1.)-sqrt(2.))/(sqrt(w0+1.)+sqrt(2.));
412 af0_0 = fplus(0.);
413#if NBGL == 3
414 af0_0 -= (af0_1 * z0 + af0_2 * z0 * z0 + af0_3 * z0 * z0 *z0) / phi_f0(z0) / ((z0 - z0p_1) / (1. - z0 * z0p_1)*(z0 - z0p_2) / (1. - z0 * z0p_2));
415#else
416 af0_0 -= (af0_1 * z0 + af0_2 * z0 * z0) / phi_f0(z0) / ((z0 - z0p_1) / (1. - z0 * z0p_1)*(z0 - z0p_2) / (1. - z0 * z0p_2));
417#endif
418 af0_0 *= phi_f0(z0)*((z0 - z0p_1) / (1. - z0 * z0p_1)*(z0 - z0p_2) / (1. - z0 * z0p_2));
419 }
420 else if (DMflag) {
421 afplus1_cache = afplus1;
422 afplus2_cache = afplus2;
423 afplus3_cache = afplus3;
424 afplus4_cache = afplus4;
425 afplus5_cache = afplus5;
426 afplus6_cache = afplus6;
427 afplus7_cache = afplus7;
428 afplus8_cache = afplus8;
429 afplus9_cache = afplus9;
430 afplus10_cache = afplus10;
431 afzero1_cache = afzero1;
432 afzero2_cache = afzero2;
433 afzero3_cache = afzero3;
434 afzero4_cache = afzero4;
435 afzero5_cache = afzero5;
436 afzero6_cache = afzero6;
437 afzero7_cache = afzero7;
438 afzero8_cache = afzero8;
439 afzero9_cache = afzero9;
440 afzero10_cache = afzero10;
441 bfplus1_cache = bfplus1;
442 bfplus2_cache = bfplus2;
443 bfplus3_cache = bfplus3;
444 bfplus4_cache = bfplus4;
445 bfplus5_cache = bfplus5;
446 bfplus6_cache = bfplus6;
447 bfplus7_cache = bfplus7;
448 bfplus8_cache = bfplus8;
449 bfplus9_cache = bfplus9;
450 bfplus10_cache = bfplus10;
451 bfzero1_cache = bfzero1;
452 bfzero2_cache = bfzero2;
453 bfzero3_cache = bfzero3;
454 bfzero4_cache = bfzero4;
455 bfzero5_cache = bfzero5;
456 bfzero6_cache = bfzero6;
457 bfzero7_cache = bfzero7;
458 bfzero8_cache = bfzero8;
459 bfzero9_cache = bfzero9;
460 bfzero10_cache = bfzero10;
461 }
462 else{ };
463 CS_cache = CS;
464 CSp_cache = CSp;
465 CP_cache = CP;
466 CPp_cache = CPp;
467 CV_cache = CV;
468 CVp_cache = CVp;
469 CA_cache = CA;
470 CAp_cache = CAp;
471 CT_cache = CT;
472 CTp_cache = CTp;
473
475
476 return;
477
478}
479
480/*******************************************************************************
481 * Kinematic functions *
482 * ****************************************************************************/
483
484double MPlnu::lambda_half(double a, double b, double c)
485{
486 return sqrt(a*a+b*b+c*c-2.*(a*b+a*c+b*c));
487}
488/*******************************************************************************
489 * Form factors *
490 * ****************************************************************************/
491
492double MPlnu::phi_fplus(double z)
493{
494 // chitildeT in GeV-2
495 double prefac = 8. * (MP / MM)*(MP / MM) / MM * sqrt(8. * nI / (3. * M_PI * chitildeT));
496 double num = (1. + z)*(1. + z) * sqrt(1. - z);
497 double den = (1. + MP / MM)*(1. - z) + 2. * sqrt(MP / MM)*(1. + z);
498 double den5 = den * den * den * den*den;
499 return prefac * num / den5;
500}
501
502double MPlnu::phi_f0(double z)
503{
504 // chiL dimensionless
505 double prefac = (MP / MM)*(1. - (MP / MM)*(MP / MM)) * sqrt(8. * nI / (M_PI * chiL));
506 double num = (1. - z * z) * sqrt((1. - z));
507 double den = (1. + MP / MM)*(1. - z) + 2. * sqrt(MP / MM)*(1. + z);
508 double den4 = den * den * den*den;
509 return prefac * num / den4;
510}
511
512double MPlnu::fplus(double q2)
513{
514 double w = w0 - q2 / (2. * MM * MP);
515 double z = (sqrt(w + 1.) - M_SQRT2) / (sqrt(w + 1.) + M_SQRT2);
516 if (CLNflag) {
517 return fplusz0 * N_0 * (1. - alpha_p*8. * rho1to2 * z + beta_p*(51. * rho1to2 - 10.) * z * z - gamma_p*(252. * rho1to2 - 84.) * z * z * z);
518 }
519 else if (BGLflag) {
520 double P_fplus = (z - z1m_1) / (1. - z * z1m_1)*(z - z1m_2) / (1. - z * z1m_2)*(z - z1m_3) / (1. - z * z1m_3);
521#if NBGL == 3
522 return (afplus_0 + afplus_1 * z + afplus_2 * z * z + afplus_3 * z * z * z) / phi_fplus(z) / P_fplus;
523#else
524 return (afplus_0 + afplus_1 * z + afplus_2 * z * z) / phi_fplus(z) / P_fplus;
525#endif
526 }
527 else if (DMflag) {
528 double fp_formfac = 0.;
529 if (w<1.06) {fp_formfac = afplus1 + bfplus1 * w;}
530 else if (w<1.12) {fp_formfac = afplus2 + bfplus2 * w;}
531 else if (w<1.18) {fp_formfac = afplus3 + bfplus3 * w;}
532 else if (w<1.24) {fp_formfac = afplus4 + bfplus4 * w;}
533 else if (w<1.30) {fp_formfac = afplus5 + bfplus5 * w;}
534 else if (w<1.36) {fp_formfac = afplus6 + bfplus6 * w;}
535 else if (w<1.42) {fp_formfac = afplus7 + bfplus7 * w;}
536 else if (w<1.48) {fp_formfac = afplus8 + bfplus8 * w;}
537 else if (w<1.54) {fp_formfac = afplus9 + bfplus9 * w;}
538 else {fp_formfac = afplus10 + bfplus10 * w;}
539 return fp_formfac;
540 }
541 else return 0.;
542}
543
544double MPlnu::f0(double q2)
545{
546 double w = w0 - q2 / (2. * MM * MP);
547 double z = (sqrt(w + 1.) - M_SQRT2) / (sqrt(w + 1.) + M_SQRT2);
548 if (CLNflag) {
549 double prefac0 = 2. * (MP / MM) / (1. + MP / MM)/ (1. + MP / MM)*(1. + w0);
550 double norm = prefac0 * (1. - alpha_0*0.0068 * (w0 - 1.) + beta_0*0.0017 * (w0 - 1.)*(w0 - 1.) - gamma_0*0.0013 * (w0 - 1.)*(w0 - 1.)*(w0 - 1.));
551 double prefac = fplus(q2)* prefac0/(1. + w0)*(1. + w);
552 // norm introduced to respect f+(q2=0)=f0(q2=0) exactly
553 return prefac/norm * (1. - alpha_0*0.0068 * (w - 1.) + beta_0*0.0017 * (w - 1.)*(w - 1.) - gamma_0*0.0013 * (w - 1.)*(w - 1.)*(w - 1.));
554 }
555 else if (BGLflag) {
556 double P_f0 = (z - z0p_1) / (1. - z * z0p_1)*(z - z0p_2) / (1. - z * z0p_2);
557#if NBGL == 3
558 return (af0_0 + af0_1 * z + af0_2 * z * z + af0_3 * z * z * z) / phi_f0(z) / P_f0;
559#else
560 return (af0_0 + af0_1 * z + af0_2 * z * z) / phi_f0(z) / P_f0;
561#endif
562 }
563 else if (DMflag) {
564 double f0_formfac = 0.;
565 if (w<1.06) {f0_formfac = afzero1 + bfzero1 * w;}
566 else if (w<1.12) {f0_formfac = afzero2 + bfzero2 * w;}
567 else if (w<1.18) {f0_formfac = afzero3 + bfzero3 * w;}
568 else if (w<1.24) {f0_formfac = afzero4 + bfzero4 * w;}
569 else if (w<1.30) {f0_formfac = afzero5 + bfzero5 * w;}
570 else if (w<1.36) {f0_formfac = afzero6 + bfzero6 * w;}
571 else if (w<1.42) {f0_formfac = afzero7 + bfzero7 * w;}
572 else if (w<1.48) {f0_formfac = afzero8 + bfzero8 * w;}
573 else if (w<1.54) {f0_formfac = afzero9 + bfzero9 * w;}
574 else {f0_formfac = afzero10 + bfzero10 * w;}
575 return f0_formfac;
576 }
577 else return 0.;
578}
579
580double MPlnu::fT(double q2)
581{
582 return 0.;
583}
584/********************************************************************************
585 * Helicity amplitudes (normalization such that all H \propto (mass scale)^-1) *
586 * *****************************************************************************/
587
588gslpp::complex MPlnu::HV(double q2)
589{
590 return lambda_half(MM*MM, MP*MP, q2) / 2. / sqrt(q2)*((CV + CVp) * fplus(q2) + 2. * Mb / (MM + MP)*(C7 + C7p) * fT(q2));
591}
592
593gslpp::complex MPlnu::HA(double q2)
594{
595 return lambda_half(MM*MM, MP*MP, q2) / 2. / sqrt(q2)*(CA + CAp) * fplus(q2);
596}
597
598gslpp::complex MPlnu::HP(double q2)
599{
600 return (MM * MM - MP * MP) / 2. * ((CP + CPp) / (Mb - Mc)+(Mlep + Mnu) / q2 * (CA + CAp)) * f0(q2);
601}
602
603gslpp::complex MPlnu::HS(double q2)
604{
605 return (MM * MM - MP * MP) / 2. * ((CS + CSp) / (Mb - Mc)+(Mlep - Mnu) / q2 * (CV + CVp)) * f0(q2);
606}
607
608gslpp::complex MPlnu::HT(double q2)
609{
610 return -gslpp::complex::i() * lambda_half(MM*MM, MP*MP, q2) / sqrt(2.) / (MM + MP)*(CT - CTp) * fT(q2);
611}
612
613gslpp::complex MPlnu::HTt(double q2)
614{
615 return -gslpp::complex::i() * lambda_half(MM*MM, MP*MP, q2) / 2. / (MM + MP)*(CT + CTp) * fT(q2);
616}
617/*******************************************************************************
618 * Generalized angular coefficients (see 1506.03970) *
619 * ****************************************************************************/
620
621gslpp::complex MPlnu::G0(double q2)
622{
623 double lambda_MM = lambda_half(MM*MM, MP*MP, q2);
624 double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
625 double lambda_lep2 = lambda_lep*lambda_lep;
626 double Elep = sqrt(Mlep * Mlep + lambda_lep2 / (4. * q2));
627 double Enu = sqrt(Mnu * Mnu + lambda_lep2 / (4. * q2));
628 double Gprefactor = lambda_MM * lambda_lep / q2;
629
630 return Gprefactor * ((4. * (Elep * Enu + Mlep * Mnu) + lambda_lep2 / (3. * q2)) * HV(q2).abs2()+
631 (4. * (Elep * Enu - Mlep * Mnu) + lambda_lep2 / (3. * q2)) * HA(q2).abs2()+
632 (4. * (Elep * Enu - Mlep * Mnu) + lambda_lep2 / q2) * HS(q2).abs2()+
633 (4. * (Elep * Enu + Mlep * Mnu) + lambda_lep2 / q2) * HP(q2).abs2()+
634 (8. * (Elep * Enu - Mlep * Mnu) - lambda_lep2 / (12. * q2)) * HT(q2).abs2()+
635 (16. * (Elep * Enu + Mlep * Mnu) - lambda_lep2 / (12. * q2)) * HTt(q2).abs2() +
636 8. * M_SQRT2 *(Enu * Mlep - Elep * Mnu)*(HA(q2) * HT(q2).conjugate()).imag() +
637 16. * (Enu * Mlep + Elep * Mnu)*(HV(q2) * HTt(q2).conjugate()).imag());
638}
639
640gslpp::complex MPlnu::G1(double q2)
641{
642 double lambda_MM = lambda_half(MM*MM, MP*MP, q2);
643 double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
644 double Gprefactor = lambda_MM * lambda_lep / q2;
645
646 return -Gprefactor * 4. * lambda_lep * ((HA(q2)*(Mlep - Mnu) * HP(q2).conjugate()
647 + HV(q2)*(Mlep + Mnu) * HS(q2).conjugate()).real() / sqrt(q2)
648 -(sqrt(2.) * HT(q2) * HP(q2).conjugate() + 2. * HTt(q2) * HS(q2).conjugate()).imag());
649}
650
651gslpp::complex MPlnu::G2(double q2)
652{
653 double lambda_MM = lambda_half(MM*MM, MP*MP, q2);
654 double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
655 double lambda_lep2 = lambda_lep*lambda_lep;
656 double Gprefactor = lambda_MM * lambda_lep / q2;
657
658 return -Gprefactor * (4. * lambda_lep2 / (3. * q2))*(HA(q2).abs2() + HV(q2).abs2() - 2. * HT(q2).abs2() - 4. * HTt(q2).abs2());
659}
660/***************************************************************************
661 * 12 independent J angular coefficients (see again 1506.03970) *
662 * ************************************************************************/
663
664double MPlnu::J1(double q2)
665{
666 if ((q2 < Mlep * Mlep) or (q2 > (MM - MP)*(MM - MP))) return 0.;
667 return amplsq_factor * (G0(q2) - G2(q2) / 2).real();
668}
669
670double MPlnu::J2(double q2)
671{
672 if ((q2 < Mlep * Mlep) or (q2 > (MM - MP)*(MM - MP))) return 0.;
673 return amplsq_factor * G1(q2).real();
674}
675
676double MPlnu::J3(double q2)
677{
678 if ((q2 < Mlep * Mlep) or (q2 > (MM - MP)*(MM - MP))) return 0.;
679 return amplsq_factor * 3. / 2. * G2(q2).real();
680}
681
682double MPlnu::dGammadq2(double q2)
683{
684 if ((q2 < q2min) or (q2 > (MM - MP)*(MM - MP))) return 0.;
685 double sqlambdaB = lambda_half(q2,MM*MM,MP*MP);
686 double prefac = (CV-CA).abs2()*MM*sqlambdaB/192./M_PI/M_PI/M_PI;
687 double coeff_fp = (1.+Mlep*Mlep/(2.*q2))*sqlambdaB*sqlambdaB/MM/MM/MM/MM;
688 double coeff_f0 = (1.-MP*MP/MM/MM)*(1.-MP*MP/MM/MM)*3.*Mlep*Mlep/(2.*q2);
689 double TotAmp2 = coeff_fp*fplus(q2)*fplus(q2)+coeff_f0*f0(q2)*f0(q2);
690 return prefac*(1.-Mlep*Mlep/q2)*(1.-Mlep*Mlep/q2)*TotAmp2;
691 }
692
693/***************************************************************************
694 * Integration of angular coefficients Js *
695 * ************************************************************************/
696
697double MPlnu::integrateJ(int i, double q2_min, double q2_max)
698{
699
700 switch (i) {
701 case 1:
702 if (lep == StandardModel::TAU) if (checkcache_int_tau && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ1_tau;
703 if (lep == StandardModel::MU) if (checkcache_int_mu && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ1_mu;
704 if (lep == StandardModel::ELECTRON) if (checkcache_int_el && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ1_el;
705 wf=ROOT::Math::Functor1D(&(*this),&MPlnu::J1);
706 ig.SetFunction(wf);
707 J_res = ig.Integral(q2_min, q2_max);
708 if (ig.Status() != 0) return std::numeric_limits<double>::quiet_NaN();
709 return J_res;
710 break;
711 case 2:
712 if (lep == StandardModel::TAU) if (checkcache_int_tau && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ2_tau;
713 if (lep == StandardModel::MU) if (checkcache_int_mu && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ2_mu;
714 if (lep == StandardModel::ELECTRON) if (checkcache_int_el && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ2_el;
715 wf=ROOT::Math::Functor1D(&(*this),&MPlnu::J2);
716 ig.SetFunction(wf);
717 J_res = ig.Integral(q2_min, q2_max);
718 if (ig.Status() != 0) return std::numeric_limits<double>::quiet_NaN();
719 return J_res;
720 break;
721 case 3:
722 if (lep == StandardModel::TAU) if (checkcache_int_tau && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ3_tau;
723 if (lep == StandardModel::MU) if (checkcache_int_mu && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ3_mu;
724 if (lep == StandardModel::ELECTRON) if (checkcache_int_el && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ3_el;
725 wf=ROOT::Math::Functor1D(&(*this),&MPlnu::J3);
726 ig.SetFunction(wf);
727 J_res = ig.Integral(q2_min, q2_max);
728 if (ig.Status() != 0) return std::numeric_limits<double>::quiet_NaN();
729 return J_res;
730 case 4:
731 wf=ROOT::Math::Functor1D(&(*this),&MPlnu::dGammadq2);
732 ig.SetFunction(wf);
733 J_res = ig.Integral(q2_min, q2_max);
734 if (ig.Status() != 0) return std::numeric_limits<double>::quiet_NaN();
735 return J_res;
736 break;
737 default:
738 std::stringstream out;
739 out << i;
740 throw std::runtime_error("MPlnu::integrateJ: index " + out.str() + " not implemented");
741 }
742}
743
744/* misleading name, needs to be changed: this is like dGammadq2 ... */
745double MPlnu::dGammadw(double q2)
746{
747 updateParameters();
748
749 return 2. * (J1(q2) + J3(q2) / 3.);
750}
751
752/* misleading name, needs to be changed: this is like DeltaGamma only ... */
753double MPlnu::getDeltaGammaDeltaw(double w_min, double w_max)
754{
755 updateParameters();
756
757 double q2_min = (2. * MM * MP)*(w0 - w_max); // min is Mlep*Mlep;
758 double q2_max = (2. * MM * MP)*(w0 - w_min); // max is (MM-MP)*(MM-MP);
759
760 double intJ1 = integrateJ(1, q2_min, q2_max);
761 double intJ3 = integrateJ(3, q2_min, q2_max);
762
763 return 2. * (intJ1 + intJ3 / 3.);
764
765 // x-check of the SM computation
766 //return integrateJ(4, q2_min, q2_max);
767}
768
770{
771 updateParameters();
772
773#if NBGL == 3
774 return afplus_0 * afplus_0 + afplus_1 * afplus_1 + afplus_2 * afplus_2 + afplus_3 * afplus_3;
775#else
776 return afplus_0 * afplus_0 + afplus_1 * afplus_1 + afplus_2 * afplus_2;
777#endif
778}
779
781{
782 updateParameters();
783
784#if NBGL == 3
785 return af0_0 * af0_0 + af0_1 * af0_1 + af0_2 * af0_2 + af0_3 * af0_3;
786#else
787 return af0_0 * af0_0 + af0_1 * af0_1 + af0_2*af0_2;
788#endif
789}
790
792{
793 updateParameters();
794
795#if NBGL == 3
796 return 1707.54 * afplus_0 * afplus_0 + 1299.57 * afplus_0 * afplus_1 + 442.82 * afplus_1 * afplus_1 - 356.01 * afplus_0 * afplus_2
797 - 101.62 * afplus_1 * afplus_2 + 34.947 * afplus_2 * afplus_2 - 206.767 * afplus_0 * afplus_3 - 127.668 * afplus_1 * afplus_3
798 + 33.234 * afplus_2 * afplus_3 + 16.475 * afplus_3 * afplus_3;
799#else
800 return 442.82 * afplus_0 * afplus_0 - 101.619 * afplus_0 * afplus_1 + 34.947* afplus_1 * afplus_1 - 127.668 * afplus_0 * afplus_2 + 33.234 * afplus_1 * afplus_2 + 16.4754 * afplus_2 * afplus_2;
801#endif
802}
803
804double MPlnu::get_fplus(double q2)
805{
806 updateParameters();
807
808 return fplus(q2);
809}
810
811double MPlnu::get_f0(double q2)
812{
813 updateParameters();
814
815 return f0(q2);
816}
817
818double MPlnu::get_fT(double q2)
819{
820 updateParameters();
821
822 return fT(q2);
823}
@ LO
Definition: OrderScheme.h:34
@ FULLNLO
Definition: OrderScheme.h:38
@ HV
Definition: OrderScheme.h:22
bool getFlagBGL() const
Definition: Flavour.h:344
bool getFlagCLN() const
Definition: Flavour.h:340
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
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
bool getFlagDM() const
Definition: Flavour.h:348
gslpp::vector< gslpp::complex > ** ComputeCoeffdiujlknu(int i, int j, int k, double mu) const
Computes the Wilson coefficient for the Hamiltonian transitions in the JMS basis ordered as CnueduVL...
Definition: Flavour.cpp:152
QCD::lepton lep
Definition: MPlnu.h:106
virtual ~MPlnu()
Destructor.
Definition: MPlnu.cpp:116
bool btocNPpmflag
Definition: MPlnu.h:113
double z0
Definition: MPlnu.h:122
bool CLNflag
Definition: MPlnu.h:110
bool NPanalysis
Definition: MPlnu.h:114
double Mc
Definition: MPlnu.h:126
double get_fT(double q2)
return fT form factor at
Definition: MPlnu.cpp:818
std::vector< std::string > mplnuParameters
Definition: MPlnu.h:109
double getDeltaGammaDeltaw(double w_min, double w_max)
The integral of from to .
Definition: MPlnu.cpp:753
QCD::meson pseudoscalarM
Definition: MPlnu.h:108
double Mnu
Definition: MPlnu.h:118
double Mb
Definition: MPlnu.h:125
QCD::meson meson
Definition: MPlnu.h:107
std::vector< std::string > initializeMPlnuParameters()
Definition: MPlnu.cpp:119
double get_fplus(double q2)
return fplus form factor at
Definition: MPlnu.cpp:804
double mu_b
Definition: MPlnu.h:124
MPlnu(const StandardModel &SM_i, QCD::meson meson_i, QCD::meson pseudoscalar_i, QCD::lepton lep_i)
Constructor.
Definition: MPlnu.cpp:21
double get_unitarity_1min_BGL()
Weak Unitarity constraint for BGL parameters related to 1- resonances.
Definition: MPlnu.cpp:769
bool BGLflag
Definition: MPlnu.h:111
double Mlep
Definition: MPlnu.h:117
double MP
Definition: MPlnu.h:120
double w0
Definition: MPlnu.h:121
double get_strong_unitarity_BGL()
Strong Unitarity constraint for BGL parameters using HQET.
Definition: MPlnu.cpp:791
double get_unitarity_0plus_BGL()
Weak Unitarity constraint for BGL parameters related to 0+ resonances.
Definition: MPlnu.cpp:780
const StandardModel & mySM
Definition: MPlnu.h:105
double get_f0(double q2)
return f0 form factor at
Definition: MPlnu.cpp:811
double width
Definition: MPlnu.h:127
double MM
Definition: MPlnu.h:119
bool DMflag
Definition: MPlnu.h:112
double RV
Definition: MPlnu.h:123
double computeWidth() const
A method to compute the width of the meson from its lifetime.
Definition: Meson.cpp:521
std::string getModelName() const
A method to fetch the name of the model.
Definition: Model.h:59
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
meson
An enum type for mesons.
Definition: QCD.h:336
@ D_P
Definition: QCD.h:342
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
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
@ BOTTOM
Definition: QCD.h:329
@ CHARM
Definition: QCD.h:326
lepton
An enum type for leptons.
Definition: QCD.h:310
@ MU
Definition: QCD.h:314
@ ELECTRON
Definition: QCD.h:312
@ TAU
Definition: QCD.h:316
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
void initializeMeson(QCD::meson meson_i) const
A method to initialize a meson.
Definition: QCD.cpp:280
A model class for the Standard Model.
virtual const double getCCC1() const
A virtual implementation for the RealWeakEFTCC class.
const Particle & getLeptons(const QCD::lepton p) const
A get method to retrieve the member object of a lepton.
virtual const double getCCC3() const
A virtual implementation for the RealWeakEFTCC class.
const double getMz() const
A get method to access the mass of the boson .
const double Ale(double mu, orders order, bool Nf_thr=true) const
The running electromagnetic coupling in the scheme.
virtual const double getCCC4() const
A virtual implementation for the RealWeakEFTCC class.
const Flavour & getFlavour() const
virtual const double getCCC2() const
A virtual implementation for the RealWeakEFTCC class.
virtual const double getCCC5() const
A virtual implementation for the RealWeakEFTCC class.