a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
EvolDB1Mll.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 <gsl/gsl_sf_zeta.h>
9#include "EvolDB1Mll.h"
10#include "StandardModel.h"
11
12EvolDB1Mll::EvolDB1Mll(unsigned int dim_i, schemes scheme, orders order, const StandardModel& model)
13: RGEvolutor(dim_i, scheme, order), model(model),
14 v(dim_i,0.), vi(dim_i,0.), js(dim_i,0.), h(dim_i,0.), gg(dim_i,0.), s_s(dim_i,0.),
15 jssv(dim_i,0.), jss(dim_i,0.), jv(dim_i,0.), vij(dim_i,0.), e(dim_i,0.), dim(dim_i)
16{
17 if (dim != 13 ) throw std::runtime_error("ERROR: EvolDB1Mll can only be of dimension 13");
18
19 /* magic numbers a & b */
20
21 for(int L=3; L>-1; L--){
22
23 /* L=3 --> u,d,s (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
25 nu = L; nd = L;
26 if(L == 3){nd = 2; nu = 1;}
27 if(L == 1){nd = 3; nu = 2;}
28 if(L == 0){nd = 3; nu = 3;}
29
30 // LO evolutor of the effective Wilson coefficients in the Chetyrkin, Misiak and Munz basis
31
32 (ToEffectiveBasis(ToRescaleBasis(LO,nu,nd))).transpose().eigensystem(v,e);
33 vi = v.inverse();
34 for(unsigned int i = 0; i < dim; i++){
35 a[L][i] = e(i).real();
36 for (unsigned int j = 0; j < dim; j++) {
37 for (unsigned int k = 0; k < dim; k++) {
38 b[L][i][j][k] = v(i, k).real() * vi(k, j).real();
39 }
40 }
41 }
42
43 // NLO evolutor of the effective Wilson coefficients in the Chetyrkin, Misiak and Munz basis
44
45 gg = vi * (ToEffectiveBasis(ToRescaleBasis(NLO,nu,nd))).transpose() * v;
46 double b0 = model.Beta0(nu+nd);
47 double b1 = model.Beta1(nu+nd);
48 for (unsigned int i = 0; i < dim; i++){
49 for (unsigned int j = 0; j < dim; j++){
50 s_s.assign( i, j, (b1 / b0) * (i==j) * e(i).real() - gg(i,j));
51 if(fabs(e(i).real() - e(j).real() + 2. * b0)>0.00000000001){
52 h.assign( i, j, s_s(i,j) / (2. * b0 + e(i) - e(j)));
53 }
54 }
55 }
56 js = v * h * vi;
57 jv = js * v;
58 vij = vi * js;
59 jss = v * s_s * vi;
60 jssv = jss * v;
61 for (unsigned int i = 0; i < dim; i++){
62 for (unsigned int j = 0; j < dim; j++){
63 if(fabs(e(i).real() - e(j).real() + 2. * b0) > 0.00000000001){
64 for(unsigned int k = 0; k < dim; k++){
65 c[L][i][j][k] = jv(i, k).real() * vi(k, j).real();
66 d[L][i][j][k] = -v(i, k).real() * vij(k, j).real();
67 }
68 }
69 else{
70 for(unsigned int k = 0; k < dim; k++){
71 c[L][i][j][k] = (1./(2. * b0)) * jssv(i, k).real() * vi(k, j).real();
72 d[L][i][j][k] = 0.;
73 }
74 }
75 }
76 }
77 }
78}
79
81{}
82
83gslpp::matrix<double> EvolDB1Mll::AnomalousDimension_M(orders order, unsigned int n_u, unsigned int n_d) const
84{
85
86 /* Delta F = 1 anomalous dimension in Misiak basis,
87 ref.: M. Misiak, Nucl. Phys. B393 (1993) 23, B439 (1995) 461 (E),
88 A.J. Buras and M. Munz, Phys. Rev. D52 (1995) 186. */
89
90 /* gamma(row, coloumn) at the LO */
91
92 unsigned int nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
93
94 gslpp::matrix<double> gammaDF1(dim, dim, 0.);
95
96 switch(order){
97
98 case LO:
99
100 gammaDF1(0,0) = -4. ;
101 gammaDF1(0,1) = 8./3. ;
102 gammaDF1(0,3) = -2./9.;
103 gammaDF1(0,8) = -32./27.;
104
105
106 gammaDF1(1,0) = 12.;
107 gammaDF1(1,3) = 4./3.;
108 gammaDF1(1,8) = -8./9.;
109
110 gammaDF1(2,3) = -52./3.;
111 gammaDF1(2,5) = 2.;
112 gammaDF1(2,8) = 8./9. + (8.*n_d)/3. - (16.*n_u)/3.;
113
114 gammaDF1(3,2) = -40./9.;
115 gammaDF1(3,3) = -160./9. + 4./3.*nf;
116 gammaDF1(3,4) = 4./9.;
117 gammaDF1(3,5) = 5./6.;
118 gammaDF1(3,8) = 32./27.;
119
120 gammaDF1(4,3) = -256./3.;
121 gammaDF1(4,5) = 20.;
122 gammaDF1(4,8) = 128./9.+(80.*n_d)/3. - (160.*n_u)/3.;
123
124 gammaDF1(5,2) = -256./9.;
125 gammaDF1(5,3) = -544./9. + (40./3.)*nf;
126 gammaDF1(5,4) = 40./9.;
127 gammaDF1(5,5) = -2./3.;
128 gammaDF1(5,8) = 512./27.;
129
130 gammaDF1(6,6) = 32./3. - 2.*model.Beta0(nf);
131
132 gammaDF1(7,6) = -32./9.;
133 gammaDF1(7,7) = 28./3. - 2.*model.Beta0(nf);
134
135 gammaDF1(8,8) = -2.*model.Beta0(nf);
136
137 gammaDF1(9,9) = -2.*model.Beta0(nf);
138
139 gammaDF1(10,10)= -2.*model.Beta0(nf);
140
141 gammaDF1(11,11)= -2.*model.Beta0(nf);
142
143 gammaDF1(12,12)= -2.*model.Beta0(nf);
144
145 break;
146 case NLO:
147
148 if (!(nf == 3 || nf == 4 || nf == 5 || nf == 6)){
149 throw std::runtime_error("EvolDF1::AnomalousDimension_M(): wrong number of flavours");
150 }
151
152 /* gamma(row, coloumn) at the NLO */
153
154 gammaDF1(0,0) = -145./3. + (16./9.)*nf;
155 gammaDF1(0,1) = -26. + (40./27.)*nf;
156 gammaDF1(0,2) = -1412./243.;
157 gammaDF1(0,3) = -1369./243.;
158 gammaDF1(0,4) = 134./243.;
159 gammaDF1(0,5) = -35./162.;
160 gammaDF1(0,6) = -232./243.;
161 gammaDF1(0,7) = +167./162.;
162 gammaDF1(0,8) = -2272./729.;
163
164 gammaDF1(1,0) = -45. + (20./3.)*nf;
165 gammaDF1(1,1) = -28./3.;
166 gammaDF1(1,2) = -416./81.;
167 gammaDF1(1,3) = 1280./81.;
168 gammaDF1(1,4) = 56./81.;
169 gammaDF1(1,5) = 35./27.;
170 gammaDF1(1,6) = 464./81.;
171 gammaDF1(1,7) = 76./27.;
172 gammaDF1(1,8) = 1952./243.;
173
174 gammaDF1(2,2) = -4468./81.;
175 gammaDF1(2,3) = -29129./81. - (52./9.)*nf;
176 gammaDF1(2,4) = 400./81.;
177 gammaDF1(2,5) = 3493./108. - (2./9.)*nf;
178 gammaDF1(2,6) = 64./81.;
179 gammaDF1(2,7) = 368./27.;
180 gammaDF1(2,8) = -4160./243. + (32.*n_d)/3. - (64.*n_u)/3.;
181
182 gammaDF1(3,2) = -13678./243. + (368.*nf)/81.;
183 gammaDF1(3,3) = -79409./243. + (1334.*nf)/81.;
184 gammaDF1(3,4) = 509./486. - (8.*nf)/81.;
185 gammaDF1(3,5) = 13499./648. - (5.*nf)/27.;
186 gammaDF1(3,6) = -680./243. + (32.*nf)/81;
187 gammaDF1(3,7) = -427./81. - (37.*nf)/54.;
188 gammaDF1(3,8) = 784./729. - (2272.*n_d)/243. + (2912.*n_u)/243.;
189
190 gammaDF1(4,2) = -244480./81. - (160./9.)*nf;
191 gammaDF1(4,3) = -29648./81. - (2200./9.)*nf;
192 gammaDF1(4,4) = 23116./81. + (16./9.)*nf;
193 gammaDF1(4,5) = 3886./27. + (148./9.)*nf;
194 gammaDF1(4,6) = -6464./81.;
195 gammaDF1(4,7) = 8192./27. + 36.*nf;
196 gammaDF1(4,8) = -58112./243. + (320.*n_d)/3. - (640.*n_u)/3.;
197
198 gammaDF1(5,2) = 77600./243. - (1264./81.)*nf;
199 gammaDF1(5,3) = -28808./243. + (164./81.)*nf;
200 gammaDF1(5,4) = -20324./243. + (400./81.)*nf;
201 gammaDF1(5,5) = -21211./162.+ (622./27.)*nf;
202 gammaDF1(5,6) = -20096./243. - (976.*n_d)/81. + (2912.*n_u)/81.;
203 gammaDF1(5,7) = -6040./81. + (220./27.)*nf;
204 gammaDF1(5,8) = -22784./729. - (20704.*n_d)/243. + (28544.*n_u)/243.;
205
206 gammaDF1(6,6) = 1936./9.-224./27.*nf-2*model.Beta1(nf);
207
208 gammaDF1(7,6) = -368./9.+224./81.*nf;
209 gammaDF1(7,7) = 1456./9.-61./27.*nf-2*model.Beta1(nf);
210
211 gammaDF1(8,8) = -2.*model.Beta1(nf);
212
213 gammaDF1(9,9) = -2.*model.Beta1(nf);
214
215 gammaDF1(10,10)= -2.*model.Beta1(nf);
216
217 gammaDF1(11,11)= -2.*model.Beta1(nf);
218
219 gammaDF1(12,12)= -2.*model.Beta1(nf);
220
221 break;
222 default:
223 throw std::runtime_error("EvolDF1bsg::AnomalousDimension_M(): order not implemented");
224 }
225 return (gammaDF1);
226}
227
228gslpp::matrix<double> EvolDB1Mll::ToRescaleBasis(orders order, unsigned int n_u, unsigned int n_d) const
229{
230
231 /* matrix entries for the anomalous dimension in the Chetyrkin, Misiak and Munz basis,
232 ref. hep-ph/9711280v1, hep-ph/0504194 */
233
234 gslpp::matrix<double> mat(dim, 0.);
235 gslpp::matrix<double> mat1(dim, 0.);
236 unsigned int nf = n_u + n_d;
237 double z3 = gsl_sf_zeta_int(3);
238
239 mat1(0,6) = - 13454./2187. + 44./2187.*nf;
240 mat1(1,6) = 20644./729. - 88./729.*nf;
241 mat1(2,6) = 119456./729. + 5440./729.*n_d -21776./729.*n_u;
242 mat1(3,6) = - 202990./2187. + 32./729.*n_d*n_d + n_d*(16888./2187. + 64./729.*n_u)
243 - 17132./2187.*n_u + 32./729.*n_u*n_u;
244 mat1(4,6) = 530240./243. + 300928./729.*n_d - 461120./729.*n_u;
245 mat1(5,6) = - 1112344./729. + 5432./729.*n_d*n_d + n_d*(419440./2187. -
246 2744./729.*n_u) + 143392./2187.*n_u - 8176./729.*n_u*n_u;
247
248 mat1(0,7) = 25759./5832. + 431./5832.*nf;
249 mat1(1,7) = 9733./486. - 917./972.*nf;
250 mat1(2,7) = 82873./243. - 3361./243.*nf;
251 mat1(3,7) = - 570773./2916. - 253./486.*n_d*n_d +n_d*(-40091./5832. -
252 253./243.*n_u) - 40091./5832.*n_u - 253./486.*n_u*n_u;
253 mat1(4,7) = 838684./81. - 14.*n_d*n_d + n_d*(129074./243. - 28.*n_u) +
254 129074./243.*n_u - 14.*n_u*n_u;
255 mat1(5,7) = - 923522./243. - 6031./486.*n_d*n_d + n_d*(-13247./1458. - 6031./243.*n_u)
256 -13247./1458.*n_u - 6031./486.*n_u*n_u;
257
258 mat1(0,8) = - 2357278./19683. + 14440./6561.*n_d + 144688./6561.*n_u + 6976./243.*z3;
259 mat1(1,8) = - 200848./6561. - 23696./2187.*n_d + 30736./2187.*n_u - 3584./81.*z3;
260 mat1(2,8) = - 1524104./6561. - 176./27.*n_d*n_d + 352./27.*n_u*n_u +
261 n_d*(257564./2187. + 176./27.*n_u - 128./3.*z3) - 256./81.*z3 +
262 n_u*(-382984./2187. + 256./3.*z3);
263 mat1(3,8) = 1535926./19683. + 1984./2187.*n_d*n_d - 5792./2187.*n_u*n_u +
264 n_d*(-256901./6561. - 3808./2187.*n_u - 2720./81.*z3) -
265 5056./243.*z3 + n_u*(34942./6561. + 1600./81.*z3);
266 mat1(4,8) = - 31433600./6561. - 2912./27.*n_d*n_d + 5824./27.*n_u*n_u +
267 n_d*(- 3786616./2187. + 2912./27.*n_u - 1280./3.*z3) -
268 4096./81.*z3 + n_u*(7525520./2187. + 2560./3.*z3);
269 mat1(5,8) = 48510784./19683. -51296./2187.*n_d*n_d + 54976./2187.*n_u*n_u +
270 n_u*(-11231648./6561. - 22016./81.*z3) + n_d*(340984./6561. +
271 3680./2187.*n_u - 8192./81.*z3) - 80896./243.*z3;
272
273
274 switch(order){
275 case(NLO):
276 mat = AnomalousDimension_M(NLO, n_u, n_d);
277 for (unsigned int i=0; i<6; i++){
278 for (unsigned int j=6; j<dim; j++){
279 mat(i,j) = mat1(i,j);
280 }
281 }
282 for (unsigned int i=6; i<dim; i++){
283 for (unsigned int j=6; j<dim; j++){
284 mat(i,j) = mat(i,j) + 2. * (i==j) * model.Beta1(nf);
285 }
286 }
287 return (mat);
288 case(LO):
289 mat = AnomalousDimension_M(LO, n_u, n_d);
290 for (unsigned int i=0; i<6; i++){
291 for (unsigned int j=6; j<dim; j++){
292 mat(i,j) = AnomalousDimension_M(NLO, n_u, n_d)(i,j);
293 }
294 }
295 for (unsigned int i=6; i<dim; i++){
296 for (unsigned int j=6; j<dim; j++){
297 mat(i,j) = mat(i,j) + 2. * (i==j) * model.Beta0(nf);
298 }
299 }
300 return (mat);
301 default:
302 throw std::runtime_error("change to rescaled operator basis: order not implemented");
303 }
304
305}
306
307gslpp::matrix<double> EvolDB1Mll::ToEffectiveBasis(gslpp::matrix<double> mat) const
308{
309
310 gslpp::matrix<double> y(dim, 0.);
311
312 y(0,0) = 1.;
313 y(1,1) = 1.;
314 y(2,2) = 1.;
315 y(3,3) = 1.;
316 y(4,4) = 1.;
317 y(5,5) = 1.;
318 y(6,6) = 1.;
319 y(7,7) = 1.;
320 y(8,8) = 1.;
321 y(9,9) = 1.;
322 y(10,10) = 1.;
323 y(11,11) = 1.;
324 y(12,12) = 1.;
325
326 y(6,2) = -1./3.;
327 y(6,3) = -4./9.;
328 y(6,4) = -20./3.;
329 y(6,5) = -80./9.;
330
331 y(7,2) = 1.;
332 y(7,3) = -1./6.;
333 y(7,4) = 20.;
334 y(7,5) = -10./3.;
335
336 y(8,2) = 4./3.;
337 y(8,4) = 64./9.;
338 y(8,5) = 64./27.; // Add terms proportional to Log(mb/mub))
339
340 return( (y.inverse()).transpose() * mat * y.transpose() );
341
342}
343
344gslpp::matrix<double>& EvolDB1Mll::Df1EvolMll(double mu, double M, orders order, schemes scheme)
345{
346
347 switch (scheme) {
348 case NDR:
349 break;
350 case LRI:
351 case HV:
352 default:
353 std::stringstream out;
354 out << scheme;
355 throw std::runtime_error("EvolDF1bsg::Df1Evolbsg(): scheme " + out.str() + " not implemented ");
356 }
357
358 double alsMZ = model.getAlsMz();
359 double Mz = model.getMz();
360 if (alsMZ == alsMZ_cache && Mz == Mz_cache) {
361 if (mu == this->mu && M == this->M && scheme == this->scheme)
362 return (*Evol(order));
363 }
364 alsMZ_cache = alsMZ;
365 Mz_cache = Mz;
366
367 if (M < mu) {
368 std::stringstream out;
369 out << "M = " << M << " < mu = " << mu;
370 throw out.str();
371 }
372
373 setScales(mu, M); // also assign evol to identity
374
375 if (M != mu) {
376 double m_down = mu;
377 double m_up = model.AboveTh(m_down);
378 double nf = model.Nf(m_down);
379
380 while (m_up < M) {
381 Df1EvolMll(m_down, m_up, nf, scheme);
382 m_down = m_up;
383 m_up = model.AboveTh(m_down);
384 nf += 1.;
385 }
386 Df1EvolMll(m_down, M, nf, scheme);
387 }
388
389 return (*Evol(order));
390
391}
392
393 void EvolDB1Mll::Df1EvolMll(double mu, double M, double nf, schemes scheme)
394 {
395
396 gslpp::matrix<double> resLO(dim, 0.), resNLO(dim, 0.), resNNLO(dim, 0.);
397
398 int L = 6 - (int) nf;
399 double alsM = model.Als(M) / 4. / M_PI;
400 double alsmu = model.Als(mu) / 4. / M_PI;
401
402 double eta = alsM / alsmu;
403
404 for (unsigned int k = 0; k < dim; k++) {
405 double etap = pow(eta, a[L][k] / 2. / model.Beta0(nf));
406 for (unsigned int i = 0; i < dim; i++){
407 for (unsigned int j = 0; j < dim; j++) {
408 resNNLO(i, j) += 0.;
409
410 if(fabs(e(i).real() - e(j).real() + 2. * model.Beta0(nf))>0.000000000001) {
411 resNLO(i, j) += c[L][i][j][k] * etap * alsmu;
412 resNLO(i, j) += d[L][i][j][k] * etap * alsM;
413 }
414 else{
415 resNLO(i, j) += - c[L][i][j][k] * etap * alsmu * log(eta);
416 }
417 resLO(i, j) += b[L][i][j][k] * etap;
418 if (fabs(resLO(i, j)) < 1.e-12) {resLO(i, j) = 0.;}
419 if (fabs(resNLO(i, j)) < 1.e-12) {resNLO(i, j) = 0.;}
420 }
421 }
422 }
423
424 switch(order) {
425 case NNLO:
426 *elem[NNLO] = 0.;
427 case NLO:
428 *elem[NLO] = (*elem[LO]) * resNLO + (*elem[NLO]) * resLO;
429 case LO:
430 *elem[LO] = (*elem[LO]) * resLO;
431 break;
432 case FULLNNLO:
433 case FULLNLO:
434 default:
435 throw std::runtime_error("Error in EvolDF1bsg::Df1Evolbsg()");
436 }
437
438 }
439
440
@ NNLO
Definition: OrderScheme.h:36
@ LO
Definition: OrderScheme.h:34
@ NLO
Definition: OrderScheme.h:35
@ FULLNNLO
Definition: OrderScheme.h:39
@ FULLNLO
Definition: OrderScheme.h:38
@ HV
Definition: OrderScheme.h:22
@ LRI
Definition: OrderScheme.h:23
@ NDR
Definition: OrderScheme.h:21
gslpp::vector< gslpp::complex > e
Definition: EvolDB1Mll.h:93
gslpp::matrix< gslpp::complex > js
Definition: EvolDB1Mll.h:92
double Mz_cache
Definition: EvolDB1Mll.h:96
double b[4][13][13][13]
Definition: EvolDB1Mll.h:82
double c[4][13][13][13]
Definition: EvolDB1Mll.h:82
gslpp::matrix< gslpp::complex > s_s
Definition: EvolDB1Mll.h:92
gslpp::matrix< gslpp::complex > jss
Definition: EvolDB1Mll.h:92
gslpp::matrix< gslpp::complex > vi
Definition: EvolDB1Mll.h:92
double alsMZ_cache
Definition: EvolDB1Mll.h:95
double d[4][13][13][13]
Definition: EvolDB1Mll.h:82
gslpp::matrix< gslpp::complex > jv
Definition: EvolDB1Mll.h:92
virtual ~EvolDB1Mll()
EvolDF1bsg destructor.
Definition: EvolDB1Mll.cpp:80
gslpp::matrix< gslpp::complex > vij
Definition: EvolDB1Mll.h:92
gslpp::matrix< gslpp::complex > gg
Definition: EvolDB1Mll.h:92
gslpp::matrix< gslpp::complex > h
Definition: EvolDB1Mll.h:92
double a[4][13]
Definition: EvolDB1Mll.h:82
gslpp::matrix< gslpp::complex > v
Definition: EvolDB1Mll.h:92
gslpp::matrix< gslpp::complex > jssv
Definition: EvolDB1Mll.h:92
gslpp::matrix< double > ToEffectiveBasis(gslpp::matrix< double > mat) const
a method returning the anomalous dimension for the evolution of the effective Wilson coefficients
Definition: EvolDB1Mll.cpp:307
gslpp::matrix< double > ToRescaleBasis(orders order, unsigned int n_u, unsigned int n_d) const
a method returning the anomalous dimension in the Chetyrkin, Misiak and Munz operator basis
Definition: EvolDB1Mll.cpp:228
gslpp::matrix< double > & Df1EvolMll(double mu, double M, orders order, schemes scheme=NDR)
a method returning the evolutor related to the high scale and the low scale
Definition: EvolDB1Mll.cpp:344
unsigned int dim
Definition: EvolDB1Mll.h:94
EvolDB1Mll(unsigned int dim, schemes scheme, orders order, const StandardModel &model)
EvolDF1bsg constructor.
Definition: EvolDB1Mll.cpp:12
gslpp::matrix< double > AnomalousDimension_M(orders order, unsigned int n_u, unsigned int n_d) const
a method returning the anomalous dimension matrix given in the Misiak basis
Definition: EvolDB1Mll.cpp:83
const StandardModel & model
Definition: EvolDB1Mll.h:83
const double Beta1(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:606
const double Beta0(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:601
const double AboveTh(const double mu) const
The active flavour threshold above the scale as defined in QCD::Thresholds().
Definition: QCD.cpp:547
const double Nf(const double mu) const
The number of active flavour at scale .
Definition: QCD.cpp:571
A class for the RG evolutor of the Wilson coefficients.
Definition: RGEvolutor.h:24
double M
Definition: RGEvolutor.h:142
void setScales(double mu, double M)
Sets the upper and lower scale for the running of the Wilson Coefficients.
Definition: RGEvolutor.cpp:85
gslpp::matrix< double > * Evol(orders order)
Evolution matrix set at a fixed order of QCD coupling.
Definition: RGEvolutor.cpp:103
A model class for the Standard Model.
const double getMz() const
A get method to access the mass of the boson .
const double getAlsMz() const
A get method to access the value of .
const double Als(const double mu, const orders order, const bool Nf_thr, const bool qed_flag) const
The running QCD coupling in the scheme including QED corrections.
gslpp::matrix< double > * elem[MAXORDER_QED+1]
orders
An enum type for orders in QCD.
Definition: OrderScheme.h:33
schemes
An enum type for regularization schemes.
Definition: OrderScheme.h:20