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.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>
19 
20 MPlnu::MPlnu(const StandardModel& SM_i, QCD::meson meson_i, QCD::meson pseudoscalar_i, QCD::lepton lep_i)
21 : mySM(SM_i)
22 {
23  lep = lep_i;
24  meson = meson_i;
25  pseudoscalarM = pseudoscalar_i;
26  CLNflag = false;
27  btocNPpmflag = false;
28 
29  w_J = gsl_integration_cquad_workspace_alloc(100);
30 
31  checkcache_int_tau = 0;
32  checkcache_int_mu = 0;
33  checkcache_int_el = 0;
34 
35  double max_double = std::numeric_limits<double>::max();
36 
37  fplusz0_cache = max_double;
38  rho1to2_cache = max_double;
39  N_0_cache = max_double;
40  alpha_0_cache = max_double;
41  alpha_p_cache = max_double;
42  beta_0_cache = max_double;
43  beta_p_cache = max_double;
44  gamma_0_cache = max_double;
45  gamma_p_cache = max_double;
46  af0_1_cache = max_double;
47  af0_2_cache = max_double;
48  afplus_0_cache = max_double;
49  afplus_1_cache = max_double;
50  afplus_2_cache = max_double;
51 
52 #if NBGL == 3
53  af0_3_cache = max_double;
54  afplus_3_cache = max_double;
55 #endif
56 
57  CS_cache = max_double;
58  CSp_cache = max_double;
59  CP_cache = max_double;
60  CPp_cache = max_double;
61  CV_cache = max_double;
62  CVp_cache = max_double;
63  CA_cache = max_double;
64  CAp_cache = max_double;
65  CT_cache = max_double;
66  CTp_cache = max_double;
67 }
68 
69 MPlnu::~MPlnu()
70 {}
71 
72 std::vector<std::string> MPlnu::initializeMPlnuParameters()
73 {
74  CLNflag = mySM.getFlavour().getFlagCLN();
75  btocNPpmflag = (mySM.getModelName().compare("RealWeakEFTCCPM") == 0);
76  NPanalysis = (mySM.getModelName().compare("RealWeakEFTCCPM") == 0 || mySM.getModelName().compare("RealWeakEFTCC") == 0);
77 
78  if (pseudoscalarM == StandardModel::D_P) mplnuParameters = make_vector<std::string>()
79  << "af0_1" << "af0_2" << "afplus_0" << "afplus_1" << "afplus_2"
80 #if NBGL == 3
81  << "af0_3" << "afplus_3"
82 #endif
83  << "mBc1m_1" << "mBc1m_2" << "mBc1m_3" << "mBc1m_4"
84  << "mBc0p_1" << "mBc0p_2" << "chitildeT" << "chiL" << "nI";
85  else {
86  std::stringstream out;
87  out << pseudoscalarM;
88  throw std::runtime_error("MPlnu: vector " + out.str() + " not implemented");
89  }
90 
91  if (CLNflag) {
92  mplnuParameters.clear();
93  if (pseudoscalarM == StandardModel::D_P) mplnuParameters = make_vector<std::string>()
94  << "fplusz0" << "rho1to2"
95  << "N_0" << "alpha_0" << "alpha_p" << "beta_0" << "beta_p" << "gamma_0" << "gamma_p";
96  }
97 
98  mySM.initializeMeson(meson);
99  mySM.initializeMeson(pseudoscalarM);
100  return mplnuParameters;
101 }
102 
103 void MPlnu::updateParameters()
104 {
105  if (!mySM.getFlavour().getUpdateFlag(meson, pseudoscalarM, lep)) return;
106 
107  GF = mySM.getGF();
108  Mlep = mySM.getLeptons(lep).getMass();
109  Mnu = 0.; // neutrinos assumed to be massless
110  MM = mySM.getMesons(meson).getMass();
111  MP = mySM.getMesons(pseudoscalarM).getMass();
112  width = mySM.getMesons(meson).computeWidth();
113  w0 = (MM*MM+MP*MP)/(2.*MM*MP);
114  RV = 2.*sqrt(MM*MP)/(MM+MP);
115  mu_b = MM; // mySM.getMub();
116  Mb = mySM.getQuarks(QCD::BOTTOM).getMass(); // add the PS b mass
117  Mc = mySM.getQuarks(QCD::CHARM).getMass(); // add the PS b mass
118  Vcb = mySM.getCKM().getV_cb(); // mySM.getOptionalParameter("AbsVcb");
119  ale_mub = mySM.Ale(mu_b,FULLNLO);
120  /* Amplitude propto 4*GF*Vij/sqrt(2) & kinematics requires 1/(2^9 pi^3 MB^3) */
121  amplsq_factor = GF*GF*Vcb.abs2()/(64.*M_PI*M_PI*M_PI*MM*MM*MM);
122  q2min = Mlep*Mlep;
123  q2max = (MM-MP)*(MM-MP);
124 
125  /* SM Wilson coefficients */
126  eta_EW = (1.+ale_mub/M_PI*log(mySM.getMz()/mu_b));
127  CV_SM = 1./2.;
128  CV = CV_SM*eta_EW;
129  CA = -CV_SM*eta_EW;
130  CVp = 0.;
131  CAp = 0.;
132  CS = 0.;
133  CSp = 0.;
134  CP = 0.;
135  CPp = 0.;
136  C7 = 0.;
137  C7p = 0.;
138  CT = 0.;
139  CTp = 0.;
140 
141  /* SM + NP Wilson coefficients */
142  if (NPanalysis) {
143  if (lep == StandardModel::TAU) {
144  if (btocNPpmflag) {
145  CV += (mySM.getCCC3() - mySM.getCCC4()) / 2. / M_SQRT2;
146  CVp = (mySM.getCCC3() + mySM.getCCC4()) / 2. / M_SQRT2;
147  CA -= (mySM.getCCC3() - mySM.getCCC4()) / 2. / M_SQRT2;
148  CAp = -(mySM.getCCC3() + mySM.getCCC4()) / 2. / M_SQRT2;
149  CS = (mySM.getCCC1() - mySM.getCCC2()) / 2. / M_SQRT2;
150  CSp = (mySM.getCCC1() + mySM.getCCC2()) / 2. / M_SQRT2;
151  CP = -(mySM.getCCC1() - mySM.getCCC2()) / 2. / M_SQRT2;
152  CPp = -(mySM.getCCC1() + mySM.getCCC2()) / 2. / M_SQRT2;
153  CTp = mySM.getCCC5();
154  } else {
155  CV += mySM.getCCC3() / 2.;
156  CVp = mySM.getCCC4() / 2.;
157  CA -= mySM.getCCC3() / 2.;
158  CAp = -mySM.getCCC4() / 2.;
159  CS = mySM.getCCC1() / 2.;
160  CSp = mySM.getCCC2() / 2.;
161  CP = -mySM.getCCC1() / 2.;
162  CPp = -mySM.getCCC2() / 2.;
163  CTp = mySM.getCCC5();
164  }
165  }
166  }
167 
168  switch (pseudoscalarM) {
169  case StandardModel::D_P:
170  if (CLNflag) {
171  fplusz0 = mySM.getOptionalParameter("fplusz0");
172  rho1to2 = mySM.getOptionalParameter("rho1to2");
173  N_0 = mySM.getOptionalParameter("N_0");
174  alpha_0 = mySM.getOptionalParameter("alpha_0");
175  alpha_p = mySM.getOptionalParameter("alpha_p");
176  beta_0 = mySM.getOptionalParameter("beta_0");
177  beta_p = mySM.getOptionalParameter("beta_p");
178  gamma_0 = mySM.getOptionalParameter("gamma_0");
179  gamma_p = mySM.getOptionalParameter("gamma_p");
180  af0_1 = 0.;
181  af0_2 = 0.;
182  afplus_1 = 0.;
183  afplus_2 = 0.;
184 #if NBGL == 3
185  af0_3 = 0.;
186  afplus_3 = 0.;
187 #endif
188  mBc1m_1 = 0.;
189  mBc1m_2 = 0.;
190  mBc1m_3 = 0.;
191  mBc1m_4 = 0.;
192  mBc0p_1 = 0.;
193  mBc0p_2 = 0.;
194  chitildeT = 0.;
195  chiL = 0.;
196  nI = 0.;
197  } else {
198  fplusz0 = 0.;
199  rho1to2 = 0.;
200  N_0 = 0.;
201  alpha_0 = 0.;
202  alpha_p = 0.;
203  beta_0 = 0.;
204  beta_p = 0.;
205  gamma_0 = 0.;
206  gamma_p = 0.;
207  af0_1 = mySM.getOptionalParameter("af0_1");
208  af0_2 = mySM.getOptionalParameter("af0_2");
209  afplus_0 = mySM.getOptionalParameter("afplus_0");
210  afplus_1 = mySM.getOptionalParameter("afplus_1");
211  afplus_2 = mySM.getOptionalParameter("afplus_2");
212 #if NBGL == 3
213  af0_3 = mySM.getOptionalParameter("af0_3");
214  afplus_3 = mySM.getOptionalParameter("afplus_3");
215 #endif
216  mBc1m_1 = mySM.getOptionalParameter("mBc1m_1");
217  mBc1m_2 = mySM.getOptionalParameter("mBc1m_2");
218  mBc1m_3 = mySM.getOptionalParameter("mBc1m_3");
219  mBc1m_4 = mySM.getOptionalParameter("mBc1m_4");
220  mBc0p_1 = mySM.getOptionalParameter("mBc0p_1");
221  mBc0p_2 = mySM.getOptionalParameter("mBc0p_2");
222  chitildeT = mySM.getOptionalParameter("chitildeT");
223  chiL = mySM.getOptionalParameter("chiL");
224  nI = mySM.getOptionalParameter("nI");
225  }
226  break;
227  default:
228  std::stringstream out;
229  out << pseudoscalarM;
230  throw std::runtime_error("MPlnu: vector " + out.str() + " not implemented");
231  }
232 
233  z1m_1 = sqrt((MM+MP)*(MM+MP)-mBc1m_1*mBc1m_1)-sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP));
234  z1m_1 /= (sqrt((MM+MP)*(MM+MP)-mBc1m_1*mBc1m_1)+sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP)));
235  z1m_2 = sqrt((MM+MP)*(MM+MP)-mBc1m_2*mBc1m_2)-sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP));
236  z1m_2 /= (sqrt((MM+MP)*(MM+MP)-mBc1m_2*mBc1m_2)+sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP)));
237  z1m_3 = sqrt((MM+MP)*(MM+MP)-mBc1m_3*mBc1m_3)-sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP));
238  z1m_3 /= (sqrt((MM+MP)*(MM+MP)-mBc1m_3*mBc1m_3)+sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP)));
239 
240  z0p_1 = sqrt((MM+MP)*(MM+MP)-mBc0p_1*mBc0p_1)-sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP));
241  z0p_1 /= (sqrt((MM+MP)*(MM+MP)-mBc0p_1*mBc0p_1)+sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP)));
242  z0p_2 = sqrt((MM+MP)*(MM+MP)-mBc0p_2*mBc0p_2)-sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP));
243  z0p_2 /= (sqrt((MM+MP)*(MM+MP)-mBc0p_2*mBc0p_2)+sqrt((MM+MP)*(MM+MP)-(MM-MP)*(MM-MP)));
244 
245  if ((fplusz0 != fplusz0_cache) || (rho1to2 != rho1to2_cache)
246  || (N_0 != N_0_cache)
247  || (alpha_0 != alpha_0_cache) || (alpha_p != alpha_p_cache)
248  || (beta_0 != beta_0_cache) || (beta_p != beta_p_cache)
249  || (gamma_0 != gamma_0_cache) || (gamma_p != gamma_p_cache)
250  || (af0_1 != af0_1_cache) || (af0_2 != af0_2_cache)
251  || (afplus_0 != afplus_0_cache) || (afplus_1 != afplus_1_cache)
252  || (afplus_2 != afplus_2_cache)
253 #if NBGL == 3
254  || (af0_3 != af0_3_cache) || (afplus_3 != afplus_3_cache)
255 #endif
256  || (CS != CS_cache) || (CSp != CSp_cache)
257  || (CP != CP_cache) || (CPp != CPp_cache)
258  || (CV != CV_cache) || (CVp != CVp_cache)
259  || (CA != CA_cache) || (CAp != CAp_cache)
260  || (CT != CT_cache) || (CTp != CTp_cache)) {
261  checkcache_int_tau = 0;
262  checkcache_int_mu = 0;
263  checkcache_int_el = 0;
264  }
265 
266  if ((checkcache_int_tau == 0) || (checkcache_int_mu == 0) || (checkcache_int_el == 0)) {
267  if (lep == StandardModel::TAU) {
268  cached_intJ1_tau = integrateJ(1, q2min, q2max);
269  cached_intJ3_tau = integrateJ(3, q2min, q2max);
270  cached_intJ2_tau = 0.;
271  // cached_intJ2_tau = integrateJ(2,q2min,q2max);
272  checkcache_int_tau = 1;
273  }
274  if (lep == StandardModel::MU) {
275  cached_intJ1_mu = integrateJ(1, q2min, q2max);
276  cached_intJ3_mu = integrateJ(3, q2min, q2max);
277  cached_intJ2_mu = 0.;
278  // cached_intJ2_mu = integrateJ(2,q2min,q2max);
279  checkcache_int_mu = 1;
280  }
281  if (lep == StandardModel::ELECTRON) {
282  cached_intJ1_el = integrateJ(1, q2min, q2max);
283  cached_intJ3_el = integrateJ(3, q2min, q2max);
284  cached_intJ2_el = 0.;
285  // cached_intJ2_el = integrateJ(2,q2min,q2max);
286  checkcache_int_el = 1;
287  }
288  }
289 
290  fplusz0_cache = fplusz0;
291  rho1to2_cache = rho1to2;
292  N_0_cache = N_0;
293  alpha_0_cache = alpha_0;
294  alpha_p_cache = alpha_p;
295  beta_0_cache = beta_0;
296  beta_p_cache = beta_p;
297  gamma_0_cache = gamma_0;
298  gamma_p_cache = gamma_p;
299 
300  af0_1_cache = af0_1;
301  af0_2_cache = af0_2;
302  afplus_0_cache = afplus_0;
303  afplus_1_cache = afplus_1;
304  afplus_2_cache = afplus_2;
305 #if NBGL == 3
306  af0_3_cache = af0_3;
307  afplus_3_cache = afplus_3;
308 #endif
309  /* f+(q2=0) = f0(q2=0) */
310  z0 = (sqrt(w0+1.)-sqrt(2.))/(sqrt(w0+1.)+sqrt(2.));
311  if (CLNflag) af0_0 = 0.;
312  else {
313  af0_0 = fplus(0.);
314 #if NBGL == 3
315  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));
316 #else
317  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));
318 #endif
319  af0_0 *= phi_f0(z0)*((z0 - z0p_1) / (1. - z0 * z0p_1)*(z0 - z0p_2) / (1. - z0 * z0p_2));
320  }
321  CS_cache = CS;
322  CSp_cache = CSp;
323  CP_cache = CP;
324  CPp_cache = CPp;
325  CV_cache = CV;
326  CVp_cache = CVp;
327  CA_cache = CA;
328  CAp_cache = CAp;
329  CT_cache = CT;
330  CTp_cache = CTp;
331 
332  mySM.getFlavour().setUpdateFlag(meson, pseudoscalarM, lep, false);
333 
334  return;
335 
336 }
337 
338 /*******************************************************************************
339  * Kinematic functions *
340  * ****************************************************************************/
341 
342 double MPlnu::lambda_half(double a, double b, double c)
343 {
344  return sqrt(a*a+b*b+c*c-2.*(a*b+a*c+b*c));
345 }
346 /*******************************************************************************
347  * Form factors *
348  * ****************************************************************************/
349 
350 double MPlnu::phi_fplus(double z)
351 {
352  // chitildeT in GeV-2
353  double prefac = 8. * (MP / MM)*(MP / MM) / MM * sqrt(8. * nI / (3. * M_PI * chitildeT));
354  double num = (1. + z)*(1. + z) * sqrt(1. - z);
355  double den = (1. + MP / MM)*(1. - z) + 2. * sqrt(MP / MM)*(1. + z);
356  double den5 = den * den * den * den*den;
357  return prefac * num / den5;
358 }
359 
360 double MPlnu::phi_f0(double z)
361 {
362  // chiL dimensionless
363  double prefac = (MP / MM)*(1. - (MP / MM)*(MP / MM)) * sqrt(8. * nI / (M_PI * chiL));
364  double num = (1. - z * z) * sqrt((1. - z));
365  double den = (1. + MP / MM)*(1. - z) + 2. * sqrt(MP / MM)*(1. + z);
366  double den4 = den * den * den*den;
367  return prefac * num / den4;
368 }
369 
370 double MPlnu::fplus(double q2)
371 {
372  double w = w0 - q2 / (2. * MM * MP);
373  double z = (sqrt(w + 1.) - M_SQRT2) / (sqrt(w + 1.) + M_SQRT2);
374  if (CLNflag) {
375  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);
376  } else {
377  double P_fplus = (z - z1m_1) / (1. - z * z1m_1)*(z - z1m_2) / (1. - z * z1m_2)*(z - z1m_3) / (1. - z * z1m_3);
378 #if NBGL == 3
379  return (afplus_0 + afplus_1 * z + afplus_2 * z * z + afplus_3 * z * z * z) / phi_fplus(z) / P_fplus;
380 #else
381  return (afplus_0 + afplus_1 * z + afplus_2 * z * z) / phi_fplus(z) / P_fplus;
382 #endif
383  }
384 }
385 
386 double MPlnu::f0(double q2)
387 {
388  double w = w0 - q2 / (2. * MM * MP);
389  double z = (sqrt(w + 1.) - M_SQRT2) / (sqrt(w + 1.) + M_SQRT2);
390  if (CLNflag) {
391  double prefac0 = 2. * (MP / MM) / (1. + MP / MM)/ (1. + MP / MM)*(1. + w0);
392  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.));
393  double prefac = fplus(q2)* prefac0/(1. + w0)*(1. + w);
394  // norm introduced to respect f+(q2=0)=f0(q2=0) exactly
395  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.));
396  } else {
397  double P_f0 = (z - z0p_1) / (1. - z * z0p_1)*(z - z0p_2) / (1. - z * z0p_2);
398 #if NBGL == 3
399  return (af0_0 + af0_1 * z + af0_2 * z * z + af0_3 * z * z * z) / phi_f0(z) / P_f0;
400 #else
401  return (af0_0 + af0_1 * z + af0_2 * z * z) / phi_f0(z) / P_f0;
402 #endif
403  }
404 }
405 
406 double MPlnu::fT(double q2)
407 {
408  return 0.;
409 }
410 /********************************************************************************
411  * Helicity amplitudes (normalization such that all H \propto (mass scale)^-1) *
412  * *****************************************************************************/
413 
414 gslpp::complex MPlnu::HV(double q2)
415 {
416  return lambda_half(MM*MM, MP*MP, q2) / 2. / sqrt(q2)*((CV + CVp) * fplus(q2) + 2. * Mb / (MM + MP)*(C7 + C7p) * fT(q2));
417 }
418 
419 gslpp::complex MPlnu::HA(double q2)
420 {
421  return lambda_half(MM*MM, MP*MP, q2) / 2. / sqrt(q2)*(CA + CAp) * fplus(q2);
422 }
423 
424 gslpp::complex MPlnu::HP(double q2)
425 {
426  return (MM * MM - MP * MP) / 2. * ((CP + CPp) / (Mb - Mc)+(Mlep + Mnu) / q2 * (CA + CAp)) * f0(q2);
427 }
428 
429 gslpp::complex MPlnu::HS(double q2)
430 {
431  return (MM * MM - MP * MP) / 2. * ((CS + CSp) / (Mb - Mc)+(Mlep - Mnu) / q2 * (CV + CVp)) * f0(q2);
432 }
433 
434 gslpp::complex MPlnu::HT(double q2)
435 {
436  return -gslpp::complex::i() * lambda_half(MM*MM, MP*MP, q2) / sqrt(2.) / (MM + MP)*(CT - CTp) * fT(q2);
437 }
438 
439 gslpp::complex MPlnu::HTt(double q2)
440 {
441  return -gslpp::complex::i() * lambda_half(MM*MM, MP*MP, q2) / 2. / (MM + MP)*(CT + CTp) * fT(q2);
442 }
443 /*******************************************************************************
444  * Generalized angular coefficients (see 1506.03970) *
445  * ****************************************************************************/
446 
447 gslpp::complex MPlnu::G0(double q2)
448 {
449  double lambda_MM = lambda_half(MM*MM, MP*MP, q2);
450  double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
451  double lambda_lep2 = lambda_lep*lambda_lep;
452  double Elep = sqrt(Mlep * Mlep + lambda_lep2 / (4. * q2));
453  double Enu = sqrt(Mnu * Mnu + lambda_lep2 / (4. * q2));
454  double Gprefactor = lambda_MM * lambda_lep / q2;
455 
456  return Gprefactor * ((4. * (Elep * Enu + Mlep * Mnu) + lambda_lep2 / (3. * q2)) * HV(q2).abs2()+
457  (4. * (Elep * Enu - Mlep * Mnu) + lambda_lep2 / (3. * q2)) * HA(q2).abs2()+
458  (4. * (Elep * Enu - Mlep * Mnu) + lambda_lep2 / q2) * HS(q2).abs2()+
459  (4. * (Elep * Enu + Mlep * Mnu) + lambda_lep2 / q2) * HP(q2).abs2()+
460  (8. * (Elep * Enu - Mlep * Mnu) - lambda_lep2 / (12. * q2)) * HT(q2).abs2()+
461  (16. * (Elep * Enu + Mlep * Mnu) - lambda_lep2 / (12. * q2)) * HTt(q2).abs2() +
462  8. * M_SQRT2 *(Enu * Mlep - Elep * Mnu)*(HA(q2) * HT(q2).conjugate()).imag() +
463  16. * (Enu * Mlep + Elep * Mnu)*(HV(q2) * HTt(q2).conjugate()).imag());
464 }
465 
466 gslpp::complex MPlnu::G1(double q2)
467 {
468  double lambda_MM = lambda_half(MM*MM, MP*MP, q2);
469  double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
470  double Gprefactor = lambda_MM * lambda_lep / q2;
471 
472  return -Gprefactor * 4. * lambda_lep * ((HA(q2)*(Mlep - Mnu) * HP(q2).conjugate()
473  + HV(q2)*(Mlep + Mnu) * HS(q2).conjugate()).real() / sqrt(q2)
474  -(sqrt(2.) * HT(q2) * HP(q2).conjugate() + 2. * HTt(q2) * HS(q2).conjugate()).imag());
475 }
476 
477 gslpp::complex MPlnu::G2(double q2)
478 {
479  double lambda_MM = lambda_half(MM*MM, MP*MP, q2);
480  double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
481  double lambda_lep2 = lambda_lep*lambda_lep;
482  double Gprefactor = lambda_MM * lambda_lep / q2;
483 
484  return -Gprefactor * (4. * lambda_lep2 / (3. * q2))*(HA(q2).abs2() + HV(q2).abs2() - 2. * HT(q2).abs2() - 4. * HTt(q2).abs2());
485 }
486 /***************************************************************************
487  * 12 independent J angular coefficients (see again 1506.03970) *
488  * ************************************************************************/
489 
490 double MPlnu::J1(double q2)
491 {
492  if ((q2 < Mlep * Mlep) or (q2 > (MM - MP)*(MM - MP))) return 0.;
493  return amplsq_factor * (G0(q2) - G2(q2) / 2).real();
494 }
495 
496 double MPlnu::J2(double q2)
497 {
498  if ((q2 < Mlep * Mlep) or (q2 > (MM - MP)*(MM - MP))) return 0.;
499  return amplsq_factor * G1(q2).real();
500 }
501 
502 double MPlnu::J3(double q2)
503 {
504  if ((q2 < Mlep * Mlep) or (q2 > (MM - MP)*(MM - MP))) return 0.;
505  return amplsq_factor * 3. / 2. * G2(q2).real();
506 }
507 
508 double MPlnu::dGammadq2(double q2)
509 {
510  if ((q2 < q2min) or (q2 > (MM - MP)*(MM - MP))) return 0.;
511  double sqlambdaB = lambda_half(q2,MM*MM,MP*MP);
512  double prefac = (CV-CA)*(CV-CA)*GF*GF*Vcb.abs2()*MM*sqlambdaB/192./M_PI/M_PI/M_PI;
513  double coeff_fp = (1.+Mlep*Mlep/(2.*q2))*sqlambdaB*sqlambdaB/MM/MM/MM/MM;
514  double coeff_f0 = (1.-MP*MP/MM/MM)*(1.-MP*MP/MM/MM)*3.*Mlep*Mlep/(2.*q2);
515  double TotAmp2 = coeff_fp*fplus(q2)*fplus(q2)+coeff_f0*f0(q2)*f0(q2);
516  return prefac*(1.-Mlep*Mlep/q2)*(1.-Mlep*Mlep/q2)*TotAmp2;
517  }
518 
519 /***************************************************************************
520  * Integration of angular coefficients Js *
521  * ************************************************************************/
522 
523 double MPlnu::integrateJ(int i, double q2_min, double q2_max)
524 {
525  old_handler = gsl_set_error_handler_off();
526 
527  switch (i) {
528  case 1:
529  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ1_tau;
530  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ1_mu;
531  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ1_el;
532  FJ = convertToGslFunction(boost::bind(&MPlnu::J1, &(*this), _1));
533  if (gsl_integration_cquad(&FJ, q2_min, q2_max, 1.e-2, 1.e-1, w_J, &J_res, &J_err, NULL) != 0) std::numeric_limits<double>::quiet_NaN();
534  gsl_set_error_handler(old_handler);
535  return J_res;
536  break;
537  case 2:
538  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ2_tau;
539  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ2_mu;
540  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ2_el;
541  FJ = convertToGslFunction(boost::bind(&MPlnu::J2, &(*this), _1));
542  if (gsl_integration_cquad(&FJ, q2_min, q2_max, 1.e-2, 1.e-1, w_J, &J_res, &J_err, NULL) != 0) std::numeric_limits<double>::quiet_NaN();
543  gsl_set_error_handler(old_handler);
544  return J_res;
545  break;
546  case 3:
547  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ3_tau;
548  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ3_mu;
549  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ3_el;
550  FJ = convertToGslFunction(boost::bind(&MPlnu::J3, &(*this), _1));
551  if (gsl_integration_cquad(&FJ, q2_min, q2_max, 1.e-2, 1.e-1, w_J, &J_res, &J_err, NULL) != 0) std::numeric_limits<double>::quiet_NaN();
552  gsl_set_error_handler(old_handler);
553  return J_res;
554  case 4:
555  FJ = convertToGslFunction(boost::bind(&MPlnu::dGammadq2, &(*this), _1));
556  if (gsl_integration_cquad(&FJ, q2_min, q2_max, 1.e-2, 1.e-1, w_J, &J_res, &J_err, NULL) != 0) std::numeric_limits<double>::quiet_NaN();
557  gsl_set_error_handler(old_handler);
558  return J_res;
559  break;
560  default:
561  gsl_set_error_handler(old_handler);
562  std::stringstream out;
563  out << i;
564  throw std::runtime_error("MPlnu::integrateJ: index " + out.str() + " not implemented");
565  }
566 }
567 
568 /* misleading name, needs to be changed: this is like dGammadq2 ... */
569 double MPlnu::dGammadw(double q2)
570 {
571  updateParameters();
572 
573  return 2. * (J1(q2) + J3(q2) / 3.);
574 }
575 
576 /* misleading name, needs to be changed: this is like DeltaGamma only ... */
577 double MPlnu::getDeltaGammaDeltaw(double w_min, double w_max)
578 {
579  updateParameters();
580 
581  double q2_min = (2. * MM * MP)*(w0 - w_max); // min is Mlep*Mlep;
582  double q2_max = (2. * MM * MP)*(w0 - w_min); // max is (MM-MP)*(MM-MP);
583 
584  double intJ1 = integrateJ(1, q2_min, q2_max);
585  double intJ3 = integrateJ(3, q2_min, q2_max);
586 
587  return 2. * (intJ1 + intJ3 / 3.);
588 
589  // x-check of the SM computation
590  //return integrateJ(4, q2_min, q2_max);
591 }
592 
593 double MPlnu::get_unitarity_1min_BGL()
594 {
595  updateParameters();
596 
597 #if NBGL == 3
598  return afplus_0 * afplus_0 + afplus_1 * afplus_1 + afplus_2 * afplus_2 + afplus_3 * afplus_3;
599 #else
600  return afplus_0 * afplus_0 + afplus_1 * afplus_1 + afplus_2 * afplus_2;
601 #endif
602 }
603 
604 double MPlnu::get_unitarity_0plus_BGL()
605 {
606  updateParameters();
607 
608 #if NBGL == 3
609  return af0_0 * af0_0 + af0_1 * af0_1 + af0_2 * af0_2 + af0_3 * af0_3;
610 #else
611  return af0_0 * af0_0 + af0_1 * af0_1 + af0_2*af0_2;
612 #endif
613 }
614 
615 double MPlnu::get_strong_unitarity_BGL()
616 {
617  updateParameters();
618 
619 #if NBGL == 3
620  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
621  - 101.62 * afplus_1 * afplus_2 + 34.947 * afplus_2 * afplus_2 - 206.767 * afplus_0 * afplus_3 - 127.668 * afplus_1 * afplus_3
622  + 33.234 * afplus_2 * afplus_3 + 16.475 * afplus_3 * afplus_3;
623 #else
624  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;
625 #endif
626 }
627 
628 double MPlnu::get_fplus(double q2)
629 {
630  updateParameters();
631 
632  return fplus(q2);
633 }
634 
635 double MPlnu::get_f0(double q2)
636 {
637  updateParameters();
638 
639  return f0(q2);
640 }
641 
642 double MPlnu::get_fT(double q2)
643 {
644  updateParameters();
645 
646  return fT(q2);
647 }
QCD::TAU
Definition: QCD.h:316
std_make_vector.h
MPlnu::initializeMPlnuParameters
std::vector< std::string > initializeMPlnuParameters()
Definition: MPlnu.cpp:72
MPlnu::get_fT
double get_fT(double q2)
return fT form factor at
Definition: MPlnu.cpp:636
MPlnu::RV
double RV
Definition: MPlnu.h:119
QCD::BOTTOM
Definition: QCD.h:329
MPlnu::mySM
const StandardModel & mySM
Definition: MPlnu.h:103
MPlnu::width
double width
Definition: MPlnu.h:123
StandardModel::getCCC5
virtual double getCCC5() const
A virtual implementation for the RealWeakEFTCC class.
Definition: StandardModel.h:1164
MPlnu::GF
double GF
Definition: MPlnu.h:112
MPlnu::~MPlnu
virtual ~MPlnu()
Destructor.
Definition: MPlnu.cpp:69
StandardModel::getCCC1
virtual double getCCC1() const
A virtual implementation for the RealWeakEFTCC class.
Definition: StandardModel.h:1144
MPlnu::btocNPpmflag
bool btocNPpmflag
Definition: MPlnu.h:109
QCD::initializeMeson
void initializeMeson(QCD::meson meson_i) const
A method to initialize a meson.
Definition: QCD.cpp:236
MPlnu::z0
double z0
Definition: MPlnu.h:118
make_vector
Definition: std_make_vector.h:15
QCD::D_P
Definition: QCD.h:342
MPlnu::get_f0
double get_f0(double q2)
return f0 form factor at
Definition: MPlnu.cpp:629
MPlnu::meson
QCD::meson meson
Definition: MPlnu.h:105
MPlnu::Mc
double Mc
Definition: MPlnu.h:122
StandardModel.h
MPlnu::MPlnu
MPlnu(const StandardModel &SM_i, QCD::meson meson_i, QCD::meson pseudoscalar_i, QCD::lepton lep_i)
Constructor.
Definition: MPlnu.cpp:20
QCD::CHARM
Definition: QCD.h:326
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:218
StandardModel::getCCC2
virtual double getCCC2() const
A virtual implementation for the RealWeakEFTCC class.
Definition: StandardModel.h:1149
gslpp::complex
A class for defining operations on and functions of complex numbers.
Definition: gslpp_complex.h:35
CKM::getV_cb
gslpp::complex getV_cb() const
A member for returning the value of the CKM element .
Definition: CKM.h:237
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:243
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
MPlnu::get_unitarity_0plus_BGL
double get_unitarity_0plus_BGL()
Weak Unitarity constraint for BGL parameters related to 0+ resonances.
Definition: MPlnu.cpp:598
QCD::ELECTRON
Definition: QCD.h:312
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:474
Flavour::getFlagCLN
bool getFlagCLN() const
Definition: Flavour.h:247
MPlnu::MM
double MM
Definition: MPlnu.h:115
Meson::computeWidth
double computeWidth() const
A method to compute the width of the meson from its lifetime.
Definition: Meson.cpp:479
MPlnu::Mb
double Mb
Definition: MPlnu.h:121
MPlnu::pseudoscalarM
QCD::meson pseudoscalarM
Definition: MPlnu.h:106
MPlnu::get_strong_unitarity_BGL
double get_strong_unitarity_BGL()
Strong Unitarity constraint for BGL parameters using HQET.
Definition: MPlnu.cpp:609
gslpp::complex::conjugate
complex conjugate() const
Definition: gslpp_complex.cpp:288
Particle::getMass
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
MPlnu::get_fplus
double get_fplus(double q2)
return fplus form factor at
Definition: MPlnu.cpp:622
StandardModel::getCCC4
virtual double getCCC4() const
A virtual implementation for the RealWeakEFTCC class.
Definition: StandardModel.h:1159
MPlnu::get_unitarity_1min_BGL
double get_unitarity_1min_BGL()
Weak Unitarity constraint for BGL parameters related to 1- resonances.
Definition: MPlnu.cpp:587
MPlnu::NPanalysis
bool NPanalysis
Definition: MPlnu.h:110
QCD::meson
meson
An enum type for mesons.
Definition: QCD.h:336
MPlnu::Mnu
double Mnu
Definition: MPlnu.h:114
StandardModel::getGF
double getGF() const
A get method to retrieve the Fermi constant .
Definition: StandardModel.h:736
gslpp::sqrt
complex sqrt(const complex &z)
Definition: gslpp_complex.cpp:385
QCD::getMesons
Meson getMesons(const QCD::meson m) const
A get method to access a meson as an object of the type Meson.
Definition: QCD.h:524
gslpp::complex::i
static const complex & i()
Definition: gslpp_complex.cpp:154
MPlnu.h
MPlnu::getDeltaGammaDeltaw
double getDeltaGammaDeltaw(double w_min, double w_max)
The integral of from to .
Definition: MPlnu.cpp:571
MPlnu::Mlep
double Mlep
Definition: MPlnu.h:113
StandardModel::Ale
double Ale(double mu, orders order, bool Nf_thr=true) const
The running electromagnetic coupling in the scheme.
Definition: StandardModel.cpp:706
gslpp_function_adapter.h
QCD::getQuarks
Particle getQuarks(const QCD::quark q) const
A get method to access a quark as an object of the type Particle.
Definition: QCD.h:534
MPlnu::w0
double w0
Definition: MPlnu.h:117
MPlnu::mu_b
double mu_b
Definition: MPlnu.h:120
StandardModel::getFlavour
const Flavour & getFlavour() const
Definition: StandardModel.h:1006
MPlnu::mplnuParameters
std::vector< std::string > mplnuParameters
Definition: MPlnu.h:107
StandardModel::getMz
double getMz() const
A get method to access the mass of the boson .
Definition: StandardModel.h:718
af0_0
Definition: MPlnuObservables.h:195
MPlnu::CLNflag
bool CLNflag
Definition: MPlnu.h:108
Model::getModelName
std::string getModelName() const
A method to fetch the name of the model.
Definition: Model.h:59
HV
Definition: OrderScheme.h:22
QCD::getOptionalParameter
double getOptionalParameter(std::string name) const
A method to get parameters that are specific to only one set of observables.
Definition: QCD.h:448
StandardModel::getCCC3
virtual double getCCC3() const
A virtual implementation for the RealWeakEFTCC class.
Definition: StandardModel.h:1154
MPlnu::MP
double MP
Definition: MPlnu.h:116
convertToGslFunction
gsl_function convertToGslFunction(const F &f)
Definition: gslpp_function_adapter.h:24
MPlnu::lep
QCD::lepton lep
Definition: MPlnu.h:104
FULLNLO
Definition: OrderScheme.h:37
StandardModel::getCKM
CKM getCKM() const
A get method to retrieve the member object of type CKM.
Definition: StandardModel.h:876
QCD::MU
Definition: QCD.h:314
QCD::lepton
lepton
An enum type for leptons.
Definition: QCD.h:310
StandardModel::getLeptons
Particle getLeptons(const QCD::lepton p) const
A get method to retrieve the member object of a lepton.
Definition: StandardModel.h:709