a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
Loading...
Searching...
No Matches
EvolBsmm.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015 HEPfit Collaboration
3 *
4 *
5 * For the licensing terms see doc/COPYING.
6 */
7
8#include "EvolBsmm.h"
9#include "StandardModel.h"
10#include <gsl/gsl_sf.h>
11
12EvolBsmm::EvolBsmm(unsigned int dim_i, schemes scheme, orders order, orders_qed order_qed, const StandardModel & model)
13: RGEvolutor(dim_i, scheme, order, order_qed), model(model), V(dim_i,0.), Vi(dim_i,0.),
14 AA(dim_i,0.), BB(dim_i,0.), CC(dim_i,0.), DD(dim_i,0.), EE(dim_i,0.), FF(dim_i,0.),
15 RR(dim_i,0.), e(dim_i,0.), vavi(0,0.), vbvi(0,0.), vcvi(0,0.), vdvi(0,0.),
16 vevi(0,0.), vfvi(0,0.), vrvi(0,0.),vaevi(0,0.), vbbvi(0,0.), vbdvi(0,0.), vbevi(0,0.),
17 vdbvi(0,0.), vdevi(0,0.), veavi(0,0.), vebvi(0,0.), vedvi(0,0.), veevi(0,0.), vbeevi(0,0.),
18 vebevi(0,0.), veebvi(0,0.), vbbevi(0,0.), vbebvi(0,0.), vebbvi(0,0.), dim(dim_i)
19{
20 unsigned int i = 0;
21 unsigned int j = 0;
22 unsigned int l = 0;
23 unsigned int m = 0;
24 unsigned int p = 0;
25 unsigned int q = 0;
26 int L = 1;
27 double b0 = 0.;
28
29
30 vavi.clear();
31 vbvi.clear();
32 vcvi.clear();
33 vdvi.clear();
34 vfvi.clear();
35 vevi.clear();
36 vrvi.clear();
37 vaevi.clear();
38 vbbvi.clear();
39 vbdvi.clear();
40 vbevi.clear();
41 vdbvi.clear();
42 vdevi.clear();
43 veavi.clear();
44 vebvi.clear();
45 vedvi.clear();
46 veevi.clear();
47 vbeevi.clear();
48 vebevi.clear();
49 veebvi.clear();
50 vbbevi.clear();
51 vbebvi.clear();
52 vebbvi.clear();
53
54
55
56/* define L, nu, nd */
57 if(L == 1){nd = 3; nu = 2;}
58 b0 = model.Beta0(6-L);
59
60 AnomalousDimension(10,nu,nd).transpose().eigensystem(V,e);
61
62 Vi = V.inverse();
63
64 /* magic numbers of U0 */
65 for(unsigned int i = 0; i < dim; i++){
66 a[L][i] = e(i).real()/2./b0;
67 for (unsigned int j = 0; j < dim; j++){
68 for (unsigned int k = 0; k < dim; k++){
69 b[L][i][j][k] = V(i, k).real() * Vi(k, j).real();
70 }
71 }
72 }
73
74
75
76 AA = BuiltB('A', nu, nd);
77 BB = BuiltB('B', nu, nd);
78 CC = BuiltB('C', nu, nd);
79 DD = BuiltB('D', nu, nd);
80 EE = BuiltB('E', nu, nd);
81 FF = BuiltB('F', nu, nd);
82 RR = BuiltB('R', nu, nd);
83
84 double cutoff = 0.000000001 ;
85
86 for(l = 0; l < dim; l++){
87 for(i = 0; i < dim; i++){
88 for(j = 0; j < dim; j++){
89 for(m = 0; m < dim; m++){
90
91 if(fabs(V(l, i).real() * AA(i, j).real() * Vi(j, m).real()) > cutoff){
92
93 vavi.push_back(l);
94 vavi.push_back(i);
95 vavi.push_back(j);
96 vavi.push_back(m);
97 vavi.push_back(V(l, i).real() * AA(i, j).real() * Vi(j, m).real());
98
99 }
100 if(fabs(V(l, i).real() * BB(i, j).real() * Vi(j, m).real()) > cutoff){
101
102 vbvi.push_back(l);
103 vbvi.push_back(i);
104 vbvi.push_back(j);
105 vbvi.push_back(m);
106 vbvi.push_back(V(l, i).real() * BB(i, j).real() * Vi(j, m).real());
107
108 }
109 if(fabs(V(l, i).real() * CC(i, j).real() * Vi(j, m).real()) > cutoff){
110
111 vcvi.push_back(l);
112 vcvi.push_back(i);
113 vcvi.push_back(j);
114 vcvi.push_back(m);
115 vcvi.push_back(V(l, i).real() * CC(i, j).real() * Vi(j, m).real());
116
117 }
118 if(fabs(V(l, i).real() * DD(i, j).real() * Vi(j, m).real()) > cutoff){
119
120 vdvi.push_back(l);
121 vdvi.push_back(i);
122 vdvi.push_back(j);
123 vdvi.push_back(m);
124 vdvi.push_back(V(l, i).real() * DD(i, j).real() * Vi(j, m).real());
125
126 }
127 if(fabs(V(l, i).real() * EE(i, j).real() * Vi(j, m).real()) > cutoff){
128
129 vevi.push_back(l);
130 vevi.push_back(i);
131 vevi.push_back(j);
132 vevi.push_back(m);
133 vevi.push_back(V(l, i).real() * EE(i, j).real() * Vi(j, m).real());
134
135 }
136 if(fabs(V(l, i).real() * FF(i, j).real() * Vi(j, m).real()) > cutoff){
137
138 vfvi.push_back(l);
139 vfvi.push_back(i);
140 vfvi.push_back(j);
141 vfvi.push_back(m);
142 vfvi.push_back(V(l, i).real() * FF(i, j).real() * Vi(j, m).real());
143
144 }
145 if(fabs(V(l, i).real() * RR(i, j).real() * Vi(j, m).real()) > cutoff){
146
147 vrvi.push_back(l);
148 vrvi.push_back(i);
149 vrvi.push_back(j);
150 vrvi.push_back(m);
151 vrvi.push_back(V(l, i).real() * RR(i, j).real() * Vi(j, m).real());
152
153 }
154
155 for(p = 0; p < dim; p++){
156
157 if(fabs(V(l, i).real() * AA(i, p).real() * EE(p, j).real() * Vi(j, m).real()) > cutoff){
158
159 vaevi.push_back(l);
160 vaevi.push_back(i);
161 vaevi.push_back(p);
162 vaevi.push_back(j);
163 vaevi.push_back(m);
164 vaevi.push_back(V(l, i).real() * AA(i, p).real() * EE(p, j).real() * Vi(j, m).real());
165
166 }
167 if(fabs(V(l, i).real() * BB(i, p).real() * BB(p, j).real() * Vi(j, m).real()) > cutoff){
168
169 vbbvi.push_back(l);
170 vbbvi.push_back(i);
171 vbbvi.push_back(p);
172 vbbvi.push_back(j);
173 vbbvi.push_back(m);
174 vbbvi.push_back(V(l, i).real() * BB(i, p).real() * BB(p, j).real() * Vi(j, m).real());
175
176 }
177 if(fabs(V(l, i).real() * BB(i, p).real() * DD(p, j).real() * Vi(j, m).real()) > cutoff){
178
179 vbdvi.push_back(l);
180 vbdvi.push_back(i);
181 vbdvi.push_back(p);
182 vbdvi.push_back(j);
183 vbdvi.push_back(m);
184 vbdvi.push_back(V(l, i).real() * BB(i, p).real() * DD(p, j).real() * Vi(j, m).real());
185
186 }
187 if(fabs(V(l, i).real() * BB(i, p).real() * EE(p, j).real() * Vi(j, m).real()) > cutoff){
188
189 vbevi.push_back(l);
190 vbevi.push_back(i);
191 vbevi.push_back(p);
192 vbevi.push_back(j);
193 vbevi.push_back(m);
194 vbevi.push_back(V(l, i).real() * BB(i, p).real() * EE(p, j).real() * Vi(j, m).real());
195
196 }
197 if(fabs(V(l, i).real() * DD(i, p).real() * BB(p, j).real() * Vi(j, m).real()) > cutoff){
198
199 vdbvi.push_back(l);
200 vdbvi.push_back(i);
201 vdbvi.push_back(p);
202 vdbvi.push_back(j);
203 vdbvi.push_back(m);
204 vdbvi.push_back(V(l, i).real() * DD(i, p).real() * BB(p, j).real() * Vi(j, m).real());
205
206 }
207 if(fabs(V(l, i).real() * DD(i, p).real() * EE(p, j).real() * Vi(j, m).real()) > cutoff){
208
209 vdevi.push_back(l);
210 vdevi.push_back(i);
211 vdevi.push_back(p);
212 vdevi.push_back(j);
213 vdevi.push_back(m);
214 vdevi.push_back(V(l, i).real() * DD(i, p).real() * EE(p, j).real() * Vi(j, m).real());
215
216 }
217 if(fabs(V(l, i).real() * EE(i, p).real() * AA(p, j).real() * Vi(j, m).real()) > cutoff){
218
219 veavi.push_back(l);
220 veavi.push_back(i);
221 veavi.push_back(p);
222 veavi.push_back(j);
223 veavi.push_back(m);
224 veavi.push_back(V(l, i).real() * EE(i, p).real() * AA(p, j).real() * Vi(j, m).real());
225
226 }
227 if(fabs(V(l, i).real() * EE(i, p).real() * BB(p, j).real() * Vi(j, m).real()) > cutoff){
228
229 vebvi.push_back(l);
230 vebvi.push_back(i);
231 vebvi.push_back(p);
232 vebvi.push_back(j);
233 vebvi.push_back(m);
234 vebvi.push_back(V(l, i).real() * EE(i, p).real() * BB(p, j).real() * Vi(j, m).real());
235
236 }
237 if(fabs(V(l, i).real() * EE(i, p).real() * DD(p, j).real() * Vi(j, m).real()) > cutoff){
238
239 vedvi.push_back(l);
240 vedvi.push_back(i);
241 vedvi.push_back(p);
242 vedvi.push_back(j);
243 vedvi.push_back(m);
244 vedvi.push_back(V(l, i).real() * EE(i, p).real() * DD(p, j).real() * Vi(j, m).real());
245
246 }
247 if(fabs(V(l, i).real() * EE(i, p).real() * EE(p, j).real() * Vi(j, m).real()) > cutoff){
248
249 veevi.push_back(l);
250 veevi.push_back(i);
251 veevi.push_back(p);
252 veevi.push_back(j);
253 veevi.push_back(m);
254 veevi.push_back(V(l, i).real() * EE(i, p).real() * EE(p, j).real() * Vi(j, m).real());
255
256 }
257
258 for(q = 0; q < dim; q++){
259
260 if(fabs(V(l, i).real() * BB(i, p).real() * EE(p, q).real() * EE(q, j).real() * Vi(j, m).real()) > cutoff){
261
262 vbeevi.push_back(l);
263 vbeevi.push_back(i);
264 vbeevi.push_back(p);
265 vbeevi.push_back(q);
266 vbeevi.push_back(j);
267 vbeevi.push_back(m);
268 vbeevi.push_back(V(l, i).real() * BB(i, p).real() * EE(p, q).real() * EE(q, j).real() * Vi(j, m).real());
269
270 }
271
272 if (fabs(V(l, i).real() * EE(i, p).real() * BB(p, q).real() * EE(q, j).real() * Vi(j, m).real()) > cutoff) {
273
274 vebevi.push_back(l);
275 vebevi.push_back(i);
276 vebevi.push_back(p);
277 vebevi.push_back(q);
278 vebevi.push_back(j);
279 vebevi.push_back(m);
280 vebevi.push_back(V(l, i).real() * EE(i, p).real() * BB(p, q).real() * EE(q, j).real() * Vi(j, m).real());
281
282 }
283 if (fabs(V(l, i).real() * EE(i, p).real() * EE(p, q).real() * BB(q, j).real() * Vi(j, m).real()) > cutoff) {
284
285 veebvi.push_back(l);
286 veebvi.push_back(i);
287 veebvi.push_back(p);
288 veebvi.push_back(q);
289 veebvi.push_back(j);
290 veebvi.push_back(m);
291 veebvi.push_back(V(l, i).real() * EE(i, p).real() * EE(p, q).real() * BB(q, j).real() * Vi(j, m).real());
292
293 }
294 if (fabs(V(l, i).real() * BB(i, p).real() * BB(p, q).real() * EE(q, j).real() * Vi(j, m).real()) > cutoff) {
295
296 vbbevi.push_back(l);
297 vbbevi.push_back(i);
298 vbbevi.push_back(p);
299 vbbevi.push_back(q);
300 vbbevi.push_back(j);
301 vbbevi.push_back(m);
302 vbbevi.push_back(V(l, i).real() * BB(i, p).real() * BB(p, q).real() * EE(q, j).real() * Vi(j, m).real());
303
304 }
305 if (fabs(V(l, i).real() * BB(i, p).real() * EE(p, q).real() * BB(q, j).real() * Vi(j, m).real()) > cutoff) {
306
307 vbebvi.push_back(l);
308 vbebvi.push_back(i);
309 vbebvi.push_back(p);
310 vbebvi.push_back(q);
311 vbebvi.push_back(j);
312 vbebvi.push_back(m);
313 vbebvi.push_back(V(l, i).real() * BB(i, p).real() * EE(p, q).real() * BB(q, j).real() * Vi(j, m).real());
314
315 }
316 if (fabs(V(l, i).real() * EE(i, p).real() * BB(p, q).real() * BB(q, j).real() * Vi(j, m).real()) > cutoff) {
317
318 vebbvi.push_back(l);
319 vebbvi.push_back(i);
320 vebbvi.push_back(p);
321 vebbvi.push_back(q);
322 vebbvi.push_back(j);
323 vebbvi.push_back(m);
324 vebbvi.push_back(V(l, i).real() * EE(i, p).real() * BB(p, q).real() * BB(q, j).real() * Vi(j, m).real());
325
326 }
327 }
328 }
329 }
330 }
331 }
332 }
333
334}
335
336
339
340gslpp::matrix<double> EvolBsmm::AnomalousDimension(int gam, unsigned int n_u, unsigned int n_d) const
341
342{
343
344 /* Delta F = 1 anomalous dimension in Misiak basis,
345 ref.: ref. hep-ph/1311.1348v2, hep-ph/0512066 */
346
347 /* gamma(row, coloumn) at the LO */
348
349 unsigned int nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
350 double zeta3 = 0;
351
352 gslpp::matrix<double> gammaDF1(dim, 0.);
353
354 zeta3 = gsl_sf_zeta_int(3);
355
356 switch(gam){
357
358 case 10:
359 if (!(nf == 3 || nf == 4 || nf == 5 || nf == 6)){
360 throw std::runtime_error("EvolBsmm::AnomalousDimension(): wrong number of flavours");
361 }
362
363 gammaDF1(0,0) = -4. ;
364 gammaDF1(0,1) = 8./3. ;
365 gammaDF1(0,3) = -2./9.;
366
367 gammaDF1(1,0) = 12.;
368 gammaDF1(1,3) = 4./3.;
369
370 gammaDF1(2,3) = -52./3.;
371 gammaDF1(2,5) = 2.;
372
373 gammaDF1(3,2) = -40./9.;
374 gammaDF1(3,3) = -100./9.;
375 gammaDF1(3,4) = 4./9.;
376 gammaDF1(3,5) = 5./6.;
377
378 gammaDF1(4,3) = -256./3.;
379 gammaDF1(4,5) = 20.;
380
381 gammaDF1(5,2) = -256./9.;
382 gammaDF1(5,3) = 56./9.;
383 gammaDF1(5,4) = 40./9.;
384 gammaDF1(5,5) = -2./3.;
385
386 break;
387 case 20:
388
389 if (!(nf == 3 || nf == 4 || nf == 5 || nf == 6)){
390 throw std::runtime_error("EvolBsmm::AnomalousDimension(): wrong number of flavours");
391 }
392
393 /* gamma(row, coloumn) at the NLO */
394
395 gammaDF1(0,0) = -355./9.;
396 gammaDF1(0,1) = -502./27.;
397 gammaDF1(0,2) = -1412./243.;
398 gammaDF1(0,3) = -1369./243.;
399 gammaDF1(0,4) = 134./243.;
400 gammaDF1(0,5) = -35./162.;
401
402 gammaDF1(1,0) = -35./3.;
403 gammaDF1(1,1) = -28./3.;
404 gammaDF1(1,2) = -416./81.;
405 gammaDF1(1,3) = 1280./81.;
406 gammaDF1(1,4) = 56./81.;
407 gammaDF1(1,5) = 35./27.;
408
409 gammaDF1(2,2) = -4468./81.;
410 gammaDF1(2,3) = -31469./81.;
411 gammaDF1(2,4) = 400./81.;
412 gammaDF1(2,5) = 3373./108.;
413
414 gammaDF1(3,2) = -8158./243.;
415 gammaDF1(3,3) = -59399./243.;
416 gammaDF1(3,4) = 269./486.;
417 gammaDF1(3,5) = 12899./648.;
418
419 gammaDF1(4,2) = -251680./81.;
420 gammaDF1(4,3) = -128648./81.;
421 gammaDF1(4,4) = 23836./81.;
422 gammaDF1(4,5) = 6106./27.;
423
424 gammaDF1(5,2) = 58640./243.;
425 gammaDF1(5,3) = -26348./243.;
426 gammaDF1(5,4) = -14324./243.;
427 gammaDF1(5,5) = -2551./162.;
428
429 break;
430
431 case 30:
432
433 if (!(nf == 3 || nf == 4 || nf == 5 || nf == 6)){
434 throw std::runtime_error("EvolBsmm::AnomalousDimension(): wrong number of flavours");
435 }
436
437 /* gamma(row, coloumn) at the NNLO */
438
439 gammaDF1(0,0) = -12773./18. + zeta3 * 1472./3.;
440 gammaDF1(0,1) = 745./9. - zeta3 *4288./9.;
441 gammaDF1(0,2) = 63187./13122. - zeta3 * 1360./81.;
442 gammaDF1(0,3) = -981796./6561. - zeta3 * 776./81.;
443 gammaDF1(0,4) = -202663./52488. + zeta3 * 124./81.;
444 gammaDF1(0,5) = -24973./69984. + zeta3 * 100./27.;
445
446 gammaDF1(1,0) = 1177./2. - zeta3 * 2144.;
447 gammaDF1(1,1) = 306. - zeta3 * 224.;
448 gammaDF1(1,2) = 110477./2187. + zeta3 * 2720./27.;
449 gammaDF1(1,3) = 133529./8748. - zeta3 * 2768./27.;
450 gammaDF1(1,4) = -42929./8748. - zeta3 * 248./27.;
451 gammaDF1(1,5) = 354319./11664. - zeta3 * 110./9.;
452
453 gammaDF1(2,2) = -3572528./2187. - zeta3 * 608./27.;
454 gammaDF1(2,3) = -58158773./8748. + zeta3 * 61424./27.;
455 gammaDF1(2,4) = 552601./4374. - zeta3 * 496./27.;
456 gammaDF1(2,5) = 6989171./11664. - zeta3 * 2821./9.;
457
458 gammaDF1(3,2) = -1651004./6561. + zeta3 * 88720./81.;
459 gammaDF1(3,3) = -155405353./52488 + zeta3 * 54272./81.;
460 gammaDF1(3,4) = 1174159./52488. - zeta3 * 9274./81.;
461 gammaDF1(3,5) = 10278809./34992. - zeta3 * 3100./27.;
462
463 gammaDF1(4,2) = -147978032./2187. + zeta3 * 87040./27.;
464 gammaDF1(4,3) = -168491372./2187. + zeta3 * 324416./27.;
465 gammaDF1(4,4) = 11213042./2187. - zeta3 * 13984./27.;
466 gammaDF1(4,5) = 17850329./2916. - zeta3 * 31420./9.;
467
468 gammaDF1(5,2) = 136797922./6561. + zeta3 * 721408./81.;
469 gammaDF1(5,3) = -72614473./13122. - zeta3 * 166432./81.;
470 gammaDF1(5,4) = -9288181./6561. - zeta3 * 95032./81.;
471 gammaDF1(5,5) = -16664027./17496. - zeta3 * 7552./27.;
472
473 break;
474 case 01:
475
476 if (!(nf == 3 || nf == 4 || nf == 5 || nf == 6)){
477 throw std::runtime_error("EvolBsmm::AnomalousDimension(): wrong number of flavours");
478 }
479
480 gammaDF1(0,0) = -8./3.;
481 gammaDF1(0,6) = -32./27.;
482
483 gammaDF1(1,1) = -8./3.;
484 gammaDF1(1,6) = -8./9.;
485
486 gammaDF1(2,6) = -16./9.;
487
488 gammaDF1(3,6) = 32./27.;
489
490 gammaDF1(4,6) = -112./9.;
491
492 gammaDF1(5,6) = 512./27.;
493
494 gammaDF1(6,6) = 8.;
495 gammaDF1(6,7) = -4.;
496
497 gammaDF1(7,6) = -4.;
498
499 break;
500 case 11:
501
502 if (!(nf == 3 || nf == 4 || nf == 5 || nf == 6)){
503 throw std::runtime_error("EvolBsmm::AnomalousDimension(): wrong number of flavours");
504 }
505
506
507 gammaDF1(0,0) = 169./9.;
508 gammaDF1(0,1) = 100./27.;
509 gammaDF1(0,3) = 254./729.;
510 gammaDF1(0,6) = -2272./729.;
511
512 gammaDF1(1,0) = 50./3.;
513 gammaDF1(1,1) = -8./3.;
514 gammaDF1(1,3) = 1076./243.;
515 gammaDF1(1,6) = 1952./243.;
516
517 gammaDF1(2,3) = 11116./243.;
518 gammaDF1(2,5) = -14./3.;
519 gammaDF1(2,6) = -6752./243.;
520
521 gammaDF1(3,2) = 280./27.;
522 gammaDF1(3,3) = 18763./729.;
523 gammaDF1(3,4) = -28./27.;
524 gammaDF1(3,5) = -35./18.;
525 gammaDF1(3,6) = -2192./729.;
526
527 gammaDF1(4,3) = 111136./243.;
528 gammaDF1(4,5) = -140./3.;
529 gammaDF1(4,6) = -84032./243.;
530
531 gammaDF1(5,2) = 2944./27.;
532 gammaDF1(5,3) = 193312./729.;
533 gammaDF1(5,4) = -280./27.;
534 gammaDF1(5,5) = -175./9.;
535 gammaDF1(5,6) = -37856./729.;
536
537 gammaDF1(6,7) = 16.;
538
539 gammaDF1(7,6) = 16.;
540
541 break;
542 case 02:
543
544 if (!(nf == 3 || nf == 4 || nf == 5 || nf == 6)){
545 throw std::runtime_error("EvolBsmm::AnomalousDimension(): wrong number of flavours");
546 }
547
548 gammaDF1(0,6)= -11680./2187.;
549 gammaDF1(0,7)= -416./81.;
550
551 gammaDF1(1,6)= -2920./729.;
552 gammaDF1(1,7)= -104./27.;
553
554 gammaDF1(2,6)= -39752./729.;
555 gammaDF1(2,7)= -136./27.;
556
557 gammaDF1(3,6)= 1024./2187.;
558 gammaDF1(3,7)= -448./81.;
559
560 gammaDF1(4,6)= -381344./729.;
561 gammaDF1(4,7)= -15616./27.;
562
563 gammaDF1(5,6)= 24832./2187.;
564 gammaDF1(5,7)= -7936./81.;
565
566
567 break;
568 case 21:
569
570 if (!(nf == 3 || nf == 4 || nf == 5 || nf == 6)){
571 throw std::runtime_error("EvolBsmm::AnomalousDimension(): wrong number of flavours");
572 }
573
574 gammaDF1(0,6)= -1359190./19683. + zeta3 * 6976./243.;
575
576 gammaDF1(1,6)= -229696./6561 - zeta3 * 3584./81.;
577
578 gammaDF1(2,6)= -1290092./6561 + zeta3 * 3200./81.;
579
580 gammaDF1(3,6)= -819971./19683. - zeta3 * 19936./243;
581
582 gammaDF1(4,6)= -16821944./6561 + zeta3 * 30464./81.;
583
584 gammaDF1(5,6)= -17787368./19683. - zeta3 * 286720./243.;
585
586
587 break;
588
589 default:
590 throw std::runtime_error("EvolBsmm::AnomalousDimension(): order not implemented");
591 }
592 return (gammaDF1);
593}
594
595
596gslpp::matrix<double> EvolBsmm::BuiltB(char letter, unsigned int n_u, unsigned int n_d)
597
598{
599
600 unsigned int nf = 5; //all the class works for nf = 5
601
602 double B00S = model.Beta0(nf), B10S = model.Beta1(nf), B20S = model.Beta2(nf),
603 B01S = -22./9., B11S = -308./27.;
604
605 double B00E = 80./9., B01E = 176./9.;
606
607 double b1 = B10S/(2. * B00S * B00S), b2 = B20S/(4. * B00S * B00S * B00S) - b1 * b1 ,
608 b3 = B01S/(2. * B00S * B00E), b4 = B11S /(4. * B00S * B00S * B00E) - 2 * b1 * b3,
609 b5 = B01E/(2. * B00S * B00E) - b1;
610
611
612
613 gslpp::matrix<double> B(dim, 0.);
614 gslpp::matrix<double> W10T(dim, 0.);
615 gslpp::matrix<double> W20T(dim, 0.);
616 gslpp::matrix<double> W30T(dim, 0.);
617 gslpp::matrix<double> W01T(dim, 0.);
618 gslpp::matrix<double> W02T(dim, 0.);
619 gslpp::matrix<double> W11T(dim, 0.);
620 gslpp::matrix<double> W21T(dim, 0.);
621
622
623 W10T = (AnomalousDimension(10, n_u, n_d).transpose())/2./B00S;
624 W20T = (AnomalousDimension(20, n_u, n_d).transpose())/4./B00S/B00S;
625 W30T = (AnomalousDimension(30, n_u, n_d).transpose())/8./B00S/B00S/B00S;
626 W01T = (AnomalousDimension(01, n_u, n_d).transpose())/2./B00E;
627 W02T = (AnomalousDimension(02, n_u, n_d).transpose())/4./B00E/B00E;
628 W11T = (AnomalousDimension(11, n_u, n_d).transpose())/4./B00S/B00E;
629 W21T = (AnomalousDimension(21, n_u, n_d).transpose())/8./B00S/B00S/B00E;
630
631
632
633 switch(letter){
634
635 case 'A':
636
637 B = Vi.real() * (W30T - b1 * W20T - b2 * W10T) * V.real();
638
639 break;
640 case 'B':
641
642 B = Vi.real() * (W20T - b1 * W10T) * V.real();
643
644 break;
645 case 'C':
646
647 B = Vi.real() * (W21T - b1 * W11T - b2 * W01T - b3 * W20T - b4 * W10T) * V.real();
648
649 break;
650 case 'D':
651
652 B = Vi.real() * (W11T - b1 * W01T - b3 * W10T) * V.real();
653
654 break;
655 case 'E':
656
657 B = Vi.real() * (W01T) * V.real();
658
659 break;
660 case 'F':
661
662 B = Vi.real() * (W02T + W11T - (b1 + b3) * W01T - b3 * W10T) * V.real();
663
664 break;
665 case 'R':
666
667 B = b5 * Vi.real() * (W01T) * V.real();
668
669 break;
670 default:
671 throw std::runtime_error("EvolBsmm::BuiltB(): order not implemented");
672 }
673 return (B);
674}
675
676
677
678gslpp::matrix< double > & EvolBsmm::Df1Evol(double mu, double M, orders order, orders_qed order_qed, schemes scheme)
679
680{
681 switch (scheme) { /* complete this method */
682 case NDR:
683 break;
684 case LRI:
685 case HV:
686 default:
687 std::stringstream out;
688 out << scheme;
689 throw std::runtime_error("EvolDF1nlep::Df1Evol(): scheme " + out.str() + " not implemented ");
690 }
691
692 if (mu == this->mu && M == this->M && scheme == this->scheme && order_qed == NO_QED)
693 return (*Evol(order));
694
695 if (mu == this->mu && M == this->M && scheme == this->scheme && order_qed != NO_QED)
696 return (*Evol(order_qed));
697
698
699 if (M < mu) {
700 std::stringstream out;
701 out << "M = " << M << " < mu = " << mu;
702 throw out.str();
703 }
704
705 setScales(mu, M); // also assign evol to identity
706 if (M != mu) {
707 double m_down = mu;
708 //double m_up = model.AboveTh(m_down);
709 //double nf = model.Nf(m_down);
710 double nf = 5; //all the process in implemented for nf = 5
711 alsM = alphatilde_s(M);
712 alsmu = alphatilde_s(mu);
713
714 eta = alsM / alsmu;
715 logeta = log(eta);
716
717 /*while (m_up < M) { //there is no thresholds
718 Df1Evol(m_down, m_up, nf, scheme);
719 //Df1threshold_nlep(m_up, nf+1.);
720 m_down = m_up;
721 m_up = model.AboveTh(m_down);
722 nf += 1.;
723 } */
724
725 Df1Evol(m_down, M, nf, scheme);
726 }
727
728 if(order_qed != NO_QED) return (*Evol(order_qed));
729 else return (*Evol(order));
730
731}
732
733void EvolBsmm::Df1Evol(double mu, double M, double nf, schemes scheme)
734
735{
736 unsigned int i = 0;
737 unsigned int j = 0;
738 unsigned int k = 0;
739
740 gslpp::matrix<double> resLO(dim, 0.);
741 gslpp::matrix<double> Ueos(dim, 0.), Ue(dim, 0.), Ues(dim,0.), Us(dim,0.), Ue2os(dim, 0.), Ue2os2(dim, 0.), Ue2(dim,0.), Us2(dim,0.);
742
743 int L = 6 - (int) nf;
744
745 double eta = alsM / alsmu;
746
747 double B00S = model.Beta0(nf);// B10S = model.Beta1(nf); /*inizializza i B*/
748
749 double B00E = 80./9.;// B01E = 176./9.; /*inizializza i B*/
750
751 double omega = 2. * B00S , lambda = B00E /B00S ; /* E ale dipende da mu?*/
752
753
754 for (k = 0; k < dim; k++) {
755 double etap = pow(eta, a[L][k]);
756 for (i = 0; i < dim; i++) {
757 for (j = 0; j < dim; j++) {
758 resLO(i, j) += b[L][i][j][k] * etap;
759 Ue2(i, j) = (i == j);
760 }
761 }
762 }
763
764 unsigned int ind = 0;
765 //double max = 0;
766
767 gslpp::vector<double> list(23, 0.);
768
769 list(0) = vbeevi.size()/7.;
770 list(1) = vebevi.size()/7.;
771 list(2) = veebvi.size()/7.;
772 list(3) = vbbevi.size()/7.;
773 list(4) = vbebvi.size()/7.;
774 list(5) = vebbvi.size()/7.;
775 list(6) = vaevi.size()/6.;
776 list(7) = vbbvi.size()/6.;
777 list(8) = vbdvi.size()/6.;
778 list(9) = vbevi.size()/6.;
779 list(10) = vdbvi.size()/6.;
780 list(11) = vdevi.size()/6.;
781 list(12) = veavi.size()/6.;
782 list(13) = vebvi.size()/6.;
783 list(14) = vedvi.size()/6.;
784 list(15) = veevi.size()/6.;
785 list(16) = vavi.size()/5.;
786 list(17) = vbvi.size()/5.;
787 list(18) = vcvi.size()/5.;
788 list(19) = vdvi.size()/5.;
789 list(20) = vevi.size()/5.;
790 list(21) = vfvi.size()/5.;
791 list(22) = vrvi.size()/5.;
792
793 double max = list.max();
794
795 for (ind = 0; ind < max; ind++) {
796
797 if (ind < list(0)) {
798
799 Ue2os(vbeevi[7 * ind], vbeevi[7 * ind + 5]) += vbeevi[7 * ind + 6] * H(vbeevi[7 * ind + 1], vbeevi[7 * ind + 2], vbeevi[7 * ind + 3], vbeevi[7 * ind + 4], 2, 4, 4, mu, M, nf);
800 }
801 if (ind < list(1)) {
802
803 Ue2os(vebevi[7 * ind], vebevi[7 * ind + 5]) += vebevi[7 * ind + 6] * H(vebevi[7 * ind + 1], vebevi[7 * ind + 2], vebevi[7 * ind + 3], vebevi[7 * ind + 4], 4, 2, 4, mu, M, nf);
804 }
805 if (ind < list(2)) {
806
807 Ue2os(veebvi[7 * ind], veebvi[7 * ind + 5]) += veebvi[7 * ind + 6] * H(veebvi[7 * ind + 1], veebvi[7 * ind + 2], veebvi[7 * ind + 3], veebvi[7 * ind + 4], 4, 4, 2, mu, M, nf);
808 }
809 if (ind < list(3)) {
810
811 Ues(vbbevi[7 * ind], vbbevi[7 * ind + 5]) += vbbevi[7 * ind + 6] * H(vbbevi[7 * ind + 1], vbbevi[7 * ind + 2], vbbevi[7 * ind + 3], vbbevi[7 * ind + 4], 2, 2, 4, mu, M, nf);
812 }
813 if (ind < list(4)) {
814
815 Ues(vbebvi[7 * ind], vbebvi[7 * ind + 5]) += vbebvi[7 * ind + 6] * H(vbebvi[7 * ind + 1], vbebvi[7 * ind + 2], vbebvi[7 * ind + 3], vbebvi[7 * ind + 4], 2, 4, 2, mu, M, nf);
816 }
817 if (ind < list(5)) {
818
819 Ues(vebbvi[7 * ind], vebbvi[7 * ind + 5]) += vebbvi[7 * ind + 6] * H(vebbvi[7 * ind + 1], vebbvi[7 * ind + 2], vebbvi[7 * ind + 3], vebbvi[7 * ind + 4], 4, 2, 2, mu, M, nf);
820 }
821
822
823 if (ind < list(6)) {
824
825 Ues(vaevi[6 * ind], vaevi[6 * ind + 4]) += vaevi[6 * ind + 5] * G(vaevi[6 * ind + 1], vaevi[6 * ind + 2], vaevi[6 * ind + 3], 1, 4, mu, M, nf);
826 }
827 if (ind < list(7)) {
828
829 Us2(vbbvi[6 * ind], vbbvi[6 * ind + 4]) += vbbvi[6 * ind + 5] * G(vbbvi[6 * ind + 1], vbbvi[6 * ind + 2], vbbvi[6 * ind + 3], 2, 2, mu, M, nf);
830 }
831 if (ind < list(8)) {
832
833 Ues(vbdvi[6 * ind], vbdvi[6 * ind + 4]) += vbdvi[6 * ind + 5] * G(vbdvi[6 * ind + 1], vbdvi[6 * ind + 2], vbdvi[6 * ind + 3], 2, 3, mu, M, nf);
834 }
835 if (ind < list(9)) {
836
837 Ue(vbevi[6 * ind], vbevi[6 * ind + 4]) += vbevi[6 * ind + 5] * G(vbevi[6 * ind + 1], vbevi[6 * ind + 2], vbevi[6 * ind + 3], 2, 4, mu, M, nf);
838 Ue2os(vbevi[6 * ind], vbevi[6 * ind + 4]) += vbevi[6 * ind + 5] * (-G(vbevi[6 * ind + 1], vbevi[6 * ind + 2], vbevi[6 * ind + 3], 2, 4, mu, M, nf)
839 + G(vbevi[6 * ind + 1], vbevi[6 * ind + 2], vbevi[6 * ind + 3], 2, 5, mu, M, nf));
840 }
841 if (ind < list(10)) {
842
843 Ues(vdbvi[6 * ind], vdbvi[6 * ind + 4]) += vdbvi[6 * ind + 5] * G(vdbvi[6 * ind + 1], vdbvi[6 * ind + 2], vdbvi[6 * ind + 3], 3, 2, mu, M, nf);
844 }
845 if (ind < list(11)) {
846
847 Ue2os(vdevi[6 * ind], vdevi[6 * ind + 4]) += vdevi[6 * ind + 5] * G(vdevi[6 * ind + 1], vdevi[6 * ind + 2], vdevi[6 * ind + 3], 3, 4, mu, M, nf);
848 }
849 if (ind < list(12)) {
850
851 Ues(veavi[6 * ind], veavi[6 * ind + 4]) += veavi[6 * ind + 5] * G(veavi[6 * ind + 1], veavi[6 * ind + 2], veavi[6 * ind + 3], 4, 1, mu, M, nf);
852
853 }
854 if (ind < list(13)) {
855
856 Ue(vebvi[6 * ind], vebvi[6 * ind + 4]) += vebvi[6 * ind + 5] * G(vebvi[6 * ind + 1], vebvi[6 * ind + 2], vebvi[6 * ind + 3], 4, 2, mu, M, nf);
857 Ue2os(vebvi[6 * ind], vebvi[6 * ind + 4]) += vebvi[6 * ind + 5] * (-G(vebvi[6 * ind + 1], vebvi[6 * ind + 2], vebvi[6 * ind + 3], 4, 2, mu, M, nf)
858 + G(vebvi[6 * ind + 1], vebvi[6 * ind + 2], vebvi[6 * ind + 3], 5, 2, mu, M, nf));
859 }
860 if (ind < list(14)) {
861
862 Ue2os(vedvi[6 * ind], vedvi[6 * ind + 4]) += vedvi[6 * ind + 5] * G(vedvi[6 * ind + 1], vedvi[6 * ind + 2], vedvi[6 * ind + 3], 4, 3, mu, M, nf);
863
864 }
865 if (ind < list(15)) {
866
867 Ue2os2(veevi[6 * ind], veevi[6 * ind + 4]) += (veevi[6 * ind + 5] * G(veevi[6 * ind + 1], veevi[6 * ind + 2], veevi[6 * ind + 3], 4, 4, mu, M, nf));
868
869 }
870
871
872 if (ind < list(16)) {
873
874 Us2(vavi[5 * ind], vavi[5 * ind + 3]) += (vavi[5 * ind + 4] * F(vavi[5 * ind + 1], vavi[5 * ind + 2], 1, mu, M, nf));
875 }
876
877
878 if (ind < list(17)) {
879
880 Us(vbvi[5 * ind], vbvi[5 * ind + 3]) += vbvi[5 * ind + 4] * F(vbvi[5 * ind + 1], vbvi[5 * ind + 2], 2, mu, M, nf);
881
882 }
883
884 if (ind < list(18)) {
885
886 Ues(vcvi[5 * ind], vcvi[5 * ind + 3]) += vcvi[5 * ind + 4] * F(vcvi[5 * ind + 1], vcvi[5 * ind + 2], 2, mu, M, nf);
887 }
888 if (ind < list(19)) {
889
890 Ue(vdvi[5 * ind], vdvi[5 * ind + 3]) += vdvi[5 * ind + 4] * F(vdvi[5 * ind + 1], vdvi[5 * ind + 2], 3, mu, M, nf);
891// Ue2os(vdvi[5 * ind], vdvi[5 * ind + 3]) += (lambda * lambda * omega)
892// * (-vdvi[5 * ind + 4] * F(vdvi[5 * ind + 1], vdvi[5 * ind + 2], 3, mu, M, nf));
893 Ue2os(vdvi[5 * ind], vdvi[5 * ind + 3]) += (-vdvi[5 * ind + 4] * F(vdvi[5 * ind + 1], vdvi[5 * ind + 2], 3, mu, M, nf));
894 }
895 if (ind < list(20)) {
896
897 Ueos(vevi[5 * ind], vevi[5 * ind + 3]) += vevi[5 * ind + 4] * F(vevi[5 * ind + 1], vevi[5 * ind + 2], 4, mu, M, nf);
898 Ue2os2(vevi[5 * ind], vevi[5 * ind + 3]) += (vevi[5 * ind + 4] * (F(vevi[5 * ind + 1], vevi[5 * ind + 2], 5, mu, M, nf)
899 - F(vevi[5 * ind + 1], vevi[5 * ind + 2], 4, mu, M, nf)));
900 }
901 if (ind < list(21)) {
902
903 Ue2os(vfvi[5 * ind], vfvi[5 * ind + 3]) += (vfvi[5 * ind + 4] * F(vfvi[5 * ind + 1], vfvi[5 * ind + 2], 4, mu, M, nf));
904 }
905 if (ind < list(22)) {
906
907 Ue2os(vrvi[5 * ind], vrvi[5 * ind + 3]) += vrvi[5 * ind + 4] * R(vrvi[5 * ind + 1], vrvi[5 * ind + 2], 4, mu, M, nf);
908
909 }
910
911 }
912
913
914 Us = omega * Us;
915 Us2 = omega * omega * Us2;
916 Ueos = lambda * Ueos;
917 Ue = lambda * omega * Ue;
918 Ue2os2 = lambda * lambda * Ue2os2;
919 Ues = (omega * omega * lambda) * Ues;
920 Ue2os = (omega * lambda * lambda) * Ue2os;
921
922 switch(order_qed) {
923
924
925 case NLO_QED22:
926
927 *elem[NLO_QED22] = (*elem[NLO_QED22]) * resLO + (*elem[NLO_QED11]) * Ue + (*elem[LO]) * Ue2 +
928 (*elem[NLO_QED21]) * Ueos + (*elem[NNLO]) * Ue2os2 + (*elem[NLO]) * Ue2os;
929
930 case NLO_QED12:
931
932 *elem[NLO_QED12] =(*elem[NLO_QED11]) * Ueos + (*elem[NLO]) * Ue2os2 + (*elem[LO]) * Ue2os;
933
934 case NLO_QED21:
935
936 *elem[NLO_QED21] = (*elem[NLO_QED21]) * resLO + (*elem[NLO_QED11]) * Us +
937 (*elem[NLO]) * Ue + (*elem[LO]) * Ues + (*elem[NNLO]) * Ueos;
938
939 case NLO_QED02:
940
941 *elem[NLO_QED02] = (*elem[LO]) * Ue2os2;
942
943 case NLO_QED11:
944
945 *elem[NLO_QED11] = (*elem[NLO_QED11]) * resLO + (*elem[LO]) * Ue + (*elem[NLO]) * Ueos;
946
947
948 case LO_QED:
949
950 *elem[LO_QED] = (*elem[LO]) * Ueos;
951 break;
952 default:
953 throw std::runtime_error("Error in EvolBsmm::Df1Evol");
954 }
955
956 switch(order) {
957 case NNLO:
958
959 *elem[NNLO] = (*elem[LO]) * Us2 + (*elem[NNLO]) * resLO + (*elem[NLO]) * Us;
960 case NLO:
961
962 *elem[NLO] = (*elem[LO]) * Us + (*elem[NLO]) * resLO;
963 case LO:
964
965 *elem[LO] = (*elem[LO]) * resLO;
966 break;
967 default:
968 throw std::runtime_error("Error in EvolBsmm::Df1Evol");
969 }
970}
971
972double EvolBsmm::F(unsigned int i, unsigned int j, int x, double mu, double M, double nf)
973
974{
975 int value = 0;
976 int L = 6 - (int) nf;
977
978 double etai = pow(eta, a[L][i]);
979 double etajx = pow(eta, a[L][j]+ x - 3.);
980
981 double cut = 0.000000001;
982 double result = 0.;
983
984 if(fabs(a[L][j] + x - 3. - a[L][i]) < cut ) value = 0;
985 else value = 1;
986
987 switch(value) {
988 case 0:
989 result = etai * logeta;
990 break;
991 case 1:
992 result = (etajx - etai)/(a[L][j] + x - 3. - a[L][i]);
993 break;
994 }
995 return (result);
996}
997
998
999double EvolBsmm::R(unsigned int i, unsigned int j, int x, double mu, double M, double nf)
1000
1001{
1002 int value = 0;
1003 int L = 6 - (int) nf;
1004
1005 double etai = pow(eta, a[L][i]);
1006 double etajx = pow(eta, a[L][j] + x - 3.);
1007
1008 double cut = 0.000000001;
1009 double result = 0.;
1010
1011 if(fabs(a[L][j] + x - 3. - a[L][i]) < cut ) value = 0;
1012 else value = 1;
1013
1014 switch(value) {
1015 case 0:
1016 result = etai * logeta * logeta * 0.5;
1017 break;
1018 case 1:
1019 result = (etajx * logeta - ((etajx - etai)/(a[L][j] + x - 3. - a[L][i])))/(a[L][j] + x - 3. - a[L][i]);
1020 break;
1021 }
1022return (result);
1023}
1024
1025double EvolBsmm::G(unsigned int i, unsigned int p, unsigned int j, int x, int y, double mu, double M, double nf)
1026
1027{
1028 int value = 0;
1029 int L = 6 - (int) nf;
1030
1031 double etai = pow(eta, a[L][i]);
1032 double etapx = pow(eta, a[L][p]+ x - 3.);
1033 double etajxy = pow(eta, a[L][j]+ x - 3. + y - 3.);
1034
1035 double cut = 0.000000001;
1036 double result = 0.;
1037
1038 if(fabs(a[L][j] + y - 3. - a[L][p]) < cut && fabs(a[L][p] + x - 3. - a[L][i]) < cut ) value = 0;
1039 else if(fabs(a[L][j] + y - 3. - a[L][p]) < cut && fabs(a[L][p] + x - 3. - a[L][i]) > cut) value = 1;
1040 else if(fabs(a[L][j] + y - 3. - a[L][p]) > cut && fabs(a[L][j] + y - 3. + x - 3. - a[L][i]) < cut && fabs(a[L][p] + x - 3. - a[L][i]) < cut ) value = 2;
1041 else if(fabs(a[L][j] + y - 3. - a[L][p]) > cut && fabs(a[L][j] + y - 3. + x - 3. - a[L][i]) > cut && fabs(a[L][p] + x - 3. - a[L][i]) > cut ) value = 3;
1042 else if(fabs(a[L][j] + y - 3. - a[L][p]) > cut && fabs(a[L][j] + y - 3. + x - 3. - a[L][i]) < cut && fabs(a[L][p] + x - 3. - a[L][i]) > cut ) value = 4;
1043 else if(fabs(a[L][j] + y - 3. - a[L][p]) > cut && fabs(a[L][j] + y - 3. + x - 3. - a[L][i]) > cut && fabs(a[L][p] + x - 3. - a[L][i]) < cut ) value = 5;
1044
1045 switch(value) {
1046 case 0:
1047 result = etai * logeta * log(eta) * 0.5;
1048 break;
1049 case 1:
1050 result = ((etapx * logeta - ((etapx - etai)/(a[L][p] + x - 3. - a[L][i])))/(a[L][p] + x - 3. - a[L][i]));
1051 break;
1052 case 2:
1053 result = (etai * logeta - etai * logeta)/(a[L][j] + y - 3. - a[L][p]);
1054 break;
1055 case 3:
1056 result = (((etajxy - etai)/(a[L][j] + x - 3. + y - 3. - a[L][i])) - ((etapx - etai)/(a[L][p] + x - 3. - a[L][i])))/(a[L][j] + y - 3. - a[L][p]);
1057 break;
1058 case 4:
1059 result = (etai * logeta - ((etapx - etai)/(a[L][p] + x - 3. - a[L][i])))/(a[L][j] + y - 3. - a[L][p]);
1060 break;
1061 case 5:
1062 result = (((etajxy - etai)/(a[L][j] + x - 3. + y - 3. - a[L][i])) - etai * log(eta))/(a[L][j] + y - 3. - a[L][p]);
1063 break;
1064 }
1065 return (result);
1066}
1067
1068double EvolBsmm::H(unsigned int i, unsigned int p, unsigned int q, unsigned int j, int x, int y, int z, double mu, double M, double nf)
1069
1070{
1071 int value = 0;
1072 int L = 6 - (int) nf;
1073
1074 double etai = pow(eta, a[L][i]);
1075 double etapx = pow(eta, a[L][p] + x - 3.);
1076 double etaqx = pow(eta, a[L][q] + x - 3.);
1077
1078 double cut = 0.000000001;
1079 double result = 0.;
1080
1081 if(fabs(a[L][p] + x - 3. - a[L][i]) < cut && fabs(a[L][q] + y - 3. - a[L][p]) < cut && fabs(a[L][j] + z - 3. - a[L][q]) < cut) value = 0;
1082 else if(fabs(a[L][p] + x - 3. - a[L][i]) > cut && fabs(a[L][q] + y - 3. - a[L][p]) < cut && fabs(a[L][j] + z - 3. - a[L][q]) < cut) value = 1;
1083 else if((a[L][q] + x - 3. + y + 3 - a[L][i] ) < cut && fabs(a[L][j] + z - 3. - a[L][q]) < cut){
1084 if(fabs(a[L][p] + x - 3. - a[L][i]) < cut) value = 2;
1085 else value = 3;
1086 }
1087 else if((a[L][q] + x - 3. + y + 3 - a[L][i]) > cut && fabs(a[L][j] + z - 3. - a[L][q]) < cut){
1088 if(fabs(a[L][p] + x - 3. - a[L][i]) > cut) value = 4;
1089 else value = 5;
1090 }
1091 else if(fabs(a[L][j] + z - 3. - a[L][q]) > cut) value = 6;
1092
1093 switch(value) {
1094 case 0:
1095 result = etai * logeta * logeta * logeta /6.;
1096 break;
1097 case 1:
1098 result = (0.5 * etapx * logeta * logeta - ((etapx * logeta
1099 - ((etapx - etai)/(a[L][p] + x - 3. - a[L][i])))/(a[L][p] + x - 3. - a[L][i])))/(a[L][p] + x - 3. - a[L][i]);
1100 break;
1101 case 2:
1102 result = (etai * logeta * logeta * 0.5 - ((etai * logeta - etai * logeta)/(a[L][q] + y - 3. - a[L][p])))/(a[L][q] + y - 3. - a[L][p]);
1103 break;
1104 case 3:
1105 result = (etai * logeta * logeta * 0.5 - ((etai * logeta- ((etapx - etai)/(a[L][p] + x - 3. - a[L][i])) )/(a[L][q] + y - 3. - a[L][p])))/(a[L][q] + y - 3. - a[L][p]);
1106 break;
1107 case 4:
1108 result = (((etaqx * logeta - ((etaqx - etai)/(a[L][q] + x - 3. - a[L][i])))/(a[L][q] + y + 3. + x - 3. - a[L][i])) - ((etai * logeta - etai * logeta)/(a[L][q] + y - 3. - a[L][p])))/(a[L][q] + y - 3. - a[L][p]);
1109 break;
1110 case 5:
1111 result = (((etaqx * logeta - ((etaqx - etai)/(a[L][q] + x - 3. - a[L][i])))/(a[L][q] + y + 3. + x - 3. - a[L][i])) - ((etai * logeta - ((etapx - etai)/(a[L][p] + x - 3. - a[L][i])) )/(a[L][q] + y - 3. - a[L][p])))/(a[L][q] + y - 3. - a[L][p]);
1112 break;
1113 case 6:
1114 result = (G(i, p, j, x , y + z - 3., mu, M, nf) - G(i, p, q, x , y, mu, M, nf))/(a[L][j] + z - 3. - a[L][q]);
1115 break;
1116 }
1117return (result);
1118}
1119
1120double EvolBsmm::alphatilde_e(double mu)
1121
1122{ // also the running is only for nf = 5
1123
1124 //double mu_0 = 91.1876;
1125 double mu_0 = model.getMz();
1126 //double alphatilde_e = 1./(127.751 * 4. * M_PI); // alpha_e at mu_0 = 91.1876 Gev
1127 double alphatilde_e = model.alphaMz()/4./M_PI;
1128 //double alphatilde_s = 0.1184/(4.* M_PI); // alpha_s at mu_0 = 91.1876 Gev
1129 double alphatilde_s = model.getAlsMz()/4./M_PI;
1130 unsigned int nf = 5;
1131
1132 double B00S = model.Beta0(nf), B10S = model.Beta1(nf);//B20S = model.Beta2(nf),
1133 double B01S = -22./9.; //B11S = -308./27.;
1134
1135 double B00E = 80./9., B01E = 176./9., B10E = 464./27.;
1136
1137 //double b1 = B10S/(2. * B00S * B00S); //b2 = B20S/(4. * B00S * B00S * B00S) - b1 * b1 ,
1138 //double b3 = B01S/(2. * B00S * B00E);// b4 = B11S /(4. * B00S * B00S * B00E) - 2 * b1 * b3,
1139 //b5 = B01E/(2. * B00S * B00E) - b1;
1140
1141 double vs= 1. + 2. * B00S * alphatilde_s * log(mu/ mu_0);
1142 double ve= 1. - 2. * B00E * alphatilde_e * log(mu/ mu_0);
1143 //double ps= B00S * alphatilde_s /(B00S * alphatilde_s + B00E * alphatilde_e);
1144 double pe= B00E * alphatilde_e /(B00S * alphatilde_s + B00E * alphatilde_e);
1145
1146 double logve = log(ve);
1147 double logvs = log(vs);
1148 double logeos = log(ve/vs);
1149 double asovs = alphatilde_s/vs;
1150 double aeove = alphatilde_e/ve;
1151
1152 double result = 0;
1153
1154 result = aeove - pow(aeove, 2) * (logve * B10E/ B00E - logvs * B01E/B00S)
1155 + pow(aeove, 2) * (asovs) * ((logvs - vs + 1.) * B01E * B10S/(B00S * B00S)
1156 + logve * vs * pe * B01E * B10E/(B00E * B00E) +(logeos * vs * pe - logve) * B01E * B01S/(B00S * B00E));
1157 return (result);
1158}
1159
1160double EvolBsmm::alphatilde_s(double mu)
1161
1162{ // also the running is only for nf = 5
1163
1164 //double mu_0 = 91.1876;
1165 double mu_0 = model.getMz();
1166 //double alphatilde_e = 1./(127.751 * 4. * M_PI); // alpha_e at mu_0 = 91.1876 Gev
1167 double alphatilde_e = model.alphaMz()/4./M_PI;
1168 //double alphatilde_s = 0.1184/(4.* M_PI); // alpha_s at mu_0 = 91.1876 Gev
1169 double alphatilde_s = model.getAlsMz()/4./M_PI;
1170 unsigned int nf = 5;
1171
1172 double B00S = model.Beta0(nf), B10S = model.Beta1(nf), B20S = model.Beta2(nf), B30S = gsl_sf_zeta_int(3) * 352864./81. - 598391./1458,
1173 B01S = -22./9., B11S = -308./27., B02S = 4945./243.;
1174
1175 double B00E = 80./9., B01E = 176./9., B10E = 464./27.;
1176
1177 //double B00S2 = B00S * B00S;
1178 double B10soB00s = B10S / B00S;
1179 double B01soB00e = B01S/B00E;
1180
1181 //double b1 = B10soB00s/(2. * B00S), b2 = B20S/(4. * B00S2 * B00S) - b1 * b1 ,
1182 // b3 = B01soB00e/(2. * B00S ), b4 = B11S /(4. * B00S2 * B00E) - 2 * b1 * b3,
1183 // b5 = B01E/(2. * B00S * B00E) - b1;
1184
1185 double vs= 1. + 2. * B00S * alphatilde_s * log(mu/ mu_0);
1186 double ve= 1. - 2. * B00E * alphatilde_e * log(mu/ mu_0);
1187 double ps= B00S * alphatilde_s /(B00S * alphatilde_s + B00E * alphatilde_e);
1188 //double pe= B00E * alphatilde_e /(B00S * alphatilde_s + B00E * alphatilde_e);
1189
1190 double logve = log(ve);
1191 double logvs = log(vs);
1192 double logeos = log(ve/vs);
1193 double logsoe = log(vs/ve);
1194 double asovs = alphatilde_s/vs;
1195 double aeove = alphatilde_e/ve;
1196
1197 double result = 0;
1198
1199 result = asovs - pow(asovs, 2) * (logvs * B10soB00s - logve * B01soB00e)
1200 + pow(asovs, 3) * ((1. - vs) * B20S / B00S + B10soB00s * B10soB00s * (logvs * logvs - logvs
1201 + vs - 1.) + B01soB00e * B01soB00e * logve * logve + (-2. * logvs * logve
1202 + ps * ve * logve) * B01S * B10S/(B00E * B00S))
1203 + pow(asovs, 4) * (0.5 * B30S *(1. - vs * vs)/ B00S + ((2. * vs - 3.) * logvs + vs * vs
1204 - vs) * B20S * B10soB00s /(B00S) + B10soB00s * B10soB00s * B10soB00s * (- pow(logvs,3)
1205 + 5. * pow(logvs,2) / 2. + 2. * (1. - vs) * logvs - (vs - 1.) * (vs - 1.)* 0.5))
1206 + pow(asovs, 2) * (aeove) * ((ve - 1.) * B02S / B00E
1207 + ps * ve * logeos * B11S /B00S +(logve - ve + 1.) * B01soB00e * B10E/(B00E)
1208 + logvs * ps * B01S * B10soB00s/(B00S) +(logsoe * ve * ps - logvs) * B01soB00e * B01E/( B00S));
1209 return (result);
1210}
orders
orders
Definition OrderScheme.h:33
NNLO
@ NNLO
Definition OrderScheme.h:36
LO
@ LO
Definition OrderScheme.h:34
NLO
@ NLO
Definition OrderScheme.h:35
schemes
schemes
Definition OrderScheme.h:20
HV
@ HV
Definition OrderScheme.h:22
LRI
@ LRI
Definition OrderScheme.h:23
NDR
@ NDR
Definition OrderScheme.h:21
orders_qed
orders_qed
Definition OrderScheme.h:56
NLO_QED11
@ NLO_QED11
Definition OrderScheme.h:59
NLO_QED02
@ NLO_QED02
Definition OrderScheme.h:61
NLO_QED12
@ NLO_QED12
Definition OrderScheme.h:62
LO_QED
@ LO_QED
Definition OrderScheme.h:58
NLO_QED22
@ NLO_QED22
Definition OrderScheme.h:63
NLO_QED21
@ NLO_QED21
Definition OrderScheme.h:60
NO_QED
@ NO_QED
Definition OrderScheme.h:57
EvolBsmm::CC
gslpp::matrix< gslpp::complex > CC
Definition EvolBsmm.h:31
EvolBsmm::vdevi
std::vector< double > vdevi
Definition EvolBsmm.h:34
EvolBsmm::veebvi
std::vector< double > veebvi
Definition EvolBsmm.h:35
EvolBsmm::vfvi
std::vector< double > vfvi
Definition EvolBsmm.h:33
EvolBsmm::EvolBsmm
EvolBsmm(unsigned int dim, schemes scheme, orders order, orders_qed order_qed, const StandardModel &model)
Definition EvolBsmm.cpp:12
EvolBsmm::AnomalousDimension
gslpp::matrix< double > AnomalousDimension(int gam, unsigned int n_u, unsigned int n_d) const
Definition EvolBsmm.cpp:340
EvolBsmm::vbvi
std::vector< double > vbvi
Definition EvolBsmm.h:33
EvolBsmm::vbdvi
std::vector< double > vbdvi
Definition EvolBsmm.h:34
EvolBsmm::nd
int nd
Definition EvolBsmm.h:27
EvolBsmm::BuiltB
gslpp::matrix< double > BuiltB(char letter, unsigned int n_u, unsigned int n_d)
Definition EvolBsmm.cpp:596
EvolBsmm::~EvolBsmm
virtual ~EvolBsmm()
Definition EvolBsmm.cpp:337
EvolBsmm::alphatilde_e
double alphatilde_e(double mu)
Definition EvolBsmm.cpp:1120
EvolBsmm::FF
gslpp::matrix< gslpp::complex > FF
Definition EvolBsmm.h:31
EvolBsmm::vdbvi
std::vector< double > vdbvi
Definition EvolBsmm.h:34
EvolBsmm::vbeevi
std::vector< double > vbeevi
Definition EvolBsmm.h:34
EvolBsmm::R
double R(unsigned int i, unsigned int j, int x, double mu, double M, double nf)
Definition EvolBsmm.cpp:999
EvolBsmm::DD
gslpp::matrix< gslpp::complex > DD
Definition EvolBsmm.h:31
EvolBsmm::vebbvi
std::vector< double > vebbvi
Definition EvolBsmm.h:35
EvolBsmm::EE
gslpp::matrix< gslpp::complex > EE
Definition EvolBsmm.h:31
EvolBsmm::vedvi
std::vector< double > vedvi
Definition EvolBsmm.h:34
EvolBsmm::vebvi
std::vector< double > vebvi
Definition EvolBsmm.h:34
EvolBsmm::a
double a[4][8]
Definition EvolBsmm.h:28
EvolBsmm::logeta
double logeta
Definition EvolBsmm.h:41
EvolBsmm::Vi
gslpp::matrix< gslpp::complex > Vi
Definition EvolBsmm.h:31
EvolBsmm::vbebvi
std::vector< double > vbebvi
Definition EvolBsmm.h:35
EvolBsmm::vcvi
std::vector< double > vcvi
Definition EvolBsmm.h:33
EvolBsmm::G
double G(unsigned int i, unsigned int p, unsigned int j, int x, int y, double mu, double M, double nf)
Definition EvolBsmm.cpp:1025
EvolBsmm::RR
gslpp::matrix< gslpp::complex > RR
Definition EvolBsmm.h:31
EvolBsmm::vevi
std::vector< double > vevi
Definition EvolBsmm.h:33
EvolBsmm::nu
int nu
Definition EvolBsmm.h:27
EvolBsmm::AA
gslpp::matrix< gslpp::complex > AA
Definition EvolBsmm.h:31
EvolBsmm::vaevi
std::vector< double > vaevi
Definition EvolBsmm.h:33
EvolBsmm::F
double F(unsigned int i, unsigned int j, int x, double mu, double M, double nf)
Definition EvolBsmm.cpp:972
EvolBsmm::H
double H(unsigned int i, unsigned int p, unsigned int q, unsigned int j, int x, int y, int z, double mu, double M, double nf)
Definition EvolBsmm.cpp:1068
EvolBsmm::vbbvi
std::vector< double > vbbvi
Definition EvolBsmm.h:33
EvolBsmm::vbbevi
std::vector< double > vbbevi
Definition EvolBsmm.h:35
EvolBsmm::alsM
double alsM
Definition EvolBsmm.h:37
EvolBsmm::veevi
std::vector< double > veevi
Definition EvolBsmm.h:34
EvolBsmm::model
const StandardModel & model
Definition EvolBsmm.h:29
EvolBsmm::b
double b[4][8][8][8]
Definition EvolBsmm.h:28
EvolBsmm::dim
unsigned int dim
Definition EvolBsmm.h:36
EvolBsmm::vbevi
std::vector< double > vbevi
Definition EvolBsmm.h:34
EvolBsmm::Df1Evol
gslpp::matrix< double > & Df1Evol(double mu, double M, orders order, orders_qed order_qed, schemes scheme=NDR)
Definition EvolBsmm.cpp:678
EvolBsmm::eta
double eta
Definition EvolBsmm.h:40
EvolBsmm::vavi
std::vector< double > vavi
Definition EvolBsmm.h:33
EvolBsmm::BB
gslpp::matrix< gslpp::complex > BB
Definition EvolBsmm.h:31
EvolBsmm::alsmu
double alsmu
Definition EvolBsmm.h:38
EvolBsmm::V
gslpp::matrix< gslpp::complex > V
Definition EvolBsmm.h:31
EvolBsmm::vebevi
std::vector< double > vebevi
Definition EvolBsmm.h:35
EvolBsmm::e
gslpp::vector< gslpp::complex > e
Definition EvolBsmm.h:32
EvolBsmm::veavi
std::vector< double > veavi
Definition EvolBsmm.h:34
EvolBsmm::vdvi
std::vector< double > vdvi
Definition EvolBsmm.h:33
EvolBsmm::alphatilde_s
double alphatilde_s(double mu)
Definition EvolBsmm.cpp:1160
EvolBsmm::vrvi
std::vector< double > vrvi
Definition EvolBsmm.h:33
QCD::Beta2
const double Beta2(const double nf) const
The coefficient for a certain number of flavours .
Definition QCD.cpp:611
QCD::Beta1
const double Beta1(const double nf) const
The coefficient for a certain number of flavours .
Definition QCD.cpp:606
QCD::Beta0
const double Beta0(const double nf) const
The coefficient for a certain number of flavours .
Definition QCD.cpp:601
RGEvolutor
A class for the RG evolutor of the Wilson coefficients.
Definition RGEvolutor.h:24
RGEvolutor::setScales
void setScales(double mu, double M)
Sets the upper and lower scale for the running of the Wilson Coefficients.
Definition RGEvolutor.cpp:85
RGEvolutor::Evol
gslpp::matrix< double > * Evol(orders order)
Evolution matrix set at a fixed order of QCD coupling.
Definition RGEvolutor.cpp:103
StandardModel
A model class for the Standard Model.
Definition StandardModel.h:521
StandardModel::getMz
const double getMz() const
A get method to access the mass of the boson .
Definition StandardModel.h:765
StandardModel::getAlsMz
const double getAlsMz() const
A get method to access the value of .
Definition StandardModel.h:783
StandardModel::alphaMz
virtual const double alphaMz() const
The electromagnetic coupling at the -mass scale, .
Definition StandardModel/src/StandardModel.cpp:944
WilsonTemplate< gslpp::matrix< double > >::elem
gslpp::matrix< double > * elem[MAXORDER_QED+1]
Definition WilsonTemplate.h:114