a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
MVlnu.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 "MVlnu.h"
10 #include "std_make_vector.h"
11 #include "gslpp_function_adapter.h"
12 #include <boost/bind.hpp>
13 #include <limits>
14 #include <gsl/gsl_sf_gegenbauer.h>
15 #include <gsl/gsl_sf_expint.h>
16 #include <limits>
17 
18 MVlnu::MVlnu(const StandardModel& SM_i, QCD::meson meson_i, QCD::meson vector_i, QCD::lepton lep_i)
19 : mySM(SM_i)
20 {
21  lep = lep_i;
22  meson = meson_i;
23  vectorM = vector_i;
24  CLNflag = false;
25  btocNPpmflag = false;
26 
27  w_J = gsl_integration_cquad_workspace_alloc (100);
28 
29  checkcache_int_tau = 0;
30  checkcache_int_mu = 0;
31  checkcache_int_el = 0;
32 
33  double max_double = std::numeric_limits<double>::max();
34 
35  hA1w1_cache = max_double;
36  rho2_cache = max_double;
37  R1w1_cache = max_double;
38  R2w1_cache = max_double;
39  N_A_cache = max_double;
40  N_1_cache =max_double;
41  N_2_cache = max_double;
42  j_A_cache= max_double;
43  j_0_cache = max_double;
44  j_1_cache = max_double;
45  j_2_cache = max_double;
46  k_A_cache = max_double;
47  k_0_cache = max_double;
48  k_1_cache = max_double;
49  k_2_cache = max_double;
50  l_A_cache = max_double;
51 
52  af0_cache = max_double;
53  af1_cache = max_double;
54  af2_cache = max_double;
55  ag0_cache = max_double;
56  ag1_cache = max_double;
57  ag2_cache = max_double;
58  aF11_cache = max_double;
59  aF12_cache = max_double;
60  aF21_cache = max_double;
61  aF22_cache = max_double;
62 
63  CS_cache = max_double;
64  CSp_cache = max_double;
65  CP_cache = max_double;
66  CPp_cache = max_double;
67  CV_cache = max_double;
68  CVp_cache = max_double;
69  CA_cache = max_double;
70  CAp_cache = max_double;
71  CT_cache = max_double;
72  CTp_cache = max_double;
73 }
74 
75 MVlnu::~MVlnu() {
76 }
77 
78 std::vector<std::string> MVlnu::initializeMVlnuParameters()
79 {
80  CLNflag = mySM.getFlavour().getFlagCLN();
81  btocNPpmflag = (mySM.getModelName().compare("RealWeakEFTCCPM") == 0);
82  NPanalysis = (mySM.getModelName().compare("RealWeakEFTCCPM") == 0 || mySM.getModelName().compare("RealWeakEFTCC") == 0);
83 
84  if (vectorM == StandardModel::D_star_P) mvlnuParameters = make_vector<std::string>()
85  << "af0" << "af1" << "af2" << "ag0" << "ag1" << "ag2"
86  << "aF11" << "aF12" << "aF21" << "aF22"
87  << "mBcstV1" << "mBcstV2" << "mBcstV3" << "mBcstV4"
88  << "mBcstA1" << "mBcstA2" << "mBcstA3" << "mBcstA4"
89  << "mBcstP1" << "mBcstP2" << "mBcstP3"
90  << "chiTV" << "chiTA" << "chiTP" << "nI";
91  else {
92  std::stringstream out;
93  out << vectorM;
94  throw std::runtime_error("MVlnu: vector " + out.str() + " not implemented");
95  }
96 
97  if (CLNflag) {
98  mvlnuParameters.clear();
99  if (vectorM == StandardModel::D_star_P) mvlnuParameters = make_vector<std::string>()
100  << "hA1w1" << "rho2" << "R1w1" << "R2w1"
101  << "N_A" << "N_1" << "N_2" << "j_A" << "j_0" << "j_1" << "j_2"
102  << "k_A" << "k_0" << "k_1" << "k_2" << "l_A";
103  }
104 
105  mySM.initializeMeson(meson);
106  mySM.initializeMeson(vectorM);
107  return mvlnuParameters;
108 }
109 
110 void MVlnu::updateParameters()
111 {
112  if (!mySM.getFlavour().getUpdateFlag(meson, vectorM, lep)) return;
113 
114  GF = mySM.getGF();
115  Mlep = mySM.getLeptons(lep).getMass();
116  Mnu = 0.; // neutrinos assumed to be massless
117  MM = mySM.getMesons(meson).getMass();
118  MV = mySM.getMesons(vectorM).getMass();
119  width = mySM.getMesons(meson).computeWidth();
120  w0 = (MM*MM+MV*MV)/(2.*MM*MV);
121  RV = 2.*sqrt(MM*MV)/(MM+MV);
122  mu_b = MM; // mySM.getMub();
123  Mb = mySM.getQuarks(QCD::BOTTOM).getMass(); // add the PS b mass
124  Mc = mySM.getQuarks(QCD::CHARM).getMass(); // add the PS b mass
125  Vcb = mySM.getCKM().getV_cb(); // mySM.getOptionalParameter("AbsVcb");
126  ale_mub = mySM.Ale(mu_b,FULLNLO);
127  /* Amplitude propto 4*GF*Vij/sqrt(2) & kinematics requires 1/(2^9 pi^3 MB^3) */
128  amplsq_factor = GF*GF*Vcb.abs2()/(64.*M_PI*M_PI*M_PI*MM*MM*MM);
129  q2min = Mlep*Mlep;
130  q2max = (MM-MV)*(MM-MV);
131 
132  MV_o_MM = MV / MM;
133  sqrtMV_o_MM = sqrt(MV_o_MM);
134 
135  /* SM Wilson coefficients */
136  CV_SM = 1./2.*(1.+ale_mub/M_PI*log(mySM.getMz()/mu_b));
137  CV = CV_SM;
138  CA = -CV_SM;
139  CVp = 0.;
140  CAp = 0.;
141  CS = 0.;
142  CSp = 0.;
143  CP = 0.;
144  CPp = 0.;
145  C7 = 0.;
146  C7p = 0.;
147  CT = 0.;
148  CTp = 0.;
149 
150  /* SM + NP Wilson coefficients */
151  if (NPanalysis) {
152  if (lep == StandardModel::TAU) {
153  if (btocNPpmflag) {
154  CV += (mySM.getCCC3() - mySM.getCCC4()) / 2. / M_SQRT2;
155  CVp = (mySM.getCCC3() + mySM.getCCC4()) / 2. / M_SQRT2;
156  CA -= (mySM.getCCC3() - mySM.getCCC4()) / 2. / M_SQRT2;
157  CAp = -(mySM.getCCC3() + mySM.getCCC4()) / 2. / M_SQRT2;
158  CS = (mySM.getCCC1() - mySM.getCCC2()) / 2. / M_SQRT2;
159  CSp = (mySM.getCCC1() + mySM.getCCC2()) / 2. / M_SQRT2;
160  CP = -(mySM.getCCC1() - mySM.getCCC2()) / 2. / M_SQRT2;
161  CPp = -(mySM.getCCC1() + mySM.getCCC2()) / 2. / M_SQRT2;
162  CTp = mySM.getOptionalParameter("CT_NP");
163  } else {
164  CV += mySM.getCCC3() / 2.;
165  CVp = mySM.getCCC4() / 2.;
166  CA -= mySM.getCCC3() / 2.;
167  CAp = -mySM.getCCC4() / 2.;
168  CS = mySM.getCCC1() / 2.;
169  CSp = mySM.getCCC2() / 2.;
170  CP = -mySM.getCCC1() / 2.;
171  CPp = -mySM.getCCC2() / 2.;
172  CTp = mySM.getCCC5();
173  }
174  }
175  }
176 
177  switch (vectorM) {
178  case StandardModel::D_star_P:
179  if (CLNflag) {
180  hA1w1 = mySM.getOptionalParameter("hA1w1");
181  rho2 = mySM.getOptionalParameter("rho2");
182  R1w1 = mySM.getOptionalParameter("R1w1");
183  R2w1 = mySM.getOptionalParameter("R2w1");
184  N_A = mySM.getOptionalParameter("N_A");
185  N_1 = mySM.getOptionalParameter("N_1");
186  N_2 = mySM.getOptionalParameter("N_2");
187  j_A = mySM.getOptionalParameter("j_A");
188  j_0 = mySM.getOptionalParameter("j_0");
189  j_1 = mySM.getOptionalParameter("j_1");
190  j_2 = mySM.getOptionalParameter("j_2");
191  k_A = mySM.getOptionalParameter("k_A");
192  k_0 = mySM.getOptionalParameter("k_0");
193  k_1 = mySM.getOptionalParameter("k_1");
194  k_2 = mySM.getOptionalParameter("k_2");
195  l_A = mySM.getOptionalParameter("l_A");
196  af0 = 0.;
197  af1 = 0.;
198  af2 = 0.;
199  ag0 = 0.;
200  ag1 = 0.;
201  ag2 = 0.;
202  aF11 = 0.;
203  aF12 = 0.;
204  aF21 = 0.;
205  aF22 = 0.;
206  mBcstV1 = 0.;
207  mBcstV2 = 0.;
208  mBcstV3 = 0.;
209  mBcstV4 = 0.;
210  mBcstA1 = 0.;
211  mBcstA2 = 0.;
212  mBcstA3 = 0.;
213  mBcstA4 = 0.;
214  mBcstP1 = 0.;
215  mBcstP2 = 0.;
216  mBcstP3 = 0.;
217  chiTV = 0.;
218  chiTA = 0.;
219  chiTP = 0.;
220  nI = 0.;
221  } else {
222  hA1w1 = 0.;
223  rho2 = 0.;
224  R1w1 = 0.;
225  R2w1 = 0.;
226  N_A = 0.;
227  N_1 = 0.;
228  N_2 = 0.;
229  j_A = 0.;
230  j_0 = 0.;
231  j_1 = 0.;
232  j_2 = 0.;
233  k_A = 0.;
234  k_0 = 0.;
235  k_1 = 0.;
236  k_2 = 0.;
237  l_A = 0.;
238  af0 = mySM.getOptionalParameter("af0");
239  af1 = mySM.getOptionalParameter("af1");
240  af2 = mySM.getOptionalParameter("af2");
241  ag0 = mySM.getOptionalParameter("ag0");
242  ag1 = mySM.getOptionalParameter("ag1");
243  ag2 = mySM.getOptionalParameter("ag2");
244  aF11 = mySM.getOptionalParameter("aF11");
245  aF12 = mySM.getOptionalParameter("aF12");
246  aF21 = mySM.getOptionalParameter("aF21");
247  aF22 = mySM.getOptionalParameter("aF22");
248  mBcstV1 = mySM.getOptionalParameter("mBcstV1");
249  mBcstV2 = mySM.getOptionalParameter("mBcstV2");
250  mBcstV3 = mySM.getOptionalParameter("mBcstV3");
251  mBcstV4 = mySM.getOptionalParameter("mBcstV4");
252  mBcstA1 = mySM.getOptionalParameter("mBcstA1");
253  mBcstA2 = mySM.getOptionalParameter("mBcstA2");
254  mBcstA3 = mySM.getOptionalParameter("mBcstA3");
255  mBcstA4 = mySM.getOptionalParameter("mBcstA4");
256  mBcstP1 = mySM.getOptionalParameter("mBcstP1");
257  mBcstP2 = mySM.getOptionalParameter("mBcstP2");
258  mBcstP3 = mySM.getOptionalParameter("mBcstP3");
259  chiTV = mySM.getOptionalParameter("chiTV");
260  chiTA = mySM.getOptionalParameter("chiTA");
261  chiTP = mySM.getOptionalParameter("chiTP");
262  nI = mySM.getOptionalParameter("nI");
263  }
264  break;
265  default:
266  std::stringstream out;
267  out << vectorM;
268  throw std::runtime_error("MVlnu: vector " + out.str() + " not implemented");
269  }
270 
271  zV1 = sqrt((MM+MV)*(MM+MV)-mBcstV1*mBcstV1)-sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV));
272  zV1 /= (sqrt((MM+MV)*(MM+MV)-mBcstV1*mBcstV1)+sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV)));
273  zV2 = sqrt((MM+MV)*(MM+MV)-mBcstV2*mBcstV2)-sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV));
274  zV2 /= (sqrt((MM+MV)*(MM+MV)-mBcstV2*mBcstV2)+sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV)));
275  zV3 = sqrt((MM+MV)*(MM+MV)-mBcstV3*mBcstV3)-sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV));
276  zV3 /= (sqrt((MM+MV)*(MM+MV)-mBcstV3*mBcstV3)+sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV)));
277  zV4 = sqrt((MM+MV)*(MM+MV)-mBcstV4*mBcstV4)-sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV));
278  zV4 /= (sqrt((MM+MV)*(MM+MV)-mBcstV4*mBcstV4)+sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV)));
279 
280  zA1 = sqrt((MM+MV)*(MM+MV)-mBcstA1*mBcstA1)-sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV));
281  zA1 /= (sqrt((MM+MV)*(MM+MV)-mBcstA1*mBcstA1)+sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV)));
282  zA2 = sqrt((MM+MV)*(MM+MV)-mBcstA2*mBcstA2)-sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV));
283  zA2 /= (sqrt((MM+MV)*(MM+MV)-mBcstA2*mBcstA2)+sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV)));
284  zA3 = sqrt((MM+MV)*(MM+MV)-mBcstA3*mBcstA3)-sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV));
285  zA3 /= (sqrt((MM+MV)*(MM+MV)-mBcstA3*mBcstA3)+sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV)));
286  zA4 = sqrt((MM+MV)*(MM+MV)-mBcstA4*mBcstA4)-sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV));
287  zA4 /= (sqrt((MM+MV)*(MM+MV)-mBcstA4*mBcstA4)+sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV)));
288 
289  zP1 = sqrt((MM+MV)*(MM+MV)-mBcstP1*mBcstP1)-sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV));
290  zP1 /= (sqrt((MM+MV)*(MM+MV)-mBcstP1*mBcstP1)+sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV)));
291  zP2 = sqrt((MM+MV)*(MM+MV)-mBcstP2*mBcstP2)-sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV));
292  zP2 /= (sqrt((MM+MV)*(MM+MV)-mBcstP2*mBcstP2)+sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV)));
293  zP3 = sqrt((MM+MV)*(MM+MV)-mBcstP3*mBcstP3)-sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV));
294  zP3 /= (sqrt((MM+MV)*(MM+MV)-mBcstP3*mBcstP3)+sqrt((MM+MV)*(MM+MV)-(MM-MV)*(MM-MV)));
295 
296  if ((hA1w1 != hA1w1_cache) || (rho2 != rho2_cache) || (R1w1 != R1w1_cache) || (R2w1 != R2w1_cache)
297  || (N_A != N_A_cache) || (N_1 != N_1_cache) || (N_2 != N_2_cache)
298  || (j_A != j_A_cache) || (j_0 != j_0_cache) || (j_1 != j_1_cache) || (j_2 != j_2_cache)
299  || (k_A != k_A_cache) || (k_0 != k_0_cache) || (k_1 != k_1_cache) || (k_2 != k_2_cache)
300  || (l_A != l_A_cache)
301  || (af0 != af0_cache) || (af1 != af1_cache) || (af2 != af2_cache)
302  || (ag0 != ag0_cache) || (ag1 != af1_cache) || (ag2 != af2_cache)
303  || (aF11 != aF11_cache) || (aF12 != aF12_cache)
304  || (aF21 != aF21_cache) || (aF22 != aF22_cache)
305  || (CS != CS_cache) || (CSp != CSp_cache)
306  || (CP != CP_cache) || (CPp != CPp_cache)
307  || (CV != CV_cache) || (CVp != CVp_cache)
308  || (CA != CA_cache) || (CAp != CAp_cache)
309  || (CT != CT_cache) || (CTp != CTp_cache)) {
310  checkcache_int_tau = 0;
311  checkcache_int_mu = 0;
312  checkcache_int_el = 0;
313  }
314 
315  if ((checkcache_int_tau == 0) || (checkcache_int_mu == 0) || (checkcache_int_el == 0)) {
316  if (lep == StandardModel::TAU) {
317  cached_intJ1s_tau = integrateJ(1, q2min, q2max);
318  cached_intJ1c_tau = integrateJ(2, q2min, q2max);
319  cached_intJ2s_tau = integrateJ(3, q2min, q2max);
320  cached_intJ2c_tau = integrateJ(4, q2min, q2max);
321  cached_intJ3_tau = integrateJ(5, q2min, q2max);
322  cached_intJ6s_tau = integrateJ(8, q2min, q2max);
323  cached_intJ6c_tau = integrateJ(9, q2min, q2max);
324  cached_intJ9_tau = integrateJ(12, q2min, q2max);
325  cached_intJ4_tau = 0.;
326  cached_intJ5_tau = 0.;
327  cached_intJ7_tau = 0.;
328  cached_intJ8_tau = 0.;
329  // not needed at present
330  /*
331  cached_intJ4_tau = integrateJ(6,q2min,q2max);
332  cached_intJ5_tau = integrateJ(7,q2min,q2max);
333  cached_intJ7_tau = integrateJ(10,q2min,q2max);
334  cached_intJ8_tau = integrateJ(11,q2min,q2max);
335  */
336  checkcache_int_tau = 1;
337  }
338  if (lep == StandardModel::MU) {
339  cached_intJ1s_mu = integrateJ(1, q2min, q2max);
340  cached_intJ1c_mu = integrateJ(2, q2min, q2max);
341  cached_intJ2s_mu = integrateJ(3, q2min, q2max);
342  cached_intJ2c_mu = integrateJ(4, q2min, q2max);
343  cached_intJ3_mu = integrateJ(5, q2min, q2max);
344  cached_intJ6s_mu = integrateJ(8, q2min, q2max);
345  cached_intJ6c_mu = integrateJ(9, q2min, q2max);
346  cached_intJ9_mu = integrateJ(12, q2min, q2max);
347  cached_intJ4_mu = 0.;
348  cached_intJ5_mu = 0.;
349  cached_intJ7_mu = 0.;
350  cached_intJ8_mu = 0.;
351  // not needed at present
352  /*
353  cached_intJ4_mu = integrateJ(6,q2min,q2max);
354  cached_intJ5_mu = integrateJ(7,q2min,q2max);
355  cached_intJ7_mu = integrateJ(10,q2min,q2max);
356  cached_intJ8_mu = integrateJ(11,q2min,q2max);
357  */
358  checkcache_int_mu = 1;
359  }
360  if (lep == StandardModel::ELECTRON) {
361  cached_intJ1s_el = integrateJ(1, q2min, q2max);
362  cached_intJ1c_el = integrateJ(2, q2min, q2max);
363  cached_intJ2s_el = integrateJ(3, q2min, q2max);
364  cached_intJ2c_el = integrateJ(4, q2min, q2max);
365  cached_intJ3_el = integrateJ(5, q2min, q2max);
366  cached_intJ6s_el = integrateJ(8, q2min, q2max);
367  cached_intJ6c_el = integrateJ(9, q2min, q2max);
368  cached_intJ9_el = integrateJ(12, q2min, q2max);
369  cached_intJ4_el = 0.;
370  cached_intJ5_el = 0.;
371  cached_intJ7_el = 0.;
372  cached_intJ8_el = 0.;
373  // not needed at present
374  /*
375  cached_intJ4_el = integrateJ(6,q2min,q2max);
376  cached_intJ5_el = integrateJ(7,q2min,q2max);
377  cached_intJ7_el = integrateJ(10,q2min,q2max);
378  cached_intJ8_el = integrateJ(11,q2min,q2max);
379  */
380  checkcache_int_el = 1;
381  }
382  }
383 
384  hA1w1_cache = hA1w1;
385  rho2_cache = rho2;
386  R1w1_cache = R1w1;
387  R2w1_cache = R2w1;
388  N_A_cache = N_A;
389  N_1_cache =N_1;
390  N_2_cache = N_2;
391  j_A_cache= j_A;
392  j_0_cache = j_0;
393  j_1_cache = j_1;
394  j_2_cache = j_2;
395  k_A_cache = k_A;
396  k_0_cache = k_0;
397  k_1_cache = k_1;
398  k_2_cache = k_2;
399  l_A_cache = l_A;
400 
401  af0_cache = af0;
402  af1_cache = af1;
403  af2_cache = af2;
404  ag0_cache = ag0;
405  ag1_cache = ag1;
406  ag2_cache = ag2;
407  aF11_cache = aF11;
408  aF12_cache = aF12;
409  aF21_cache = aF21;
410  aF22_cache = aF22;
411 
412  CS_cache = CS;
413  CSp_cache = CSp;
414  CP_cache = CP;
415  CPp_cache = CPp;
416  CV_cache = CV;
417  CVp_cache = CVp;
418  CA_cache = CA;
419  CAp_cache = CAp;
420  CT_cache = CT;
421  CTp_cache = CTp;
422 
423  mySM.getFlavour().setUpdateFlag(meson, vectorM, lep, false);
424 
425  return;
426 
427 }
428 
429 /*******************************************************************************
430  * Kinematic functions *
431  * ****************************************************************************/
432 
433 double MVlnu::lambda_half(double a, double b, double c)
434 {
435  return sqrt(a*a+b*b+c*c-2.*(a*b+a*c+b*c));
436 }
437 /*******************************************************************************
438  * Form factors *
439  * ****************************************************************************/
440 
441 double MVlnu::phi_f(double z)
442 {
443  double prefac = 4. * (MV_o_MM) / MM / MM * sqrt(nI / (3. * M_PI * chiTA));
444  double num = (1. + z) * sqrt((1. - z)*(1. - z)*(1. - z));
445  double den = (1. + MV_o_MM)*(1. - z) + 2. * sqrtMV_o_MM * (1. + z);
446  double den4 = den * den * den*den;
447  return prefac * num / den4;
448 }
449 
450 double MVlnu::f_BGL(double q2)
451 {
452  double w = w0 - q2 / (2. * MM * MV);
453  double z = (sqrt(w + 1.) - M_SQRT2) / (sqrt(w + 1.) + M_SQRT2);
454  double Pfacf = (z - zA1) / (1. - z * zA1)*(z - zA2) / (1. - z * zA2)*(z - zA3) / (1. - z * zA3)*(z - zA4) / (1. - z * zA4);
455  double phif = phi_f(z);
456  return (af0 + af1 * z + af2 * z * z) / phif / Pfacf;
457 }
458 
459 double MVlnu::phi_g(double z)
460 {
461  double prefac = sqrt(nI / (3. * M_PI * chiTV));
462  double num = 16. * (MV_o_MM)*(MV_o_MM)*(1. + z)*(1. + z) / sqrt(1. - z);
463  double den = (1. + MV_o_MM)*(1. - z) + 2. * sqrtMV_o_MM * (1. + z);
464  double den4 = den * den * den*den;
465  return prefac * num / den4;
466 }
467 
468 double MVlnu::g_BGL(double q2)
469 {
470  double w = w0 - q2 / (2. * MM * MV);
471  double z = (sqrt(w + 1.) - M_SQRT2) / (sqrt(w + 1.) + M_SQRT2);
472  double Pfacg = (z - zV1) / (1. - z * zV1)*(z - zV2) / (1. - z * zV2)*(z - zV3) / (1. - z * zV3)*(z - zV4) / (1. - z * zV4);
473  double phig = phi_g(z);
474  return (ag0 + ag1 * z + ag2 * z * z) / phig / Pfacg;
475 }
476 
477 double MVlnu::phi_F1(double z)
478 {
479  double prefac = 4. * (MV_o_MM) / MM / MM / MM * sqrt(nI / (6. * M_PI * chiTA));
480  double num = (1. + z) * sqrt((1. - z)*(1. - z)*(1. - z)*(1. - z)*(1. - z));
481  double den = (1. + MV_o_MM)*(1. - z) + 2. * sqrtMV_o_MM * (1. + z);
482  double den5 = den * den * den * den*den;
483  return prefac * num / den5;
484 }
485 
486 double MVlnu::F1_BGL(double q2)
487 {
488  double w = w0 - q2 / (2. * MM * MV);
489  double z = (sqrt(w + 1.) - M_SQRT2) / (sqrt(w + 1.) + M_SQRT2);
490  double PfacF1 = (z - zA1) / (1. - z * zA1)*(z - zA2) / (1. - z * zA2)*(z - zA3) / (1. - z * zA3)*(z - zA4) / (1. - z * zA4);
491  double phiF1 = phi_F1(z);
492  double aF10 = (MM - MV)*(phi_F1(0.) / phi_f(0.)) * af0; // F1(z=0) = (MM-MV)*f(z=0)
493  return (aF10 + aF11 * z + aF12 * z * z) / phiF1 / PfacF1;
494 }
495 
496 double MVlnu::phi_F2(double z)
497 {
498  double prefac = 8. * sqrt(2.)*(MV_o_MM)*(MV_o_MM) * sqrt(nI / (M_PI * chiTP));
499  double num = (1. + z)*(1. + z) / sqrt(1. - z);
500  double den = (1. + MV_o_MM)*(1. - z) + 2. * sqrtMV_o_MM * (1. + z);
501  double den4 = den * den * den*den;
502  return prefac * num / den4;
503 }
504 
505 double MVlnu::F2_BGL(double q2)
506 {
507  double w = w0 - q2 / (2. * MM * MV);
508  double z = (sqrt(w + 1.) - M_SQRT2) / (sqrt(w + 1.) + M_SQRT2);
509  double z0 = (sqrt(w0 + 1.) - M_SQRT2) / (sqrt(w0 + 1.) + M_SQRT2);
510  double PfacF2 = (z - zP1) / (1. - z * zP1)*(z - zP2) / (1. - z * zP2)*(z - zP3) / (1. - z * zP3);
511  double PfacF2z0 = (z0 - zP1) / (1. - z0 * zP1)*(z0 - zP2) / (1. - z0 * zP2)*(z0 - zP3) / (1. - z0 * zP3);
512  double phiF2 = phi_F2(z);
513  double phiF2z0 = phi_F2(z0);
514  double aF20 = PfacF2z0 * phiF2z0 * 2. * F1_BGL(0.) / (MM * MM - MV * MV) - aF21 * z0 - aF22 * z0*z0; // F2(q2=0) = 2.*F1(q2=0)/(MM*MM-MV*MV)
515  return (aF20 + aF21 * z + aF22 * z * z) / phiF2 / PfacF2;
516 }
517 
518 double MVlnu::hA1(double q2)
519 {
520  double w = w0 - q2 / (2. * MM * MV);
521  double z = (sqrt(w + 1.) - M_SQRT2) / (sqrt(w + 1.) + M_SQRT2);
522  if (CLNflag) return hA1w1 * N_A * (1. - j_A * 8. * rho2 * z + k_A * (53. * rho2 - 15.) * z * z - l_A * (231. * rho2 - 91.) * z * z * z);
523  else return f_BGL(q2) / sqrt(MM * MV) / (1. + w);
524 }
525 
526 double MVlnu::R1(double q2)
527 {
528  double w = w0 - q2 / (2. * MM * MV);
529  return N_1 * R1w1 - j_1 * 0.12 * (w - 1.) + k_1 * 0.05 * (w - 1.)*(w - 1.);
530 }
531 
532 double MVlnu::R2(double q2)
533 {
534  double w = w0 - q2 / (2. * MM * MV);
535  return N_2 * R2w1 + j_2 * 0.11 * (w - 1.) - k_2 * 0.06 * (w - 1.)*(w - 1.);
536 }
537 
538 double MVlnu::R0(double q2)
539 {
540  double w = w0 - q2 / (2. * MM * MV);
541  /* form factor relation among A0, A1 and A2 at q2=0 */
542  double R2q2at0 = R2(0.);
543  double R0q2at0 = (MM + MV - (MM - MV) * R2q2at0) / (2. * MV);
544  // caveat: HQET rel at the kinematic endpoint, q2 = 0 ...
545  double R0w1 = R0q2at0 + j_0 * 0.11 * (w0 - 1.) - k_0 * 0.01 * (w0 - 1.)*(w0 - 1.);
546  // one may consider "lattice" R0w1 = 1.14 +- O(10%) + consistency rel at q2 = 0 ...
547  return R0w1 - j_0 * 0.11 * (w - 1.) + k_0 * 0.01 * (w - 1.)*(w - 1.);
548 }
549 
550 double MVlnu::V(double q2)
551 {
552  if (CLNflag) return R1(q2) / RV * hA1(q2);
553  else return (MM + MV) * g_BGL(q2) / 2.;
554 }
555 
556 double MVlnu::A0(double q2)
557 {
558  double w = w0 - q2 / (2. * MM * MV);
559  if (CLNflag) return R0(q2) / RV * hA1(q2);
560  else return F2_BGL(q2) / RV / (1. + w);
561 }
562 
563 double MVlnu::A1(double q2)
564 {
565  double w = w0 - q2 / (2. * MM * MV);
566  return (w + 1.)*RV / 2. * hA1(q2);
567 }
568 
569 double MVlnu::A2(double q2)
570 {
571  double w = w0 - q2 / (2. * MM * MV);
572  if (CLNflag) return R2(q2) / RV * hA1(q2);
573  else return (MM + MV) / 2. / (w * w - 1.) / MM / MV * ((w - MV_o_MM) * f_BGL(q2) - F1_BGL(q2) / MM);
574 }
575 
576 double MVlnu::T1(double q2)
577 {
578  double delta_T1 = 0.;
579  return (Mb + Mc) / (MM + MV) * V(q2)*(1. + delta_T1);
580 }
581 
582 double MVlnu::T2(double q2)
583 {
584  double delta_T2 = 0.;
585  return (Mb - Mc) / (MM - MV) * A1(q2)*(1. + delta_T2);
586 }
587 
588 double MVlnu::A12(double q2)
589 {
590  return (A1(q2)*(MM + MV)*(MM + MV)*(MM * MM - MV * MV - q2) -
591  A2(q2)*(MM * MM * MM * MM + (MV * MV - q2)*(MV * MV - q2) -
592  2. * MM * MM * (MV * MV + q2))) / (16. * MM * MV * MV * (MM + MV));
593 }
594 
595 double MVlnu::T23(double q2)
596 {
597  double delta_T23 = 0.;
598  return ((Mb - Mc)*((MM - MV)*(MM - MV) - q2)*((MM + MV)*(MM + MV) - q2) * A0(q2) +
599  8 * MM * MV * (MV * MV - MM * MM) * A12(q2)) / (4. * MM * (MV - MM) * MV * q2)*(1. + delta_T23);
600 }
601 /********************************************************************************
602  * Helicity amplitudes (normalization such that all H \propto (mass scale)^-1) *
603  * *****************************************************************************/
604 
605 gslpp::complex MVlnu::HV0(double q2)
606 {
607  return 4. * gslpp::complex::i() * MM * MV / (sqrt(q2)*(MM + MV))*((CV - CVp)*(MM + MV) * A12(q2) + Mb * (C7 - C7p) * T23(q2));
608 }
609 
610 gslpp::complex MVlnu::HVp(double q2)
611 {
612  return gslpp::complex::i()*((((CV + CVp) * lambda_half(MM*MM, MV*MV, q2) * V(q2)-(MM + MV)*(MM + MV)*(CV - CVp) * A1(q2))) / (2. * (MM + MV))
613  + (Mb / q2)*((C7 + C7p) * lambda_half(MM*MM, MV*MV, q2) * T1(q2)-(C7 - C7p)*(MM * MM - MV * MV) * T2(q2)));
614 }
615 
616 gslpp::complex MVlnu::HVm(double q2)
617 {
618  return gslpp::complex::i()*(((-(CV + CVp) * lambda_half(MM*MM, MV*MV, q2) * V(q2)-(MM + MV)*(MM + MV)*(CV - CVp) * A1(q2))) / (2. * (MM + MV))
619  + (Mb / q2)*(-(C7 + C7p) * lambda_half(MM*MM, MV*MV, q2) * T1(q2)-(C7 - C7p)*(MM * MM - MV * MV) * T2(q2)));
620 }
621 
622 gslpp::complex MVlnu::HAp(double q2)
623 {
624  return gslpp::complex::i()*((CA + CAp) * lambda_half(MM*MM, MV*MV, q2) * V(q2)-(MM + MV)*(MM + MV)*(CA - CAp) * A1(q2)) / (2. * (MM + MV));
625 }
626 
627 gslpp::complex MVlnu::HAm(double q2)
628 {
629  return gslpp::complex::i()*(-(CA + CAp) * lambda_half(MM*MM, MV*MV, q2) * V(q2)-(MM + MV)*(MM + MV)*(CA - CAp) * A1(q2)) / (2. * (MM + MV));
630 }
631 
632 gslpp::complex MVlnu::HA0(double q2)
633 {
634  return 4. * gslpp::complex::i() * MV * MM / (sqrt(q2))*(CA - CAp) * A12(q2);
635 }
636 
637 gslpp::complex MVlnu::HP(double q2)
638 {
639  return gslpp::complex::i() * lambda_half(MM*MM, MV*MV, q2) / 2. * ((CP - CPp) / (Mb + Mc)+(Mlep + Mnu) / q2 * (CA - CAp)) * A0(q2);
640 }
641 
642 gslpp::complex MVlnu::HS(double q2)
643 {
644  return gslpp::complex::i() * lambda_half(MM*MM, MV*MV, q2) / 2. * ((CS - CSp) / (Mb + Mc)+(Mlep - Mnu) / q2 * (CV - CVp)) * A0(q2);
645 }
646 
647 gslpp::complex MVlnu::HT0(double q2)
648 {
649  return 2. * M_SQRT2 * (MM * MV) / (MM + MV)*(CT + CTp) * T23(q2);
650 }
651 
652 gslpp::complex MVlnu::HT0t(double q2)
653 {
654  return 2. * (MM * MV) / (MM + MV)*(CT - CTp) * T23(q2);
655 }
656 
657 gslpp::complex MVlnu::HTp(double q2)
658 {
659  return ((CT - CTp) * lambda_half(MM*MM, MV*MV, q2) * T1(q2)-(CT + CTp)*(MM * MM - MV * MV) * T2(q2)) / (sqrt(2. * q2));
660 }
661 
662 gslpp::complex MVlnu::HTpt(double q2)
663 {
664  return ((CT + CTp) * lambda_half(MM*MM, MV*MV, q2) * T1(q2)-(CT - CTp)*(MM * MM - MV * MV) * T2(q2)) / (sqrt(2. * q2));
665 }
666 
667 gslpp::complex MVlnu::HTm(double q2)
668 {
669  return (-(CT - CTp) * lambda_half(MM*MM, MV*MV, q2) * T1(q2)-(CT + CTp)*(MM * MM - MV * MV) * T2(q2)) / (sqrt(2. * q2));
670 }
671 
672 gslpp::complex MVlnu::HTmt(double q2)
673 {
674  return (-(CT + CTp) * lambda_half(MM*MM, MV*MV, q2) * T1(q2)-(CT - CTp)*(MM * MM - MV * MV) * T2(q2)) / (sqrt(2. * q2));
675 }
676 /*******************************************************************************
677  * Generalized angular coefficients (see 1506.03970) *
678  * ****************************************************************************/
679 
680 gslpp::complex MVlnu::G000(double q2)
681 {
682  double lambda_MM = lambda_half(MM*MM, MV*MV, q2);
683  double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
684  double lambda_lep2 = lambda_lep*lambda_lep;
685  double Elep = sqrt(Mlep * Mlep + lambda_lep2 / (4. * q2));
686  double Enu = sqrt(Mnu * Mnu + lambda_lep2 / (4. * q2));
687  double Gprefactor = lambda_MM * lambda_lep / q2;
688 
689  return Gprefactor * (4. / 9. * (3. * Elep * Enu + lambda_lep2 / (4. * q2))*(HVp(q2).abs2() + HVm(q2).abs2() + HV0(q2).abs2() + HAp(q2).abs2() + HAm(q2).abs2() + HA0(q2).abs2()) +
690  4. * Mlep * Mnu / 3. * (HVp(q2).abs2() + HVm(q2).abs2() + HV0(q2).abs2() - HAp(q2).abs2() - HAm(q2).abs2() - HA0(q2).abs2()) +
691  4. / 3. * ((Elep * Enu - Mlep * Mnu + lambda_lep2 / (4. * q2)) * HS(q2).abs2()+(Elep * Enu + Mlep * Mnu + lambda_lep2 / (4. * q2)) * HP(q2).abs2()) +
692  16. / 9. * (3. * (Elep * Enu + Mlep * Mnu) - lambda_lep2 / (4. * q2))*(HTpt(q2).abs2() + HTmt(q2).abs2() + HT0t(q2).abs2()) +
693  8. / 9. * (3. * (Elep * Enu - Mlep * Mnu) - lambda_lep2 / (4. * q2))*(HTpt(q2).abs2() + HTmt(q2).abs2() + HT0t(q2).abs2()) +
694  16. / 3. * (Mlep * Enu + Mnu * Elep)*(HVp(q2) * HTpt(q2).conjugate() + HVm(q2) * HTmt(q2).conjugate() + HV0(q2) * HT0t(q2).conjugate()).imag() +
695  8. * M_SQRT2 / 3. * (Mlep * Enu - Mnu * Elep)*(HAp(q2) * HTp(q2).conjugate() + HAm(q2) * HTm(q2).conjugate() + HA0(q2) * HT0(q2).conjugate()).imag());
696 }
697 
698 gslpp::complex MVlnu::G010(double q2)
699 {
700  double lambda_MM = lambda_half(MM*MM, MV*MV, q2);
701  double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
702  double Gprefactor = lambda_MM * lambda_lep / q2;
703 
704  return Gprefactor * 4. / 3. * lambda_lep * ((HVp(q2) * HAp(q2).conjugate() - HVm(q2) * HAm(q2).conjugate()).real()+
705  (2. * M_SQRT2) / q2 * (Mlep * Mlep - Mnu * Mnu)*(HTp(q2) * HTpt(q2).conjugate() - HTm(q2) * HTmt(q2).conjugate()).real() +
706  2. * (Mlep + Mnu) / sqrt(q2)*(HAp(q2) * HTpt(q2).conjugate() - HAm(q2) * HTmt(q2).conjugate()).imag() +
707  M_SQRT2 * (Mlep - Mnu) / sqrt(q2)*(HVp(q2) * HTp(q2).conjugate() - HVm(q2) * HTm(q2).conjugate()).imag()-
708  (Mlep - Mnu) / sqrt(q2)*(HA0(q2) * HP(q2).conjugate()).real()-(Mlep + Mnu) / sqrt(q2)*(HV0(q2) * HS(q2).conjugate()).real()+
709  (M_SQRT2 * HT0(q2) * HP(q2).conjugate() + 2. * HT0t(q2) * HS(q2).conjugate()).imag());
710 
711 }
712 
713 gslpp::complex MVlnu::G020(double q2)
714 {
715  double lambda_MM = lambda_half(MM*MM, MV*MV, q2);
716  double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
717  double lambda_lep2 = lambda_lep*lambda_lep;
718  double Gprefactor = lambda_MM * lambda_lep / q2;
719 
720  return -Gprefactor * 2. / 9. * lambda_lep2 / q2 * ((-HVp(q2).abs2() - HVm(q2).abs2() + 2. * HV0(q2).abs2() - HAp(q2).abs2() - HAm(q2).abs2() + 2. * HA0(q2).abs2()) -
721  2. * (2. * HT0(q2).abs2() - HTp(q2).abs2() - HTm(q2).abs2()) - 4. * (2. * HT0t(q2).abs2() - HTpt(q2).abs2() - HTmt(q2).abs2()));
722 }
723 
724 gslpp::complex MVlnu::G200(double q2)
725 {
726  double lambda_MM = lambda_half(MM*MM, MV*MV, q2);
727  double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
728  double lambda_lep2 = lambda_lep*lambda_lep;
729  double Elep = sqrt(Mlep * Mlep + lambda_lep2 / (4. * q2));
730  double Enu = sqrt(Mnu * Mnu + lambda_lep2 / (4. * q2));
731  double Gprefactor = lambda_MM * lambda_lep / q2;
732 
733  return Gprefactor * (-4. / 9. * (3. * Elep * Enu + lambda_lep2 / (4. * q2))*(HVp(q2).abs2() + HVm(q2).abs2() - 2. * HV0(q2).abs2() + HAp(q2).abs2() + HAm(q2).abs2() - 2. * HA0(q2).abs2()) -
734  4. / 3. * Mlep * Mnu * (HVp(q2).abs2() + HVm(q2).abs2() - 2. * HV0(q2).abs2() - HAp(q2).abs2() - HAm(q2).abs2() + 2. * HA0(q2).abs2()) +
735  8. / 3. * (Elep * Enu - Mlep * Mnu + lambda_lep2 / (4. * q2)) * HS(q2).abs2() + 8. / 3. * (Elep * Enu + Mlep * Mnu + lambda_lep2 / (4. * q2)) * HP(q2).abs2() -
736  16. / 9. * (3. * (Elep * Enu + Mlep * Mnu) - lambda_lep2 / (4. * q2))*(HTpt(q2).abs2() + HTmt(q2).abs2() - 2. * HT0t(q2).abs2()) - 8. / 9. * (3. * (Elep * Enu - Mlep * Mnu) - lambda_lep2 / (4. * q2))*
737  (HTp(q2).abs2() + HTm(q2).abs2() - 2 * HT0(q2).abs2()) - 16. / 3. * (Mlep * Enu + Mnu * Elep)*(HVp(q2) * HTpt(q2).conjugate() + HVm(q2) * HTmt(q2).conjugate() - 2. * HV0(q2) * HT0t(q2).conjugate()).imag() -
738  8. * M_SQRT2 / 3. * (Mlep * Enu - Mnu * Elep)*(HAp(q2) * HTp(q2).conjugate() + HAm(q2) * HTm(q2).conjugate() - 2. * HA0(q2) * HT0(q2).conjugate()).imag());
739 }
740 
741 gslpp::complex MVlnu::G210(double q2)
742 {
743  double lambda_MM = lambda_half(MM*MM, MV*MV, q2);
744  double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
745  double Gprefactor = lambda_MM * lambda_lep / q2;
746 
747  return -Gprefactor * 4. * lambda_lep / 3. * ((HVp(q2) * HAp(q2).conjugate() - HVm(q2) * HAm(q2).conjugate()).real() +
748  2. * M_SQRT2 * (Mlep * Mlep - Mnu * Mnu) / q2 * (HTp(q2) * HTpt(q2).conjugate() - HTm(q2) * HTmt(q2).conjugate()).real() +
749  2. * (Mlep + Mnu) / sqrt(q2)*(HAp(q2) * HTpt(q2).conjugate() - HAm(q2) * HTmt(q2).conjugate()).imag() +
750  M_SQRT2 * (Mlep - Mnu) / sqrt(q2)*(HVp(q2) * HTp(q2).conjugate() - HVm(q2) * HTm(q2).conjugate()).imag() +
751  2. * (Mlep - Mnu) / sqrt(q2)*(HA0(q2) * HP(q2).conjugate()).real() + 2. * (Mlep + Mnu) / sqrt(q2)*(HV0(q2) * HS(q2).conjugate()).real() -
752  2. * (M_SQRT2 * HT0(q2) * HP(q2).conjugate() + 2. * HT0t(q2) * HS(q2).conjugate()).imag());
753 }
754 
755 gslpp::complex MVlnu::G220(double q2)
756 {
757  double lambda_MM = lambda_half(MM*MM, MV*MV, q2);
758  double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
759  double lambda_lep2 = lambda_lep*lambda_lep;
760  double Gprefactor = lambda_MM * lambda_lep / q2;
761 
762  return -Gprefactor * 2. / 9. * lambda_lep2 / q2 * (HVp(q2).abs2() + HVm(q2).abs2() + 4. * HV0(q2).abs2() + HAp(q2).abs2() +
763  HAm(q2).abs2() + 4. * HA0(q2).abs2() - 2. * (HTp(q2).abs2() + HTm(q2).abs2() + 4. * HT0(q2).abs2()) -
764  4. * (HTpt(q2).abs2() + HTmt(q2).abs2() + 4. * HT0t(q2).abs2()));
765 }
766 
767 gslpp::complex MVlnu::G211(double q2)
768 {
769  double lambda_MM = lambda_half(MM*MM, MV*MV, q2);
770  double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
771  double Gprefactor = lambda_MM * lambda_lep / q2;
772 
773  return Gprefactor * 4. / sqrt(3.) * lambda_lep * (HVp(q2) * HA0(q2).conjugate() + HAp(q2) * HV0(q2).conjugate() -
774  HV0(q2) * HAm(q2).conjugate() - HA0(q2) * HVm(q2).conjugate()+(Mlep + Mnu) / sqrt(q2)*(HVp(q2) * HS(q2).conjugate() + HS(q2) * HVm(q2).conjugate()) -
775  gslpp::complex::i() * M_SQRT2 * (HP(q2) * HTm(q2).conjugate() - HTp(q2) * HP(q2).conjugate() + M_SQRT2 * (HS(q2) * HTmt(q2).conjugate() - HTpt(q2) * HS(q2).conjugate()))+
776  (Mlep - Mnu) / sqrt(q2)*(HAp(q2) * HP(q2).conjugate() + HP(q2) * HAm(q2).conjugate()) -
777  gslpp::complex::i()*2. * (Mlep + Mnu) / sqrt(q2)*(HAp(q2) * HT0t(q2).conjugate() + HT0t(q2) * HAm(q2).conjugate() - HTpt(q2) * HA0(q2).conjugate() - HA0(q2) * HTmt(q2).conjugate()) -
778  gslpp::complex::i() * M_SQRT2 * (Mlep - Mnu) / sqrt(q2)*(HVp(q2) * HT0(q2).conjugate() + HT0(q2) * HVm(q2).conjugate() - HTp(q2) * HV0(q2).conjugate() - HV0(q2) * HTm(q2).conjugate()) +
779  2. * M_SQRT2 * (Mlep * Mlep - Mnu * Mnu) / q2 * (HTp(q2) * HT0t(q2).conjugate() + HTpt(q2) * HT0(q2).conjugate() - HT0(q2) * HTmt(q2).conjugate() - HT0t(q2) * HTm(q2).conjugate()));
780 }
781 
782 gslpp::complex MVlnu::G221(double q2)
783 {
784  double lambda_MM = lambda_half(MM*MM, MV*MV, q2);
785  double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
786  double lambda_lep2 = lambda_lep*lambda_lep;
787  double Gprefactor = lambda_MM * lambda_lep / q2;
788 
789  return Gprefactor * 4. / 3. * lambda_lep2 / q2 * (HVp(q2) * HV0(q2).conjugate() + HV0(q2) * HVm(q2).conjugate() +
790  HAp(q2) * HA0(q2).conjugate() + HA0(q2) * HAm(q2).conjugate() - 2. * (HTp(q2) * HT0(q2).conjugate() +
791  HT0(q2) * HTm(q2).conjugate() + 2. * (HTpt(q2) * HT0t(q2).conjugate() + HT0t(q2) * HTmt(q2).conjugate())));
792 }
793 
794 gslpp::complex MVlnu::G222(double q2)
795 {
796  double lambda_MM = lambda_half(MM*MM, MV*MV, q2);
797  double lambda_lep = lambda_half(Mlep*Mlep, Mnu*Mnu, q2);
798  double lambda_lep2 = lambda_lep*lambda_lep;
799  double Gprefactor = lambda_MM * lambda_lep / q2;
800 
801  return -Gprefactor * 8. / 3. * lambda_lep2 / q2 * (HVp(q2) * HVm(q2).conjugate() + HAp(q2) * HAm(q2).conjugate() -
802  2. * (HTp(q2) * HTm(q2).conjugate() + 2. * HTpt(q2) * HTmt(q2).conjugate()));
803 }
804 /***************************************************************************
805  * 12 independent J angular coefficients (see again 1506.03970) *
806  * ************************************************************************/
807 
808 double MVlnu::J1s(double q2)
809 {
810  if ((q2 < Mlep * Mlep) or (q2 > (MM - MV)*(MM - MV))) return 0.;
811  return amplsq_factor * (8. * G000(q2) + 2. * G020(q2) - 4. * G200(q2) - G220(q2)).real() / 3.;
812 }
813 
814 double MVlnu::J1c(double q2)
815 {
816  if ((q2 < Mlep * Mlep) or (q2 > (MM - MV)*(MM - MV))) return 0.;
817  return amplsq_factor * (8. * G000(q2) + 2. * G020(q2) + 8. * G200(q2) + 2. * G220(q2)).real() / 3.;
818 }
819 
820 double MVlnu::J2s(double q2)
821 {
822  if ((q2 < Mlep * Mlep) or (q2 > (MM - MV)*(MM - MV))) return 0.;
823  return amplsq_factor * (2. * G020(q2) - G220(q2)).real();
824 }
825 
826 double MVlnu::J2c(double q2)
827 {
828  if ((q2 < Mlep * Mlep) or (q2 > (MM - MV)*(MM - MV))) return 0.;
829  return amplsq_factor * (2. * (G020(q2) + G220(q2))).real();
830 }
831 
832 double MVlnu::J3(double q2)
833 {
834  if ((q2 < Mlep * Mlep) or (q2 > (MM - MV)*(MM - MV))) return 0.;
835  return amplsq_factor * (G222(q2).real());
836 }
837 
838 double MVlnu::J4(double q2)
839 {
840  if ((q2 < Mlep * Mlep) or (q2 > (MM - MV)*(MM - MV))) return 0.;
841  return -amplsq_factor * (G221(q2).real());
842 }
843 
844 double MVlnu::J5(double q2)
845 {
846  if ((q2 < Mlep * Mlep) or (q2 > (MM - MV)*(MM - MV))) return 0.;
847  return amplsq_factor * (2. * G211(q2).real() / sqrt(3.));
848 }
849 
850 double MVlnu::J6s(double q2)
851 {
852  if ((q2 < Mlep * Mlep) or (q2 > (MM - MV)*(MM - MV))) return 0.;
853  return -amplsq_factor * (4. * (2. * G010(q2) - G210(q2)).real() / 3.);
854 }
855 
856 double MVlnu::J6c(double q2)
857 {
858  if (q2 < Mlep * Mlep) return 0.;
859  if (q2 > (MM - MV)*(MM - MV)) return 0.;
860  return -amplsq_factor * (8. * (G010(q2) + G210(q2)).real() / 3.);
861 }
862 
863 double MVlnu::J7(double q2)
864 {
865  if ((q2 < Mlep * Mlep) or (q2 > (MM - MV)*(MM - MV))) return 0.;
866  return -amplsq_factor * (2. * sqrt(3.)*(G211(q2).imag()) / 3.);
867 }
868 
869 double MVlnu::J8(double q2)
870 {
871  if ((q2 < Mlep * Mlep) or (q2 > (MM - MV)*(MM - MV))) return 0.;
872  return amplsq_factor * (G221(q2).imag());
873 }
874 
875 double MVlnu::J9(double q2)
876 {
877  if ((q2 < Mlep * Mlep) or (q2 > (MM - MV)*(MM - MV))) return 0.;
878  return -amplsq_factor * (G222(q2).imag());
879 }
880 /***************************************************************************
881  * Integration of angular coefficients Js *
882  * ************************************************************************/
883 
884 double MVlnu::integrateJ(int i, double q2_min, double q2_max)
885 {
886  old_handler = gsl_set_error_handler_off();
887 
888  switch (i) {
889  case 1:
890  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ1s_tau;
891  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ1s_mu;
892  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ1s_el;
893  FJ = convertToGslFunction(boost::bind(&MVlnu::J1s, &(*this), _1));
894  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();
895  gsl_set_error_handler(old_handler);
896  return J_res;
897  break;
898  case 2:
899  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ1c_tau;
900  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ1c_mu;
901  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ1c_el;
902  FJ = convertToGslFunction(boost::bind(&MVlnu::J1c, &(*this), _1));
903  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();
904  gsl_set_error_handler(old_handler);
905  return J_res;
906  break;
907  case 3:
908  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ2s_tau;
909  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ2s_mu;
910  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ2s_el;
911  FJ = convertToGslFunction(boost::bind(&MVlnu::J2s, &(*this), _1));
912  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();
913  gsl_set_error_handler(old_handler);
914  return J_res;
915  break;
916  case 4:
917  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ2c_tau;
918  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ2c_mu;
919  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ2c_el;
920  FJ = convertToGslFunction(boost::bind(&MVlnu::J2c, &(*this), _1));
921  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();
922  gsl_set_error_handler(old_handler);
923  return J_res;
924  break;
925  case 5:
926  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ3_tau;
927  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ3_mu;
928  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ3_el;
929  FJ = convertToGslFunction(boost::bind(&MVlnu::J3, &(*this), _1));
930  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();
931  gsl_set_error_handler(old_handler);
932  gsl_set_error_handler(old_handler);
933  return J_res;
934  break;
935  case 6:
936  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ4_tau;
937  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ4_mu;
938  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ4_el;
939  FJ = convertToGslFunction(boost::bind(&MVlnu::J4, &(*this), _1));
940  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();
941  gsl_set_error_handler(old_handler);
942  return J_res;
943  break;
944  case 7:
945  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ5_tau;
946  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ5_mu;
947  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ5_el;
948  FJ = convertToGslFunction(boost::bind(&MVlnu::J5, &(*this), _1));
949  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();
950  gsl_set_error_handler(old_handler);
951  return J_res;
952  break;
953  case 8:
954  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ6s_tau;
955  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ6s_mu;
956  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ6s_el;
957  FJ = convertToGslFunction(boost::bind(&MVlnu::J6s, &(*this), _1));
958  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();
959  gsl_set_error_handler(old_handler);
960  return J_res;
961  break;
962  case 9:
963  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ6c_mu;
964  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ6c_mu;
965  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ6c_el;
966  FJ = convertToGslFunction(boost::bind(&MVlnu::J6c, &(*this), _1));
967  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();
968  return J_res;
969  break;
970  case 10:
971  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ7_tau;
972  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ7_mu;
973  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ7_el;
974  FJ = convertToGslFunction(boost::bind(&MVlnu::J7, &(*this), _1));
975  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();
976  gsl_set_error_handler(old_handler);
977  return J_res;
978  break;
979  case 11:
980  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ8_tau;
981  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ8_mu;
982  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ8_el;
983  FJ = convertToGslFunction(boost::bind(&MVlnu::J8, &(*this), _1));
984  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();
985  gsl_set_error_handler(old_handler);
986  return J_res;
987  break;
988  case 12:
989  if (lep == StandardModel::TAU) if ((checkcache_int_tau == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ9_tau;
990  if (lep == StandardModel::MU) if ((checkcache_int_mu == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ9_mu;
991  if (lep == StandardModel::ELECTRON) if ((checkcache_int_el == 1) && (q2_min == q2min) && (q2_max == q2max)) return cached_intJ9_el;
992  FJ = convertToGslFunction(boost::bind(&MVlnu::J9, &(*this), _1));
993  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();
994  gsl_set_error_handler(old_handler);
995  return J_res;
996  break;
997  default:
998  gsl_set_error_handler(old_handler);
999  std::stringstream out;
1000  out << i;
1001  throw std::runtime_error("MVlnu::integrateJ: index " + out.str() + " not implemented");
1002  }
1003 }
1004 
1005 double MVlnu::dGammadw(double q2)
1006 {
1007  updateParameters();
1008 
1009  return 3. / 4. * (2. * J1s(q2) + J1c(q2)) - 1. / 4. * (2. * J2s(q2) + J2c(q2));
1010 }
1011 
1012 double MVlnu::getDeltaGammaDeltaw(double w_min, double w_max)
1013 {
1014  updateParameters();
1015 
1016  double q2_min = (2. * MM * MV)*(w0 - w_max); // min is Mlep*Mlep;
1017  double q2_max = (2. * MM * MV)*(w0 - w_min); // max is (MM-MV)*(MM-MV);
1018 
1019  double intJ1s = integrateJ(1, q2_min, q2_max);
1020  double intJ1c = integrateJ(2, q2_min, q2_max);
1021  double intJ2s = integrateJ(3, q2_min, q2_max);
1022  double intJ2c = integrateJ(4, q2_min, q2_max);
1023 
1024  return 3. / 4. * (2. * intJ1s + intJ1c) - 1. / 4. * (2. * intJ2s + intJ2c);
1025 }
1026 
1027 double MVlnu::dGammadcldq2(double q2, double cl)
1028 {
1029  updateParameters();
1030 
1031  return 3. / 8. * ((J1s(q2) + 2. * J1c(q2)) + cl * (J6s(q2) + 2. * J6c(q2))+(2. * cl * cl - 1.)*(J2s(q2) + 2. * J2c(q2)));
1032 }
1033 
1034 double MVlnu::dGammadcl(double cl)
1035 {
1036  updateParameters();
1037 
1038  double intJ1s = integrateJ(1, q2min, q2max);
1039  double intJ1c = integrateJ(2, q2min, q2max);
1040  double intJ2s = integrateJ(3, q2min, q2max);
1041  double intJ2c = integrateJ(4, q2min, q2max);
1042  double intJ6s = integrateJ(8, q2min, q2max);
1043  double intJ6c = integrateJ(9, q2min, q2max);
1044 
1045  return 3. / 8. * ((intJ1s + 2. * intJ1c) + cl * (intJ6s + 2. * intJ6c)+(2. * cl * cl - 1.)*(intJ2s + 2. * intJ2c));
1046 }
1047 
1048 double MVlnu::getDeltaGammaDeltacl(double cl_min, double cl_max)
1049 {
1050  updateParameters();
1051 
1052  double intJ1s = integrateJ(1, q2min, q2max);
1053  double intJ1c = integrateJ(2, q2min, q2max);
1054  double intJ2s = integrateJ(3, q2min, q2max);
1055  double intJ2c = integrateJ(4, q2min, q2max);
1056  double intJ6s = integrateJ(8, q2min, q2max);
1057  double intJ6c = integrateJ(9, q2min, q2max);
1058 
1059  return 3. / 8. * ((cl_max - cl_min)*(intJ1c + 2. * intJ1s)+
1060  (cl_max * cl_max - cl_min * cl_min) / 2. * (intJ6c + 2. * intJ6s)+
1061  (2. / 3. * (cl_max * cl_max * cl_max - cl_min * cl_min * cl_min)-(cl_max - cl_min))*(intJ2c + 2. * intJ2s));
1062 }
1063 
1064 double MVlnu::dGammadcVdq2(double q2, double cl)
1065 {
1066  updateParameters();
1067 
1068  return 3. / 8. * ((J1s(q2) + 2. * J1c(q2)) + cl * (J6s(q2) + 2. * J6c(q2))+(2. * cl * cl - 1.)*(J2s(q2) + 2. * J2c(q2)));
1069 }
1070 
1071 double MVlnu::dGammadcV(double cV)
1072 {
1073  updateParameters();
1074 
1075  double intJ1s = integrateJ(1, q2min, q2max);
1076  double intJ1c = integrateJ(2, q2min, q2max);
1077  double intJ2s = integrateJ(3, q2min, q2max);
1078  double intJ2c = integrateJ(4, q2min, q2max);
1079 
1080  return 3. / 8. * (cV * cV * (3. * intJ1c - intJ2c)+(1. - cV * cV)*(3. * intJ1s - intJ2s));
1081 }
1082 
1083 double MVlnu::getDeltaGammaDeltacV(double cV_min, double cV_max)
1084 {
1085  updateParameters();
1086 
1087  double intJ1s = integrateJ(1, q2min, q2max);
1088  double intJ1c = integrateJ(2, q2min, q2max);
1089  double intJ2s = integrateJ(3, q2min, q2max);
1090  double intJ2c = integrateJ(4, q2min, q2max);
1091 
1092  return 3. / 8. * ((cV_max * cV_max * cV_max - cV_min * cV_min * cV_min) / 3. * (3. * intJ1c - intJ2c)+
1093  ((cV_max - cV_min)-(cV_max * cV_max * cV_max - cV_min * cV_min * cV_min) / 3.)*(3. * intJ1s - intJ2s));
1094 }
1095 
1096 double MVlnu::dGammadchidq2(double q2, double chi)
1097 {
1098  updateParameters();
1099 
1100  return (3. * J1c(q2) + 6. * J1s(q2) - J2c(q2) - 2. * J2s(q2)) / 8. / M_PI +
1101  cos(2. * chi) / 2. / M_PI * J3(q2) + sin(2. * chi) / 2. / M_PI * J9(q2);
1102 }
1103 
1104 double MVlnu::dGammadchi(double chi)
1105 {
1106  updateParameters();
1107 
1108  double intJ1s = integrateJ(1, q2min, q2max);
1109  double intJ1c = integrateJ(2, q2min, q2max);
1110  double intJ2s = integrateJ(3, q2min, q2max);
1111  double intJ2c = integrateJ(4, q2min, q2max);
1112  double intJ3 = integrateJ(5, q2min, q2max);
1113  double intJ9 = integrateJ(12, q2min, q2max);
1114 
1115  return ((3. * intJ1c + 6. * intJ1s - intJ2c - 2. * intJ2s) / 4. +
1116  cos(2. * chi) * intJ3 + sin(2. * chi) * intJ9) / 2. / M_PI;
1117 }
1118 
1119 double MVlnu::getDeltaGammaDeltachi(double chi_min, double chi_max)
1120 {
1121  updateParameters();
1122 
1123  double intJ1s = integrateJ(1, q2min, q2max);
1124  double intJ1c = integrateJ(2, q2min, q2max);
1125  double intJ2s = integrateJ(3, q2min, q2max);
1126  double intJ2c = integrateJ(4, q2min, q2max);
1127  double intJ3 = integrateJ(5, q2min, q2max);
1128  double intJ9 = integrateJ(12, q2min, q2max);
1129 
1130  return ((chi_max - chi_min)*(3. * intJ1c + 6. * intJ1s - intJ2c - 2. * intJ2s) / 4. +
1131  (sin(2. * chi_max) - sin(2. * chi_min)) / 2. * intJ3 -
1132  (cos(2. * chi_max) - cos(2. * chi_min)) / 2. * intJ9) / (2. * M_PI);
1133 }
1134 
1135 double MVlnu::getFL()
1136 {
1137  updateParameters();
1138 
1139  double intJ1s = integrateJ(1, q2min, q2max);
1140  double intJ1c = integrateJ(2, q2min, q2max);
1141  double intJ2s = integrateJ(3, q2min, q2max);
1142  double intJ2c = integrateJ(4, q2min, q2max);
1143 
1144  double DeltaJL = (3. * intJ1c - intJ2c) / 4.;
1145  double DeltaJ = 3. / 4. * (2. * intJ1s + intJ1c) - 1. / 4. * (2. * intJ2s + intJ2c);
1146  return DeltaJL/DeltaJ;
1147 
1148 }
1149 
1150 double MVlnu::get_unitarity_V_BGL()
1151 {
1152  updateParameters();
1153 
1154  return ag0 * ag0 + ag1 * ag1 + ag2*ag2;
1155 
1156 }
1157 
1158 double MVlnu::get_unitarity_A_BGL()
1159 {
1160  updateParameters();
1161 
1162  double aF10 = (MM - MV)*(phi_F1(0.) / phi_f(0.)) * af0;
1163  return af0 * af0 + af1 * af1 + af2 * af2 + aF10 * aF10 + aF11 * aF11 + aF12*aF12;
1164 }
1165 
1166 double MVlnu::get_unitarity_P_BGL()
1167 {
1168  updateParameters();
1169 
1170  double z0 = (sqrt(w0 + 1.) - M_SQRT2) / (sqrt(w0 + 1.) + M_SQRT2);
1171  double PfacF2z0 = (z0 - zP1) / (1. - z0 * zP1)*(z0 - zP2) / (1. - z0 * zP2)*(z0 - zP3) / (1. - z0 * zP3);
1172  double phiF2z0 = phi_F2(z0);
1173  double aF20 = PfacF2z0 * phiF2z0 * 2. * F1_BGL(0.) / (MM * MM - MV * MV) - aF21 * z0 - aF22 * z0*z0;
1174  return aF20 * aF20 + aF21 * aF21 + aF22*aF22;
1175 }
1176 
1177 double MVlnu::get_hA1w1()
1178 {
1179  updateParameters();
1180 
1181  return hA1(q2max);
1182 }
1183 
1184 double MVlnu::get_hA1(double w)
1185 {
1186  updateParameters();
1187  double q2 = (2. * MM * MV)*(w0 - w);
1188 
1189  return hA1(q2);
1190 }
1191 
1192 double MVlnu::get_R1(double w)
1193 {
1194  updateParameters();
1195  double q2 = (2. * MM * MV)*(w0 - w);
1196 
1197  if (CLNflag) return R1(q2);
1198  else return V(q2) * RV / hA1(q2);
1199 }
1200 
1201 double MVlnu::get_R2(double w)
1202 {
1203  updateParameters();
1204  double q2 = (2. * MM * MV)*(w0 - w);
1205 
1206  if (CLNflag) return R2(q2);
1207  else return A2(q2) * RV / hA1(q2);
1208 }
1209 
1210 double MVlnu::get_R0(double w)
1211 {
1212  updateParameters();
1213  double q2 = (2. * MM * MV)*(w0 - w);
1214 
1215  if (CLNflag) return R0(q2);
1216  else return A0(q2) * RV / hA1(q2);
1217 }
1218 
1219 
1220 /***************************************************************************
1221  * SM computation ... lep hel asymmetry, see 1203.2654, 1707.09509 *
1222  * ************************************************************************/
1223 
1224 double MVlnu::Hplus(double q2){
1225  double abs_p = lambda_half(MM*MM,MV*MV,q2);
1226  return (MM+MV)*A1(q2)-2.*MM/(MM+MV)*abs_p*V(q2);
1227 }
1228 
1229 double MVlnu::Hminus(double q2){
1230  double abs_p = lambda_half(MM*MM,MV*MV,q2);
1231  return (MM+MV)*A1(q2)+2.*MM/(MM+MV)*abs_p*V(q2);
1232 }
1233 
1234 double MVlnu::H0(double q2){
1235  double abs_p = lambda_half(MM*MM,MV*MV,q2);
1236  return ((MM*MM-MV*MV-q2)*(MM+MV)*A1(q2)-4.*MM*MM*abs_p*abs_p/(MM+MV)*A2(q2))/(2.*MV*sqrt(q2));
1237 }
1238 
1239 double MVlnu::H0t(double q2){
1240  double abs_p = lambda_half(MM*MM,MV*MV,q2);
1241  return 2.*MM*abs_p*A0(q2)/sqrt(q2);
1242 }
1243 
1244 double MVlnu::dGpdq2(double q2){
1245 
1246  updateParameters();
1247 
1248  double abs_p = lambda_half(MM*MM,MV*MV,q2);
1249  double lep_factor = (1.-Mlep*Mlep/q2)*(1.-Mlep*Mlep/q2)*Mlep*Mlep/(2.*q2);
1250  return 2./3.*amplsq_factor*abs_p*q2*lep_factor*(Hplus(q2)*Hplus(q2)+Hminus(q2)*Hminus(q2)+H0(q2)*H0(q2)
1251  +3.*H0t(q2)*H0t(q2));
1252 }
1253 
1254 double MVlnu::dGmdq2(double q2){
1255 
1256  updateParameters();
1257 
1258  double abs_p = lambda_half(MM*MM,MV*MV,q2);
1259  double lep_factor = (1.-Mlep*Mlep/q2)*(1.-Mlep*Mlep/q2);
1260  return 2./3.*amplsq_factor*abs_p*q2*lep_factor*(Hplus(q2)*Hplus(q2)+Hminus(q2)*Hminus(q2)+H0(q2)*H0(q2));
1261 }
1262 
1263 double MVlnu::integrateGpm(int i, double q2_min, double q2_max)
1264 {
1265  old_handler = gsl_set_error_handler_off();
1266 
1267  switch (i) {
1268  case 1:
1269  FJ = convertToGslFunction(boost::bind(&MVlnu::dGpdq2, &(*this), _1));
1270  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();
1271  gsl_set_error_handler(old_handler);
1272  return J_res;
1273  break;
1274  case 2:
1275  FJ = convertToGslFunction(boost::bind(&MVlnu::dGmdq2, &(*this), _1));
1276  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();
1277  gsl_set_error_handler(old_handler);
1278  return J_res;
1279  break;
1280  default:
1281  gsl_set_error_handler(old_handler);
1282  std::stringstream out;
1283  out << i;
1284  throw std::runtime_error("MVlnu::integrateGpm: index " + out.str() + " not implemented");
1285  }
1286 }
1287 
1288 double MVlnu::getPlep()
1289 {
1290  updateParameters();
1291 
1292  double DeltaGammaPlus = integrateGpm(1,q2min,q2max);
1293  double DeltaGammaMinus = integrateGpm(2,q2min,q2max);
1294 
1295  return (DeltaGammaPlus-DeltaGammaMinus)/(DeltaGammaPlus+DeltaGammaMinus);
1296 
1297 }
QCD::TAU
Definition: QCD.h:316
MVlnu::getDeltaGammaDeltachi
double getDeltaGammaDeltachi(double chi_min, double chi_max)
The integral of from to .
Definition: MVlnu.cpp:1113
std_make_vector.h
gslpp::cos
complex cos(const complex &z)
Definition: gslpp_complex.cpp:429
MVlnu::getDeltaGammaDeltaw
double getDeltaGammaDeltaw(double w_min, double w_max)
The integral of from to .
Definition: MVlnu.cpp:1006
MVlnu::get_R1
double get_R1(double w)
return at
Definition: MVlnu.cpp:1186
MVlnu::GF
double GF
Definition: MVlnu.h:147
MVlnu::sqrtMV_o_MM
double sqrtMV_o_MM
Definition: MVlnu.h:145
MVlnu::meson
QCD::meson meson
Definition: MVlnu.h:137
MVlnu::CLNflag
bool CLNflag
Definition: MVlnu.h:140
QCD::BOTTOM
Definition: QCD.h:329
StandardModel::getCCC5
virtual double getCCC5() const
A virtual implementation for the RealWeakEFTCC class.
Definition: StandardModel.h:1164
MVlnu::getDeltaGammaDeltacl
double getDeltaGammaDeltacl(double cl_min, double cl_max)
The integral of from to .
Definition: MVlnu.cpp:1042
MVlnu::~MVlnu
virtual ~MVlnu()
Destructor.
Definition: MVlnu.cpp:75
MVlnu::mySM
const StandardModel & mySM
Definition: MVlnu.h:135
StandardModel::getCCC1
virtual double getCCC1() const
A virtual implementation for the RealWeakEFTCC class.
Definition: StandardModel.h:1144
gslpp::sin
complex sin(const complex &z)
Definition: gslpp_complex.cpp:420
QCD::initializeMeson
void initializeMeson(QCD::meson meson_i) const
A method to initialize a meson.
Definition: QCD.cpp:236
MVlnu::getPlep
double getPlep()
Binned lepton helicity asymmetry .
Definition: MVlnu.cpp:1281
MVlnu::btocNPpmflag
bool btocNPpmflag
Definition: MVlnu.h:141
make_vector
Definition: std_make_vector.h:15
MVlnu::MV
double MV
Definition: MVlnu.h:151
MVlnu::Mnu
double Mnu
Definition: MVlnu.h:149
StandardModel.h
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
MVlnu::getDeltaGammaDeltacV
double getDeltaGammaDeltacV(double cV_min, double cV_max)
The integral of from to .
Definition: MVlnu.cpp:1077
MVlnu::MM
double MM
Definition: MVlnu.h:150
MVlnu::initializeMVlnuParameters
std::vector< std::string > initializeMVlnuParameters()
Definition: MVlnu.cpp:78
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
MVlnu::vectorM
QCD::meson vectorM
Definition: MVlnu.h:138
QCD::ELECTRON
Definition: QCD.h:312
MVlnu::MVlnu
MVlnu(const StandardModel &SM_i, QCD::meson meson_i, QCD::meson vector_i, QCD::lepton lep_i)
Constructor.
Definition: MVlnu.cpp:18
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:474
MVlnu::Mc
double Mc
Definition: MVlnu.h:156
Flavour::getFlagCLN
bool getFlagCLN() const
Definition: Flavour.h:247
MVlnu::w0
double w0
Definition: MVlnu.h:152
Meson::computeWidth
double computeWidth() const
A method to compute the width of the meson from its lifetime.
Definition: Meson.cpp:479
MVlnu::mu_b
double mu_b
Definition: MVlnu.h:154
MVlnu::Mlep
double Mlep
Definition: MVlnu.h:148
MVlnu::get_unitarity_P_BGL
double get_unitarity_P_BGL()
Pseudoscalar unitarity constraint for BGL parameters.
Definition: MVlnu.cpp:1160
MVlnu::lep
QCD::lepton lep
Definition: MVlnu.h:136
Particle::getMass
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
StandardModel::getCCC4
virtual double getCCC4() const
A virtual implementation for the RealWeakEFTCC class.
Definition: StandardModel.h:1159
QCD::meson
meson
An enum type for mesons.
Definition: QCD.h:336
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
MVlnu::get_R0
double get_R0(double w)
return at
Definition: MVlnu.cpp:1204
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
StandardModel::getFlavour
const Flavour & getFlavour() const
Definition: StandardModel.h:1006
MVlnu::width
double width
Definition: MVlnu.h:157
QCD::D_star_P
Definition: QCD.h:350
MVlnu.h
StandardModel::getMz
double getMz() const
A get method to access the mass of the boson .
Definition: StandardModel.h:718
MVlnu::RV
double RV
Definition: MVlnu.h:153
MVlnu::Mb
double Mb
Definition: MVlnu.h:155
MVlnu::get_hA1w1
double get_hA1w1()
return A1 form factor at
Definition: MVlnu.cpp:1171
Model::getModelName
std::string getModelName() const
A method to fetch the name of the model.
Definition: Model.h:59
MVlnu::get_unitarity_V_BGL
double get_unitarity_V_BGL()
Vector unitarity constraint for BGL parameters.
Definition: MVlnu.cpp:1144
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
MVlnu::get_R2
double get_R2(double w)
return at
Definition: MVlnu.cpp:1195
convertToGslFunction
gsl_function convertToGslFunction(const F &f)
Definition: gslpp_function_adapter.h:24
MVlnu::NPanalysis
bool NPanalysis
Definition: MVlnu.h:142
MVlnu::get_hA1
double get_hA1(double w)
return at
Definition: MVlnu.cpp:1178
MVlnu::MV_o_MM
double MV_o_MM
Definition: MVlnu.h:144
MVlnu::mvlnuParameters
std::vector< std::string > mvlnuParameters
Definition: MVlnu.h:139
MVlnu::get_unitarity_A_BGL
double get_unitarity_A_BGL()
Axial unitarity constraint for BGL parameters.
Definition: MVlnu.cpp:1152
FULLNLO
Definition: OrderScheme.h:37
MVlnu::getFL
double getFL()
Binned D* polarization fraction .
Definition: MVlnu.cpp:1129
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