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
16std::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
29EvolDF1::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),
31evec(blocks_nops.at(reqblocks), 0.), evec_i(blocks_nops.at(reqblocks), 0.), js(blocks_nops.at(reqblocks), 0.),
32h(blocks_nops.at(reqblocks), 0.), gg(blocks_nops.at(reqblocks), 0.), s_s(blocks_nops.at(reqblocks), 0.),
33jssv(blocks_nops.at(reqblocks), 0.), jss(blocks_nops.at(reqblocks), 0.), jv(blocks_nops.at(reqblocks), 0.),
34vij(blocks_nops.at(reqblocks), 0.), eval(blocks_nops.at(reqblocks), 0.), evecc(blocks_nops.at(reqblocks), 0.),
35evalc(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
192EvolDF1::~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
206double 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
233double 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
243double 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
253double 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
269void 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
279gslpp::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
330gslpp::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
384gslpp::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
442gslpp::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
484gslpp::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
521gslpp::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
611gslpp::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
689gslpp::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
742gslpp::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
796gslpp::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
850gslpp::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
883gslpp::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
959gslpp::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
994gslpp::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
1030gslpp::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
1116gslpp::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
1154gslpp::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
1184gslpp::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
1215gslpp::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
1251gslpp::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
1359const 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
1453void 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, true, order_qed == QED0 ? false : true);
1477 eta = alsM / model.Als(mu, FULLNNNLO, true, 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, 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}
blocks_nops
std::map< std::string, uint > blocks_nops
Definition: EvolDF1.cpp:16
uint
unsigned int uint
Definition: EvolDF1.h:20
indices
unsigned int indices
Definition: EvolDF1.h:21
FULLNNNLO
@ FULLNNNLO
Definition: OrderScheme.h:40
FULLNLO
@ FULLNLO
Definition: OrderScheme.h:38
HV
@ HV
Definition: OrderScheme.h:22
LRI
@ LRI
Definition: OrderScheme.h:23
NDR
@ NDR
Definition: OrderScheme.h:21
QED1
@ QED1
Definition: OrderScheme.h:92
QED0
@ QED0
Definition: OrderScheme.h:91
QED2
@ QED2
Definition: OrderScheme.h:93
QCD1
@ QCD1
Definition: OrderScheme.h:76
QCD0
@ QCD0
Definition: OrderScheme.h:75
QCD2
@ QCD2
Definition: OrderScheme.h:77
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::vM131vi
std::map< std::vector< uint >, double > vM131vi[NF]
Definition: EvolDF1.h:331
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
EvolDF1::vM14vi
std::map< std::vector< uint >, double > vM14vi[NF]
Definition: EvolDF1.h:331
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::vM34vi
std::map< std::vector< uint >, double > vM34vi[NF]
Definition: EvolDF1.h:330
EvolDF1::vM2vi
std::map< std::vector< uint >, double > vM2vi[NF]
Definition: EvolDF1.h:329
EvolDF1::nfmin
uint nfmin
Definition: EvolDF1.h:336
EvolDF1::vM3vi
std::map< std::vector< uint >, double > vM3vi[NF]
Definition: EvolDF1.h:329
EvolDF1::vM5vi
std::map< std::vector< uint >, double > vM5vi[NF]
Definition: EvolDF1.h:329
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
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:331
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
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::evec
gslpp::matrix< double > evec
Definition: EvolDF1.h:339
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
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::vM133vi
std::map< std::vector< uint >, double > vM133vi[NF]
Definition: EvolDF1.h:331
EvolDF1::f_f_d
double f_f_d[2][F_iCacheSize]
Definition: EvolDF1.h:348
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::vM331vi
std::map< std::vector< uint >, double > vM331vi[NF]
Definition: EvolDF1.h:331
EvolDF1::nops
uint nops
Definition: EvolDF1.h:336
EvolDF1::MAls_cache
double MAls_cache
Definition: EvolDF1.h:343
EvolDF1::evecc
gslpp::matrix< gslpp::complex > evecc
Definition: EvolDF1.h:341
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::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::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::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::alsM_cache
double alsM_cache
Definition: EvolDF1.h:343
EvolDF1::~EvolDF1
virtual ~EvolDF1()
EvolDF1 destructor.
Definition: EvolDF1.cpp:192
EvolDF1::vM33vi
std::map< std::vector< uint >, double > vM33vi[NF]
Definition: EvolDF1.h:330
EvolDF1::vM43vi
std::map< std::vector< uint >, double > vM43vi[NF]
Definition: EvolDF1.h:330
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::evalc
gslpp::vector< gslpp::complex > evalc
Definition: EvolDF1.h:342
EvolDF1::vM11vi
std::map< std::vector< uint >, double > vM11vi[NF]
Definition: EvolDF1.h:330
EvolDF1::model
const StandardModel & model
Definition: EvolDF1.h:333
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::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::zeta
friend double gslpp_special_functions::zeta(int i)
EvolDF1::vM13vi
std::map< std::vector< uint >, double > vM13vi[NF]
Definition: EvolDF1.h:330
EvolDF1::f_f_c
int f_f_c[4][F_iCacheSize]
Definition: EvolDF1.h:347
EvolDF1::vM313vi
std::map< std::vector< uint >, double > vM313vi[NF]
Definition: EvolDF1.h:331
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
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::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
EvolDF1::vM1vi
std::map< std::vector< uint >, double > vM1vi[NF]
Definition: EvolDF1.h:329
EvolDF1::vM32vi
std::map< std::vector< uint >, double > vM32vi[NF]
Definition: EvolDF1.h:330
EvolDF1::vM41vi
std::map< std::vector< uint >, double > vM41vi[NF]
Definition: EvolDF1.h:331
EvolDF1::blocks
std::string blocks
Definition: EvolDF1.h:337
EvolDF1::vM6vi
std::map< std::vector< uint >, double > vM6vi[NF]
Definition: EvolDF1.h:330
EvolDF1::vM0vi
std::map< std::vector< uint >, double > vM0vi[NF]
Definition: EvolDF1.h:329
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
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
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
EvolDF1::vM23vi
std::map< std::vector< uint >, double > vM23vi[NF]
Definition: EvolDF1.h:330
EvolDF1::nfmax
uint nfmax
Definition: EvolDF1.h:336
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
EvolDF1::vM4vi
std::map< std::vector< uint >, double > vM4vi[NF]
Definition: EvolDF1.h:329
EvolDF1::vM113vi
std::map< std::vector< uint >, double > vM113vi[NF]
Definition: EvolDF1.h:331
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
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::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::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
EvolDF1::evec_i
gslpp::matrix< double > evec_i
Definition: EvolDF1.h:339
EvolDF1::vM31vi
std::map< std::vector< uint >, double > vM31vi[NF]
Definition: EvolDF1.h:330
EvolDF1::ai
std::map< uint, double > ai[NF]
Definition: EvolDF1.h:328
EvolDF1::eval
gslpp::vector< double > eval
Definition: EvolDF1.h:340
QCD::getAlsM
const double getAlsM() const
A get method to access the value of .
Definition: QCD.h:555
QCD::AboveTh
const double AboveTh(const double mu) const
The active flavour threshold above the scale as defined in QCD::Thresholds().
Definition: QCD.cpp:547
QCD::Nf
const double Nf(const double mu) const
The number of active flavour at scale .
Definition: QCD.cpp:571
QCD::getMAls
const double getMAls() const
A get method to access the mass scale at which the strong coupling constant measurement is provided.
Definition: QCD.h:564
QCD::CacheShift
void CacheShift(double cache[][5], int n) const
A member used to manage the caching for this class.
RGEvolutorNew::getEvol
const Expanded< gslpp::matrix< double > > & getEvol() const
Definition: RGEvolutorNew.cpp:37
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
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:509
StandardModel::Ale
const double Ale(double mu, orders order, bool Nf_thr=true) const
The running electromagnetic coupling in the scheme.
Definition: StandardModel/src/StandardModel.cpp:777
StandardModel::Beta_s
const double Beta_s(int nm, unsigned int nf) const
QCD beta function coefficients including QED corrections - eq. (36) hep-ph/0512066.
Definition: StandardModel/src/StandardModel.cpp:626
StandardModel::Als
const double Als(const double mu, const orders order, const bool Nf_thr, const bool qed_flag) const
The running QCD coupling in the scheme including QED corrections.
Definition: StandardModel.h:1088
StandardModel::Beta_e
const double Beta_e(int nm, unsigned int nf) const
QED beta function coefficients - eq. (36) hep-ph/0512066.
Definition: StandardModel/src/StandardModel.cpp:654
WilsonTemplateNew< gslpp::matrix< double > >::wilson
Expanded< gslpp::matrix< double > > wilson
Definition: WilsonTemplateNew.h:116
schemes
schemes
An enum type for regularization schemes.
Definition: OrderScheme.h:20
qed_orders
qed_orders
An enum type for qed_orders in electroweak.
Definition: OrderScheme.h:90
qcd_orders
qcd_orders
An enum type for qcd_orders in QCD.
Definition: OrderScheme.h:74