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