a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
EvolDF1nlep.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 "EvolDF1nlep.h"
9#include "StandardModel.h"
10
11EvolDF1nlep::EvolDF1nlep(unsigned int dim_i, schemes scheme, orders order, orders_qed order_qed, const StandardModel& model)
12: RGEvolutor(dim_i, scheme, order, order_qed), model(model), V(dim_i, 0.), Vi(dim_i, 0.),
13gs(dim_i, 0.), Js(dim_i, 0.), ge0(dim_i, 0.), K0(dim_i, 0.), ge11(dim_i, 0.), K11(dim_i, 0.),
14JsK0V(dim_i, 0.), ViK0Js(dim_i, 0.), Gamma_s0T(dim_i, 0.), Gamma_s1T(dim_i, 0.),
15Gamma_eT(dim_i, 0.), Gamma_seT(dim_i, 0.), JsV(dim_i, 0.), ViJs(dim_i, 0.), K0V(dim_i, 0.),
16ViK0(dim_i, 0.), ge0sing(dim_i, 0.), K0sing(dim_i, 0.), K0singV(dim_i, 0.), K11V(dim_i, 0.),
17ViK11(dim_i, 0.), ge11sing(dim_i, 0.), K11sing(dim_i, 0.),K11singV(dim_i, 0.),
18JsK0singV(dim_i, 0.), e(dim_i, 0.), dim(dim_i)
19{
20
21 int nu = 0, nd = 0;
22 double b0 = 0., b1 = 0.;
23
24 /* L=3 --> u,d,s,c (nf=3) L=2 --> u,d,s,c (nf=4) L=1 --> u,d,s,c,b (nf=5) L=0 --> u,d,s,c,b,t (nf=6)*/
25 for (int L = 3; L>-1; L--) {
26
27 b0 = model.Beta0(6 - L);
28 b1 = model.Beta1(6 - L);
29
30 if (L == 3) {
31 nd = 2;
32 nu = 1;
33 }
34 if (L == 2) {
35 nd = 2;
36 nu = 2;
37 }
38 if (L == 1) {
39 nd = 3;
40 nu = 2;
41 }
42 if (L == 0) {
43 nd = 3;
44 nu = 3;
45 }
46
47 Gamma_s0T = AnomalousDimension_nlep_S(LO, nu, nd).transpose();
48 Gamma_s1T = AnomalousDimension_nlep_S(NLO, nu, nd).transpose();
49 Gamma_eT = AnomalousDimension_nlep_EM(LO, nu, nd).transpose();
50 Gamma_seT = AnomalousDimension_nlep_EM(NLO, nu, nd).transpose();
51
52 AnomalousDimension_nlep_S(LO, nu, nd).transpose().eigensystem(V, e);
53 Vi = V.inverse();
54
55 /* magic numbers of U0 */
56 for (unsigned int i = 0; i < dim; i++) {
57 a[L][i] = e(i).real() / 2. / b0;
58 for (unsigned int j = 0; j < dim; j++) {
59 for (unsigned int k = 0; k < dim; k++) {
60 b[L][i][j][k] = V(i, k).real() * Vi(k, j).real();
61 }
62 }
63 }
64
65 gs = (b1 / 2. / b0 / b0) * Vi * Gamma_s0T * V - (1. / 2. / b0) * Vi * Gamma_s1T * V;
66 for (unsigned int i = 0; i < dim; i++) {
67 for (unsigned int j = 0; j < dim; j++) {
68 gs.assign(i, j, gs(i, j) / (1. + a[L][i] - a[L][j]));
69 }
70 }
71 Js = V * gs * Vi;
72
73 /*magic numbers related to Js*/
74 JsV = Js*V;
75 ViJs = Vi*Js;
76 for (unsigned int i = 0; i < dim; i++) {
77 for (unsigned int j = 0; j < dim; j++) {
78 for (unsigned int k = 0; k < dim; k++) {
79 c[L][i][j][k] = JsV(i, k).real() * Vi(k, j).real();
80 d[L][i][j][k] = -V(i, k).real() * ViJs(k, j).real();
81 }
82 }
83 }
84
85 ge0 = (1. / 2. / b0) * Vi * Gamma_eT * V;
86 for (unsigned int i = 0; i < dim; i++) {
87 for (unsigned int j = 0; j < dim; j++) {
88 if (fabs(a[L][j] + 1. - a[L][i]) > 0.00000000001) {
89 ge0.assign(i, j, ge0(i, j) / (1. - a[L][i] + a[L][j]));
90 } else {
91 ge0sing.assign(i, j, ge0(i, j) / 2. / b0);
92 ge0.assign(i, j, 0.);
93 }
94 }
95 }
96 K0 = V * ge0 * Vi;
97 K0sing = V * ge0sing * Vi;
98
99 /*magic numbers related to K0*/
100 K0V = K0*V;
101 ViK0 = Vi * K0;
102 K0singV = K0sing*V;
103 for (unsigned int i = 0; i < dim; i++) {
104 for (unsigned int j = 0; j < dim; j++) {
105 for (unsigned int k = 0; k < dim; k++) {
106 m[L][i][j][k] = K0V(i, k).real() * Vi(k, j).real();
107 n[L][i][j][k] = -V(i, k).real() * ViK0(k, j).real();
108 mn[L][i][j][k] = K0singV(i, k).real() * Vi(k, j).real();
109 }
110 }
111 }
112
113 ge11 = Gamma_seT - (b1 / b0) * Gamma_eT + Gamma_eT * Js - Js * Gamma_eT;
114 ge11 = Vi * ge11;
115 ge11 = ge11 * V;
116 for (unsigned int i = 0; i < dim; i++) {
117 for (unsigned int j = 0; j < dim; j++) {
118 if (fabs(a[L][j] - a[L][i]) > 0.00000000001) {
119 ge11.assign(i, j, ge11(i, j) / (2. * b0 * (a[L][j] - a[L][i])));
120 } else {
121 ge11sing.assign(i, j, ge11(i, j) / 2. / b0);
122 ge11.assign(i, j, 0.);
123 }
124 }
125 }
126 K11 = V * ge11 * Vi;
127 K11sing = V * ge11sing * Vi;
128 /*magic numbers related to K11*/
129 K11V = K11 * V;
130 ViK11 = Vi * K11;
131 K11singV = K11sing * V;
132 for (unsigned int i = 0; i < dim; i++) {
133 for (unsigned int j = 0; j < dim; j++) {
134 for (unsigned int k = 0; k < dim; k++) {
135 o[L][i][j][k] = K11V(i, k).real() * Vi(k, j).real();
136 p[L][i][j][k] = -V(i, k).real() * ViK11(k, j).real();
137 op[L][i][j][k] = K11singV(i, k).real() * Vi(k, j).real();
138 }
139 }
140 }
141
142 /*magic numbers related to K12 and K13*/
143 JsK0V = Js * K0 * V;
144 ViK0Js = Vi * K0 * Js;
145 JsK0singV = Js * K0sing * V;
146 for (unsigned int i = 0; i < dim; i++) {
147 for (unsigned int j = 0; j < dim; j++) {
148 for (unsigned int k = 0; k < dim; k++) {
149 q[L][i][j][k] = JsK0V(i, k).real() * Vi(k, j).real();
150 r[L][i][j][k] = V(i, k).real() * ViK0Js(k, j).real();
151 s[L][i][j][k] = -JsV(i, k).real() * ViK0(k, j).real();
152 t[L][i][j][k] = -K0V(i, k).real() * ViJs(k, j).real();
153 qq[L][i][j][k] = JsK0singV(i, k).real() * Vi(k, j).real();
154 rr[L][i][j][k] = -K0singV(i, k).real() * ViJs(k, j).real();
155 }
156 }
157 }
158 }
159}
160
162{
163}
164
165gslpp::matrix<double> EvolDF1nlep::AnomalousDimension_nlep_S(orders order, unsigned int n_u, unsigned int n_d) const
166{
167
168 /* anomalous dimension related to Delta F = 1 operators in Buras basis, hep-ph/9512380v1 */
169
170 /*gamma(row, column) leading order*/
171
172 unsigned int nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
173 gslpp::matrix<double> gammaDF1(dim, 0.);
174
175 switch (order) {
176
177 case LO:
178
179 gammaDF1(0, 0) = -2.;
180 gammaDF1(0, 1) = 6.;
181
182
183 gammaDF1(1, 0) = 6.;
184 gammaDF1(1, 1) = -2.;
185 gammaDF1(1, 2) = -2. / 9.;
186 gammaDF1(1, 3) = 2. / 3.;
187 gammaDF1(1, 4) = -2. / 9.;
188 gammaDF1(1, 5) = 2. / 3.;
189
190 gammaDF1(2, 2) = -22. / 9.;
191 gammaDF1(2, 3) = 22. / 3.;
192 gammaDF1(2, 4) = -4. / 9.;
193 gammaDF1(2, 5) = 4. / 3.;
194
195 gammaDF1(3, 2) = 6. - 2. / 9. * nf;
196 gammaDF1(3, 3) = -2. + 2. / 3. * nf;
197 gammaDF1(3, 4) = -2. / 9. * nf;
198 gammaDF1(3, 5) = 2. / 3. * nf;
199
200 gammaDF1(4, 4) = 2.;
201 gammaDF1(4, 5) = -6.;
202
203 gammaDF1(5, 2) = -2. / 9. * nf;
204 gammaDF1(5, 3) = 2. / 3. * nf;
205 gammaDF1(5, 4) = -2. / 9. * nf;
206 gammaDF1(5, 5) = -16. + 2. / 3. * nf;
207
208 gammaDF1(6, 6) = 2.;
209 gammaDF1(6, 7) = -6.;
210
211 gammaDF1(7, 2) = -2. / 9. * (n_u - n_d / 2.);
212 gammaDF1(7, 3) = 2. / 3. * (n_u - n_d / 2.);
213 gammaDF1(7, 4) = -2. / 9. * (n_u - n_d / 2.);
214 gammaDF1(7, 5) = 2. / 3. * (n_u - n_d / 2.);
215 gammaDF1(7, 7) = -16.;
216
217 gammaDF1(8, 2) = 2. / 9.;
218 gammaDF1(8, 3) = -2. / 3.;
219 gammaDF1(8, 4) = 2. / 9.;
220 gammaDF1(8, 5) = -2. / 3.;
221 gammaDF1(8, 8) = -2.;
222 gammaDF1(8, 9) = 6.;
223
224 gammaDF1(9, 2) = -2. / 9. * (n_u - n_d / 2.);
225 gammaDF1(9, 3) = 2. / 3. * (n_u - n_d / 2.);
226 gammaDF1(9, 4) = -2. / 9. * (n_u - n_d / 2.);
227 gammaDF1(9, 5) = 2. / 3. * (n_u - n_d / 2.);
228 gammaDF1(9, 8) = 6.;
229 gammaDF1(9, 9) = -2.;
230
231 break;
232
233 case NLO:
234
235 if (!(nf == 3 || nf == 4 || nf == 5 || nf == 6)) {
236 throw std::runtime_error("EvolDF1nlep::AnomalousDimension_nlep_S("
237 "orders order, unsigned int n_u, unsigned int n_d) " " wrong number of flavour");
238 }
239
240 /*gamma(riga, colonna) next to leading order*/
241
242 gammaDF1(0, 0) = -21. / 2. - 2. / 9. * nf;
243 gammaDF1(0, 1) = 7. / 2. + 2. / 3. * nf;
244 gammaDF1(0, 2) = 79. / 9.;
245 gammaDF1(0, 3) = -7. / 3.;
246 gammaDF1(0, 4) = -65. / 9.;
247 gammaDF1(0, 5) = -7. / 3.;
248
249
250 gammaDF1(1, 0) = 7. / 2. + 2. / 3. * nf;
251 gammaDF1(1, 1) = -21. / 2. - 2. / 9. * nf;
252 gammaDF1(1, 2) = -202. / 243.;
253 gammaDF1(1, 3) = 1354. / 81.;
254 gammaDF1(1, 4) = -1192. / 243.;
255 gammaDF1(1, 5) = 904. / 81.;
256
257 gammaDF1(2, 2) = -5911. / 486. + 71. / 9. * nf;
258 gammaDF1(2, 3) = 5983. / 162. + 1. / 3. * nf;
259 gammaDF1(2, 4) = -2384. / 243. - 71. / 9. * nf;
260 gammaDF1(2, 5) = 1808. / 81. - 1. / 3. * nf;
261
262 gammaDF1(3, 2) = 379. / 18. + 56. / 243. * nf;
263 gammaDF1(3, 3) = -91. / 6. + 808. / 81. * nf;
264 gammaDF1(3, 4) = -130. / 9. - 502. / 243. * nf;
265 gammaDF1(3, 5) = -14. / 3. + 646. / 81. * nf;
266
267 gammaDF1(4, 2) = -61. / 9. * nf;
268 gammaDF1(4, 3) = -11. / 3. * nf;
269 gammaDF1(4, 4) = 71. / 3. + 61. / 9. * nf;
270 gammaDF1(4, 5) = -99. + 11. / 3. * nf;
271
272 gammaDF1(5, 2) = -682. / 243. * nf;
273 gammaDF1(5, 3) = 106. / 81. * nf;
274 gammaDF1(5, 4) = -225. / 2. + 1676. / 243. * nf;
275 gammaDF1(5, 5) = -1343. / 6. + 1348. / 81. * nf;
276
277 gammaDF1(6, 2) = -61. / 9. * (n_u - n_d / 2.);
278 gammaDF1(6, 3) = -11. / 3. * (n_u - n_d / 2.);
279 gammaDF1(6, 4) = 83. / 9. * (n_u - n_d / 2.);
280 gammaDF1(6, 5) = -11. / 3. * (n_u - n_d / 2.);
281 gammaDF1(6, 6) = 71. / 3. - 22. / 9. * nf;
282 gammaDF1(6, 7) = -99. + 22. / 3. * nf;
283
284 gammaDF1(7, 2) = -682. / 243. * (n_u - n_d / 2.);
285 gammaDF1(7, 3) = 106. / 81. * (n_u - n_d / 2.);
286 gammaDF1(7, 4) = 704. / 243. * (n_u - n_d / 2.);
287 gammaDF1(7, 5) = 736. / 81. * (n_u - n_d / 2.);
288 gammaDF1(7, 6) = -225. / 2. + 4 * nf;
289 gammaDF1(7, 7) = -1343. / 6. + 68. / 9. * nf;
290
291 gammaDF1(8, 2) = 202. / 243. + 73. / 9. * (n_u - n_d / 2.);
292 gammaDF1(8, 3) = -1354. / 81. - 1. / 3. * (n_u - n_d / 2.);
293 gammaDF1(8, 4) = 1192. / 243. - 71. / 9. * (n_u - n_d / 2.);
294 gammaDF1(8, 5) = -904. / 81. - 1. / 3. * (n_u - n_d / 2.);
295 gammaDF1(8, 8) = -21. / 2. - 2. / 9. * nf;
296 gammaDF1(8, 9) = 7. / 2. + 2. / 3. * nf;
297
298 gammaDF1(9, 2) = -79. / 9. - 106. / 243. * (n_u - n_d / 2.);
299 gammaDF1(9, 3) = 7. / 3. + 826. / 81. * (n_u - n_d / 2.);
300 gammaDF1(9, 4) = 65. / 9. - 502. / 243. * (n_u - n_d / 2.);
301 gammaDF1(9, 5) = 7. / 3. + 646. / 81. * (n_u - n_d / 2.);
302 gammaDF1(9, 8) = 7. / 2. + 2. / 3. * nf;
303 gammaDF1(9, 9) = -21. / 2. - 2. / 9. * nf;
304
305 break;
306
307 default:
308 std::stringstream out;
309 out << order;
310 throw std::runtime_error("EvolDF1nlep::AnomalousDimension_nlep_S("
311 "orders order, unsigned int n_u, unsigned int n_d) "
312 + out.str() + " not implemented");
313
314 }
315
316 return (gammaDF1);
317
318}
319
320gslpp::matrix<double> EvolDF1nlep::AnomalousDimension_nlep_EM(orders order, unsigned int n_u, unsigned int n_d) const
321{
322
323 /* anomalous dimension related to Buras operators hep-ph/9512380v1 */
324 /*gamma(riga, colonna) leading order*/
325 unsigned int nf = n_u + n_d; /*n_u\d = active type up/down flavor d.o.f.*/
326 gslpp::matrix<double> gammaDF1(dim, 0.);
327
328 switch (order) {
329
330 case LO:
331
332 gammaDF1(0, 0) = -8. / 3.;
333 gammaDF1(0, 6) = 16. / 9.;
334 gammaDF1(0, 8) = 16. / 9.;
335
336 gammaDF1(1, 1) = -8. / 3.;
337 gammaDF1(1, 6) = 16. / 27.;
338 gammaDF1(1, 8) = 16. / 27.;
339
340 gammaDF1(2, 6) = -16. / 27. + 16. / 9. * (n_u - n_d / 2.);
341 gammaDF1(2, 8) = -88. / 27. + 16. / 9. * (n_u - n_d / 2.);
342
343 gammaDF1(3, 6) = -16. / 9. + 16. / 27. * (n_u - n_d / 2.);
344 gammaDF1(3, 8) = -16. / 9. + 16. / 27. * (n_u - n_d / 2.);
345 gammaDF1(3, 9) = -8. / 3.;
346
347 gammaDF1(4, 6) = 8. / 3. + 16. / 9. * (n_u - n_d / 2.);
348 gammaDF1(4, 8) = 16. / 9. * (n_u - n_d / 2.);
349
350 gammaDF1(5, 6) = 16. / 27. * (n_u - n_d / 2.);
351 gammaDF1(5, 7) = 8. / 3.;
352 gammaDF1(5, 8) = 16. / 27. * (n_u - n_d / 2.);
353
354 gammaDF1(6, 4) = 4. / 3.;
355 gammaDF1(6, 6) = 4. / 3. + 16. / 9. * (n_u + n_d / 4.);
356 gammaDF1(6, 8) = 16. / 9. * (n_u + n_d / 4.);
357
358 gammaDF1(7, 5) = 4. / 3.;
359 gammaDF1(7, 6) = 16. / 27. * (n_u + n_d / 4.);
360 gammaDF1(7, 7) = 4. / 3.;
361 gammaDF1(7, 8) = 16. / 27. * (n_u + n_d / 4.);
362
363 gammaDF1(8, 2) = -4. / 3.;
364 gammaDF1(8, 6) = 8. / 27. + 16. / 9. * (n_u + n_d / 4.);
365 gammaDF1(8, 8) = -28. / 27. + 16. / 9. * (n_u + n_d / 4.);
366
367 gammaDF1(9, 3) = -4. / 3.;
368 gammaDF1(9, 6) = 8. / 9. + 16. / 27. * (n_u + n_d / 4.);
369 gammaDF1(9, 8) = 8. / 9. + 16. / 27. * (n_u + n_d / 4.);
370 gammaDF1(9, 9) = -4. / 3.;
371
372 break;
373
374 case NLO:
375
376 if (!(nf == 3 || nf == 4 || nf == 5 || nf == 6)) {
377 throw std::runtime_error("EvolDF1nlep::AnomalousDimension_nlep_EM("
378 "orders order, unsigned int n_u, unsigned int n_d) " " wrong number of flavour");
379 }
380
381 /*gamma(riga, colonna) next to leading order*/
382
383 gammaDF1(0, 0) = 194. / 9.;
384 gammaDF1(0, 1) = -2. / 3.;
385 gammaDF1(0, 2) = -88. / 243.;
386 gammaDF1(0, 3) = 88. / 81.;
387 gammaDF1(0, 4) = -88. / 243.;
388 gammaDF1(0, 5) = 88. / 81.;
389 gammaDF1(0, 6) = 152. / 27.;
390 gammaDF1(0, 7) = 40. / 9.;
391 gammaDF1(0, 8) = 136. / 27.;
392 gammaDF1(0, 9) = 56. / 9.;
393
394 gammaDF1(1, 0) = 25. / 3.;
395 gammaDF1(1, 1) = -49. / 9.;
396 gammaDF1(1, 2) = -556. / 729.;
397 gammaDF1(1, 3) = 556. / 243.;
398 gammaDF1(1, 4) = -556. / 729.;
399 gammaDF1(1, 5) = 556. / 243.;
400 gammaDF1(1, 6) = -484. / 729.;
401 gammaDF1(1, 7) = -124. / 27.;
402 gammaDF1(1, 8) = -3148. / 729.;
403 gammaDF1(1, 9) = 172. / 27.;
404
405 gammaDF1(2, 2) = 1690. / 729. - 136. / 243. * (n_u - n_d / 2.);
406 gammaDF1(2, 3) = -1690. / 243. + 136. / 81. * (n_u - n_d / 2.);
407 gammaDF1(2, 4) = 232. / 729. - 136. / 243. * (n_u - n_d / 2.);
408 gammaDF1(2, 5) = -232. / 243. + 136. / 81. * (n_u - n_d / 2.);
409 gammaDF1(2, 6) = 3136. / 729. + 104. / 27. * (n_u - n_d / 2.);
410 gammaDF1(2, 7) = 64. / 27. + 88. / 9. * (n_u - n_d / 2.);
411 gammaDF1(2, 8) = 20272. / 729. + 184. / 27. * (n_u - n_d / 2.);
412 gammaDF1(2, 9) = -112. / 27. + 8. / 9. * (n_u - n_d / 2.);
413
414 gammaDF1(3, 2) = -641. / 243. - 388. / 729. * n_u + 32. / 729. * n_d;
415 gammaDF1(3, 3) = -655. / 81. + 388. / 243. * n_u - 32. / 243. * n_d;
416 gammaDF1(3, 4) = 88. / 243. - 388. / 729 * n_u + 32. / 729. * n_d;
417 gammaDF1(3, 5) = -88. / 81. + 388. / 243. * n_u - 32. / 243. * n_d;
418 gammaDF1(3, 6) = -152. / 27. + 3140. / 729. * n_u + 656. / 729. * n_d;
419 gammaDF1(3, 7) = -40. / 9. - 100. / 27. * n_u - 16. / 27. * n_d;
420 gammaDF1(3, 8) = 170. / 27. + 908. / 729. * n_u + 1232. / 729. * n_d;
421 gammaDF1(3, 9) = -14. / 3. + 148. / 27. * n_u - 80. / 27 * n_d;
422
423 gammaDF1(4, 2) = -136. / 243. * (n_u - n_d / 2.);
424 gammaDF1(4, 3) = 136. / 81. * (n_u - n_d / 2.);
425 gammaDF1(4, 4) = -2. - 136. / 243. * (n_u - n_d / 2.);
426 gammaDF1(4, 5) = 6. + 136. / 81. * (n_u - n_d / 2.);
427 gammaDF1(4, 6) = -232. / 9. + 104. / 27. * (n_u - n_d / 2.);
428 gammaDF1(4, 7) = 40. / 3. + 88. / 9. * (n_u - n_d / 2.);
429 gammaDF1(4, 8) = 184. / 27. * (n_u - n_d / 2.);
430 gammaDF1(4, 9) = 8. / 9. * (n_u - n_d / 2.);
431
432 gammaDF1(5, 2) = -748. / 729. * n_u + 212. / 729. * n_d;
433 gammaDF1(5, 3) = 748. / 243. * n_u - 212. / 243. * n_d;
434 gammaDF1(5, 4) = 3. - 748. / 729. * n_u + 212. / 729. * n_d;
435 gammaDF1(5, 5) = 7. + 748. / 243. * n_u - 212. / 243. * n_d;
436 gammaDF1(5, 6) = -2. - 5212. / 729. * n_u + 4832. / 729. * n_d;
437 gammaDF1(5, 7) = 182. / 9. + 188. / 27. * n_u - 160. / 27. * n_d;
438 gammaDF1(5, 8) = -2260. / 729. * n_u + 2816. / 729. * n_d;
439 gammaDF1(5, 9) = -140. / 27. * n_u + 64. / 27. * n_d;
440
441 gammaDF1(6, 2) = -136. / 243. * (n_u + n_d / 4.);
442 gammaDF1(6, 3) = 136. / 81. * (n_u + n_d / 4.);
443 gammaDF1(6, 4) = -116. / 9. - 136. / 243. * (n_u + n_d / 4.);
444 gammaDF1(6, 5) = 20. / 3. + 136. / 81. * (n_u + n_d / 4.);
445 gammaDF1(6, 6) = -134. / 9. + 104. / 27. * (n_u + n_d / 4.);
446 gammaDF1(6, 7) = 38. / 3. + 88. / 9. * (n_u + n_d / 4.);
447 gammaDF1(6, 8) = 184. / 27. * (n_u + n_d / 4.);
448 gammaDF1(6, 9) = 8. / 9. * (n_u + n_d / 4.);
449
450 gammaDF1(7, 2) = -748. / 729. * n_u - 106. / 729. * n_d;
451 gammaDF1(7, 3) = 748. / 243. * n_u + 106. / 243. * n_d;
452 gammaDF1(7, 4) = -1. - 748. / 729. * n_u - 106. / 729. * n_d;
453 gammaDF1(7, 5) = 91. / 9. + 748. / 243. * n_u + 106. / 243. * n_d;
454 gammaDF1(7, 6) = 2. - 5212. / 729. * n_u - 2416. / 729. * n_d;
455 gammaDF1(7, 7) = 154. / 9. + 188. / 27. * n_u + 80. / 27. * n_d;
456 gammaDF1(7, 8) = -2260. / 729. * n_u - 1408. / 729. * n_d;
457 gammaDF1(7, 9) = -140. / 27. * n_u - 32. / 27. * n_d;
458
459 gammaDF1(8, 2) = 7012. / 729. - 136. / 243. * (n_u + n_d / 4.);
460 gammaDF1(8, 3) = 764. / 243. + 136. / 81. * (n_u + n_d / 4.);
461 gammaDF1(8, 4) = -116. / 729. - 136. / 243. * (n_u + n_d / 4.);
462 gammaDF1(8, 5) = 116. / 243. + 136. / 81. * (n_u + n_d / 4.);
463 gammaDF1(8, 6) = -1568. / 729. + 104. / 27. * (n_u + n_d / 4.);
464 gammaDF1(8, 7) = -32. / 27. + 88. / 9. * (n_u + n_d / 4.);
465 gammaDF1(8, 8) = 5578. / 729. + 184. / 27. * (n_u + n_d / 4.);
466 gammaDF1(8, 9) = 38. / 27. + 8. / 9. * (n_u + n_d / 4.);
467
468 gammaDF1(9, 2) = 1333. / 243. - 388. / 729. * n_u - 16. / 729. * n_d;
469 gammaDF1(9, 3) = 107. / 81. + 388. / 243. * n_u + 16. / 243. * n_d;
470 gammaDF1(9, 4) = -44. / 243. - 388. / 729. * n_u - 16. / 729. * n_d;
471 gammaDF1(9, 5) = 44. / 81. + 388. / 243. * n_u + 16. / 243. * n_d;
472 gammaDF1(9, 6) = 76. / 27. + 3140. / 729. * n_u - 328. / 729. * n_d;
473 gammaDF1(9, 7) = 20. / 9. - 100. / 27. * n_u + 8. / 27. * n_d;
474 gammaDF1(9, 8) = 140. / 27. + 908. / 729. * n_u - 616. / 729. * n_d;
475 gammaDF1(9, 9) = -28. / 9. + 148. / 27. * n_u + 40. / 27. * n_d;
476
477 break;
478
479 default:
480 std::stringstream out;
481 out << order;
482 throw std::runtime_error("EvolDF1nlep::AnomalousDimension_nlep_EM("
483 "orders order, unsigned int n_u, unsigned int n_d) "
484 + out.str() + " not implemented");
485
486 }
487
488 return (gammaDF1);
489
490}
491
492gslpp::matrix<double> EvolDF1nlep::Df1threshold_deltarsT(double nf) const
493{
494
495 gslpp::matrix <double> delta_rsT(dim, 0.);
496
497 delta_rsT(2, 3) = 5. / 27.;
498 delta_rsT(2, 5) = 5. / 27.;
499 delta_rsT(3, 3) = -5. / 9.;
500 delta_rsT(4, 5) = -5. / 9.;
501 delta_rsT(4, 3) = 5. / 27.;
502 delta_rsT(4, 5) = 5. / 27.;
503 delta_rsT(5, 3) = -5. / 9.;
504 delta_rsT(5, 5) = -5. / 9.;
505
506 if (nf == 3. || nf == 5.) {
507
508 delta_rsT(2, 7) = -5. / 54.;
509 delta_rsT(2, 9) = -5. / 54.;
510 delta_rsT(3, 7) = 5. / 18.;
511 delta_rsT(3, 9) = 5. / 18.;
512 delta_rsT(4, 7) = -5. / 54.;
513 delta_rsT(4, 9) = -5. / 54.;
514 delta_rsT(5, 7) = 5. / 18.;
515 delta_rsT(5, 9) = 5. / 18.;
516 }
517
518
519
520 else {
521
522 delta_rsT(2, 7) = 5. / 27.;
523 delta_rsT(2, 9) = 5. / 27.;
524 delta_rsT(3, 7) = -5. / 9.;
525 delta_rsT(3, 9) = -5. / 9.;
526 delta_rsT(4, 7) = 5. / 27.;
527 delta_rsT(4, 9) = 5. / 27.;
528 delta_rsT(5, 7) = -5. / 9.;
529 delta_rsT(5, 9) = -5. / 9.;
530
531 }
532
533 return (delta_rsT);
534
535}
536
537gslpp::matrix<double> EvolDF1nlep::Df1threshold_deltareT(double nf) const
538{
539
540 gslpp::matrix<double> delta_reT(dim, 0.);
541
542 if (nf == 3. || nf == 5.) {
543
544 delta_reT(6, 2) = 20. / 27.;
545 delta_reT(6, 4) = 20. / 81.;
546 delta_reT(6, 4) = 20. / 27.;
547 delta_reT(6, 5) = 20. / 81.;
548 delta_reT(6, 6) = -10. / 27.;
549 delta_reT(6, 7) = -10. / 81.;
550 delta_reT(6, 8) = -10. / 27.;
551 delta_reT(6, 9) = -10. / 81.;
552 delta_reT(8, 2) = 20. / 27.;
553 delta_reT(8, 3) = 20. / 81.;
554 delta_reT(8, 4) = 20. / 27.;
555 delta_reT(8, 5) = 20. / 81.;
556 delta_reT(8, 6) = -10. / 27.;
557 delta_reT(8, 7) = -10. / 81.;
558 delta_reT(8, 8) = -10. / 27.;
559 delta_reT(8, 9) = -10. / 81.;
560
561 }
562 else {
563
564 delta_reT(6, 2) = -40. / 27.;
565 delta_reT(6, 3) = -40. / 81.;
566 delta_reT(6, 4) = -40. / 27.;
567 delta_reT(6, 5) = -40. / 81.;
568 delta_reT(6, 5) = -40. / 27.;
569 delta_reT(6, 6) = -40. / 81.;
570 delta_reT(6, 7) = -40. / 27.;
571 delta_reT(6, 8) = -40. / 81.;
572 delta_reT(8, 2) = -40. / 27.;
573 delta_reT(8, 3) = -40. / 81.;
574 delta_reT(8, 4) = -40. / 27.;
575 delta_reT(8, 5) = -40. / 81.;
576 delta_reT(8, 6) = -40. / 27.;
577 delta_reT(8, 7) = -40. / 81.;
578 delta_reT(8, 8) = -40. / 27.;
579 delta_reT(8, 9) = -40. / 81.;
580
581 }
582
583 return (delta_reT);
584
585}
586
587// Check if one could use a single order argument since effectively only one of the two is used
588gslpp::matrix<double>& EvolDF1nlep::Df1Evolnlep(double mu, double M, orders order, orders_qed order_qed, schemes scheme)
589{
590 switch (scheme) {
591 case NDR:
592 break;
593 case LRI:
594 case HV:
595 default:
596 std::stringstream out;
597 out << scheme;
598 throw std::runtime_error("EvolDF1nlep::Df1Evolnlep_EM(): scheme " + out.str()
599 + " not implemented ");
600 }
601
602 double alsMZ = model.getAlsMz();
603 double Mz = model.getMz();
604 double Ale = model.getAle();
605 if (alsMZ == alsMZ_cache && Mz == Mz_cache && Ale == Ale_cache) {
606 if (mu == this->mu && M == this->M && scheme == this->scheme && order_qed == NO_QED)
607 return (*Evol(order));
608
609 if (mu == this->mu && M == this->M && scheme == this->scheme && order_qed == NLO_QED11)
610 return (*Evol(order_qed));
611 }
612 alsMZ_cache = alsMZ;
613 Mz_cache = Mz;
614 Ale_cache = Ale;
615
616 if (M < mu) {
617 std::stringstream out;
618 out << "M = " << M << " < mu = " << mu;
619 throw out.str();
620 }
621
622 setScales(mu, M); // also assign evol to identity
623
624 double m_down = mu;
625 double m_up = model.AboveTh(m_down);
626 double nf = model.Nf(m_down);
627
628 while (m_up < M) {
629 Df1Evolnlep(m_down, m_up, nf, scheme);
630 Df1threshold_nlep(m_up, nf + 1.);
631 m_down = m_up;
632 m_up = model.AboveTh(m_down);
633 nf += 1.;
634 }
635
636 Df1Evolnlep(m_down, M, nf, scheme);
637
638 if (order_qed != NO_QED) {
639 return (*Evol(order_qed));
640 }
641 else {
642 return (*Evol(order));
643 }
644
645}
646
647void EvolDF1nlep::Df1Evolnlep(double mu, double M, double nf, schemes scheme)
648{
649
650 gslpp::matrix<double> resLO(dim, 0.), resNLO(dim, 0.), resLO_ew(dim, 0.), resNLO_QED(dim, 0.);
651
652 int L = 6 - (int) nf;
653 double alsM = model.Als(M) / 4. / M_PI;
654 double alsmu = model.Als(mu) / 4. / M_PI;
655 double ale = model.getAle() / 4. / M_PI;
656
657 double eta = alsM / alsmu;
658
659 //The calculation of matrix entries for the various orders should be moved inside the switch to avoid wasting time on unnecessary calculations
660
661 for (unsigned int k = 0; k < dim; k++) {
662 double etap = pow(eta, a[L][k]);
663 for (unsigned int i = 0; i < dim; i++) {
664 for (unsigned int j = 0; j < dim; j++) {
665
666 resLO(i, j) += b[L][i][j][k] * etap;
667
668 resNLO(i, j) += c[L][i][j][k] * etap * alsmu;
669 resNLO(i, j) += d[L][i][j][k] * etap * alsM;
670
671 resLO_ew(i, j) += m[L][i][j][k] * etap * ale / alsmu;
672 resLO_ew(i, j) += n[L][i][j][k] * etap * ale / alsM;
673 resLO_ew(i, j) += mn[L][i][j][k] * etap * ale / alsM * log(eta);
674
675 resNLO_QED(i, j) += o[L][i][j][k] * etap * ale;
676 resNLO_QED(i, j) += p[L][i][j][k] * etap * ale;
677 resNLO_QED(i, j) += op[L][i][j][k] * etap * ale * log(eta);
678
679 resNLO_QED(i, j) += q[L][i][j][k] * etap * ale;
680 resNLO_QED(i, j) += r[L][i][j][k] * etap * ale;
681 resNLO_QED(i, j) += s[L][i][j][k] * etap * ale / eta;
682 resNLO_QED(i, j) += t[L][i][j][k] * etap * ale * eta;
683 resNLO_QED(i, j) += qq[L][i][j][k] * etap * ale * log(eta);
684 resNLO_QED(i, j) += rr[L][i][j][k] * etap * ale * log(eta);
685
686 // unreasonable large entries: this fixes the issue, weird numerical instability needs investigation
687 if(L==3){
688 if((i==6) and (j==6)){
689 resNLO_QED(i, j) -= op[L][i][j][k] * etap * ale * log(eta);
690 resNLO_QED(i, j) -= qq[L][i][j][k] * etap * ale * log(eta);
691 resNLO_QED(i, j) -= rr[L][i][j][k] * etap * ale * log(eta);
692 }
693 if((i==7) and ((j==6) or (j==7))){
694 resNLO_QED(i, j) -= op[L][i][j][k] * etap * ale * log(eta);
695 resNLO_QED(i, j) -= qq[L][i][j][k] * etap * ale * log(eta);
696 resNLO_QED(i, j) -= rr[L][i][j][k] * etap * ale * log(eta);
697 }
698 }
699 }
700 }
701 }
702
703 switch (order_qed) {
704 case NLO_QED11:
705 *elem[NLO_QED11] = (*elem[NLO]) * resLO_ew +
706 (*elem[NLO_QED11]) * resLO + (*elem[LO]) * resNLO_QED +
707 (*elem[LO_QED]) * resNLO;
708 case LO_QED:
709 *elem[LO_QED] = (*elem[LO]) * resLO_ew + (*elem[LO_QED]) * resLO;
710 break;
711 default:
712 throw std::runtime_error("Error in EvolDF1nlep::Df1Evolnlep()");
713 }
714
715 switch (order) {
716 case NNLO:
717 *elem[NNLO] = 0.;
718 case NLO:
719 *elem[NLO] = (*elem[LO]) * resNLO + (*elem[NLO]) * resLO;
720 case LO:
721 *elem[LO] = (*elem[LO]) * resLO;
722 break;
723 default:
724 throw std::runtime_error("Error in EvolDF1nlep::Df1Evolnlep()");
725 }
726}
727
728void EvolDF1nlep::Df1threshold_nlep(double M, double nf)
729{
730
731 gslpp::matrix<double> drsT(dim, 0.), dreT(dim, 0.);
732
733 double alsM = model.Als(M) / 4. / M_PI;
734 double ale = model.getAle() / 4. / M_PI;
735
736 drsT = alsM * Df1threshold_deltarsT(nf);
737 dreT = ale * Df1threshold_deltareT(nf);
738
739 switch (order_qed) {
740 case NLO_QED11:
741 *elem[NLO_QED11] += (*elem[LO]) * dreT + (*elem[LO_QED]) * drsT;
742 break;
743 default:
744 throw std::runtime_error("Error in EvolDF1nlep::Df1threshold_nlep()");
745 }
746
747 switch (order) {
748 case NNLO:
749 *elem[NNLO] = 0.;
750 case NLO:
751 *elem[NLO] += (*elem[LO]) * drsT;
752 break;
753 default:
754 throw std::runtime_error("Error in EvolDF1nlep::Df1threshold_nlep()");
755 }
756
757
758}
759
760gslpp::matrix<double>& EvolDF1nlep::Df1Evolnlep3flav(double mu, double M, orders order, orders_qed order_qed, schemes scheme)
761{
762 switch (scheme) {
763 case NDR:
764 break;
765 case LRI:
766 case HV:
767 default:
768 std::stringstream out;
769 out << scheme;
770 throw std::runtime_error("EvolDF1nlep::Df1Evolnlep_EM(): scheme " + out.str()
771 + " not implemented ");
772 }
773
774 double alsMZ = model.getAlsMz();
775 double Mz = model.getMz();
776 double Ale = model.getAle();
777 if (alsMZ == alsMZ_cache && Mz == Mz_cache && Ale == Ale_cache) {
778 if (mu == this->mu && M == this->M && scheme == this->scheme && order_qed == NO_QED)
779 return (*Evol(order));
780
781 if (mu == this->mu && M == this->M && scheme == this->scheme && order_qed == NLO_QED11)
782 return (*Evol(order_qed));
783 }
784
785 alsMZ_cache = alsMZ;
786 Mz_cache = Mz;
787 Ale_cache = Ale;
788
789 setScales(mu, M); // also assign evol to identity
790
791 Df1Evolnlep(mu, M, 3, scheme);
792
793 if (order_qed != NO_QED) {
794 return (*Evol(order_qed));
795 }
796 else {
797 return (*Evol(order));
798 }
799
800}
@ NNLO
Definition: OrderScheme.h:36
@ LO
Definition: OrderScheme.h:34
@ NLO
Definition: OrderScheme.h:35
@ HV
Definition: OrderScheme.h:22
@ LRI
Definition: OrderScheme.h:23
@ NDR
Definition: OrderScheme.h:21
@ NLO_QED11
Definition: OrderScheme.h:59
@ LO_QED
Definition: OrderScheme.h:58
@ NO_QED
Definition: OrderScheme.h:57
gslpp::matrix< gslpp::complex > K0singV
Definition: EvolDF1nlep.h:124
gslpp::matrix< gslpp::complex > K0sing
Definition: EvolDF1nlep.h:124
double o[4][10][10][10]
Definition: EvolDF1nlep.h:104
gslpp::matrix< gslpp::complex > Gamma_eT
Definition: EvolDF1nlep.h:123
double d[4][10][10][10]
Definition: EvolDF1nlep.h:103
gslpp::matrix< gslpp::complex > K0V
Definition: EvolDF1nlep.h:124
gslpp::matrix< gslpp::complex > ge11sing
Definition: EvolDF1nlep.h:124
const StandardModel & model
Definition: EvolDF1nlep.h:107
gslpp::matrix< gslpp::complex > JsK0singV
Definition: EvolDF1nlep.h:125
gslpp::matrix< gslpp::complex > Gamma_s1T
Definition: EvolDF1nlep.h:123
gslpp::matrix< gslpp::complex > Vi
Definition: EvolDF1nlep.h:122
double alsMZ_cache
Definition: EvolDF1nlep.h:128
gslpp::matrix< gslpp::complex > ViJs
Definition: EvolDF1nlep.h:123
gslpp::matrix< double > Df1threshold_deltareT(double nf) const
a method returning the matrix threshold for the QED penguins at the NLO
double m[4][10][10][10]
Definition: EvolDF1nlep.h:104
gslpp::matrix< gslpp::complex > V
Definition: EvolDF1nlep.h:122
double a[4][10]
Definition: EvolDF1nlep.h:103
double mn[4][10][10][10]
Definition: EvolDF1nlep.h:104
gslpp::matrix< gslpp::complex > Gamma_seT
Definition: EvolDF1nlep.h:123
gslpp::matrix< double > & Df1Evolnlep3flav(double mu, double M, orders order, orders_qed order_qed, schemes scheme=NDR)
gslpp::matrix< gslpp::complex > K11V
Definition: EvolDF1nlep.h:124
gslpp::matrix< double > & Df1Evolnlep(double mu, double M, orders order, orders_qed order_qed, schemes scheme=NDR)
a method returning the evolutor related to the high scale and the low scale
double Ale_cache
Definition: EvolDF1nlep.h:130
gslpp::matrix< gslpp::complex > gs
Definition: EvolDF1nlep.h:122
gslpp::matrix< double > AnomalousDimension_nlep_EM(orders order, unsigned int n_u, unsigned int n_d) const
a method returning the anomalous dimension matrix given in the standard basis
gslpp::matrix< gslpp::complex > K11
Definition: EvolDF1nlep.h:122
double r[4][10][10][10]
Definition: EvolDF1nlep.h:106
gslpp::matrix< gslpp::complex > ge0
Definition: EvolDF1nlep.h:122
double p[4][10][10][10]
Definition: EvolDF1nlep.h:105
gslpp::matrix< gslpp::complex > JsK0V
Definition: EvolDF1nlep.h:122
gslpp::matrix< gslpp::complex > ViK0
Definition: EvolDF1nlep.h:124
EvolDF1nlep(unsigned int dim, schemes scheme, orders order, orders_qed order_qed, const StandardModel &model)
EvolDF1nlep constructor.
Definition: EvolDF1nlep.cpp:11
gslpp::matrix< gslpp::complex > ge11
Definition: EvolDF1nlep.h:122
double Mz_cache
Definition: EvolDF1nlep.h:129
virtual ~EvolDF1nlep()
EvolDF1nlep destructor.
double op[4][10][10][10]
Definition: EvolDF1nlep.h:105
gslpp::matrix< gslpp::complex > JsV
Definition: EvolDF1nlep.h:123
double b[4][10][10][10]
Definition: EvolDF1nlep.h:103
gslpp::matrix< gslpp::complex > K11singV
Definition: EvolDF1nlep.h:125
gslpp::matrix< gslpp::complex > Js
Definition: EvolDF1nlep.h:122
double n[4][10][10][10]
Definition: EvolDF1nlep.h:104
gslpp::matrix< gslpp::complex > ViK11
Definition: EvolDF1nlep.h:124
double c[4][10][10][10]
Definition: EvolDF1nlep.h:103
double rr[4][10][10][10]
Definition: EvolDF1nlep.h:106
unsigned int dim
Definition: EvolDF1nlep.h:127
gslpp::matrix< gslpp::complex > ViK0Js
Definition: EvolDF1nlep.h:122
double q[4][10][10][10]
Definition: EvolDF1nlep.h:105
double qq[4][10][10][10]
Definition: EvolDF1nlep.h:105
gslpp::vector< gslpp::complex > e
Definition: EvolDF1nlep.h:126
void Df1threshold_nlep(double M, double nf)
a void type method for the implementation of the NLO threshold effects in the evolutor
gslpp::matrix< double > AnomalousDimension_nlep_S(orders order, unsigned int n_u, unsigned int n_d) const
a method returning the anomalous dimension matrix given in the standard basis
gslpp::matrix< gslpp::complex > K0
Definition: EvolDF1nlep.h:122
gslpp::matrix< gslpp::complex > ge0sing
Definition: EvolDF1nlep.h:124
gslpp::matrix< gslpp::complex > K11sing
Definition: EvolDF1nlep.h:125
gslpp::matrix< gslpp::complex > Gamma_s0T
Definition: EvolDF1nlep.h:123
gslpp::matrix< double > Df1threshold_deltarsT(double nf) const
a method returning the matrix threshold for the QCD penguins at the NLO
const double Beta1(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:606
const double Beta0(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:601
const double AboveTh(const double mu) const
The active flavour threshold above the scale as defined in QCD::Thresholds().
Definition: QCD.cpp:547
const double Nf(const double mu) const
The number of active flavour at scale .
Definition: QCD.cpp:571
A class for the RG evolutor of the Wilson coefficients.
Definition: RGEvolutor.h:24
double M
Definition: RGEvolutor.h:142
void setScales(double mu, double M)
Sets the upper and lower scale for the running of the Wilson Coefficients.
Definition: RGEvolutor.cpp:85
gslpp::matrix< double > * Evol(orders order)
Evolution matrix set at a fixed order of QCD coupling.
Definition: RGEvolutor.cpp:103
A model class for the Standard Model.
const double getMz() const
A get method to access the mass of the boson .
const double getAlsMz() const
A get method to access the value of .
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.
const double getAle() const
A get method to retrieve the fine-structure constant .
gslpp::matrix< double > * elem[MAXORDER_QED+1]
Test Observable.
Test Observable.
orders
An enum type for orders in QCD.
Definition: OrderScheme.h:33
schemes
An enum type for regularization schemes.
Definition: OrderScheme.h:20
orders_qed
An enum type for orders in electroweak.
Definition: OrderScheme.h:56
StdVectorFiller< int > Vi