a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
EvolDF2.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2012 HEPfit Collaboration
3  *
4  *
5  * For the licensing terms see doc/COPYING.
6  */
7 
8 #include <cstring>
9 #include "EvolDF2.h"
10 #include "StandardModel.h"
11 
12 EvolDF2::EvolDF2(unsigned int dim_i, schemes scheme, orders order, const StandardModel& model)
13 : RGEvolutor(dim_i, scheme, order),
14  model(model),
15  dim(dim_i)
16 {
17  //double Nc = model.getNc();
18  int basis = 0; //0: Gabbiani, 1: Buras
19  gslpp::matrix<double> g0t(AnomalousDimension(LO, 3, basis).transpose());
20 
21  gslpp::matrix<gslpp::complex>vv(dim, dim, 0.);
22  gslpp::vector<gslpp::complex>ee(dim, 0.);
23 
24  g0t.eigensystem(vv, ee);
25 
26  gslpp::matrix<double> v(vv.real());
27  gslpp::vector<double> e(ee.real());
28 
29  gslpp::matrix<double> vi = v.inverse();
30  for (unsigned int k = 0; k < dim; k++) {
31  a[k] = e(k);
32 // std::cout << "a[" << k << "] = " << a[k] << std::endl;
33  for (unsigned int j = 0; j < dim; j++) {
34  for (unsigned int i = 0; i < dim; i++) {
35  b[i][j][k] = v(i, k) * vi(k, j);
36 // if (b[i][j][k] != 0.)
37 // std::cout << "b[" << i << "][" << j << "][" << k << "] = " << b[i][j][k] << std::endl;
38  }
39  }
40  }
41 
42  gslpp::matrix<double> h(dim, dim, 0.);
43  for (int l = 0; l < 3; l++) {
44  gslpp::matrix<double> gg = vi * (AnomalousDimension(NLO, 6 - l, basis).transpose()) * v;
45  double b0 = model.Beta0(6 - l);
46  for (unsigned int i = 0; i < dim; i++)
47  for (unsigned int j = 0; j < dim; j++)
48  h(i, j) = (i == j) * e(i) * model.Beta1(6 - l) / (2. * b0 * b0) -
49  gg(i, j) / (2. * b0 + e(i) - e(j));
50  gslpp::matrix<double> j = v * h * vi;
51  gslpp::matrix<double> jv = j*v;
52  gslpp::matrix<double> vij = vi*j;
53  for (unsigned int i = 0; i < dim; i++)
54  for (unsigned int j = 0; j < dim; j++)
55  for (unsigned int k = 0; k < dim; k++) {
56  c[l][i][j][k] = jv(i, k) * vi(k, j);
57 // if (c[l][i][j][k] != 0.)
58 // std::cout << "c[" << l << "][" << i << "][" << j << "][" << k << "] = " << c[l][i][j][k] << std::endl;
59  d[l][i][j][k] = -v(i, k) * vij(k, j);
60 // if (d[l][i][j][k] != 0.)
61 // std::cout << "d[" << l << "][" << i << "][" << j << "][" << k << "] = " << d[l][i][j][k] << std::endl;
62  }
63  }
64 }
65 
66 EvolDF2::~EvolDF2()
67 {}
68 
69 gslpp::matrix<double> EvolDF2::AnomalousDimension(orders order, unsigned int nf, int basis) const
70 {
71  gslpp::matrix<double> ad(dim, dim, 0.);
72  double Nc = model.getNc();
73  switch (basis) {
74  case 0:
75  switch (order) {
76  case LO:
77  ad(0, 0) = (6. * Nc - 6.) / Nc;
78  ad(1, 1) = -(6. * Nc * Nc - 8. * Nc - 2.) / Nc;
79  ad(1, 2) = (4. * Nc - 8.) / Nc;
80  ad(2, 1) = (4. * Nc * Nc - 4. * Nc - 8.) / Nc;
81  ad(2, 2) = (2. * Nc * Nc + 4. * Nc + 2.) / Nc;
82  ad(3, 3) = -(6. * Nc * Nc - 6.) / Nc;
83  ad(4, 3) = -6.;
84  ad(4, 4) = 6. / Nc;
85  break;
86  case NLO:
87 
88  // MSbar-NDR scheme with evanescent operators of Buras, Misiak & Urban
89 
90  if (!(nf == 3 || nf == 4 || nf == 5 || nf == 6))
91  throw std::runtime_error("EvolDF2::AnomalousDimension(): wrong number of flavours");
92  ad(0, 0) = -(-1. + Nc)*(-171. + 19. * Nc * Nc + Nc * (63. - 4. * nf)) / (6. * Nc * Nc);
93  ad(1, 1) = (-1251. - 609. * Nc * Nc * Nc * Nc + Nc * (432. - 52. * nf) - 8. * Nc * Nc * (-71 + 2. * nf) +
94  20. * Nc * Nc * Nc * (32. + 3. * nf)) / (18. * Nc * Nc);
95  ad(1, 2) = -(2. * (-2. + Nc)*(-72. + Nc * Nc + 2. * Nc * (63. + nf))) / (9. * Nc * Nc);
96  ad(2, 1) = 2. * (119. * Nc * Nc * Nc + 8. * (-9. + 5. * nf) + 2. * Nc * (-125. + 7. * nf) -
97  Nc * Nc * (101. + 14. * nf)) / (9. * Nc);
98  ad(2, 2) = (477. + 343. * Nc * Nc * Nc * Nc + Nc * Nc * Nc * (380. - 52. * nf) - 4. * Nc * (-36. + nf) -
99  8. * Nc * Nc * (16. + 13. * nf)) / (18. * Nc * Nc);
100  ad(3, 3) = (45. + 479. * Nc * Nc - 203. * Nc * Nc * Nc * Nc - 44. * Nc * nf + 20. * Nc * Nc * Nc * nf) /
101  (6. * Nc * Nc);
102  ad(3, 4) = -(18. / Nc)-(71. * Nc) / 2. + 4. * nf;
103  ad(4, 3) = 3. / Nc - (100. * Nc) / 3. + (22. * nf) / 3.;
104  ad(4, 4) = (45. + 137. * Nc * Nc - 44. * Nc * nf) / (6. * Nc * Nc);
105  break;
106  default:
107  throw std::runtime_error("EvolDF2::AnomalousDimension(): order not implemented");
108  }
109  break;
110  case 1:
111  switch (order) {
112  case LO:
113  ad(0, 0) = 6. - 6. / Nc;
114  ad(1, 1) = 6. / Nc;
115  ad(1, 2) = 12.;
116  ad(2, 2) = 6. / Nc - 6. * Nc;
117  ad(3, 3) = 6. + 6. / Nc - 6. * Nc;
118  ad(3, 4) = 1. / 2. - 1. / Nc;
119  ad(4, 3) = -24. - 48. / Nc;
120  ad(4, 4) = 6. - 2. / Nc + 2. * Nc;
121  break;
122  case NLO:
123 
124  // MSbar-NDR scheme with evanescent operators of Buras, Misiak & Urban
125 
126  if (!(nf == 3 || nf == 4 || nf == 5 || nf == 6))
127  throw std::runtime_error("EvolDF2::AnomalousDimension(): wrong number of flavours");
128  ad(0, 0) = -22. / 3. - 57. / (2. * Nc * Nc) + 39. / Nc - (19. * Nc) / 6. + (2. * nf) / 3. - (2. * nf) / 3. / Nc;
129  ad(1, 1) = 137. / 6. + 15. / (2 * Nc * Nc) - (22 * nf) / (3 * Nc);
130  ad(1, 2) = -(6. / Nc) + (200. * Nc) / 3. - (44 * nf) / 3.;
131  ad(2, 1) = 9. / Nc + (71. * Nc) / 4. - 2. * nf;
132  ad(2, 2) = 479. / 6. + 15. / (2. * Nc * Nc) - (203. * Nc * Nc) / 6. - (22. * nf) / (3. * Nc) + (10. * Nc * nf) / 3.;
133  ad(3, 3) = 136. / 3. - 107. / (2. * Nc * Nc) - 12. / Nc + (107. * Nc) / 3. - (203. * Nc * Nc) / 6. - (2. * nf) / 3. - (10. * nf) / (3. * Nc) + (10. * Nc * nf) / 3.;
134  ad(3, 4) = -31. / 9. - 4. / (Nc * Nc) + 9. / Nc - Nc / 36. - nf / 18. + nf / (9 * Nc);
135  ad(4, 3) = -704. / 3. - 320. / (Nc * Nc) - 208. / Nc - (364. * Nc) / 3. + (136. * nf) / 3. + (176. * nf) / (3. * Nc);
136  ad(4, 4) = -188. / 9. + 21. / (2. * Nc * Nc) + 44. / Nc + 21. * Nc + (343. * Nc * Nc) / 18. - 6. * nf + (2. * nf) / (9. * Nc) - (26. * Nc * nf) / 9.;
137  break;
138  default:
139  throw std::runtime_error("EvolDF2::AnomalousDimension(): order not implemented");
140  }
141  break;
142  default:
143  throw std::runtime_error("EvolDF2::AnomalousDimension(): basis not implemented");
144  }
145  return (ad);
146 }
147 
148 gslpp::matrix<double>& EvolDF2::Df2Evol(double mu, double M, orders order, schemes scheme)
149 {
150  switch (scheme) {
151  case NDR:
152  break;
153  case LRI:
154  case HV:
155  default:
156  std::stringstream out;
157  out << scheme;
158  throw std::runtime_error("EvolDF2::Df2Evol(): scheme " + out.str() + " not implemented ");
159  }
160 
161  double alsMZ = model.getAlsMz();
162  double Mz = model.getMz();
163  if(alsMZ == alsMZ_cache && Mz == Mz_cache) {
164  if (mu == this->mu && M == this->M && scheme == this->scheme)
165  return (*Evol(order));
166  }
167  alsMZ_cache = alsMZ;
168  Mz_cache = Mz;
169 
170  if (M < mu) {
171  std::stringstream out;
172  out << "M = " << M << " < mu = " << mu;
173  throw out.str();
174  }
175 
176  setScales(mu, M); // also assign evol to identity
177  if (M != mu) {
178  double m_down = mu;
179  double m_up = model.AboveTh(m_down);
180  double nf = model.Nf(m_down);
181 
182  while (m_up < M) {
183  Df2Evol(m_down, m_up, nf, scheme);
184  m_down = m_up;
185  m_up = model.AboveTh(m_down);
186  nf += 1.;
187  }
188  Df2Evol(m_down, M, nf, scheme);
189  }
190  return (*Evol(order));
191 }
192 
193 void EvolDF2::Df2Evol(double mu, double M, double nf, schemes scheme)
194 {
195 
196  gslpp::matrix<double> resLO(dim, 0.), resNLO(dim, 0.), resNNLO(dim, 0.);
197 
198  int l = 6 - (int) nf;
199  double alsM = model.Als(M, FULLNLO) / 4. / M_PI;
200  double alsmu = model.Als(mu, FULLNLO) / 4. / M_PI;
201 
202  double eta = alsM / alsmu;
203 
204  for (unsigned int k = 0; k < dim; k++) {
205  double etap = pow(eta, a[k] / 2. / model.Beta0(nf));
206  for (unsigned int i = 0; i < dim; i++)
207  for (unsigned int j = 0; j < dim; j++) {
208  resNNLO(i, j) += 0.;
209  resNLO(i, j) += c[l][i][j][k] * etap * alsmu;
210  resNLO(i, j) += d[l][i][j][k] * etap * alsM;
211  resLO(i, j) += b[i][j][k] * etap;
212  }
213  }
214  switch (order) {
215  case NNLO:
216  *elem[NNLO] = 0.; // Marco can implement it if he wishes to!
217  case NLO:
218  *elem[NLO] = (*elem[LO]) * resNLO + (*elem[NLO]) * resLO;
219  case LO:
220  *elem[LO] = (*elem[LO]) * resLO;
221  break;
222  case FULLNNLO:
223  case FULLNLO:
224  default:
225  throw std::runtime_error("Error in EvolDF2::Df2Evol()");
226  }
227 }
228 
229 double EvolDF2::etacc(double mu) const
230 {
231  double N = model.getNc();
232  double N2 = N * N;
233  double gamma0[2] = {0.}, gamma1[2][4] = {
234  {0.}
235  }, d[2][4] = {
236  {0.}
237  },
238  J[2][4] = {
239  {0.}
240  }, dd[2][2][4] = {
241  {
242  {0.}
243  }
244  }, JJ[2][2][4] = {
245  {
246  {0.}
247  }
248  },
249  B[2] = {0.}; // 0 = + ; 1 = - ;
250  double tau[2][2] = {
251  {0.}
252  }, K[2][2] = {
253  {0.}
254  }, beta[2][2] = {
255  {0.}
256  };
257 
258  gamma0[0] = 6. * (N - 1.) / N;
259  gamma0[1] = -6. * (N + 1.) / N;
260 
261  gamma1[0][0] = ((N - 1.) / (2. * N)) * (-21. + 57. / N - 19. / 3. * N + 4. / 3. * 3.);
262  gamma1[1][0] = ((N + 1.) / (2. * N)) * (-21. - 57. / N + 19. / 3. * N - 4. / 3. * 3.);
263  gamma1[0][1] = ((N - 1.) / (2. * N)) * (-21. + 57. / N - 19. / 3. * N + 4. / 3. * 4.);
264  gamma1[1][1] = ((N + 1.) / (2. * N)) * (-21. - 57. / N + 19. / 3. * N - 4. / 3. * 4.);
265  gamma1[0][2] = ((N - 1.) / (2. * N)) * (-21. + 57. / N - 19. / 3. * N + 4. / 3. * 5.);
266  gamma1[1][2] = ((N + 1.) / (2. * N)) * (-21. - 57. / N + 19. / 3. * N - 4. / 3. * 5.);
267  gamma1[0][3] = ((N - 1.) / (2. * N)) * (-21. + 57. / N - 19. / 3. * N + 4. / 3. * 6.);
268  gamma1[1][3] = ((N + 1.) / (2. * N)) * (-21. - 57. / N + 19. / 3. * N - 4. / 3. * 6.);
269 
270  for (int i = 0; i < 2; i++) { // 0 = + ; 1 = - ;
271  for (int j = 0; j < 4; j++) { // j = nf - 3;
272  d[i][j] = gamma0[i] / 2. / model.Beta0(j + 3);
273  J[i][j] = d[i][j] / model.Beta0(j + 3) * model.Beta1(j + 3) -
274  gamma1[i][j] / 2. / model.Beta0(j + 3);
275  }
276 
277  }
278 
279  for (int i = 0; i < 2; i++) { // 0 = + ; 1 = - ;
280  for (int j = 0; j < 2; j++) { // 0 = + ; 1 = - ;
281  for (int k = 0; k < 4; k++) { // k = nf - 3;
282  dd[i][j][k] = d[i][k] + d[j][k];
283  JJ[i][j][k] = J[i][k] + J[j][k];
284  }
285  }
286  }
287 
288  tau[0][0] = (N + 3.) / 4.;
289  tau[1][0] = -(N - 1.) / 4.;
290  tau[0][1] = tau[1][0];
291  tau[1][1] = (N - 1.) / 4.;
292 
293  K[0][0] = 3. * (N - 1.) * tau[0][0];
294  K[1][0] = 3. * (N + 1.) * tau[0][1];
295  K[0][1] = K[1][0];
296  K[1][1] = 3. * (N + 3.) * tau[1][1];
297 
298  beta[0][0] = (1. - N) * (M_PI * M_PI * (N2 - 6.) / (12. * N) + 3. * (-N2 + 2. * N + 13.) / (4. * N));
299  beta[1][0] = (1. - N) * (M_PI * M_PI * (-N2 + 2. * N - 2.) / (12. * N) + (3. * N2 + 13.) / (4. * N));
300  beta[0][1] = beta[1][0];
301  beta[1][1] = (1. - N) * (M_PI * M_PI * (N2 - 4. * N + 2) / (12. * N) - (3. * N2 + 10. * N + 13.) / (4. * N));
302 
303  B[0] = 11. * (N - 1.) / (2. * N);
304  B[1] = -11. * (N + 1.) / (2. * N);
305 
306  double eta = 0.;
307  double as_mb__as_muc = model.Als(model.getMub(), FULLNLO) / model.Als(model.getMuc(), FULLNLO);
308  double as_muw__as_mb = model.Als(model.getMuw(), FULLNLO) / model.Als(model.getMub(), FULLNLO);
309  double as_muc__4pi = model.Als(model.getMuc(), FULLNLO) / 4. / M_PI;
310  double log_muc__mc = log(model.getMuc() / model.getQuarks(QCD::CHARM).getMass());
311  double as_mb__4pi = model.Als(model.getMub(), FULLNLO) / 4. / M_PI;
312  double as_muw__4pi = model.Als(model.getMuw(), FULLNLO) / 4. / M_PI;
313 
314 
315  for (int i = 0; i < 2; i++) {
316  for (int j = 0; j < 2; j++) {
317  eta += pow(as_mb__as_muc, dd[i][j][1]) *
318  pow(as_muw__as_mb, dd[i][j][2]) *
319  (tau[i][j] + as_muc__4pi * (2. * K[i][j] *
320  log_muc__mc + tau[i][j] * (JJ[i][j][1] - J[0][0]) + beta[i][j]) +
321  tau[i][j]*(as_mb__4pi * (JJ[i][j][2] - JJ[i][j][1]) +
322  as_muw__4pi * (B[i] + B[j] - JJ[i][j][2])));
323  }
324  }
325 
326  eta *= pow(model.Als(model.getMuc()), d[0][0]);
327  eta *= pow(model.Als(mu, FULLNLO), -2. / 9.);
328 
329  return (eta * (1. + model.Als(mu, FULLNLO) / 4. / M_PI * J[0][0]));
330 }
331 
332 double EvolDF2::etact(double mu) const
333 {
334 
335  //temporary fix waiting for NNLO
336 
337  double K = model.Als4(model.getMuw()) / model.Als4(model.getMuc());
338 #if SUSYFIT_DEBUG & 2
339  std::cout << "K = " << K << std::endl;
340 #endif
341  double Kpp = pow(K, 12. / 25.);
342  double Kpm = pow(K, -6. / 25.);
343  double Kmm = pow(K, -24. / 25.);
344  double K7 = pow(K, 1. / 5.);
345  double xc = model.Mrun4(model.getMuc(), model.getQuarks(QCD::CHARM).getMass_scale(), model.getQuarks(QCD::CHARM).getMass()) / model.Mw();
346  xc *= xc;
347  double xt = model.Mrun4(model.getMuw(), model.getQuarks(QCD::TOP).getMass_scale(), model.getQuarks(QCD::TOP).getMass()) / model.Mw();
348  xt *= xt;
349 
350  double J3 = 6. * (3. - 1.) / 3. / 2. / model.Beta0(3) / model.Beta0(3) * model.Beta1(3) -
351  ((3. - 1.) / (2. * 3.)) * (-21. + 57. / 3. - 19. + 4.) / 2. / model.Beta0(3);
352 
353 
354  double eta = 0.;
355  double AlsC = model.Als4(model.getMuc());
356 
357  eta = (M_PI / AlsC * (-18. / 7. * Kpp - 12. / 11. * Kpm + 6. / 29. * Kmm + 7716. / 2233. * K7) *
358  (1. - AlsC / (4. * M_PI) * 307. / 162.) + (log(model.getMuc() /
359  model.getQuarks(QCD::CHARM).getMass()) - 0.25) * (3. * Kpp - 2. * Kpm + Kmm) +
360  262497. / 35000. * Kpp - 123. / 625. * Kpm + 1108657. / 1305000. * Kmm - 277133. / 50750. * K7 +
361  K * (-21093. / 8750. * Kpp + 13331. / 13750. * Kpm - 10181. / 18125. * Kmm - 1731104. / 2512125. * K7) +
362  (log(xt) - (3. * xt) / (4. - 4. * xt) - log(xt) * (3. * xt * xt) / (4. * (1. - xt) * (1. - xt)) + 0.5) * K * K7) * xc / (model.getMatching().S0(xc, xt)) * pow(AlsC, 2. / 9.);
363 
364  return (eta * (1. + model.Als(mu, FULLNLO) / 4. / M_PI * J3) * pow(model.Als(mu, FULLNLO), -2. / 9.));
365 }
366 
367 double EvolDF2::etatt(double m) const
368 {
369  double N = model.getNc();
370  double x = model.getMatching().x_t(model.getMut());
371  double x2 = x * x;
372  double x3 = x2 * x;
373  double x4 = x3 * x;
374  double xm2 = (1 - x) * (1 - x);
375  double xm3 = xm2 * (1 - x);
376  //double xm4 = xm3 * (1 - x);
377  double logx = log(x);
378  //double Li2 = gsl_sf_dilog(1-x);
379 
380  double S0tt = (4. * x - 11. * x2 + x3) / 4. / xm2 -
381  3. * x3 / (2. * xm3) * logx;
382 
383  double Bt = 5. * (N - 1.) / 2. / N + 3. * (N * N - 1.) / 2. / N;
384 
385  double gamma0 = 6. * (N - 1.) / N;
386 
387  double gamma1[4] = {0.}, J[4] = {0.};
388 
389  for (int i = 0; i < 4; ++i) {
390  gamma1[i] = (N - 1.) / (2. * N) * (-21. + 57. / N - 19. / 3. * N + 4. / 3. * (i + 3.));
391  J[i] = gamma0 * model.Beta1(3 + i) / 2. / model.Beta0(3 + i) / model.Beta0(3 + i)
392  - gamma1[i] / 2. / model.Beta0(3 + i);
393  }
394 
395  double b = (4. - 22. * x + 15. * x2 + 2. * x3 + x4 - 18. * x2 * logx)
396  / ((-1. + x) * (-4. + 15. * x - 12. * x2 + x3 + 6. * x2 * logx));
397 
398  double AlsT = model.Als(model.getMut());
399  double AlsB = model.Als(model.getMub());
400  double AlsC = model.Als(model.getMuc());
401  double eta = pow(AlsC, 6. / 27.) *
402  pow(AlsB / AlsC, 6. / 25.) *
403  pow(AlsT / AlsB, 6. / 23.) *
404  (1. + AlsC / 4. / M_PI * (J[1] - J[0]) +
405  AlsB / 4. / M_PI * (J[2] - J[1])
406  + AlsT / 4. / M_PI * (model.getMatching().S1(x) / S0tt
407  + Bt - J[2] + gamma0 * log(model.getMut() / model.getMuw())
408  + 6 * (N * N - 1) / N * log(model.getMut() / model.getMuw()) * b));
409  /* double J3 = 6. * (N - 1.) / N * (model.Beta1(3) / 2. / model.Beta0(3) / model.Beta0(3)) -
410  (N - 1.) / (2. * N) * (-21. + 57. / N - 19./3. * N + 4.) / 2. / model.Beta0(3);*/
411 #if SUSYFIT_DEBUG & 2
412  std::cout << "AlsT = " << AlsT << " AlsC = " << AlsC << " AlsB = " << AlsB << " mub = " << model.getMub() << std::endl;
413  std::cout << "etatt = " << eta*(1. + model.Als(m, FULLNLO) / 4. / M_PI * J[0]) * pow(model.Als(m, FULLNLO), -2. / 9.)
414  << " mut = " << model.getMut() << " Muw = " << model.getMuw() << std::endl;
415 #endif
416  return (eta * (1. + model.Als(m, FULLNLO) / 4. / M_PI * J[0]) * pow(model.Als(m, FULLNLO), -2. / 9.));
417 }
418 
419 double EvolDF2::etabS0(double m) const
420 {
421  double N = model.getNc();
422  double x = model.getMatching().x_t(model.getMut());
423  double x2 = x * x;
424  double x3 = x2 * x;
425  double x4 = x3 * x;
426  double xm2 = (1 - x) * (1 - x);
427  double xm3 = xm2 * (1 - x);
428  //double xm4 = xm3 * (1 - x);
429  double logx = log(x);
430  //double Li2 = gsl_sf_dilog(1-x);
431 
432  double S0tt = (4. * x - 11. * x2 + x3) / 4. / xm2 -
433  3. * x3 / (2. * xm3) * logx;
434 
435  double Bt = 5. * (N - 1.) / 2. / N + 3. * (N * N - 1.) / 2. / N;
436 
437  double gamma0 = 6. * (N - 1.) / N;
438 
439  double gamma1[4] = {0.}, J[4] = {0.};
440 
441  for (int i = 0; i < 4; ++i) {
442  gamma1[i] = (N - 1.) / (2. * N) * (-21. + 57. / N - 19. / 3. * N + 4. / 3. * (i + 3.));
443  J[i] = gamma0 * model.Beta1(3 + i) / 2. / model.Beta0(3 + i) / model.Beta0(3 + i)
444  - gamma1[i] / 2. / model.Beta0(3 + i);
445  }
446 
447  double b = (4. - 22. * x + 15. * x2 + 2. * x3 + x4 - 18. * x2 * logx)
448  / ((-1. + x) * (-4. + 15. * x - 12. * x2 + x3 + 6. * x2 * logx));
449 
450  double AlsT = model.Als(model.getMut());
451  double Alsm = model.Als(m);
452  double eta =
453  pow(AlsT / Alsm, 6. / 23.) *
454  (1.
455  + AlsT / 4. / M_PI * (model.getMatching().S1(x) / S0tt
456  + Bt - J[2] + gamma0 * log(model.getMut() / model.getMuw())
457  + 6 * (N * N - 1) / N * log(model.getMut() / model.getMuw()) * b)
458  + Alsm / 4. / M_PI * J[2]);
459  return (eta * S0tt);
460 }
461 
462 /*double EvolDF2::S1tt() const {
463  double N = model.getNc();
464  double x = model.getMyMatching()->x_t(model.getMut());
465  double x2 = x * x;
466  double x3 = x2 * x;
467  double x4 = x3 * x;
468  double xm2 = (1 - x) * (1 - x);
469  double xm3 = xm2 * (1 - x);
470  double xm4 = xm3 * (1 - x);
471  double logx = log(x);
472  double Li2 = gsl_sf_dilog(1-x);
473 
474  double S0tt = (4. * x - 11. * x2 + x3) / 4. / xm2 -
475  3. * x3 / (2. * xm3) * logx;
476 
477  double Bt = 5. * (N - 1.) / 2. / N + 3. * (N * N - 1.) / 2. / N;
478 
479  double gamma0 = 6. * (N - 1.) / N;
480 
481  double gamma1[4] = {0.}, J[4] = {0.};
482 
483  for(int i = 0; i < 4; ++i){
484  gamma1[i] = (N - 1.)/(2. * N) * (-21. + 57./N - 19./3. * N + 4./3. * (i + 3.));
485  J[i] = gamma0 * model.Beta1(3 + i) / 2. / model.Beta0(3 + i) / model.Beta0(3 + i)
486  - gamma1[i] / 2. / model.Beta0(3 + i);
487  }
488 
489  double b = (4. - 22. * x + 15. * x2 + 2. * x3 + x4 - 18. * x2 * logx)
490  / ((-1. + x) * (-4. + 15. * x - 12. * x2 + x3 + 6. * x2 * logx));
491  double AlsT = model.Als(model.getMut());
492  double AlsB = model.Als(model.getMub());
493  double AlsC = model.Als(model.getMuc());
494 
495  return (pow(model.Als(model.getMuc()), 6./27.) *
496  pow(model.Als(model.getMub())/model.Als(model.getMuc()), 6./25.) *
497  pow(model.Als(model.getMut())/model.Als(model.getMub()), 6./23.) *
498  (1. + model.Als(model.getMuc())/4./M_PI * (J[1]-J[0]) +
499  model.Als(model.getMub())/4./M_PI * (J[2]-J[1])
500  + model.Als(model.getMut())/4./M_PI * (model.getMyMatching()->S1(x)/S0tt
501  + Bt - J[2] + gamma0 * log(model.getMut() / model.getMuw())
502  + 6 * (N * N - 1) / N * log(model.getMut() / model.getMuw()) * b)));
503 }*/
WilsonTemplate< gslpp::matrix< double > >::scheme
schemes scheme
Definition: WilsonTemplate.h:117
gslpp::matrix< double >
A class for constructing and defining operations on real matrices.
Definition: gslpp_matrix_double.h:48
EvolDF2::alsMZ_cache
double alsMZ_cache
Definition: EvolDF2.h:99
QCD::Beta1
double Beta1(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:471
RGEvolutor
A class for the RG evolutor of the Wilson coefficients.
Definition: RGEvolutor.h:24
WilsonTemplate< gslpp::matrix< double > >::mu
double mu
Definition: WilsonTemplate.h:116
RGEvolutor::M
double M
Definition: RGEvolutor.h:142
QCD::Nf
double Nf(const double mu) const
The number of active flavour at scale .
Definition: QCD.cpp:438
EvolDF2::etacc
double etacc(double mu) const
Buras et al, hep-ph/9512380.
Definition: EvolDF2.cpp:229
StandardModel::getAlsMz
double getAlsMz() const
A get method to access the value of .
Definition: StandardModel.h:727
WilsonTemplate< gslpp::matrix< double > >::order
orders order
Definition: WilsonTemplate.h:118
EvolDF2::d
double d[3][5][5][5]
Definition: EvolDF2.h:96
LO
Definition: OrderScheme.h:33
StandardModel.h
QCD::CHARM
Definition: QCD.h:326
NDR
Definition: OrderScheme.h:21
gslpp::matrix< double >::transpose
matrix< double > transpose() const
Definition: gslpp_matrix_double.cpp:166
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
EvolDF2::AnomalousDimension
gslpp::matrix< double > AnomalousDimension(orders order, unsigned int nf, int basis=0) const
ADM in the basis used in Ciuchini et.al. hep-ph/9711402 (basis = 0, default) or in the basis (QVLL,...
Definition: EvolDF2.cpp:69
gslpp::matrix< gslpp::complex >
Particle::getMass_scale
double getMass_scale() const
A get method to access the scale at which the particle mass is defined.
Definition: Particle.h:133
EvolDF2::etabS0
double etabS0(double mu) const
Buras et al, hep-ph/9512380.
Definition: EvolDF2.cpp:419
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:474
EvolDF2::b
double b[5][5][5]
Definition: EvolDF2.h:94
QCD::Beta0
double Beta0(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:466
EvolDF2::a
double a[5]
Definition: EvolDF2.h:93
QCD::Als4
double Als4(const double mu) const
The value of at any scale with the number of flavours .
Definition: QCD.cpp:536
Particle::getMass
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
schemes
schemes
An enum type for regularization schemes.
Definition: OrderScheme.h:19
gslpp::matrix< double >::inverse
matrix< double > inverse()
Definition: gslpp_matrix_double.cpp:178
EvolDF2::Mz_cache
double Mz_cache
Definition: EvolDF2.h:100
QCD::TOP
Definition: QCD.h:328
EvolDF2::etatt
double etatt(double mu) const
Buras et al, hep-ph/9512380.
Definition: EvolDF2.cpp:367
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
QCD::getMut
double getMut() const
A get method to access the threshold between six- and five-flavour theory in GeV.
Definition: QCD.h:561
QCD::getMuc
double getMuc() const
A get method to access the threshold between four- and three-flavour theory in GeV.
Definition: QCD.h:579
NNLO
Definition: OrderScheme.h:35
EvolDF2::EvolDF2
EvolDF2(unsigned int dim_i, schemes scheme, orders order, const StandardModel &model)
constructor
Definition: EvolDF2.cpp:12
StandardModel::Als
double Als(double mu, orders order=FULLNLO, bool qed_flag=false, bool Nf_thr=true) const
The running QCD coupling in the scheme including QED corrections.
Definition: StandardModel.cpp:576
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
gslpp::vector< double >
A class for constructing and defining operations on real vectors.
Definition: gslpp_vector_double.h:33
orders
orders
An enum type for orders in QCD.
Definition: OrderScheme.h:31
EvolDF2.h
LRI
Definition: OrderScheme.h:23
StandardModel::getMz
double getMz() const
A get method to access the mass of the boson .
Definition: StandardModel.h:718
EvolDF2::c
double c[3][5][5][5]
Definition: EvolDF2.h:95
EvolDF2::etact
double etact(double mu) const
Buras et al, hep-ph/9512380.
Definition: EvolDF2.cpp:332
QCD::getMub
double getMub() const
A get method to access the threshold between five- and four-flavour theory in GeV.
Definition: QCD.h:570
RGEvolutor::Evol
gslpp::matrix< double > * Evol(orders order)
Evolution matrix set at a fixed order of QCD coupling.
Definition: RGEvolutor.cpp:103
StandardModel::getMatching
virtual StandardModelMatching & getMatching() const
A get method to access the member reference of type StandardModelMatching.
Definition: StandardModel.h:949
EvolDF2::dim
unsigned int dim
Definition: EvolDF2.h:98
HV
Definition: OrderScheme.h:22
EvolDF2::model
const StandardModel & model
Definition: EvolDF2.h:97
QCD::AboveTh
double AboveTh(const double mu) const
The active flavour threshold above the scale as defined in QCD::Thresholds().
Definition: QCD.cpp:420
QCD::getNc
double getNc() const
A get method to access the number of colours .
Definition: QCD.h:505
RGEvolutor::setScales
void setScales(double mu, double M)
Sets the upper and lower scale for the running of the Wilson Coefficients.
Definition: RGEvolutor.cpp:85
NLO
Definition: OrderScheme.h:34
StandardModel::Mw
virtual double Mw() const
The SM prediction for the -boson mass in the on-shell scheme, .
Definition: StandardModel.cpp:944
EvolDF2::Df2Evol
gslpp::matrix< double > & Df2Evol(double mu, double M, orders order, schemes scheme=NDR)
Definition: EvolDF2.cpp:148
FULLNNLO
Definition: OrderScheme.h:38
QCD::Mrun4
double Mrun4(const double mu_f, const double mu_i, const double m) const
The running of a mass with the number of flavours .
Definition: QCD.cpp:1199
gslpp::matrix< double >::eigensystem
void eigensystem(matrix< complex > &U, vector< complex > &S)
Definition: gslpp_matrix_double.cpp:280
FULLNLO
Definition: OrderScheme.h:37
WilsonTemplate< gslpp::matrix< double > >::elem
gslpp::matrix< double > * elem[MAXORDER_QED+1]
Definition: WilsonTemplate.h:114
gslpp::vector< gslpp::complex >
EvolDF2::~EvolDF2
virtual ~EvolDF2()
destructor
Definition: EvolDF2.cpp:66
StandardModel::getMuw
double getMuw() const
A get method to retrieve the matching scale around the weak scale.
Definition: StandardModel.h:935