a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
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 
12 EvolBsmm::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 
337 EvolBsmm::~EvolBsmm()
338 {}
339 
340 gslpp::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 
596 gslpp::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 
678 gslpp::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 
733 void 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 
972 double 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 
999 double 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  }
1022 return (result);
1023 }
1024 
1025 double 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 
1068 double 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  }
1117 return (result);
1118 }
1119 
1120 double 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 
1160 double 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 }
WilsonTemplate< gslpp::matrix< double > >::scheme
schemes scheme
Definition: WilsonTemplate.h:117
EvolBsmm::RR
gslpp::matrix< gslpp::complex > RR
Definition: EvolBsmm.h:31
EvolBsmm::vbbevi
std::vector< double > vbbevi
Definition: EvolBsmm.h:33
EvolBsmm::eta
double eta
Definition: EvolBsmm.h:40
WilsonTemplate< gslpp::matrix< double > >::order_qed
orders_qed order_qed
Definition: WilsonTemplate.h:119
NO_QED
Definition: OrderScheme.h:49
EvolBsmm::vbvi
std::vector< double > vbvi
Definition: EvolBsmm.h:33
gslpp::matrix< double >
A class for constructing and defining operations on real matrices.
Definition: gslpp_matrix_double.h:48
EvolBsmm::AA
gslpp::matrix< gslpp::complex > AA
Definition: EvolBsmm.h:31
EvolBsmm::a
double a[4][8]
Definition: EvolBsmm.h:28
NLO_QED11
Definition: OrderScheme.h:51
QCD::Beta1
double Beta1(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:471
RGEvolutor
A class for the RG evolutor of the Wilson coefficients.
Definition: RGEvolutor.h:24
WilsonTemplate< gslpp::matrix< double > >::mu
double mu
Definition: WilsonTemplate.h:116
EvolBsmm::vedvi
std::vector< double > vedvi
Definition: EvolBsmm.h:33
EvolBsmm::vebbvi
std::vector< double > vebbvi
Definition: EvolBsmm.h:33
RGEvolutor::M
double M
Definition: RGEvolutor.h:142
EvolBsmm::vaevi
std::vector< double > vaevi
Definition: EvolBsmm.h:33
StandardModel::getAlsMz
double getAlsMz() const
A get method to access the value of .
Definition: StandardModel.h:730
WilsonTemplate< gslpp::matrix< double > >::order
orders order
Definition: WilsonTemplate.h:118
EvolBsmm::Vi
gslpp::matrix< gslpp::complex > Vi
Definition: EvolBsmm.h:31
EvolBsmm::AnomalousDimension
gslpp::matrix< double > AnomalousDimension(int gam, unsigned int n_u, unsigned int n_d) const
Definition: EvolBsmm.cpp:340
LO
Definition: OrderScheme.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
StandardModel.h
StandardModel::alphaMz
double alphaMz() const
The electromagnetic coupling at the -mass scale, .
Definition: StandardModel.cpp:893
EvolBsmm::vdvi
std::vector< double > vdvi
Definition: EvolBsmm.h:33
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
NDR
Definition: OrderScheme.h:21
EvolBsmm::EE
gslpp::matrix< gslpp::complex > EE
Definition: EvolBsmm.h:31
EvolBsmm::CC
gslpp::matrix< gslpp::complex > CC
Definition: EvolBsmm.h:31
EvolBsmm::alphatilde_s
double alphatilde_s(double mu)
Definition: EvolBsmm.cpp:1160
EvolBsmm::alphatilde_e
double alphatilde_e(double mu)
Definition: EvolBsmm.cpp:1120
EvolBsmm::vavi
std::vector< double > vavi
Definition: EvolBsmm.h:33
gslpp::matrix< double >::transpose
matrix< double > transpose() const
Definition: gslpp_matrix_double.cpp:166
EvolBsmm::BB
gslpp::matrix< gslpp::complex > BB
Definition: EvolBsmm.h:31
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:477
EvolBsmm::vbeevi
std::vector< double > vbeevi
Definition: EvolBsmm.h:33
EvolBsmm::DD
gslpp::matrix< gslpp::complex > DD
Definition: EvolBsmm.h:31
QCD::Beta0
double Beta0(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:466
QCD::Beta2
double Beta2(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:476
EvolBsmm::V
gslpp::matrix< gslpp::complex > V
Definition: EvolBsmm.h:31
EvolBsmm::vdbvi
std::vector< double > vdbvi
Definition: EvolBsmm.h:33
EvolBsmm::vrvi
std::vector< double > vrvi
Definition: EvolBsmm.h:33
EvolBsmm::vevi
std::vector< double > vevi
Definition: EvolBsmm.h:33
EvolBsmm::veavi
std::vector< double > veavi
Definition: EvolBsmm.h:33
schemes
schemes
An enum type for regularization schemes.
Definition: OrderScheme.h:19
EvolBsmm::veebvi
std::vector< double > veebvi
Definition: EvolBsmm.h:33
gslpp::vector< double >::max
double max() const
Definition: gslpp_vector_double.cpp:100
EvolBsmm::R
double R(unsigned int i, unsigned int j, int x, double mu, double M, double nf)
Definition: EvolBsmm.cpp:999
EvolBsmm::BuiltB
gslpp::matrix< double > BuiltB(char letter, unsigned int n_u, unsigned int n_d)
Definition: EvolBsmm.cpp:596
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
EvolBsmm::FF
gslpp::matrix< gslpp::complex > FF
Definition: EvolBsmm.h:31
EvolBsmm::vbbvi
std::vector< double > vbbvi
Definition: EvolBsmm.h:33
EvolBsmm::vdevi
std::vector< double > vdevi
Definition: EvolBsmm.h:33
EvolBsmm::nu
int nu
Definition: EvolBsmm.h:27
EvolBsmm::veevi
std::vector< double > veevi
Definition: EvolBsmm.h:33
LO_QED
Definition: OrderScheme.h:50
EvolBsmm::vbevi
std::vector< double > vbevi
Definition: EvolBsmm.h:33
NNLO
Definition: OrderScheme.h:35
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
NLO_QED22
Definition: OrderScheme.h:55
orders_qed
orders_qed
An enum type for orders in electroweak.
Definition: OrderScheme.h:47
EvolBsmm::b
double b[4][8][8][8]
Definition: EvolBsmm.h:28
EvolBsmm::Df1Evol
gslpp::matrix< double > & Df1Evol(double mu, double M, orders order, orders_qed order_qed, schemes scheme=NDR)
Definition: EvolBsmm.cpp:678
gslpp::vector< double >
A class for constructing and defining operations on real vectors.
Definition: gslpp_vector_double.h:33
EvolBsmm::EvolBsmm
EvolBsmm(unsigned int dim, schemes scheme, orders order, orders_qed order_qed, const StandardModel &model)
Definition: EvolBsmm.cpp:12
EvolBsmm::alsmu
double alsmu
Definition: EvolBsmm.h:38
EvolBsmm::alsM
double alsM
Definition: EvolBsmm.h:37
orders
orders
An enum type for orders in QCD.
Definition: OrderScheme.h:31
NLO_QED12
Definition: OrderScheme.h:54
LRI
Definition: OrderScheme.h:23
StandardModel::getMz
double getMz() const
A get method to access the mass of the boson .
Definition: StandardModel.h:721
EvolBsmm::vbdvi
std::vector< double > vbdvi
Definition: EvolBsmm.h:33
EvolBsmm::nd
int nd
Definition: EvolBsmm.h:27
RGEvolutor::Evol
gslpp::matrix< double > * Evol(orders order)
Evolution matrix set at a fixed order of QCD coupling.
Definition: RGEvolutor.cpp:103
HV
Definition: OrderScheme.h:22
EvolBsmm::vbebvi
std::vector< double > vbebvi
Definition: EvolBsmm.h:33
EvolBsmm::model
const StandardModel & model
Definition: EvolBsmm.h:29
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
EvolBsmm::vebvi
std::vector< double > vebvi
Definition: EvolBsmm.h:33
NLO
Definition: OrderScheme.h:34
NLO_QED21
Definition: OrderScheme.h:52
NLO_QED02
Definition: OrderScheme.h:53
EvolBsmm::logeta
double logeta
Definition: EvolBsmm.h:41
EvolBsmm::vebevi
std::vector< double > vebevi
Definition: EvolBsmm.h:33
gslpp::matrix< double >::eigensystem
void eigensystem(matrix< complex > &U, vector< complex > &S)
Definition: gslpp_matrix_double.cpp:280
EvolBsmm::vfvi
std::vector< double > vfvi
Definition: EvolBsmm.h:33
EvolBsmm::dim
unsigned int dim
Definition: EvolBsmm.h:36
EvolBsmm::e
gslpp::vector< gslpp::complex > e
Definition: EvolBsmm.h:32
WilsonTemplate< gslpp::matrix< double > >::elem
gslpp::matrix< double > * elem[MAXORDER_QED+1]
Definition: WilsonTemplate.h:114
EvolBsmm::~EvolBsmm
virtual ~EvolBsmm()
Definition: EvolBsmm.cpp:337
EvolBsmm::vcvi
std::vector< double > vcvi
Definition: EvolBsmm.h:33