a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
EvolDF1.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 <map>
9 #include <vector>
10 #include <initializer_list>
11 #include "QCD.h"
12 #include "EvolDF1.h"
13 
14 #define EPS 1.e-15
15 
16 std::map<std::string, uint> blocks_nops = {
17  {"C", 2},
18  {"CP", 6},
19  {"CPM", 8},
20  {"L", 2},
21  {"CPL", 8},
22  {"CPQ", 10},
23  {"CPML", 10},
24  {"CPQB", 11},
25  {"CPMQB", 13},
26  {"CPMLQB", 15}
27 };
28 
29 EvolDF1::EvolDF1(std::string reqblocks, schemes scheme, const StandardModel& model_i, qcd_orders ord_qcd, qed_orders ord_qed)
30 : RGEvolutorNew(blocks_nops.at(reqblocks), scheme, ord_qcd, ord_qed), model(model_i), blocks(reqblocks),
31 evec(blocks_nops.at(reqblocks), 0.), evec_i(blocks_nops.at(reqblocks), 0.), js(blocks_nops.at(reqblocks), 0.),
32 h(blocks_nops.at(reqblocks), 0.), gg(blocks_nops.at(reqblocks), 0.), s_s(blocks_nops.at(reqblocks), 0.),
33 jssv(blocks_nops.at(reqblocks), 0.), jss(blocks_nops.at(reqblocks), 0.), jv(blocks_nops.at(reqblocks), 0.),
34 vij(blocks_nops.at(reqblocks), 0.), eval(blocks_nops.at(reqblocks), 0.), evecc(blocks_nops.at(reqblocks), 0.),
35 evalc(blocks_nops.at(reqblocks), 0.)
36 {
37  // blocks_ord = {
38  // {"C", NNLO},
39  // {"CP", NNLO},
40  // {"CPM", NNLO},
41  // {"L", NNLO},
42  // {"CPML", NNLO},
43  // {"CPQB", NLO},
44  // {"CPMQB", NLO},
45  // {"CPMLQB", NLO}};
46 
47  // if (blocks_nops[blocks] != nops)
48  // throw std::runtime_error("EvolDF1(): number of operators does not match block specification");
49 
50  this->nops = blocks_nops.at(reqblocks);
51  uint nf, nnf, nu, nd, a, b, i, j, p, q;
52  double b0, b0e, b1, b2, b3, b4, term;
53 
54  alsM_cache = 0.;
55  MAls_cache = 0.;
56 
57  gslpp::matrix<double> W10(nops, nops, 0.), W20(nops, nops, 0.), W30(nops, nops, 0.),
58  W01(nops, nops, 0.), W02(nops, nops, 0.), W11(nops, nops, 0.), W21(nops, nops, 0.);
59  gslpp::matrix<double> M1(nops, nops, 0.), M2(nops, nops, 0.), M3(nops, nops, 0.), M4(nops, nops, 0.),
60  M5(nops, nops, 0.), M6(nops, nops, 0.);
61 
62  if (order_qed == QED0 && blocks.find("L") == std::string::npos &&
63  blocks.find("Q") == std::string::npos && blocks.find("B") == std::string::npos)
64  {
65  nfmin = 3;
66  nfmax = 6;
67  } else
68  nfmin = nfmax = 5;
69 
70 
71  for (nf = nfmin; nf <= nfmax; nf++)
72  {
73  nnf = nf - nfmin;
74  nu = nf / 2;
75  nd = nf % 2 == 0 ? nf / 2 : nf / 2 + 1;
76 
77  b0 = model.Beta_s(00, nf);
78  b1 = model.Beta_s(10, nf) / 2. / b0 / b0;
79  b2 = model.Beta_s(20, nf) / 4. / b0 / b0 / b0 - b1 * b1;
80 
81  W10 = AnomalousDimension(10, nu, nd).transpose() / 2. / b0;
82  W20 = AnomalousDimension(20, nu, nd).transpose() / 4. / b0 / b0;
83  W30 = AnomalousDimension(30, nu, nd).transpose() / 8. / b0 / b0 / b0;
84 
85  // std::cout << AnomalousDimension(01, nu, nd).transpose() << std::endl;
86  // std::cout << W10 << std::endl;
87 
88  // Misiak-Munz basis, defined as in T. Huber et al., hep-ph/0512066
89  W10.eigensystem(evecc, evalc);
90  for (size_t jj = 0; jj < evalc.size(); jj++)
91  if (fabs(evalc(jj).imag()) > EPS) throw ("check the imaginary part of eigenvalues of W10");
92  eval = evalc.real();
93  evec = evecc.real();
94  evec_i = (evecc.inverse()).real();
95 
96  // QCD magic numbers
97  // M2: B(-2), M1: B(-1)_10
98  M2 = evec_i * (W30 - b1 * W20 - b2 * W10) * evec;
99  M1 = evec_i * (W20 - b1 * W10) * evec;
100 
101  for (a = 0; a < nops; a++)
102  {
103  ai[nnf].insert(std::pair<uint, double > (a, eval(a)));
104  for (b = 0; b < nops; b++)
105  for (i = 0; i < nops; i++)
106  {
107  if (fabs(term = evec(a, i) * evec_i(i, b)) > EPS)
108  vM0vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i}, term)); // QCD LO evolutor
109  for (j = 0; j < nops; j++)
110  {
111  if (fabs(term = evec(a, i) * M1(i, j) * evec_i(j, b)) > EPS)
112  vM1vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j}, term)); // QCD NLO evolutor
113  if (fabs(term = evec(a, i) * M2(i, j) * evec_i(j, b)) > EPS)
114  vM2vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j}, term)); // QCD NNLO evolutor
115  for (p = 0; p < nops; p++)
116  if (fabs(term = evec(a, i) * M1(i, p) * M1(p, j) * evec_i(j, b)) > EPS)
117  vM11vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p}, term)); // QCD NNLO evolutor
118  }
119  }
120  }
121 
122  if (order_qed != QED0)
123  {
124  b0e = model.Beta_e(00, nf);
125  b3 = model.Beta_s(01, nf) / 2. / b0 / b0e;
126  b4 = model.Beta_s(11, nf) / 4. / b0 / b0 / b0e - 2. * b1*b3;
127  W01 = AnomalousDimension(01, nu, nd).transpose() / 2. / b0e;
128  W02 = AnomalousDimension(02, nu, nd).transpose() / 4. / b0e / b0e;
129  W11 = AnomalousDimension(11, nu, nd).transpose() / 4. / b0 / b0e;
130  W21 = AnomalousDimension(21, nu, nd).transpose() / 8. / b0 / b0 / b0e;
131 
132  // QED magic numbers
133  // M3: B(1)_01, B(2)_02, B(1)_02, R5_12, M4: B(0)_11, B(0)_12, M5: B(-1)_21, M6: B(1)_12 - proper powers of omega and lambda added in DF1Evol
134  M3 = evec_i * W01 * evec;
135  M4 = evec_i * (W11 - b1 * W01 - b3 * W10) * evec;
136  M5 = evec_i * (W21 - b1 * W11 - b2 * W01 - b3 * W20 - b4 * W10) * evec;
137  M6 = evec_i * (W02 + W11 - (b1 + b3) * W01 - b3 * W10) * evec;
138  for (a = 0; a < nops; a++)
139  for (b = 0; b < nops; b++)
140  for (i = 0; i < nops; i++)
141  for (j = 0; j < nops; j++)
142  {
143  if (fabs(term = evec(a, i) * M3(i, j) * evec_i(j, b)) > EPS)
144  vM3vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j}, term));
145  if (fabs(term = evec(a, i) * M4(i, j) * evec_i(j, b)) > EPS)
146  vM4vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j}, term));
147  if (fabs(term = evec(a, i) * M5(i, j) * evec_i(j, b)) > EPS)
148  vM5vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j}, term));
149  if (fabs(term = evec(a, i) * M6(i, j) * evec_i(j, b)) > EPS)
150  vM6vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j}, term));
151  for (p = 0; p < nops; p++)
152  {
153  if (fabs(term = evec(a, i) * M3(i, p) * M3(p, j) * evec_i(j, b)) > EPS)
154  vM33vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p}, term));
155  if (fabs(term = evec(a, i) * M3(i, p) * M1(p, j) * evec_i(j, b)) > EPS)
156  vM31vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p}, term));
157  if (fabs(term = evec(a, i) * M1(i, p) * M3(p, j) * evec_i(j, b)) > EPS)
158  vM13vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p}, term));
159  if (fabs(term = evec(a, i) * M3(i, p) * M4(p, j) * evec_i(j, b)) > EPS)
160  vM34vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p}, term));
161  if (fabs(term = evec(a, i) * M4(i, p) * M3(p, j) * evec_i(j, b)) > EPS)
162  vM43vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p}, term));
163  if (fabs(term = evec(a, i) * M2(i, p) * M3(p, j) * evec_i(j, b)) > EPS)
164  vM23vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p}, term));
165  if (fabs(term = evec(a, i) * M3(i, p) * M2(p, j) * evec_i(j, b)) > EPS)
166  vM32vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p}, term));
167  if (fabs(term = evec(a, i) * M1(i, p) * M4(p, j) * evec_i(j, b)) > EPS)
168  vM14vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p}, term));
169  if (fabs(term = evec(a, i) * M4(i, p) * M1(p, j) * evec_i(j, b)) > EPS)
170  vM41vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p}, term));
171  for (q = 0; q < nops; q++)
172  {
173  if (fabs(term = evec(a, i) * M1(i, p) * M1(p, q) * M3(q, j) * evec_i(j, b)) > EPS)
174  vM113vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p, q}, term));
175  if (fabs(term = evec(a, i) * M1(i, p) * M3(p, q) * M1(q, j) * evec_i(j, b)) > EPS)
176  vM131vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p, q}, term));
177  if (fabs(term = evec(a, i) * M3(i, p) * M1(p, q) * M1(q, j) * evec_i(j, b)) > EPS)
178  vM311vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p, q}, term));
179  if (fabs(term = evec(a, i) * M3(i, p) * M3(p, q) * M1(q, j) * evec_i(j, b)) > EPS)
180  vM331vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p, q}, term));
181  if (fabs(term = evec(a, i) * M3(i, p) * M1(p, q) * M3(q, j) * evec_i(j, b)) > EPS)
182  vM313vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p, q}, term));
183  if (fabs(term = evec(a, i) * M1(i, p) * M3(p, q) * M3(q, j) * evec_i(j, b)) > EPS)
184  vM133vi[nnf].insert(std::pair<std::vector<uint>, double > ({a, b, i, j, p, q}, term));
185  }
186  }
187  }
188  }
189  }
190 }
191 
192 EvolDF1::~EvolDF1()
193 {
194 }
195 
196 /* Delta F = 1 anomalous dimension in Misiak basis, in the effective basis (C7eff, C8eff)
197  ref. for CC,CP,PP QED NLO, QP,QQ,BP,BB QCD NLO, CL,PL NNLO, PQ,QL,LL,BL NLO: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
198  ref. for CC,CP,PP QCD NNLO with nf: Gorbahn, Haisch, Nucl. Phys. B 713, 291, hep-ph/0411071
199  ref. for CM,PM QCD NNLO with nf: Czakon, Haisch, Misiak, JHEP 0703, 008, hep-ph/0612329
200  ref. for MM QCD NNLO with nf: Gorbahn, Haisch, Misiak, Phys. Rev. Lett. 95, 102004, hep-ph/0504194
201  ref. for CM,PM QED LO, QM QCD LO: Baranowski, Misiak, Phys. Lett. B 483, 410, hep-ph/9907427
202  ref. for MM,QP,QQ,BP,BB QED NLO, BQ NLO: Bobeth, Gambino, Gorbahn, Haisch, JHEP 0404, 071, hep-ph/0312090
203  QM QED,BM??? (in 031209 C7,8 are normalised with alpha_s and not in the effective basis)
204  */
205 
206 double EvolDF1::f_f(uint nnf, uint i, uint j, int k, double eta)
207 {
208  for (int il = 0; il < F_iCacheSize; il++)
209  if (f_f_c[0][il] == (int) nnf && f_f_c[1][il] == (int) i && f_f_c[2][il] == (int) j &&
210  f_f_c[3][il] == k && f_f_d[0][il] == eta)
211  return f_f_d[1][il];
212 
213  double aii = ai[nnf].at(i), aij = ai[nnf].at(j);
214  double den = aij + k - aii, ret;
215 
216  if (fabs(den) < EPS)
217  ret = pow(eta, aii) * log(eta);
218  else
219  ret = (pow(eta, aij + k) - pow(eta, aii)) / den;
220 
221  model.CacheShift(f_f_c, 4);
222  model.CacheShift(f_f_d, 2);
223  f_f_c[0][0] = (int) nnf;
224  f_f_c[1][0] = (int) i;
225  f_f_c[2][0] = (int) j;
226  f_f_c[3][0] = k;
227  f_f_d[0][0] = eta;
228  f_f_d[1][0] = ret;
229 
230  return ret;
231 }
232 
233 double EvolDF1::f_r(uint nnf, uint i, uint j, int k, double eta)
234 {
235  double ll = log(eta), den = ai[nnf].at(j) + k - ai[nnf].at(i);
236 
237  if (fabs(den) < EPS)
238  return (0.5 * pow(eta, ai[nnf].at(i)) * ll * ll);
239  else
240  return ((pow(eta, ai[nnf].at(j) + k) * ll - f_f(nnf, i, j, k, eta)) / den);
241 }
242 
243 double EvolDF1::f_g(uint nnf, uint i, uint p, uint j, int k, int l, double eta)
244 {
245  double den = ai[nnf].at(j) + l - ai[nnf].at(p);
246 
247  if (fabs(den) < EPS)
248  return (f_r(nnf, i, p, k, eta));
249  else
250  return ((f_f(nnf, i, j, k + l, eta) - f_f(nnf, i, p, k, eta)) / den);
251 }
252 
253 double EvolDF1::f_h(uint nnf, uint i, uint p, uint q, uint j, int k, int l, int m, double eta)
254 {
255  double ll = log(eta), den1 = ai[nnf].at(j) + m - ai[nnf].at(q),
256  den2 = ai[nnf].at(q) + l - ai[nnf].at(p),
257  den3 = ai[nnf].at(p) + k - ai[nnf].at(i);
258 
259  if (fabs(den1) < EPS && fabs(den2) < EPS && fabs(den3) < EPS)
260  return (pow(eta, ai[nnf].at(i)) * ll * ll * ll / 6.);
261  else if (fabs(den1) < EPS && fabs(den2) < EPS)
262  return ((0.5 * pow(eta, ai[nnf].at(p) + k) * ll * ll - f_r(nnf, i, p, k, eta)) / den3);
263  else if (fabs(den1) < EPS)
264  return ((f_r(nnf, i, q, k + l, eta) - f_g(nnf, i, p, q, k, l, eta)) / den2);
265  else
266  return ((f_g(nnf, i, p, j, k, l + m, eta) - f_g(nnf, i, p, q, k, l, eta)) / den1);
267 }
268 
269 void EvolDF1::CheckNf(indices nm, uint nf) const
270 {
271  if (nm % 10 == 0)
272  {
273  if (!(nf == 3 || nf == 4 || nf == 5 || nf == 6))
274  throw std::runtime_error("EvolDF1::CheckNf(): Wrong number of flavours in anoumalous dimensions");
275  } else if (nf != 5)
276  throw std::runtime_error("EvolDF1::CheckNf(): Wrong number of flavours in anoumalous dimensions");
277 }
278 
279 gslpp::matrix<double> EvolDF1::GammaCC(indices nm, uint n_u, uint n_d) const
280 {
281  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
282  CheckNf(nm, nf);
283 
284  gslpp::matrix<double> gammaDF1(2, 2, 0.);
285  double z3 = gslpp_special_functions::zeta(3);
286 
287  switch (nm)
288  {
289  // QCD
290  // ref.: Gorbahn, Haisch, Nucl. Phys. B 713, 291, hep-ph/0411071
291  case 10:
292  gammaDF1(0, 0) = -4.;
293  gammaDF1(0, 1) = 8. / 3.;
294  gammaDF1(1, 0) = 12.;
295  break;
296  case 20:
297  gammaDF1(0, 0) = -145. / 3. + nf * 16. / 9.; // -355./9.
298  gammaDF1(0, 1) = -26. + nf * 40. / 27.; // -502./27.
299  gammaDF1(1, 0) = -45. + nf * 20. / 3.; // -35./3.
300  gammaDF1(1, 1) = -28. / 3.;
301  break;
302  case 30:
303  gammaDF1(0, 0) = -1927. / 2. + nf * 257. / 9. + nf * nf * 40. / 9. + z3 * (224. + nf * 160. / 3.); // -12773./18. + z3*1472./3.
304  gammaDF1(0, 1) = 475. / 9. + nf * 362. / 27. - nf * nf * 40. / 27. - z3 * (896. / 3. + nf * 320. / 9.); // 745./9. - z3*4288./9.
305  gammaDF1(1, 0) = 307. / 2. + nf * 361. / 3. - nf * nf * 20. / 3. - z3 * (1344. + nf * 160.); // 1177./2. - z3*2144.
306  gammaDF1(1, 1) = 1298. / 3. - nf * 76. / 3. - z3 * 224.; // 306. + z3*224.
307  break;
308  // QED
309  // only available for nf = 5
310  // ref.: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
311  case 01:
312  gammaDF1(0, 0) = -8. / 3.;
313  gammaDF1(1, 1) = -8. / 3.;
314  break;
315  case 11:
316  gammaDF1(0, 0) = 169. / 9.;
317  gammaDF1(0, 1) = 100. / 27.;
318  gammaDF1(1, 0) = 50. / 3.;
319  gammaDF1(1, 1) = -8. / 3.;
320  break;
321  case 21:
322  case 02:
323  break;
324  default:
325  throw std::runtime_error("EvolDF1::GammaCC(): order not implemented");
326  }
327  return (gammaDF1);
328 }
329 
330 gslpp::matrix<double> EvolDF1::GammaCP(indices nm, uint n_u, uint n_d) const
331 {
332  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
333  CheckNf(nm, nf);
334 
335  gslpp::matrix<double> gammaDF1(2, 4, 0.);
336  double z3 = gslpp_special_functions::zeta(3);
337 
338  switch (nm)
339  {
340  // QCD
341  // ref.: Gorbahn, Haisch, Nucl. Phys. B 713, 291, hep-ph/0411071
342  case 10:
343  gammaDF1(0, 1) = -2. / 9.;
344  gammaDF1(1, 1) = 4. / 3.;
345  break;
346  case 20:
347  gammaDF1(0, 0) = -1412. / 243.;
348  gammaDF1(0, 1) = -1369. / 243.;
349  gammaDF1(0, 2) = 134. / 243.;
350  gammaDF1(0, 3) = -35. / 162.;
351  gammaDF1(1, 0) = -416. / 81.;
352  gammaDF1(1, 1) = 1280. / 81.;
353  gammaDF1(1, 2) = 56. / 81.;
354  gammaDF1(1, 3) = 35. / 27.;
355  break;
356  case 30:
357  gammaDF1(0, 0) = 269107. / 13122. - nf * 2288. / 729. - z3 * 1360. / 81.; // 63187./13122. - z3*1360./81.
358  gammaDF1(0, 1) = -2425817. / 13122. + nf * 30815. / 4374. - z3 * 776. / 81.; // -981796./6561. - z3*776./81.
359  gammaDF1(0, 2) = -343783. / 52488. + nf * 392. / 729. + z3 * 124. / 81.; // -202663./52488. + z3*124./81.
360  gammaDF1(0, 3) = -37573. / 69984. + nf * 35. / 972. + z3 * 100. / 27.; // -24973./69984. + z3*100./27.
361  gammaDF1(1, 0) = 69797. / 2187. + nf * 904. / 243. + z3 * 2720. / 27.; // 110477./2187. + z3*2720./.27.
362  gammaDF1(1, 1) = 1457549. / 8748. - nf * 22067. / 729. - z3 * 2768. / 27.; // 133529./8748. - z3*2768./27.
363  gammaDF1(1, 2) = -37889. / 8748. - nf * 28. / 243. - z3 * 248. / 27.; // -42929./8748. - z3*248./27.
364  gammaDF1(1, 3) = 366919. / 11664. - nf * 35. / 162. - z3 * 110. / 9.; // 354319./11664. - z3*110./9.
365  break;
366  // QED
367  // only available for nf = 5
368  // ref.: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
369  case 01:
370  break;
371  case 11:
372  gammaDF1(0, 1) = 254. / 729.;
373  gammaDF1(1, 1) = 1076. / 243.;
374  break;
375  case 21:
376  case 02:
377  break;
378  default:
379  throw std::runtime_error("EvolDF1::GammaCP(): order not implemented");
380  }
381  return (gammaDF1);
382 }
383 
384 gslpp::matrix<double> EvolDF1::GammaCM(indices nm, uint n_u, uint n_d) const
385 {
386  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
387  CheckNf(nm, nf);
388 
389  gslpp::matrix<double> gammaDF1(2, 2, 0.);
390  double Qu = 2. / 3.;
391  double Qd = -1. / 3.;
392  double Qbar = n_u * Qu + n_d*Qd;
393  double z3 = gslpp_special_functions::zeta(3);
394 
395  switch (nm)
396  {
397  // QCD
398  // ref.: Czakon, Haisch, Misiak, JHEP 0703, 008, hep-ph/0612329
399  case 10:
400  gammaDF1(0, 0) = 8. / 243. - Qu * 4. / 3.;
401  gammaDF1(0, 1) = 173. / 162.;
402  gammaDF1(1, 0) = -16. / 81. + Qu * 8.;
403  gammaDF1(1, 1) = 70. / 27.;
404  break;
405  case 20:
406  gammaDF1(0, 0) = 12614. / 2187. - nf * 64. / 2187. - Qu * 374. / 27. + nf * Qu * 2. / 27.;
407  gammaDF1(0, 1) = 65867. / 5832. + nf * 431. / 5832.;
408  gammaDF1(1, 0) = -2332. / 729. + nf * 128. / 729. + Qu * 136. / 9. - nf * Qu * 4. / 9.;
409  gammaDF1(1, 1) = 10577. / 486. - nf * 917. / 972.;
410  break;
411  case 30:
412  gammaDF1(0, 0) = 77506102. / 531441. - nf * 875374. / 177147. + nf * nf * 560. / 19683. - Qu * 9731. / 162. +
413  nf * Qu * 11045. / 729. + nf * nf * Qu * 316. / 729. + Qbar * 3695. / 486. + z3 * (-112216. / 6561. + nf * 728. / 729. +
414  Qu * 25508. / 81. - nf * Qu * 64. / 81. - Qbar * 100. / 27.);
415  gammaDF1(0, 1) = -421272953. / 1417176. - nf * 8210077. / 472392. - nf * nf * 1955. / 6561. + z3 * (-953042. / 2187. -
416  nf * 10381. / 486.);
417  gammaDF1(1, 0) = -15463055. / 177147. + nf * 242204. / 59049. - nf * nf * 1120. / 6561. + Qu * 55748. / 27. -
418  nf * Qu * 33970. / 243. - nf * nf * Qu * 632. / 243. - Qbar * 3695. / 81. + z3 * (365696. / 2187. - nf * 1168. / 243. -
419  Qu * 51232. / 27. - nf * Qu * 1024. / 27. + Qbar * 200. / 9.);
420  gammaDF1(1, 1) = 98548513. / 472392. - nf * 5615165. / 78732. - nf * nf * 2489. / 2187. + z3 * (-607103. / 729. -
421  nf * 1679. / 81.);
422  break;
423  // QED
424  // only available for nf = 5
425  // ref.: Baranowski, Misiak, Phys. Lett. B 483, 410, hep-ph/9907427
426  case 01:
427  gammaDF1(0, 0) = -832. / 729.;
428  gammaDF1(0, 1) = 22. / 243.;
429  gammaDF1(1, 0) = -208. / 243.;
430  gammaDF1(1, 1) = -116. / 81.;
431  break;
432  case 11:
433  case 21:
434  case 02:
435  break;
436  default:
437  throw std::runtime_error("EvolDF1::GammaCM(): order not implemented");
438  }
439  return (gammaDF1);
440 }
441 
442 gslpp::matrix<double> EvolDF1::GammaCL(indices nm, uint n_u, uint n_d) const
443 {
444  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
445  if (nf != 5)
446  throw std::runtime_error("EvolDF1::GammaCL(): Wrong number of flavours in anoumalous dimensions");
447 
448 
449  gslpp::matrix<double> gammaDF1(2, 2, 0.);
450  double z3 = gslpp_special_functions::zeta(3);
451 
452  switch (nm)
453  {
454  // QED
455  // ref.: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
456  case 01:
457  gammaDF1(0, 0) = -32. / 27.;
458  gammaDF1(1, 0) = -8. / 9.;
459  break;
460  case 11:
461  gammaDF1(0, 0) = -2272. / 729.;
462  gammaDF1(1, 0) = 1952. / 243.;
463  break;
464  case 21:
465  gammaDF1(0, 0) = -1359190. / 19683. + z3 * 6976. / 243.;
466  gammaDF1(1, 0) = -229696. / 6561. - z3 * 3584. / 81.;
467  break;
468  case 02:
469  gammaDF1(0, 0) = -11680. / 2187.;
470  gammaDF1(1, 0) = -2920. / 729.;
471  gammaDF1(0, 1) = -416. / 81.;
472  gammaDF1(1, 1) = -104. / 27.;
473  break;
474  case 10:
475  case 20:
476  case 30:
477  break;
478  default:
479  throw std::runtime_error("EvolDF1::GammaCL(): order not implemented");
480  }
481  return (gammaDF1);
482 }
483 
484 gslpp::matrix<double> EvolDF1::GammaCQ(indices nm, uint n_u, uint n_d) const
485 {
486  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
487  if (nf != 5)
488  throw std::runtime_error("EvolDF1::GammaCQ(): Wrong number of flavours in anoumalous dimensions");
489 
490 
491  gslpp::matrix<double> gammaDF1(2, 4, 0.);
492 
493  switch (nm)
494  {
495  // QED
496  // ref.: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
497  case 01:
498  gammaDF1(0, 0) = 32. / 27.;
499  gammaDF1(1, 0) = 8. / 9.;
500  break;
501  case 11:
502  gammaDF1(0, 0) = 2272. / 729.;
503  gammaDF1(0, 1) = 122. / 81.;
504  gammaDF1(0, 3) = 49. / 81.;
505  gammaDF1(1, 0) = -1952. / 243.;
506  gammaDF1(1, 1) = -748. / 27.;
507  gammaDF1(1, 3) = 82. / 27.;
508  break;
509  case 10:
510  case 20:
511  case 30:
512  case 21:
513  case 02:
514  break;
515  default:
516  throw std::runtime_error("EvolDF1::GammaCQ(): order not implemented");
517  }
518  return (gammaDF1);
519 }
520 
521 gslpp::matrix<double> EvolDF1::GammaPP(indices nm, uint n_u, uint n_d) const
522 {
523  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
524  CheckNf(nm, nf);
525 
526  gslpp::matrix<double> gammaDF1(4, 4, 0.);
527  double z3 = gslpp_special_functions::zeta(3);
528 
529  switch (nm)
530  {
531  // QCD
532  // ref.: Gorbahn, Haisch, Nucl. Phys. B 713, 291, hep-ph/0411071
533  case 10:
534  gammaDF1(0, 1) = -52. / 3.;
535  gammaDF1(0, 3) = 2.;
536  gammaDF1(1, 0) = -40. / 9.;
537  gammaDF1(1, 1) = -160. / 9. + (double) nf * 4. / 3.; // -100./9.
538  gammaDF1(1, 2) = 4. / 9.;
539  gammaDF1(1, 3) = 5. / 6.;
540  gammaDF1(2, 1) = -256. / 3.;
541  gammaDF1(2, 3) = 20.;
542  gammaDF1(3, 0) = -256. / 9.;
543  gammaDF1(3, 1) = -544. / 9. + (double) nf * 40. / 3.; // 56./9.
544  gammaDF1(3, 2) = 40. / 9.;
545  gammaDF1(3, 3) = -2. / 3.;
546  break;
547  case 20:
548  gammaDF1(0, 0) = -4468. / 81.;
549  gammaDF1(0, 1) = -29129. / 81. - nf * 52. / 9.; // -31469./81.
550  gammaDF1(0, 2) = 400. / 81.;
551  gammaDF1(0, 3) = 3493. / 108. - nf * 2. / 9.; // 3373./108.
552  gammaDF1(1, 0) = -13678. / 243. + nf * 368. / 81.; // -8158./243.
553  gammaDF1(1, 1) = -79409. / 243. + nf * 1334. / 81.; // -59399./243.
554  gammaDF1(1, 2) = 509. / 486. - nf * 8. / 81.; // 269./486.
555  gammaDF1(1, 3) = 13499. / 648. - nf * 5. / 27.; // 12899./648.
556  gammaDF1(2, 0) = -244480. / 81. - nf * 160. / 9.; // -251680./81.
557  gammaDF1(2, 1) = -29648. / 81. - nf * 2200. / 9.; // -128648./81.
558  gammaDF1(2, 2) = 23116. / 81. + nf * 16. / 9.; // 23836./81.
559  gammaDF1(2, 3) = 3886. / 27. + nf * 148. / 9.; // 6106./27.
560  gammaDF1(3, 0) = 77600. / 243. - nf * 1264. / 81.; // 58640./243.
561  gammaDF1(3, 1) = -28808. / 243. + nf * 164. / 81.; // -26348./243.
562  gammaDF1(3, 2) = -20324. / 243. + nf * 400. / 81.; // -14324./243.
563  gammaDF1(3, 3) = -21211. / 162. + nf * 622. / 27.; // -2551./162.
564  break;
565  case 30:
566  gammaDF1(0, 0) = -4203068. / 2187. + nf * 14012. / 243. - z3 * 608. / 27.; // -3572528./2187. - z3*608./27.
567  gammaDF1(0, 1) = -18422762. / 2187. + nf * 888605. / 2916. + nf * nf * 272. / 27. + z3 * (39824. / 27. + nf * 160.); // -58158773./8748. + z3*61424./27.
568  gammaDF1(0, 2) = 674281. / 4374. - nf * 1352. / 243. - z3 * 496. / 27.; // 552601./4374. - z3*496./27.
569  gammaDF1(0, 3) = 9284531. / 11664. - nf * 2798. / 81. - nf * nf * 26. / 27. - z3 * (1921. / 9. + nf * 20.); // 6989171./11664. - z3*2821./9.
570  gammaDF1(1, 0) = -5875184. / 6561. + nf * 217892. / 2187. + nf * nf * 472. / 81. + z3 * (27520. / 81. + nf * 1360. / 9.); // -1651004./6561. + z3*88720./81.
571  gammaDF1(1, 1) = -70274587. / 13122. + nf * 8860733. / 17496. - nf * nf * 4010. / 729. + z3 * (16592. / 81. + nf * 2512. / 27.); // -155405353./52488. + z3*54272./81.
572  gammaDF1(1, 2) = 2951809. / 52488. - nf * 31175. / 8748. - nf * nf * 52. / 81. - z3 * (3154. / 81. + nf * 136. / 9.); // 1174159./52488. - z3*9274./81.
573  gammaDF1(1, 3) = 3227801. / 8748. - nf * 105293. / 11664. - nf * nf * 65. / 54. + z3 * (200. / 27. - nf * 220. / 9.); // 10278809./34992. - z3*3100./27.
574  gammaDF1(2, 0) = -194951552. / 2187. + nf * 358672. / 81. - nf * nf * 2144. / 81. + z3 * 87040. / 27.; // -147978032./2187. + z3*87040./27.
575  gammaDF1(2, 1) = -130500332. / 2187. - nf * 2949616. / 729. + nf * nf * 3088. / 27. + z3 * (238016. / 27. + nf * 640.); // -168491372./2187. + z3*324416./27.
576  gammaDF1(2, 2) = 14732222. / 2187. - nf * 27428. / 81. + nf * nf * 272. / 81. - z3 * 13984. / 27.; // 11213042./2187. - z3*13984./27.
577  gammaDF1(2, 3) = 16521659. / 2916. + nf * 8081. / 54. - nf * nf * 316. / 27. - z3 * (22420. / 9. + nf * 200.); // 17850329./2916. - z3*31420./9.
578  gammaDF1(3, 0) = 162733912. / 6561. - nf * 2535466. / 2187. + nf * nf * 17920. / 243. + z3 * (174208. / 81. + nf * 12160. / 9.); // 136797922./6561. + z3*721408./81.
579  gammaDF1(3, 1) = 13286236. / 6561. - nf * 1826023. / 4374. - nf * nf * 159548. / 729. - z3 * (24832. / 81. + nf * 9440. / 27.); // -72614473./13122. - z3*166432./81.
580  gammaDF1(3, 2) = -22191107. / 13122. + nf * 395783. / 4374. - nf * nf * 1720. / 243. - z3 * (33832. / 81. + nf * 1360. / 9.); // -9288181./6561. - z3*95032./81.
581  gammaDF1(3, 3) = -32043361. / 8748. + nf * 3353393. / 5832. - nf * nf * 533. / 81. + z3 * (9248. / 27. - nf * 1120. / 9.); // -16664027./17496. - z3*7552./27.
582  break;
583  // QED
584  // only available for nf = 5
585  // ref.: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
586  case 01:
587  break;
588  case 11:
589  gammaDF1(0, 1) = 11116. / 243.;
590  gammaDF1(0, 3) = -14. / 3.;
591  gammaDF1(1, 0) = 280. / 27.;
592  gammaDF1(1, 1) = 18763. / 729.;
593  gammaDF1(1, 2) = -28. / 27.;
594  gammaDF1(1, 3) = -35. / 18.;
595  gammaDF1(2, 1) = 111136. / 243.;
596  gammaDF1(2, 3) = -140. / 3.;
597  gammaDF1(3, 0) = 2944. / 27.;
598  gammaDF1(3, 1) = 193312. / 729.;
599  gammaDF1(3, 2) = -280. / 27.;
600  gammaDF1(3, 3) = -175. / 9.;
601  break;
602  case 21:
603  case 02:
604  break;
605  default:
606  throw std::runtime_error("EvolDF1::GammaPP(): order not implemented");
607  }
608  return (gammaDF1);
609 }
610 
611 gslpp::matrix<double> EvolDF1::GammaPM(indices nm, uint n_u, uint n_d) const
612 {
613  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
614  CheckNf(nm, nf);
615 
616  gslpp::matrix<double> gammaDF1(4, 2, 0.);
617  double Qu = 2. / 3.;
618  double Qd = -1. / 3.;
619  double Qbar = n_u * Qu + n_d*Qd;
620  double z3 = gslpp_special_functions::zeta(3);
621 
622  switch (nm)
623  {
624  // QCD
625  // ref.: Czakon, Haisch, Misiak, JHEP 0703, 008, hep-ph/0612329
626  case 10:
627  gammaDF1(0, 0) = -176. / 81.;
628  gammaDF1(0, 1) = 14. / 27.;
629  gammaDF1(1, 0) = 88. / 243. - nf * 16. / 81.;
630  gammaDF1(1, 1) = 74. / 81. - nf * 49. / 54.;
631  gammaDF1(2, 0) = -6272. / 81.;
632  gammaDF1(2, 1) = 1736. / 27. + nf * 36.;
633  gammaDF1(3, 0) = 3136. / 243. - nf * 160. / 81. + Qbar * 48.;
634  gammaDF1(3, 1) = 2372. / 81. + nf * 160. / 27.;
635  break;
636  case 20:
637  gammaDF1(0, 0) = 97876. / 729. - nf * 4352. / 729. - Qbar * 112. / 3.;
638  gammaDF1(0, 1) = 42524. / 243. - nf * 2398. / 243.;
639  gammaDF1(1, 0) = -70376. / 2187. - nf * 15788. / 2187. + nf * nf * 32. / 729. - Qbar * 140. / 9.;
640  gammaDF1(1, 1) = -159718. / 729. - nf * 39719. / 5832. - nf * nf * 253. / 486.;
641  gammaDF1(2, 0) = 1764752. / 729. - nf * 65408. / 729. - Qbar * 3136. / 3.;
642  gammaDF1(2, 1) = 2281576. / 243. + nf * 140954. / 243. - nf * nf * 14.;
643  gammaDF1(3, 0) = 4193840. / 2187. - nf * 324128. / 2187. + nf * nf * 896. / 729. - Qbar * 1136. / 9. - nf * Qbar * 56. / 3.;
644  gammaDF1(3, 1) = -3031517. / 729. - nf * 15431. / 1458. - nf * nf * 6031. / 486.;
645  break;
646  case 30:
647  gammaDF1(0, 0) = 102439553. / 177147. - nf * 12273398 / 59049. + nf * nf * 5824. / 6561. + Qbar * 26639. / 81. - nf * Qbar * 8. / 27. +
648  z3 * (3508864. / 2187. - nf * 1904. / 243. - Qbar * 1984. / 9. - nf * Qbar * 64. / 9.);
649  gammaDF1(0, 1) = 3205172129. / 472392. - nf * 108963529. / 314928. + nf * nf * 58903. / 4374. + z3 * (-1597588. / 729. +
650  nf * 13028. / 81. - nf * nf * 20. / 9.);
651  gammaDF1(1, 0) = -2493414077. / 1062882. - nf * 9901031. / 354294. + nf * nf * 243872. / 59049. - nf * nf * nf * 1184. / 6561. -
652  Qbar * 49993. / 972. + nf * Qbar * 305. / 27. + z3 * (-1922264. / 6561. + nf * 308648. / 2187. - nf * nf * 1280. / 243. +
653  Qbar * 1010. / 9. - nf * Qbar * 200. / 27.);
654  gammaDF1(1, 1) = -6678822461. / 2834352. + nf * 127999025. / 1889568. + nf * nf * 1699073. / 157464. + nf * nf * nf * 505. / 4374. +
655  z3 * (2312684. / 2187. + nf * 128347. / 729. + nf * nf * 920. / 81.);
656  gammaDF1(2, 0) = 8808397748. / 177147. - nf * 174839456. / 59049. + nf * nf * 1600. / 729. - Qbar * 669694. / 81. + nf * Qbar * 10672. / 27. +
657  z3 * (123543040. / 2187. - nf * 207712. / 243. + nf * nf * 128. / 27. - Qbar * 24880. / 9. - nf * Qbar * 640. / 9.);
658  gammaDF1(2, 1) = 29013624461. / 118098. - nf * 64260772. / 19683. - nf * nf * 230962. / 243. - nf * nf * nf * 148. / 27. +
659  z3 * (-69359224. / 729. - nf * 885356. / 81. - nf * nf * 5080. / 9.);
660  gammaDF1(3, 0) = 7684242746. / 531441. - nf * 351775414. / 177147. - nf * nf * 479776. / 59049. - nf * nf * nf * 11456. / 6561. +
661  Qbar * 3950201. / 243. - nf * Qbar * 130538. / 81. - nf * nf * Qbar * 592. / 81. + z3 * (7699264. / 6561. + nf * 2854976. / 2187. -
662  nf * nf * 12320. / 243. - Qbar * 108584. / 9. - nf * Qbar * 1136. / 27.);
663  gammaDF1(3, 1) = -72810260309. / 708588. + nf * 2545824851. / 472392. - nf * nf * 33778271. / 78732. - nf * nf * nf * 3988. / 2187. +
664  z3 * (-61384768. / 2187. - nf * 685472. / 729. + nf * nf * 350. / 81.);
665  break;
666  // QED
667  // only available for nf = 5
668  // ref.: Baranowski, Misiak, Phys. Lett. B 483, 410, hep-ph/9907427
669  case 01:
670  gammaDF1(0, 0) = -20. / 243.;
671  gammaDF1(0, 1) = 20. / 81.;
672  gammaDF1(1, 0) = -176. / 729.;
673  gammaDF1(1, 1) = 14. / 243.;
674  gammaDF1(2, 0) = -22712. / 243.;
675  gammaDF1(2, 1) = 1328. / 81.;
676  gammaDF1(3, 0) = -6272. / 729.;
677  gammaDF1(3, 1) = -1180. / 243.;
678  break;
679  case 11:
680  case 21:
681  case 02:
682  break;
683  default:
684  throw std::runtime_error("EvolDF1::GammaPM(): order not implemented");
685  }
686  return (gammaDF1);
687 }
688 
689 gslpp::matrix<double> EvolDF1::GammaPL(indices nm, uint n_u, uint n_d) const
690 {
691  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
692  if (nf != 5)
693  throw std::runtime_error("EvolDF1::GammaPL(): Wrong number of flavours in anoumalous dimensions");
694 
695 
696  gslpp::matrix<double> gammaDF1(4, 2, 0.);
697  double z3 = gslpp_special_functions::zeta(3);
698 
699  switch (nm)
700  {
701  // QED
702  // ref.: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
703 
704  case 01:
705  gammaDF1(0, 0) = -16. / 9.;
706  gammaDF1(1, 0) = 32. / 27.;
707  gammaDF1(2, 0) = -112. / 9.;
708  gammaDF1(3, 0) = 512. / 27.;
709  break;
710  case 11:
711  gammaDF1(0, 0) = -6752. / 243.;
712  gammaDF1(1, 0) = -2192. / 729.;
713  gammaDF1(2, 0) = -84032. / 243.;
714  gammaDF1(3, 0) = -37856. / 729.;
715  break;
716  case 21:
717  gammaDF1(0, 0) = -1290092. / 6561. + z3 * 3200. / 81.;
718  gammaDF1(1, 0) = -819971. / 19683. - z3 * 19936. / 243.;
719  gammaDF1(2, 0) = -16821944. / 6561. + z3 * 30464. / 81.;
720  gammaDF1(3, 0) = -17787368. / 19683. - z3 * 286720. / 243.;
721  break;
722  case 02:
723  gammaDF1(0, 0) = -39752. / 729.;
724  gammaDF1(1, 0) = 1024. / 2187.;
725  gammaDF1(2, 0) = -381344. / 729.;
726  gammaDF1(3, 0) = 24832. / 2187.;
727  gammaDF1(0, 1) = -136. / 27.;
728  gammaDF1(1, 1) = -448. / 81.;
729  gammaDF1(2, 1) = -15616. / 27.;
730  gammaDF1(3, 1) = -7936. / 81.;
731  break;
732  case 10:
733  case 20:
734  case 30:
735  break;
736  default:
737  throw std::runtime_error("EvolDF1::GammaPL(): order not implemented");
738  }
739  return (gammaDF1);
740 }
741 
742 gslpp::matrix<double> EvolDF1::GammaPQ(indices nm, uint n_u, uint n_d) const
743 {
744  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
745  if (nf != 5)
746  throw std::runtime_error("EvolDF1::GammaPQ(): Wrong number of flavours in anoumalous dimensions");
747 
748  gslpp::matrix<double> gammaDF1(4, 4, 0.);
749 
750  switch (nm)
751  {
752  // QED
753  // ref.: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
754  case 01:
755  gammaDF1(0, 0) = 76. / 9.;
756  gammaDF1(0, 2) = -2. / 3.;
757  gammaDF1(1, 0) = -32. / 27.;
758  gammaDF1(1, 1) = 20. / 3.;
759  gammaDF1(1, 3) = -2. / 3.;
760  gammaDF1(2, 0) = 496. / 9.;
761  gammaDF1(2, 2) = -20. / 3.;
762  gammaDF1(3, 0) = -512. / 27.;
763  gammaDF1(3, 1) = 128. / 3.;
764  gammaDF1(3, 3) = -20. / 3.;
765  break;
766  case 11:
767  gammaDF1(0, 0) = -23488. / 243.;
768  gammaDF1(0, 1) = 6280. / 27.;
769  gammaDF1(0, 2) = 112. / 9.;
770  gammaDF1(0, 3) = -538. / 27.;
771  gammaDF1(1, 0) = 31568. / 729.;
772  gammaDF1(1, 1) = 9481. / 81.;
773  gammaDF1(1, 2) = -92. / 27.;
774  gammaDF1(1, 3) = -1012. / 81.;
775  gammaDF1(2, 0) = -233920. / 243.;
776  gammaDF1(2, 1) = 68848. / 27.;
777  gammaDF1(2, 2) = 1120. / 9.;
778  gammaDF1(2, 3) = -5056. / 27.;
779  gammaDF1(3, 0) = 352352. / 729.;
780  gammaDF1(3, 1) = 116680. / 81.;
781  gammaDF1(3, 2) = -752. / 27.;
782  gammaDF1(3, 3) = -10147. / 81.;
783  break;
784  case 10:
785  case 20:
786  case 30:
787  case 21:
788  case 02:
789  break;
790  default:
791  throw std::runtime_error("EvolDF1::GammaPQ(): order not implemented");
792  }
793  return (gammaDF1);
794 }
795 
796 gslpp::matrix<double> EvolDF1::GammaMM(indices nm, uint n_u, uint n_d) const
797 {
798  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
799  CheckNf(nm, nf);
800 
801  gslpp::matrix<double> gammaDF1(2, 2, 0.);
802  double Qu = 2. / 3.;
803  double Qd = -1. / 3.;
804  double Qbar = n_u * Qu + n_d*Qd;
805  double z3 = gslpp_special_functions::zeta(3);
806 
807  switch (nm)
808  {
809  // QCD
810  //ref.: Gorbahn, Haisch, Misiak, Phys. Rev. Lett. 95, 102004, hep-ph/0504194
811  case 10:
812  gammaDF1(0, 0) = 32. / 3.;
813  gammaDF1(1, 0) = Qd * 32. / 3.;
814  gammaDF1(1, 1) = 28. / 3.;
815  break;
816  case 20:
817  gammaDF1(0, 0) = 1936. / 9. - nf * 224. / 27.;
818  gammaDF1(1, 0) = Qd * 368. / 3. - nf * Qd * 224. / 27.;
819  gammaDF1(1, 1) = 1456. / 9. - nf * 61. / 27.;
820  break;
821  case 30:
822  gammaDF1(0, 0) = 307448. / 81. - nf * 23776. / 81. - nf * nf * 352. / 81. + z3 * (-1856. / 27. - nf * 1280. / 9.);
823  gammaDF1(1, 0) = -Qbar * 1600. / 27. + Qd * 159872. / 81. - nf * Qd * 17108. / 81. - nf * nf * Qd * 352. / 81. + z3 * (Qbar * 640. / 9. -
824  Qd * 1856. / 27. - nf * Qd * 1280. / 9.);
825  gammaDF1(1, 1) = 268807. / 81. - nf * 4343. / 27. - nf * nf * 461. / 81. + z3 * (-28624. / 27. - nf * 1312. / 9.);
826  break;
827  // QED
828  // only available for nf = 5
829  // ref.: Bobeth, Gambino, Gorbahn, Haisch, JHEP 0404, 071, hep-ph/0312090
830  case 01:
831  gammaDF1(0, 0) = 16. / 9.;
832  gammaDF1(0, 1) = -8. / 3.;
833  gammaDF1(1, 1) = 8. / 9.;
834  break;
835  case 11:
836  gammaDF1(0, 0) = -256. / 27.;
837  gammaDF1(0, 1) = -52. / 9.;
838  gammaDF1(1, 0) = 128. / 81.;
839  gammaDF1(1, 1) = -40. / 27.;
840  break;
841  case 21:
842  case 02:
843  break;
844  default:
845  throw std::runtime_error("EvolDF1::GammaMM(): order not implemented");
846  }
847  return (gammaDF1);
848 }
849 
850 gslpp::matrix<double> EvolDF1::GammaLL(indices nm, uint n_u, uint n_d) const
851 {
852  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
853  if (nf != 5)
854  throw std::runtime_error("EvolDF1::GammaLL(): Wrong number of flavours in anoumalous dimensions");
855 
856  gslpp::matrix<double> gammaDF1(2, 2, 0.);
857 
858  switch (nm)
859  {
860  // QED
861  // ref.: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
862  case 01:
863  gammaDF1(0, 0) = 8.;
864  gammaDF1(0, 1) = -4.;
865  gammaDF1(1, 0) = -4.;
866  break;
867  case 11:
868  gammaDF1(0, 1) = 16.;
869  gammaDF1(1, 0) = 16.;
870  break;
871  case 10:
872  case 20:
873  case 30:
874  case 21:
875  case 02:
876  break;
877  default:
878  throw std::runtime_error("EvolDF1::GammaLL(): order not implemented");
879  }
880  return (gammaDF1);
881 }
882 
883 gslpp::matrix<double> EvolDF1::GammaQP(indices nm, uint n_u, uint n_d) const
884 {
885  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
886  if (nf != 5)
887  throw std::runtime_error("EvolDF1::GammaQP(): Wrong number of flavours in anoumalous dimensions");
888 
889  gslpp::matrix<double> gammaDF1(4, 4, 0.);
890 
891  switch (nm)
892  {
893  // QCD
894  //ref.: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
895  case 10:
896  gammaDF1(0, 1) = -8. / 9.;
897  gammaDF1(1, 1) = 16. / 27.;
898  gammaDF1(2, 1) = -128. / 9.;
899  gammaDF1(3, 1) = 184. / 27.;
900  break;
901  case 20:
902  gammaDF1(0, 0) = 832. / 243.;
903  gammaDF1(0, 1) = -4000. / 243.;
904  gammaDF1(0, 2) = -112. / 243.;
905  gammaDF1(0, 3) = -70. / 81.;
906  gammaDF1(1, 0) = 3376. / 729.;
907  gammaDF1(1, 1) = 6344. / 729.;
908  gammaDF1(1, 2) = -280. / 729.;
909  gammaDF1(1, 3) = 55. / 486.;
910  gammaDF1(2, 0) = 2272. / 243.;
911  gammaDF1(2, 1) = -72088. / 243.;
912  gammaDF1(2, 2) = -688. / 243.;
913  gammaDF1(2, 3) = -1240. / 81.;
914  gammaDF1(3, 0) = 45424. / 729.;
915  gammaDF1(3, 1) = 84236. / 729.;
916  gammaDF1(3, 2) = -3880. / 729.;
917  gammaDF1(3, 3) = 1220. / 243.;
918  break;
919  // QED
920  // ref.: Bobeth, Gambino, Gorbahn, Haisch, JHEP 0404, 071, hep-ph/0312090
921  case 01:
922  gammaDF1(0, 0) = 40. / 27.;
923  gammaDF1(0, 2) = -4. / 27.;
924  gammaDF1(1, 1) = 40. / 27.;
925  gammaDF1(1, 3) = -4. / 27.;
926  gammaDF1(2, 0) = 256. / 27.;
927  gammaDF1(2, 2) = -40. / 27.;
928  gammaDF1(3, 1) = 256. / 27.;
929  gammaDF1(3, 3) = -40. / 27.;
930  break;
931  case 11:
932  gammaDF1(0, 0) = -2240. / 81.;
933  gammaDF1(0, 1) = 39392. / 729.;
934  gammaDF1(0, 2) = 224. / 81.;
935  gammaDF1(0, 3) = -92. / 27.;
936  gammaDF1(1, 0) = 2176. / 243.;
937  gammaDF1(1, 1) = 84890. / 2187.;
938  gammaDF1(1, 2) = -184. / 243.;
939  gammaDF1(1, 3) = -224. / 81.;
940  gammaDF1(2, 0) = -23552. / 81.;
941  gammaDF1(2, 1) = 399776. / 729.;
942  gammaDF1(2, 2) = 2240. / 81.;
943  gammaDF1(2, 3) = -752. / 27.;
944  gammaDF1(3, 0) = 23296. / 243.;
945  gammaDF1(3, 1) = 933776. / 2187.;
946  gammaDF1(3, 2) = -1504. / 243.;
947  gammaDF1(3, 3) = -2030. / 81.;
948  break;
949  case 30:
950  case 21:
951  case 02:
952  break;
953  default:
954  throw std::runtime_error("EvolDF1::GammaQP(): order not implemented");
955  }
956  return (gammaDF1);
957 }
958 
959 gslpp::matrix<double> EvolDF1::GammaQM(indices nm, uint n_u, uint n_d) const
960 {
961  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
962  if (nf != 5)
963  throw std::runtime_error("EvolDF1::GammaQM(): Wrong number of flavours in anoumalous dimensions");
964 
965  gslpp::matrix<double> gammaDF1(4, 2, 0.);
966 
967  switch (nm)
968  {
969  // QCD
970  // ref.: Baranowski, Misiak, Phys. Lett. B 483, 410, hep-ph/9907427
971  case 10:
972  gammaDF1(0, 0) = 176. / 243.;
973  gammaDF1(0, 1) = -14. / 81.;
974  gammaDF1(1, 0) = -136. / 729.;
975  gammaDF1(1, 1) = -295. / 486.;
976  gammaDF1(2, 0) = 6272. / 243.;
977  gammaDF1(2, 1) = -764. / 81.;
978  gammaDF1(3, 0) = 39152. / 729.;
979  gammaDF1(3, 1) = -1892. / 243.;
980  break;
981  case 20:
982  case 30:
983  case 01:
984  case 11:
985  case 21:
986  case 02:
987  break;
988  default:
989  throw std::runtime_error("EvolDF1::GammaQM(): order not implemented");
990  }
991  return (gammaDF1);
992 }
993 
994 gslpp::matrix<double> EvolDF1::GammaQL(indices nm, uint n_u, uint n_d) const
995 {
996  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
997  if (nf != 5)
998  throw std::runtime_error("EvolDF1::GammaQL(): Wrong number of flavours in anoumalous dimensions");
999 
1000  gslpp::matrix<double> gammaDF1(4, 2, 0.);
1001 
1002  switch (nm)
1003  {
1004  // QED
1005  // ref.: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
1006  case 01:
1007  gammaDF1(0, 0) = -272. / 27.;
1008  gammaDF1(1, 0) = -32. / 81.;
1009  gammaDF1(2, 0) = -2768. / 27.;
1010  gammaDF1(3, 0) = -512. / 81.;
1011  break;
1012  case 11:
1013  gammaDF1(0, 0) = -24352. / 729.;
1014  gammaDF1(1, 0) = 54608. / 2187.;
1015  gammaDF1(2, 0) = -227008. / 729.;
1016  gammaDF1(3, 0) = 551648. / 2187.;
1017  break;
1018  case 10:
1019  case 20:
1020  case 30:
1021  case 21:
1022  case 02:
1023  break;
1024  default:
1025  throw std::runtime_error("EvolDF1::GammaQL(): order not implemented");
1026  }
1027  return (gammaDF1);
1028 }
1029 
1030 gslpp::matrix<double> EvolDF1::GammaQQ(indices nm, uint n_u, uint n_d) const
1031 {
1032  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
1033  if (nf != 5)
1034  throw std::runtime_error("EvolDF1::GammaQQ(): Wrong number of flavours in anoumalous dimensions");
1035 
1036  gslpp::matrix<double> gammaDF1(4, 4, 0.);
1037 
1038  switch (nm)
1039  {
1040  // QCD
1041  // ref.: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
1042  case 10:
1043  gammaDF1(0, 1) = -20.;
1044  gammaDF1(0, 3) = 2.;
1045  gammaDF1(1, 0) = -40. / 9.;
1046  gammaDF1(1, 1) = -52. / 3.;
1047  gammaDF1(1, 2) = 4. / 9.;
1048  gammaDF1(1, 3) = 5. / 6.;
1049  gammaDF1(2, 1) = -128.;
1050  gammaDF1(2, 3) = 20.;
1051  gammaDF1(3, 0) = -256. / 9.;
1052  gammaDF1(3, 1) = -160. / 3.;
1053  gammaDF1(3, 2) = 40. / 9.;
1054  gammaDF1(3, 3) = -2. / 3.;
1055  break;
1056  case 20:
1057  gammaDF1(0, 0) = -404. / 9.;
1058  gammaDF1(0, 1) = -3077. / 9.;
1059  gammaDF1(0, 2) = 32. / 9.;
1060  gammaDF1(0, 3) = 1031. / 36.;
1061  gammaDF1(1, 0) = -2698. / 81.;
1062  gammaDF1(1, 1) = -8035. / 27.;
1063  gammaDF1(1, 2) = -49. / 162.;
1064  gammaDF1(1, 3) = 4493. / 216.;
1065  gammaDF1(2, 0) = -19072. / 9.;
1066  gammaDF1(2, 1) = -14096. / 9.;
1067  gammaDF1(2, 2) = 1708. / 9.;
1068  gammaDF1(2, 3) = 1622. / 9.;
1069  gammaDF1(3, 0) = 32288. / 81.;
1070  gammaDF1(3, 1) = -15976. / 27.;
1071  gammaDF1(3, 2) = -6692. / 81.;
1072  gammaDF1(3, 3) = -2437. / 54.;
1073  break;
1074  // QED
1075  // ref.: Bobeth, Gambino, Gorbahn, Haisch, JHEP 0404, 071, hep-ph/0312090
1076  case 01:
1077  gammaDF1(0, 0) = 332. / 27.;
1078  gammaDF1(0, 2) = -2. / 9.;
1079  gammaDF1(1, 0) = 32. / 81.;
1080  gammaDF1(1, 1) = 20. / 9.;
1081  gammaDF1(1, 3) = -2. / 9.;
1082  gammaDF1(2, 0) = 3152. / 27.;
1083  gammaDF1(2, 2) = -20. / 9.;
1084  gammaDF1(3, 0) = 512. / 81.;
1085  gammaDF1(3, 1) = 128. / 9.;
1086  gammaDF1(3, 3) = -20. / 9.;
1087  break;
1088  case 11:
1089  gammaDF1(0, 0) = -5888. / 729.;
1090  gammaDF1(0, 1) = 13916. / 81.;
1091  gammaDF1(0, 2) = 112. / 27.;
1092  gammaDF1(0, 3) = -812. / 81.;
1093  gammaDF1(1, 0) = -2552. / 2187.;
1094  gammaDF1(1, 1) = 15638. / 243.;
1095  gammaDF1(1, 2) = -176. / 81.;
1096  gammaDF1(1, 3) = -2881. / 486.;
1097  gammaDF1(2, 0) = -90944. / 729.;
1098  gammaDF1(2, 1) = 90128. / 81.;
1099  gammaDF1(2, 2) = 1120. / 27.;
1100  gammaDF1(2, 3) = -1748. / 81.;
1101  gammaDF1(3, 0) = 1312. / 2187.;
1102  gammaDF1(3, 1) = 102488. / 243.;
1103  gammaDF1(3, 2) = -1592. / 81.;
1104  gammaDF1(3, 3) = -6008. / 243.;
1105  break;
1106  case 30:
1107  case 21:
1108  case 02:
1109  break;
1110  default:
1111  throw std::runtime_error("EvolDF1::GammaQQ(): order not implemented");
1112  }
1113  return (gammaDF1);
1114 }
1115 
1116 gslpp::matrix<double> EvolDF1::GammaBP(indices nm, uint n_u, uint n_d) const
1117 {
1118  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
1119  if (nf != 5)
1120  throw std::runtime_error("EvolDF1::GammaBP(): Wrong number of flavours in anoumalous dimensions");
1121 
1122  gslpp::matrix<double> gammaDF1(1, 4, 0.);
1123 
1124  switch (nm)
1125  {
1126  // QCD
1127  // ref.: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
1128  case 10:
1129  gammaDF1(0, 1) = 4. / 3.;
1130  break;
1131  case 20:
1132  gammaDF1(0, 0) = -1576. / 81.;
1133  gammaDF1(0, 1) = 446. / 27.;
1134  gammaDF1(0, 2) = 172. / 81.;
1135  gammaDF1(0, 3) = 40. / 27.;
1136  break;
1137  // QED
1138  // ref.: Bobeth, Gambino, Gorbahn, Haisch, JHEP 0404, 071, hep-ph/0312090
1139  case 01:
1140  break;
1141  case 11:
1142  gammaDF1(0, 1) = -232. / 81.;
1143  break;
1144  case 30:
1145  case 21:
1146  case 02:
1147  break;
1148  default:
1149  throw std::runtime_error("EvolDF1::GammaBP(): order not implemented");
1150  }
1151  return (gammaDF1);
1152 }
1153 
1154 gslpp::matrix<double> EvolDF1::GammaBL(indices nm, uint n_u, uint n_d) const
1155 {
1156  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
1157  if (nf != 5)
1158  throw std::runtime_error("EvolDF1::GammaBL(): Wrong number of flavours in anoumalous dimensions");
1159 
1160  gslpp::matrix<double> gammaDF1(1, 2, 0.);
1161 
1162  switch (nm)
1163  {
1164  // QED
1165  // ref.: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
1166  case 01:
1167  gammaDF1(0, 0) = 16. / 9.;
1168  break;
1169  case 11:
1170  gammaDF1(0, 0) = -8. / 9.;
1171  break;
1172  case 10:
1173  case 20:
1174  case 30:
1175  case 21:
1176  case 02:
1177  break;
1178  default:
1179  throw std::runtime_error("EvolDF1::GammaBL(): order not implemented");
1180  }
1181  return (gammaDF1);
1182 }
1183 
1184 gslpp::matrix<double> EvolDF1::GammaBQ(indices nm, uint n_u, uint n_d) const
1185 {
1186  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
1187  if (nf != 5)
1188  throw std::runtime_error("EvolDF1::GammaBQ(): Wrong number of flavours in anoumalous dimensions");
1189 
1190  gslpp::matrix<double> gammaDF1(1, 4, 0.);
1191 
1192  switch (nm)
1193  {
1194  // QED
1195  // ref.: Bobeth, Gambino, Gorbahn, Haisch, JHEP 0404, 071, hep-ph/0312090
1196  case 01:
1197  gammaDF1(0, 0) = -16. / 9.;
1198  break;
1199  case 11:
1200  gammaDF1(0, 1) = 580. / 27.;
1201  gammaDF1(0, 3) = -94. / 27.;
1202  break;
1203  case 10:
1204  case 20:
1205  case 30:
1206  case 21:
1207  case 02:
1208  break;
1209  default:
1210  throw std::runtime_error("EvolDF1::GammaBQ(): order not implemented");
1211  }
1212  return (gammaDF1);
1213 }
1214 
1215 gslpp::matrix<double> EvolDF1::GammaBB(indices nm, uint n_u, uint n_d) const
1216 {
1217  uint nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
1218  if (nf != 5)
1219  throw std::runtime_error("EvolDF1::GammaBB(): Wrong number of flavours in anoumalous dimensions");
1220 
1221  gslpp::matrix<double> gammaDF1(1, 1, 0.);
1222 
1223  switch (nm)
1224  {
1225  // QCD
1226  // ref.: Huber, Lunghi, Misiak, Wyler, Nucl. Phys. B 740, 105, hep-ph/0512066
1227  case 10:
1228  gammaDF1(0, 0) = 4.;
1229  break;
1230  case 20:
1231  gammaDF1(0, 0) = 325. / 9.;
1232  break;
1233  // QED
1234  // ref.: Bobeth, Gambino, Gorbahn, Haisch, JHEP 0404, 071, hep-ph/0312090
1235  case 01:
1236  gammaDF1(0, 0) = 4. / 3.;
1237  break;
1238  case 11:
1239  gammaDF1(0, 0) = -388. / 9.;
1240  break;
1241  case 30:
1242  case 21:
1243  case 02:
1244  break;
1245  default:
1246  throw std::runtime_error("EvolDF1::GammaBB(): order not implemented");
1247  }
1248  return (gammaDF1);
1249 }
1250 
1251 gslpp::matrix<double> EvolDF1::AnomalousDimension(indices nm, uint n_u, uint n_d) const
1252 {
1253  gslpp::matrix<double> gammaDF1(nops, nops, 0.);
1254 
1255  // assign blocks according to user request: "C", "CP", "CPM", "L", "CPML", "CPQB", "CPMQB", "CPMLQB"
1256 
1257  if (blocks.compare("C") == 0)
1258  gammaDF1.assign(0, 0, GammaCC(nm, n_u, n_d));
1259  else if (blocks.compare("CP") == 0)
1260  {
1261  gammaDF1.assign(0, 0, GammaCC(nm, n_u, n_d));
1262  gammaDF1.assign(0, 2, GammaCP(nm, n_u, n_d));
1263  gammaDF1.assign(2, 2, GammaPP(nm, n_u, n_d));
1264  } else if (blocks.compare("CPM") == 0)
1265  {
1266  gammaDF1.assign(0, 0, GammaCC(nm, n_u, n_d));
1267  gammaDF1.assign(0, 2, GammaCP(nm, n_u, n_d));
1268  gammaDF1.assign(2, 2, GammaPP(nm, n_u, n_d));
1269  gammaDF1.assign(0, 6, GammaCM(nm, n_u, n_d));
1270  gammaDF1.assign(2, 6, GammaPM(nm, n_u, n_d));
1271  gammaDF1.assign(6, 6, GammaMM(nm, n_u, n_d));
1272  } else if (blocks.compare("CPQ") == 0)
1273  {
1274  gammaDF1.assign(0, 0, GammaCC(nm, n_u, n_d));
1275  gammaDF1.assign(0, 2, GammaCP(nm, n_u, n_d));
1276  gammaDF1.assign(2, 2, GammaPP(nm, n_u, n_d));
1277  gammaDF1.assign(0, 6, GammaCQ(nm, n_u, n_d));
1278  gammaDF1.assign(2, 6, GammaPQ(nm, n_u, n_d));
1279  gammaDF1.assign(6, 2, GammaQP(nm, n_u, n_d));
1280  gammaDF1.assign(6, 6, GammaQQ(nm, n_u, n_d));
1281  } else if (blocks.compare("L") == 0)
1282  {
1283  gammaDF1.assign(0, 0, GammaLL(nm, n_u, n_d));
1284  } else if (blocks.compare("CPL") == 0)
1285  {
1286  gammaDF1.assign(0, 0, GammaCC(nm, n_u, n_d));
1287  gammaDF1.assign(0, 2, GammaCP(nm, n_u, n_d));
1288  gammaDF1.assign(2, 2, GammaPP(nm, n_u, n_d));
1289  gammaDF1.assign(0, 6, GammaCL(nm, n_u, n_d));
1290  gammaDF1.assign(2, 6, GammaPL(nm, n_u, n_d));
1291  gammaDF1.assign(6, 6, GammaLL(nm, n_u, n_d));
1292  } else if (blocks.compare("CPML") == 0)
1293  {
1294  gammaDF1.assign(0, 0, GammaCC(nm, n_u, n_d));
1295  gammaDF1.assign(0, 2, GammaCP(nm, n_u, n_d));
1296  gammaDF1.assign(2, 2, GammaPP(nm, n_u, n_d));
1297  gammaDF1.assign(0, 6, GammaCM(nm, n_u, n_d));
1298  gammaDF1.assign(2, 6, GammaPM(nm, n_u, n_d));
1299  gammaDF1.assign(6, 6, GammaMM(nm, n_u, n_d));
1300  gammaDF1.assign(0, 8, GammaCL(nm, n_u, n_d));
1301  gammaDF1.assign(2, 8, GammaPL(nm, n_u, n_d));
1302  gammaDF1.assign(8, 8, GammaLL(nm, n_u, n_d));
1303  } else if (blocks.compare("CPQB") == 0)
1304  {
1305  gammaDF1.assign(0, 0, GammaCC(nm, n_u, n_d));
1306  gammaDF1.assign(0, 2, GammaCP(nm, n_u, n_d));
1307  gammaDF1.assign(2, 2, GammaPP(nm, n_u, n_d));
1308  gammaDF1.assign(0, 6, GammaCQ(nm, n_u, n_d));
1309  gammaDF1.assign(2, 6, GammaPQ(nm, n_u, n_d));
1310  gammaDF1.assign(6, 2, GammaQP(nm, n_u, n_d));
1311  gammaDF1.assign(6, 6, GammaQQ(nm, n_u, n_d));
1312  gammaDF1.assign(10, 2, GammaBP(nm, n_u, n_d));
1313  gammaDF1.assign(10, 6, GammaBQ(nm, n_u, n_d));
1314  gammaDF1.assign(10, 10, GammaBB(nm, n_u, n_d));
1315  } else if (blocks.compare("CPMQB") == 0)
1316  {
1317  gammaDF1.assign(0, 0, GammaCC(nm, n_u, n_d));
1318  gammaDF1.assign(0, 2, GammaCP(nm, n_u, n_d));
1319  gammaDF1.assign(2, 2, GammaPP(nm, n_u, n_d));
1320  gammaDF1.assign(0, 6, GammaCM(nm, n_u, n_d));
1321  gammaDF1.assign(2, 6, GammaPM(nm, n_u, n_d));
1322  gammaDF1.assign(6, 6, GammaMM(nm, n_u, n_d));
1323  gammaDF1.assign(0, 8, GammaCQ(nm, n_u, n_d));
1324  gammaDF1.assign(2, 8, GammaPQ(nm, n_u, n_d));
1325  gammaDF1.assign(8, 2, GammaQP(nm, n_u, n_d));
1326  gammaDF1.assign(8, 6, GammaQM(nm, n_u, n_d));
1327  gammaDF1.assign(8, 8, GammaQQ(nm, n_u, n_d));
1328  gammaDF1.assign(12, 2, GammaBP(nm, n_u, n_d));
1329  gammaDF1.assign(12, 8, GammaBQ(nm, n_u, n_d));
1330  gammaDF1.assign(12, 12, GammaBB(nm, n_u, n_d)); // *** does BM exists?
1331  } else if (blocks.compare("CPMLQB") == 0)
1332  {
1333  gammaDF1.assign(0, 0, GammaCC(nm, n_u, n_d));
1334  gammaDF1.assign(0, 2, GammaCP(nm, n_u, n_d));
1335  gammaDF1.assign(2, 2, GammaPP(nm, n_u, n_d));
1336  gammaDF1.assign(0, 6, GammaCM(nm, n_u, n_d));
1337  gammaDF1.assign(2, 6, GammaPM(nm, n_u, n_d));
1338  gammaDF1.assign(6, 6, GammaMM(nm, n_u, n_d));
1339  gammaDF1.assign(0, 8, GammaCL(nm, n_u, n_d));
1340  gammaDF1.assign(2, 8, GammaPL(nm, n_u, n_d));
1341  gammaDF1.assign(8, 8, GammaLL(nm, n_u, n_d));
1342 
1343  gammaDF1.assign(0, 10, GammaCQ(nm, n_u, n_d));
1344  gammaDF1.assign(2, 10, GammaPQ(nm, n_u, n_d));
1345  gammaDF1.assign(10, 2, GammaQP(nm, n_u, n_d));
1346  gammaDF1.assign(10, 6, GammaQM(nm, n_u, n_d));
1347  gammaDF1.assign(10, 8, GammaQL(nm, n_u, n_d));
1348  gammaDF1.assign(10, 10, GammaQQ(nm, n_u, n_d));
1349  gammaDF1.assign(14, 2, GammaBP(nm, n_u, n_d));
1350  gammaDF1.assign(14, 8, GammaBL(nm, n_u, n_d));
1351  gammaDF1.assign(14, 10, GammaBQ(nm, n_u, n_d));
1352  gammaDF1.assign(14, 14, GammaBB(nm, n_u, n_d)); // *** does BM exists?
1353  } else
1354  throw std::runtime_error("EvolDF1::AnomalousDimension(): block not implemented");
1355 
1356  return (gammaDF1);
1357 }
1358 
1359 const Expanded<gslpp::matrix<double> >& EvolDF1::DF1Evol(double mu, double M, schemes scheme)
1360 {
1361  if (nfmin == 5 && nfmax == 5 && (model.Nf(mu) != 5. || model.Nf(M) != 5.))
1362  throw std::runtime_error("EvolDF1::Df1Evol(): only nf = 5 available.");
1363 
1364  switch (scheme)
1365  {
1366  case NDR:
1367  break;
1368  case LRI:
1369  case HV:
1370  default:
1371  std::stringstream out;
1372  out << scheme;
1373  throw std::runtime_error("EvolDF1::Df1Evol(): scheme " + out.str() + " not implemented ");
1374  }
1375 
1376  double alsM = model.getAlsM();
1377  double MAls = model.getMAls();
1378  if (alsM == alsM_cache && MAls == MAls_cache)
1379  {
1380  if (mu == this->mu && M == this->M && scheme == this->scheme)
1381  return (getEvol());
1382  }
1383  alsM_cache = alsM;
1384  MAls_cache = MAls;
1385 
1386  if (M < mu)
1387  {
1388  std::stringstream out;
1389  out << "M = " << M << " < mu = " << mu;
1390  throw std::runtime_error("EvolDF1::Df1Evol(): " + out.str() + ".");
1391  }
1392 
1393  setScales(mu, M); // also assign evol to identity
1394 
1395  double m_down = mu;
1396  double m_up = model.AboveTh(m_down);
1397  double nf = model.Nf(m_down);
1398 
1399  while (m_up < M)
1400  { // where are the nf thresholds? ???? <<<<<<<<<<
1401  DF1Ev(m_down, m_up, (int) nf, scheme);
1402  m_down = m_up;
1403  m_up = model.AboveTh(m_down);
1404  nf += 1.;
1405  }
1406  DF1Ev(m_down, M, (int) nf, scheme);
1407 
1408  return (getEvol());
1409 }
1410 
1411 //gslpp::matrix<double>& EvolDF1::DF1Evol(double mu, double M, orders_qed ord, schemes scheme)
1412 //{
1413 // if(ord > order_qed)
1414 // throw std::runtime_error("EvolDF1::Df1Evol(): order not present in this Hamiltonian.");
1415 // double MAls = model.getMAls();
1416 // if(model.Nf(mu) != 5. || model.Nf(M) != 5. || model.Nf(MAls) != 5.)
1417 // throw std::runtime_error("EvolDF1::Df1Evol(): only nf = 5 available.");
1418 //
1419 // switch (scheme) {
1420 // case NDR:
1421 // break;
1422 // case LRI:
1423 // case HV:
1424 // default:
1425 // std::stringstream out;
1426 // out << scheme;
1427 // throw std::runtime_error("EvolDF1::Df1Evol(): scheme " + out.str() + " not implemented ");
1428 // }
1429 //
1431 //
1432 // double alsM = model.getAlsM();
1433 // if(alsM == alsM_cache && MAls == MAls_cache) {
1434 // if (mu == this->mu && M == this->M && scheme == this->scheme)
1435 // return (*Evol(ord));
1436 // }
1437 // alsM_cache = alsM;
1438 // MAls_cache = MAls;
1439 //
1440 // if (M < mu) {
1441 // std::stringstream out;
1442 // out << "M = " << M << " < mu = " << mu;
1443 // throw std::runtime_error("EvolDF1::Df1Evol(): " + out.str() + ".");
1444 // }
1445 //
1446 // setScales(mu, M); // also assign evol to identity
1447 //
1448 // DF1Ev(mu, M, 5, scheme); // only 5 flavour *******************
1449 //
1450 // return (*Evol(ord));
1451 //}
1452 
1453 void EvolDF1::DF1Ev(double mu, double M, int nf, schemes scheme)
1454 {
1455  gslpp::matrix<double> mtmp(nops, 0.);
1456  std::vector<std::vector<gslpp::matrix<double> > > vtmp2;
1457  std::vector<gslpp::matrix<double> > vtmp;
1458  for (int j = 0; j <= order_qed; j++)
1459  vtmp.push_back(mtmp);
1460  for (int i = 0; i <= order_qcd; i++)
1461  vtmp2.push_back(vtmp);
1462  Expanded<gslpp::matrix<double> > res(vtmp2);
1463  // gslpp::matrix<double> res01(nops, 0.), res02(nops, 0.), res11(nops, 0.), res12(nops, 0.),
1464  // res21(nops, 0.), resLO(nops, 0.), resNLO(nops, 0.), resNNLO(nops, 0.);
1465 
1466  // uint a, b, i, j, p, q
1467  uint nnf = nf - nfmin;
1468  double b0, b0e, b5, alsM, eta, omega, lambda; //, term;
1469  std::map< std::vector<uint>, double >::iterator itr;
1470  // std::vector<uint> v;
1471 
1472 
1473  // alsM = model.Als(M) / 4. / M_PI;
1474  // double alsmu = model.Als(mu) / 4. / M_PI;
1475  b0 = model.Beta_s(00, nf);
1476  alsM = model.Als(M, FULLNNNLO, order_qed == QED0 ? false : true);
1477  eta = alsM / model.Als(mu, FULLNNNLO, order_qed == QED0 ? false : true);
1478  // eta = alsM / model.Als(mu);
1479  omega = 2. * b0 * alsM / 4. / M_PI;
1480 
1481  for (itr = vM0vi[nnf].begin(); itr != vM0vi[nnf].end(); ++itr)
1482  {
1483  const std::vector<uint> &v = itr->first;
1484  const uint &a = v[0];
1485  const uint &b = v[1];
1486  const uint &i = v[2];
1487 
1488  res.setMatrixElement(QCD0, QED0, a, b, res.getOrd(QCD0, QED0)(a, b) + itr->second * pow(eta, ai[nnf].at(i)));
1489  }
1490 
1491  for (itr = vM1vi[nnf].begin(); itr != vM1vi[nnf].end(); ++itr)
1492  {
1493  const std::vector<uint> &v = itr->first;
1494  const uint &a = v[0];
1495  const uint &b = v[1];
1496  const uint &i = v[2];
1497  const uint &j = v[3];
1498 
1499  res.setMatrixElement(QCD1, QED0, a, b, res.getOrd(QCD1, QED0)(a, b) + omega * itr->second * f_f(nnf, i, j, -1, eta));
1500  }
1501 
1502  for (itr = vM2vi[nnf].begin(); itr != vM2vi[nnf].end(); ++itr)
1503  {
1504  const std::vector<uint> &v = itr->first;
1505  const uint &a = v[0];
1506  const uint &b = v[1];
1507  const uint &i = v[2];
1508  const uint &j = v[3];
1509 
1510  res.setMatrixElement(QCD2, QED0, a, b, res.getOrd(QCD2, QED0)(a, b) + omega * omega * itr->second * f_f(nnf, i, j, -2, eta));
1511  }
1512 
1513  for (itr = vM11vi[nnf].begin(); itr != vM11vi[nnf].end(); ++itr)
1514  {
1515  const std::vector<uint> &v = itr->first;
1516  const uint &a = v[0];
1517  const uint &b = v[1];
1518  const uint &i = v[2];
1519  const uint &j = v[3];
1520  const uint &p = v[4];
1521 
1522  res.setMatrixElement(QCD2, QED0, a, b, res.getOrd(QCD2, QED0)(a, b) + omega * omega * itr->second * f_g(nnf, i, p, j, -1, -1, eta));
1523  }
1524 
1525  if (order_qed != QED0)
1526  {
1527  b0e = model.Beta_e(00, nf);
1528  b5 = model.Beta_e(01, nf) / 2. / b0 / b0e - model.Beta_s(10, nf) / 2. / b0 / b0;
1529  lambda = b0e * model.Ale(M, FULLNLO) / b0 / model.Als(M, FULLNNNLO, true); // WARNING: CHANGE ME!!!
1530 
1531  for (itr = vM3vi[nnf].begin(); itr != vM3vi[nnf].end(); ++itr)
1532  {
1533  const std::vector<uint> &v = itr->first;
1534  const uint &a = v[0];
1535  const uint &b = v[1];
1536  const uint &i = v[2];
1537  const uint &j = v[3];
1538  const double &term = itr->second;
1539 
1540  res.setMatrixElement(QCD0, QED1, a, b, res.getOrd(QCD0, QED1)(a, b) + lambda * term * f_f(nnf, i, j, 1, eta));
1541  res.setMatrixElement(QCD0, QED2, a, b, res.getOrd(QCD0, QED2)(a, b) + lambda * lambda * term * (f_f(nnf, i, j, 2, eta) - f_f(nnf, i, j, 1, eta)));
1542  res.setMatrixElement(QCD1, QED2, a, b, res.getOrd(QCD1, QED2)(a, b) + omega * lambda * lambda * b5 * term * f_r(nnf, i, j, 1, eta));
1543  }
1544 
1545  for (itr = vM4vi[nnf].begin(); itr != vM4vi[nnf].end(); ++itr)
1546  {
1547  const std::vector<uint> &v = itr->first;
1548  const uint &a = v[0];
1549  const uint &b = v[1];
1550  const uint &i = v[2];
1551  const uint &j = v[3];
1552  const double &term = itr->second;
1553 
1554  res.setMatrixElement(QCD1, QED1, a, b, res.getOrd(QCD1, QED1)(a, b) + omega * lambda * term * f_f(nnf, i, j, 0, eta));
1555  res.setMatrixElement(QCD1, QED2, a, b, res.getOrd(QCD1, QED2)(a, b) - omega * lambda * lambda * term * f_f(nnf, i, j, 0, eta));
1556  }
1557 
1558  for (itr = vM5vi[nnf].begin(); itr != vM5vi[nnf].end(); ++itr)
1559  {
1560  const std::vector<uint> &v = itr->first;
1561  const uint &a = v[0];
1562  const uint &b = v[1];
1563  const uint &i = v[2];
1564  const uint &j = v[3];
1565 
1566  res.setMatrixElement(QCD2, QED1, a, b, res.getOrd(QCD2, QED1)(a, b) + omega * omega * lambda * itr->second * f_f(nnf, i, j, -1, eta));
1567  }
1568 
1569  for (itr = vM6vi[nnf].begin(); itr != vM6vi[nnf].end(); ++itr)
1570  {
1571  const std::vector<uint> &v = itr->first;
1572  const uint &a = v[0];
1573  const uint &b = v[1];
1574  const uint &i = v[2];
1575  const uint &j = v[3];
1576 
1577  res.setMatrixElement(QCD1, QED2, a, b, res.getOrd(QCD1, QED2)(a, b) + omega * lambda * lambda * itr->second * f_f(nnf, i, j, 1, eta));
1578  }
1579 
1580  for (itr = vM33vi[nnf].begin(); itr != vM33vi[nnf].end(); ++itr)
1581  {
1582  const std::vector<uint> &v = itr->first;
1583  const uint &a = v[0];
1584  const uint &b = v[1];
1585  const uint &i = v[2];
1586  const uint &j = v[3];
1587  const uint &p = v[4];
1588 
1589  res.setMatrixElement(QCD0, QED2, a, b, res.getOrd(QCD0, QED2)(a, b) + lambda * lambda * itr->second * f_g(nnf, i, p, j, 1, 1, eta));
1590  }
1591 
1592  for (itr = vM13vi[nnf].begin(); itr != vM13vi[nnf].end(); ++itr)
1593  {
1594  const std::vector<uint> &v = itr->first;
1595  const uint &a = v[0];
1596  const uint &b = v[1];
1597  const uint &i = v[2];
1598  const uint &j = v[3];
1599  const uint &p = v[4];
1600  const double &term = itr->second;
1601 
1602  res.setMatrixElement(QCD1, QED1, a, b, res.getOrd(QCD1, QED1)(a, b) + omega * lambda * term * f_g(nnf, i, p, j, -1, 1, eta));
1603  res.setMatrixElement(QCD1, QED2, a, b, res.getOrd(QCD1, QED2)(a, b) + omega * lambda * lambda * term * (f_g(nnf, i, p, j, -1, 2, eta) - f_g(nnf, i, p, j, -1, 1, eta)));
1604  }
1605 
1606  for (itr = vM31vi[nnf].begin(); itr != vM31vi[nnf].end(); ++itr)
1607  {
1608  const std::vector<uint> &v = itr->first;
1609  const uint &a = v[0];
1610  const uint &b = v[1];
1611  const uint &i = v[2];
1612  const uint &j = v[3];
1613  const uint &p = v[4];
1614  const double &term = itr->second;
1615 
1616  res.setMatrixElement(QCD1, QED1, a, b, res.getOrd(QCD1, QED1)(a, b) + omega * lambda * term * f_g(nnf, i, p, j, 1, -1, eta));
1617  res.setMatrixElement(QCD1, QED2, a, b, res.getOrd(QCD1, QED2)(a, b) + omega * lambda * lambda * term * (f_g(nnf, i, p, j, 2, -1, eta) - f_g(nnf, i, p, j, 1, -1, eta)));
1618  }
1619 
1620  for (itr = vM34vi[nnf].begin(); itr != vM34vi[nnf].end(); ++itr)
1621  {
1622  const std::vector<uint> &v = itr->first;
1623  const uint &a = v[0];
1624  const uint &b = v[1];
1625  const uint &i = v[2];
1626  const uint &j = v[3];
1627  const uint &p = v[4];
1628 
1629  res.setMatrixElement(QCD1, QED2, a, b, res.getOrd(QCD1, QED2)(a, b) + omega * lambda * lambda * itr->second * f_g(nnf, i, p, j, 1, 0, eta));
1630  }
1631 
1632  for (itr = vM43vi[nnf].begin(); itr != vM43vi[nnf].end(); ++itr)
1633  {
1634  const std::vector<uint> &v = itr->first;
1635  const uint &a = v[0];
1636  const uint &b = v[1];
1637  const uint &i = v[2];
1638  const uint &j = v[3];
1639  const uint &p = v[4];
1640 
1641  res.setMatrixElement(QCD1, QED2, a, b, res.getOrd(QCD1, QED2)(a, b) + omega * lambda * lambda * itr->second * f_g(nnf, i, p, j, 0, 1, eta));
1642  }
1643 
1644  for (itr = vM23vi[nnf].begin(); itr != vM23vi[nnf].end(); ++itr)
1645  {
1646  const std::vector<uint> &v = itr->first;
1647  const uint &a = v[0];
1648  const uint &b = v[1];
1649  const uint &i = v[2];
1650  const uint &j = v[3];
1651  const uint &p = v[4];
1652 
1653  res.setMatrixElement(QCD2, QED1, a, b, res.getOrd(QCD2, QED1)(a, b) + omega * omega * lambda * itr->second * f_g(nnf, i, p, j, -2, 1, eta));
1654  }
1655 
1656  for (itr = vM32vi[nnf].begin(); itr != vM32vi[nnf].end(); ++itr)
1657  {
1658  const std::vector<uint> &v = itr->first;
1659  const uint &a = v[0];
1660  const uint &b = v[1];
1661  const uint &i = v[2];
1662  const uint &j = v[3];
1663  const uint &p = v[4];
1664 
1665  res.setMatrixElement(QCD2, QED1, a, b, res.getOrd(QCD2, QED1)(a, b) + omega * omega * lambda * itr->second * f_g(nnf, i, p, j, 1, -2, eta));
1666  }
1667 
1668  for (itr = vM14vi[nnf].begin(); itr != vM14vi[nnf].end(); ++itr)
1669  {
1670  const std::vector<uint> &v = itr->first;
1671  const uint &a = v[0];
1672  const uint &b = v[1];
1673  const uint &i = v[2];
1674  const uint &j = v[3];
1675  const uint &p = v[4];
1676 
1677  res.setMatrixElement(QCD2, QED1, a, b, res.getOrd(QCD2, QED1)(a, b) + omega * omega * lambda * itr->second * f_g(nnf, i, p, j, -1, 0, eta));
1678  }
1679 
1680  for (itr = vM41vi[nnf].begin(); itr != vM41vi[nnf].end(); ++itr)
1681  {
1682  const std::vector<uint> &v = itr->first;
1683  const uint &a = v[0];
1684  const uint &b = v[1];
1685  const uint &i = v[2];
1686  const uint &j = v[3];
1687  const uint &p = v[4];
1688 
1689  res.setMatrixElement(QCD2, QED1, a, b, res.getOrd(QCD2, QED1)(a, b) + omega * omega * lambda * itr->second * f_g(nnf, i, p, j, 0, -1, eta));
1690  }
1691 
1692  for (itr = vM113vi[nnf].begin(); itr != vM113vi[nnf].end(); ++itr)
1693  {
1694  const std::vector<uint> &v = itr->first;
1695  const uint &a = v[0];
1696  const uint &b = v[1];
1697  const uint &i = v[2];
1698  const uint &j = v[3];
1699  const uint &p = v[4];
1700  const uint &q = v[5];
1701 
1702  res.setMatrixElement(QCD2, QED1, a, b, res.getOrd(QCD2, QED1)(a, b) + omega * omega * lambda * itr->second * f_h(nnf, i, p, q, j, -1, -1, 1, eta));
1703  }
1704 
1705  for (itr = vM131vi[nnf].begin(); itr != vM131vi[nnf].end(); ++itr)
1706  {
1707  const std::vector<uint> &v = itr->first;
1708  const uint &a = v[0];
1709  const uint &b = v[1];
1710  const uint &i = v[2];
1711  const uint &j = v[3];
1712  const uint &p = v[4];
1713  const uint &q = v[5];
1714 
1715  res.setMatrixElement(QCD2, QED1, a, b, res.getOrd(QCD2, QED1)(a, b) + omega * omega * lambda * itr->second * f_h(nnf, i, p, q, j, -1, 1, -1, eta));
1716  }
1717 
1718  for (itr = vM311vi[nnf].begin(); itr != vM311vi[nnf].end(); ++itr)
1719  {
1720  const std::vector<uint> &v = itr->first;
1721  const uint &a = v[0];
1722  const uint &b = v[1];
1723  const uint &i = v[2];
1724  const uint &j = v[3];
1725  const uint &p = v[4];
1726  const uint &q = v[5];
1727 
1728  res.setMatrixElement(QCD2, QED1, a, b, res.getOrd(QCD2, QED1)(a, b) + omega * omega * lambda * itr->second * f_h(nnf, i, p, q, j, 1, -1, -1, eta));
1729  }
1730 
1731  for (itr = vM133vi[nnf].begin(); itr != vM133vi[nnf].end(); ++itr)
1732  {
1733  const std::vector<uint> &v = itr->first;
1734  const uint &a = v[0];
1735  const uint &b = v[1];
1736  const uint &i = v[2];
1737  const uint &j = v[3];
1738  const uint &p = v[4];
1739  const uint &q = v[5];
1740 
1741  res.setMatrixElement(QCD1, QED2, a, b, res.getOrd(QCD1, QED2)(a, b) + omega * lambda * lambda * itr->second * f_h(nnf, i, p, q, j, -1, 1, 1, eta));
1742  }
1743 
1744  for (itr = vM313vi[nnf].begin(); itr != vM313vi[nnf].end(); ++itr)
1745  {
1746  const std::vector<uint> &v = itr->first;
1747  const uint &a = v[0];
1748  const uint &b = v[1];
1749  const uint &i = v[2];
1750  const uint &j = v[3];
1751  const uint &p = v[4];
1752  const uint &q = v[5];
1753 
1754  res.setMatrixElement(QCD1, QED2, a, b, res.getOrd(QCD1, QED2)(a, b) + omega * lambda * lambda * itr->second * f_h(nnf, i, p, q, j, 1, -1, 1, eta));
1755  }
1756 
1757  for (itr = vM331vi[nnf].begin(); itr != vM331vi[nnf].end(); ++itr)
1758  {
1759  const std::vector<uint> &v = itr->first;
1760  const uint &a = v[0];
1761  const uint &b = v[1];
1762  const uint &i = v[2];
1763  const uint &j = v[3];
1764  const uint &p = v[4];
1765  const uint &q = v[5];
1766 
1767  res.setMatrixElement(QCD1, QED2, a, b, res.getOrd(QCD1, QED2)(a, b) + omega * lambda * lambda * itr->second * f_h(nnf, i, p, q, j, 1, 1, -1, eta));
1768  }
1769  /*
1770  for (a = 0; a < nops; a++)
1771  for (b = 0; b < nops; b++)
1772  {
1773  for (i = 0; i < nops; i++)
1774  for (j = 0; j < nops; j++)
1775  {
1776  res01(a, b) += (lambda * vM3vi.at(idx(nf, a, b, i, j)) * f_f(nf, i, j, 1, eta)).real();
1777  res02(a, b) += (lambda * lambda * vM3vi.at(idx(nf, a, b, i, j)) * (f_f(nf, i, j, 2, eta) - f_f(nf, i, j, 1, eta))).real();
1778  res11(a, b) += (omega * lambda * vM4vi.at(idx(nf, a, b, i, j)) * f_f(nf, i, j, 0, eta)).real();
1779  res12(a, b) += (omega * lambda * lambda * (-vM4vi.at(idx(nf, a, b, i, j)) * f_f(nf, i, j, 0, eta) + vM6vi.at(idx(nf, a, b, i, j)) * f_f(nf, i, j, 1, eta)) +
1780  b5 * vM3vi.at(idx(nf, a, b, i, j)) * f_r(nf, i, j, 1, eta)).real();
1781  res21(a, b) += (omega * omega * lambda * vM5vi.at(idx(nf, a, b, i, j)) * f_f(nf, i, j, -1, eta)).real();
1782  for (p = 0; p < nops; p++)
1783  {
1784  res02(a, b) += (lambda * lambda * vM33vi.at(idx(nf, a, b, i, j, p)) * f_g(nf, i, p, j, 1, 1, eta)).real();
1785  res11(a, b) += (omega * lambda * (vM13vi.at(idx(nf, a, b, i, j, p)) * f_g(nf, i, p, j, -1, 1, eta) + vM31vi.at(idx(nf, a, b, i, j, p)) * f_g(nf, i, p, j, 1, -1, eta))).real();
1786  res12(a, b) += (omega * lambda * lambda * (vM13vi.at(idx(nf, a, b, i, j, p)) * (f_g(nf, i, p, j, -1, 2, eta) - f_g(nf, i, p, j, -1, 1, eta)) +
1787  vM31vi.at(idx(nf, a, b, i, j, p)) * (f_g(nf, i, p, j, 2, -1, eta) - f_g(nf, i, p, j, 1, -1, eta)) + vM34vi.at(idx(nf, a, b, i, j, p)) * f_g(nf, i, p, j, 1, 0, eta) +
1788  vM43vi.at(idx(nf, a, b, i, j, p)) * f_g(nf, i, p, j, 0, 1, eta))).real();
1789  res21(a, b) += (omega * omega * lambda * (vM14vi.at(idx(nf, a, b, i, j, p)) * f_g(nf, i, p, j, -1, 0, eta) + vM41vi.at(idx(nf, a, b, i, j, p)) * f_g(nf, i, p, j, 0, -1, eta) +
1790  vM23vi.at(idx(nf, a, b, i, j, p)) * f_g(nf, i, p, j, -2, 1, eta) + vM32vi.at(idx(nf, a, b, i, j, p)) * f_g(nf, i, p, j, 1, -2, eta))).real();
1791 
1792  for (q = 0; q < nops; q++)
1793  {
1794  res12(a, b) += (omega * lambda * lambda * (vM133vi.at(idx(nf, a, b, i, j, p, q)) * f_h(nf, i, p, q, j, -1, 1, 1, eta) + vM313vi.at(idx(nf, a, b, i, j, p, q)) *
1795  f_h(nf, i, p, q, j, 1, -1, 1, eta) + vM331vi.at(idx(nf, a, b, i, j, p, q)) * f_h(nf, i, p, q, j, 1, 1, -1, eta))).real();
1796  res21(a, b) += (omega * omega * lambda * (vM113vi.at(idx(nf, a, b, i, j, p, q)) * f_h(nf, i, p, q, j, -1, -1, 1, eta) + vM131vi.at(idx(nf, a, b, i, j, p, q)) *
1797  f_h(nf, i, p, q, j, -1, 1, -1, eta) + vM311vi.at(idx(nf, a, b, i, j, p, q)) * f_h(nf, i, p, q, j, 1, -1, -1, eta))).real();
1798  }
1799  }
1800  }
1801  if (fabs(res01(a, b)) < EPS) res01(a, b) = 0.;
1802  if (fabs(res02(a, b)) < EPS) res02(a, b) = 0.;
1803  if (fabs(res11(a, b)) < EPS) res11(a, b) = 0.;
1804  if (fabs(res12(a, b)) < EPS) res12(a, b) = 0.;
1805  if (fabs(res21(a, b)) < EPS) res21(a, b) = 0.;
1806  }
1807  */
1808  for (uint a = 0; a < nops; a++)
1809  for (uint b = 0; b < nops; b++)
1810  {
1811  if (fabs(res.getOrd(QCD0, QED0)(a, b)) < EPS) res.setMatrixElement(QCD0, QED0, a, b, 0.);
1812  if (fabs(res.getOrd(QCD1, QED0)(a, b)) < EPS) res.setMatrixElement(QCD1, QED0, a, b, 0.);
1813  if (fabs(res.getOrd(QCD2, QED0)(a, b)) < EPS) res.setMatrixElement(QCD2, QED0, a, b, 0.);
1814  if (fabs(res.getOrd(QCD0, QED1)(a, b)) < EPS) res.setMatrixElement(QCD0, QED1, a, b, 0.);
1815  if (fabs(res.getOrd(QCD0, QED2)(a, b)) < EPS) res.setMatrixElement(QCD0, QED2, a, b, 0.);
1816  if (fabs(res.getOrd(QCD1, QED1)(a, b)) < EPS) res.setMatrixElement(QCD1, QED1, a, b, 0.);
1817  if (fabs(res.getOrd(QCD1, QED2)(a, b)) < EPS) res.setMatrixElement(QCD1, QED2, a, b, 0.);
1818  if (fabs(res.getOrd(QCD2, QED1)(a, b)) < EPS) res.setMatrixElement(QCD2, QED1, a, b, 0.);
1819  }
1820  } else
1821  for (uint a = 0; a < nops; a++)
1822  for (uint b = 0; b < nops; b++)
1823  {
1824  if (fabs(res.getOrd(QCD0, QED0)(a, b)) < EPS) res.setMatrixElement(QCD0, QED0, a, b, 0.);
1825  if (fabs(res.getOrd(QCD1, QED0)(a, b)) < EPS) res.setMatrixElement(QCD1, QED0, a, b, 0.);
1826  if (fabs(res.getOrd(QCD2, QED0)(a, b)) < EPS) res.setMatrixElement(QCD2, QED0, a, b, 0.);
1827  }
1828 
1829  wilson = res * wilson;
1830 }
EvolDF1::vM13vi
std::map< std::vector< uint >, double > vM13vi[NF]
Definition: EvolDF1.h:329
gslpp_special_functions::zeta
double zeta(int i)
Definition: gslpp_special_functions.cpp:20
EvolDF1::GammaBP
gslpp::matrix< double > GammaBP(indices nm, uint n_u, uint n_d) const
BP block of the QCD+QED anomalous dimension.
Definition: EvolDF1.cpp:1116
EvolDF1::vM43vi
std::map< std::vector< uint >, double > vM43vi[NF]
Definition: EvolDF1.h:329
EvolDF1::f_f_c
int f_f_c[4][F_iCacheSize]
Definition: EvolDF1.h:347
Expanded::getOrd
const T & getOrd(int j) const
Get an element of a single expansion.
Definition: Expanded.h:694
RGEvolutorNew::getEvol
const Expanded< gslpp::matrix< double > > & getEvol() const
Definition: RGEvolutorNew.cpp:37
EvolDF1::vM1vi
std::map< std::vector< uint >, double > vM1vi[NF]
Definition: EvolDF1.h:329
EvolDF1::GammaCC
gslpp::matrix< double > GammaCC(indices nm, uint n_u, uint n_d) const
CC block of the QCD+QED anomalous dimension.
Definition: EvolDF1.cpp:279
Expanded::setMatrixElement
std::enable_if< isMat(T, Q), void >::type setMatrixElement(int j, int i, int h, int k, Q x)
Set an element of a matrix in a double Expanded<matrix<*> >
Definition: Expanded.h:741
EvolDF1::GammaQP
gslpp::matrix< double > GammaQP(indices nm, uint n_u, uint n_d) const
QP block of the QCD+QED anomalous dimension.
Definition: EvolDF1.cpp:883
gslpp::matrix< double >::assign
void assign(const size_t &i, const size_t &j, const double &a)
Definition: gslpp_matrix_double.cpp:108
EvolDF1::vM34vi
std::map< std::vector< uint >, double > vM34vi[NF]
Definition: EvolDF1.h:329
gslpp::matrix< double >
A class for constructing and defining operations on real matrices.
Definition: gslpp_matrix_double.h:48
EvolDF1::vM33vi
std::map< std::vector< uint >, double > vM33vi[NF]
Definition: EvolDF1.h:329
EvolDF1::GammaQM
gslpp::matrix< double > GammaQM(indices nm, uint n_u, uint n_d) const
QM block of the QCD anomalous dimension.
Definition: EvolDF1.cpp:959
QED0
Definition: OrderScheme.h:83
EvolDF1::vM31vi
std::map< std::vector< uint >, double > vM31vi[NF]
Definition: EvolDF1.h:329
EvolDF1::f_r
double f_r(uint nf, uint i, uint j, int k, double eta)
auxiliary function r - eq. (51) of Huber, Lunghi, Misiak, Wyler, hep-ph/0512066
Definition: EvolDF1.cpp:233
EvolDF1::vM2vi
std::map< std::vector< uint >, double > vM2vi[NF]
Definition: EvolDF1.h:329
EvolDF1::evec
gslpp::matrix< double > evec
Definition: EvolDF1.h:339
EvolDF1::vM41vi
std::map< std::vector< uint >, double > vM41vi[NF]
Definition: EvolDF1.h:329
EvolDF1::GammaMM
gslpp::matrix< double > GammaMM(indices nm, uint n_u, uint n_d) const
MM block of the QCD+QED anomalous dimension.
Definition: EvolDF1.cpp:796
QCD::Nf
double Nf(const double mu) const
The number of active flavour at scale .
Definition: QCD.cpp:438
EvolDF1::DF1Evol
const Expanded< gslpp::matrix< double > > & DF1Evol(double mu, double M, schemes scheme=NDR)
a method returning the evolutor related to the high scale and the low scale
Definition: EvolDF1.cpp:1359
EvolDF1::GammaCQ
gslpp::matrix< double > GammaCQ(indices nm, uint n_u, uint n_d) const
CQ block of the QED anomalous dimension.
Definition: EvolDF1.cpp:484
FULLNNNLO
Definition: OrderScheme.h:39
qcd_orders
qcd_orders
Definition: OrderScheme.h:65
RGEvolutorNew::M
double M
Definition: RGEvolutorNew.h:107
StandardModel::Beta_e
double Beta_e(int nm, unsigned int nf) const
QED beta function coefficients - eq. (36) hep-ph/0512066.
Definition: StandardModel.cpp:582
EvolDF1::vM133vi
std::map< std::vector< uint >, double > vM133vi[NF]
Definition: EvolDF1.h:329
NDR
Definition: OrderScheme.h:21
QCD::getAlsM
double getAlsM() const
A get method to access the value of .
Definition: QCD.h:543
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
EvolDF1::AnomalousDimension
gslpp::matrix< double > AnomalousDimension(indices nm, uint n_u, uint n_d) const
a method returning the anomalous dimension matrix given in the Misiak basis
Definition: EvolDF1.cpp:1251
EvolDF1::blocks
std::string blocks
Definition: EvolDF1.h:337
EvolDF1::vM11vi
std::map< std::vector< uint >, double > vM11vi[NF]
Definition: EvolDF1.h:329
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:477
EvolDF1::eval
gslpp::vector< double > eval
Definition: EvolDF1.h:340
EvolDF1::f_g
double f_g(uint nf, uint i, uint p, uint j, int k, int l, double eta)
auxiliary function g - eq. (52) of Huber, Lunghi, Misiak, Wyler, hep-ph/0512066
Definition: EvolDF1.cpp:243
EvolDF1::GammaPQ
gslpp::matrix< double > GammaPQ(indices nm, uint n_u, uint n_d) const
PQ block of the QED anomalous dimension.
Definition: EvolDF1.cpp:742
WilsonTemplateNew< gslpp::matrix< double > >::order_qed
qed_orders order_qed
Definition: WilsonTemplateNew.h:123
EvolDF1::f_h
double f_h(uint nf, uint i, uint p, uint q, uint j, int k, int l, int m, double eta)
auxiliary function h - eq. (53) of Huber, Lunghi, Misiak, Wyler, hep-ph/0512066
Definition: EvolDF1.cpp:253
WilsonTemplateNew< gslpp::matrix< double > >::wilson
Expanded< gslpp::matrix< double > > wilson
Definition: WilsonTemplateNew.h:116
QCD1
Definition: OrderScheme.h:68
WilsonTemplateNew< gslpp::matrix< double > >::order_qcd
qcd_orders order_qcd
Definition: WilsonTemplateNew.h:122
EvolDF1::vM313vi
std::map< std::vector< uint >, double > vM313vi[NF]
Definition: EvolDF1.h:329
EvolDF1::evalc
gslpp::vector< gslpp::complex > evalc
Definition: EvolDF1.h:342
EvolDF1::vM4vi
std::map< std::vector< uint >, double > vM4vi[NF]
Definition: EvolDF1.h:329
QCD2
Definition: OrderScheme.h:69
RGEvolutorNew
Definition: RGEvolutorNew.h:24
schemes
schemes
An enum type for regularization schemes.
Definition: OrderScheme.h:19
EvolDF1::vM32vi
std::map< std::vector< uint >, double > vM32vi[NF]
Definition: EvolDF1.h:329
EvolDF1::GammaPP
gslpp::matrix< double > GammaPP(indices nm, uint n_u, uint n_d) const
PP block of the QCD+QED anomalous dimension.
Definition: EvolDF1.cpp:521
EvolDF1::vM6vi
std::map< std::vector< uint >, double > vM6vi[NF]
Definition: EvolDF1.h:329
EvolDF1::GammaQQ
gslpp::matrix< double > GammaQQ(indices nm, uint n_u, uint n_d) const
QQ block of the QCD+QED anomalous dimension.
Definition: EvolDF1.cpp:1030
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
EvolDF1::evecc
gslpp::matrix< gslpp::complex > evecc
Definition: EvolDF1.h:341
EvolDF1::MAls_cache
double MAls_cache
Definition: EvolDF1.h:343
EvolDF1::model
const StandardModel & model
Definition: EvolDF1.h:333
qed_orders
qed_orders
Definition: OrderScheme.h:81
EvolDF1::EvolDF1
EvolDF1(std::string reqblocks, schemes scheme, const StandardModel &model_i, qcd_orders order_qcd, qed_orders order_qed)
EvolDF1 constructor.
Definition: EvolDF1.cpp:29
Expanded
A template class for Taylor double expansion of several objects.
Definition: Expanded.h:55
EvolDF1::vM131vi
std::map< std::vector< uint >, double > vM131vi[NF]
Definition: EvolDF1.h:329
indices
unsigned int indices
Definition: EvolDF1.h:21
StandardModel::Ale
double Ale(double mu, orders order, bool Nf_thr=true) const
The running electromagnetic coupling in the scheme.
Definition: StandardModel.cpp:732
uint
unsigned int uint
Definition: EvolDF1.h:20
EvolDF1::nops
uint nops
Definition: EvolDF1.h:336
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:602
QCD::CacheShift
void CacheShift(double cache[][5], int n) const
A member used to manage the caching for this class.
EvolDF1::vM331vi
std::map< std::vector< uint >, double > vM331vi[NF]
Definition: EvolDF1.h:329
QCD::getMAls
double getMAls() const
A get method to access the mass scale at which the strong coupling constant measurement is provided.
Definition: QCD.h:552
QED2
Definition: OrderScheme.h:85
QCD.h
EvolDF1::vM0vi
std::map< std::vector< uint >, double > vM0vi[NF]
Definition: EvolDF1.h:329
WilsonTemplateNew< gslpp::matrix< double > >::scheme
schemes scheme
Definition: WilsonTemplateNew.h:121
StandardModel::Beta_s
double Beta_s(int nm, unsigned int nf) const
QCD beta function coefficients including QED corrections - eq. (36) hep-ph/0512066.
Definition: StandardModel.cpp:554
QCD0
Definition: OrderScheme.h:67
EvolDF1::GammaBB
gslpp::matrix< double > GammaBB(indices nm, uint n_u, uint n_d) const
BB block of the QCD+QED anomalous dimension.
Definition: EvolDF1.cpp:1215
EvolDF1::vM311vi
std::map< std::vector< uint >, double > vM311vi[NF]
Definition: EvolDF1.h:329
LRI
Definition: OrderScheme.h:23
WilsonTemplateNew< gslpp::matrix< double > >::mu
double mu
Definition: WilsonTemplateNew.h:120
EvolDF1::ai
std::map< uint, double > ai[NF]
Definition: EvolDF1.h:328
EvolDF1::vM113vi
std::map< std::vector< uint >, double > vM113vi[NF]
Definition: EvolDF1.h:329
EvolDF1::CheckNf
void CheckNf(indices nm, uint nf) const
a method returning the anomalous dimension in the Chetyrkin, Misiak and Munz operator basis
Definition: EvolDF1.cpp:269
EvolDF1.h
blocks_nops
std::map< std::string, uint > blocks_nops
Definition: EvolDF1.cpp:16
EvolDF1::nfmin
uint nfmin
Definition: EvolDF1.h:336
EvolDF1::vM5vi
std::map< std::vector< uint >, double > vM5vi[NF]
Definition: EvolDF1.h:329
HV
Definition: OrderScheme.h:22
EvolDF1::vM14vi
std::map< std::vector< uint >, double > vM14vi[NF]
Definition: EvolDF1.h:329
QCD::AboveTh
double AboveTh(const double mu) const
The active flavour threshold above the scale as defined in QCD::Thresholds().
Definition: QCD.cpp:420
EvolDF1::f_f_d
double f_f_d[2][F_iCacheSize]
Definition: EvolDF1.h:348
EvolDF1::GammaLL
gslpp::matrix< double > GammaLL(indices nm, uint n_u, uint n_d) const
LL block of the QED anomalous dimension.
Definition: EvolDF1.cpp:850
EvolDF1::vM23vi
std::map< std::vector< uint >, double > vM23vi[NF]
Definition: EvolDF1.h:329
EvolDF1::~EvolDF1
virtual ~EvolDF1()
EvolDF1 destructor.
Definition: EvolDF1.cpp:192
EvolDF1::alsM_cache
double alsM_cache
Definition: EvolDF1.h:343
EvolDF1::GammaPL
gslpp::matrix< double > GammaPL(indices nm, uint n_u, uint n_d) const
PL block of the QED anomalous dimension.
Definition: EvolDF1.cpp:689
EvolDF1::evec_i
gslpp::matrix< double > evec_i
Definition: EvolDF1.h:339
EvolDF1::GammaCP
gslpp::matrix< double > GammaCP(indices nm, uint n_u, uint n_d) const
CP block of the QCD+QED anomalous dimension.
Definition: EvolDF1.cpp:330
EvolDF1::vM3vi
std::map< std::vector< uint >, double > vM3vi[NF]
Definition: EvolDF1.h:329
RGEvolutorNew::setScales
void setScales(double mu, double M)
Sets the upper and lower scale for the running of the Wilson Coefficients.
Definition: RGEvolutorNew.cpp:47
EvolDF1::GammaCL
gslpp::matrix< double > GammaCL(indices nm, uint n_u, uint n_d) const
CL block of the QED anomalous dimension.
Definition: EvolDF1.cpp:442
EvolDF1::GammaBL
gslpp::matrix< double > GammaBL(indices nm, uint n_u, uint n_d) const
BL block of the QED anomalous dimension.
Definition: EvolDF1.cpp:1154
gslpp::matrix< double >::eigensystem
void eigensystem(matrix< complex > &U, vector< complex > &S)
Definition: gslpp_matrix_double.cpp:280
EvolDF1::f_f
double f_f(uint nf, uint i, uint j, int k, double eta)
auxiliary function f - eq. (50) of Huber, Lunghi, Misiak, Wyler, hep-ph/0512066
Definition: EvolDF1.cpp:206
EvolDF1::GammaCM
gslpp::matrix< double > GammaCM(indices nm, uint n_u, uint n_d) const
CM block of the QCD+QED anomalous dimension.
Definition: EvolDF1.cpp:384
FULLNLO
Definition: OrderScheme.h:37
QED1
Definition: OrderScheme.h:84
EvolDF1::nfmax
uint nfmax
Definition: EvolDF1.h:336
EvolDF1::GammaBQ
gslpp::matrix< double > GammaBQ(indices nm, uint n_u, uint n_d) const
BQ block of the QED anomalous dimension.
Definition: EvolDF1.cpp:1184
EvolDF1::GammaQL
gslpp::matrix< double > GammaQL(indices nm, uint n_u, uint n_d) const
QL block of the QED anomalous dimension.
Definition: EvolDF1.cpp:994
EvolDF1::DF1Ev
void DF1Ev(double mu, double M, int nf, schemes scheme)
a void type method storing properly the magic numbers for the implementation of the evolutor
Definition: EvolDF1.cpp:1453
EvolDF1::GammaPM
gslpp::matrix< double > GammaPM(indices nm, uint n_u, uint n_d) const
PM block of the QCD+QED anomalous dimension.
Definition: EvolDF1.cpp:611