a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
EvolDC1Buras.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 "EvolDC1Buras.h"
10 #include "StandardModel.h"
11 #include <stdexcept>
12 
13 
14 EvolDC1Buras::EvolDC1Buras(unsigned int dim_i, schemes scheme, orders order, const StandardModel& model)
15 : RGEvolutor(dim_i, scheme, order), model(model),
16  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.),
17  jssv(dim_i,0.), jss(dim_i,0.), jv(dim_i,0.), vij(dim_i,0.), e(dim_i,0.), dim(dim_i)
18 {
19 
20  /*magic numbers a & b */
21 
22  for(int L=2; L>-1; L--){
23 
24  /* L=2 --> u,d,s,c (nf=4) L=1 --> u,d,s,c,b (nf=5) L=0 --> u,d,s,c,b,t (nf=6)*/
25 
26  nu = L; nd = L;
27  if(L == 1){nd = 3; nu = 2;}
28  if(L == 0){nd = 3; nu = 3;}
29 
30  AnomalousDimension_DC1_Buras(LO,nu,nd).transpose().eigensystem(v,e);
31  vi = v.inverse();
32  for(unsigned int i = 0; i < dim; i++){
33  a[L][i] = e(i).real();
34  for (unsigned int j = 0; j < dim; j++) {
35  for (unsigned int k = 0; k < dim; k++) {
36  b[L][i][j][k] = v(i, k).real() * vi(k, j).real();
37  }
38  }
39  }
40 
41  // LO evolutor in the standard basis
42 
43  gg = vi * AnomalousDimension_DC1_Buras(NLO,nu,nd).transpose() * v;
44  double b0 = model.Beta0(6-L);
45  double b1 = model.Beta1(6-L);
46  for (unsigned int i = 0; i < dim; i++){
47  for (unsigned int j = 0; j < dim; j++){
48  s_s.assign( i, j, (b1 / b0) * (i==j) * e(i).real() - gg(i,j));
49  if(fabs(e(i).real() - e(j).real() + 2. * b0)>0.00000000001){
50  h.assign( i, j, s_s(i,j) / (2. * b0 + e(i) - e(j)));
51  }
52  }
53  }
54  js = v * h * vi;
55  jv = js * v;
56  vij = vi * js;
57  jss = v * s_s * vi;
58  jssv = jss * v;
59  for (unsigned int i = 0; i < dim; i++){
60  for (unsigned int j = 0; j < dim; j++){
61  if(fabs(e(i).real() - e(j).real() + 2. * b0) > 0.00000000001){
62  for(unsigned int k = 0; k < dim; k++){
63  c[L][i][j][k] = jv(i, k).real() * vi(k, j).real();
64  d[L][i][j][k] = -v(i, k).real() * vij(k, j).real();
65  }
66  }
67  else{
68  for(unsigned int k = 0; k < dim; k++){
69  c[L][i][j][k] = (1./(2. * b0)) * jssv(i, k).real() * vi(k, j).real();
70  d[L][i][j][k] = 0.;
71  }
72  }
73  }
74  }
75  }
76 }
77 
78 EvolDC1Buras::~EvolDC1Buras()
79 {}
80 
81 gslpp::matrix<double> EvolDC1Buras::AnomalousDimension_DC1_Buras(orders order, unsigned int n_u, unsigned int n_d) const
82 {
83 
84  /* anomalous dimension related to Delta F = 1 operators in Buras basis, hep-ph/9512380v1 */
85 
86  /* gamma(row, column) at the LO */
87 
88  unsigned int nf = n_u + n_d; /*n_u/d = active type up/down flavor d.o.f.*/
89  gslpp::matrix<double> gammaDF1(dim, 0.);
90 
91  switch(order){
92 
93  case LO:
94 
95  gammaDF1(0,0) = -2.;
96  gammaDF1(0,1) = 6. ;
97 
98  gammaDF1(1,0) = 6.;
99  gammaDF1(1,1) = -2.;
100 
101  if(nf<5){
102  gammaDF1(1,2) = -2./9.;
103  gammaDF1(1,3) = 2./3.;
104  gammaDF1(1,4) = -2./9.;
105  gammaDF1(1,5) = 2./3.;
106 
107  }
108 
109  gammaDF1(2,2) = -22./9.;
110  gammaDF1(2,3) = 22./3.;
111  gammaDF1(2,4) = -4./9.;
112  gammaDF1(2,5) = 4./3.;
113 
114  gammaDF1(3,2) = 6.-2./9.*nf;
115  gammaDF1(3,3) = -2.+2./3.*nf;
116  gammaDF1(3,4) = -2./9.*nf;
117  gammaDF1(3,5) = 2./3.*nf;
118 
119  gammaDF1(4,4) = 2.;
120  gammaDF1(4,5) = -6.;
121 
122  gammaDF1(5,2) = -2./9.*nf;
123  gammaDF1(5,3) = 2./3.*nf;
124  gammaDF1(5,4) = -2./9.*nf;
125  gammaDF1(5,5) = -16.+2./3.*nf;
126 
127  break;
128 
129  case NLO:
130 
131  if (!(nf == 3 || nf == 4 || nf == 5 || nf == 6)){
132  throw std::runtime_error("EvolDF1nlep::AnomalousDimension_B(): wrong number of flavours");
133  }
134 
135  /* gamma(row, column) at the NLO */
136 
137  gammaDF1(0,0) = -21./2.-2./9.*nf;
138  gammaDF1(0,1) = 7./2.+2./3.*nf;
139 
140  gammaDF1(1,0) = 7./2.+2./3.*nf;
141  gammaDF1(1,1) = -21./2.-2./9.*nf;
142 
143  if(nf<5){
144 
145  gammaDF1(0,2) = 79./9.;
146  gammaDF1(0,3) = -7./3.;
147  gammaDF1(0,4) = -65./9.;
148  gammaDF1(0,5) = -7./3.;
149 
150  gammaDF1(1,2) = -202./243.;
151  gammaDF1(1,3) = 1354./81.;
152  gammaDF1(1,4) = -1192./243.;
153  gammaDF1(1,5) = 904./81.;
154 
155  }
156 
157  gammaDF1(2,2) = -5911./486.+71./9.*nf;
158  gammaDF1(2,3) = 5983./162.+1./3.*nf;
159  gammaDF1(2,4) = -2384./243.-71./9.*nf;
160  gammaDF1(2,5) = 1808./81.-1./3.*nf;
161 
162 
163  gammaDF1(3,2) = 379./18.+56./243.*nf;
164  gammaDF1(3,3) = -91./6.+808./81.*nf;
165  gammaDF1(3,4) = -130./9.-502./243.*nf;
166  gammaDF1(3,5) = -14./3.+646./81.*nf;
167 
168  gammaDF1(4,2) = -61./9.*nf;
169  gammaDF1(4,3) = -11./3.*nf;
170  gammaDF1(4,4) = 71./3.+61./9.*nf;
171  gammaDF1(4,5) = -99.+11./3.*nf;
172 
173  gammaDF1(5,2) = -682./243.*nf;
174  gammaDF1(5,3) = 106./81.*nf;
175  gammaDF1(5,4) = -225./2.+1676./243.*nf;
176  gammaDF1(5,5) = -1343./6.+1348./81.*nf;
177 
178  break;
179 
180  default:
181  throw std::runtime_error("EvolDF1nlep::AnomalousDimensio_B_S(): order not implemented");
182  }
183 
184  return (gammaDF1);
185 
186  }
187 
188 gslpp::matrix<double>& EvolDC1Buras::DC1EvolBuras(double mu, double M, orders order, schemes scheme)
189 {
190  switch (scheme) {
191  case NDR:
192  break;
193  case LRI:
194  case HV:
195  default:
196  std::stringstream out;
197  out << scheme;
198  throw std::runtime_error("EvolDC1::Df1EvolDC1(): scheme " + out.str() + " not implemented ");
199  }
200 
201  double alsMZ = model.getAlsMz();
202  double Mz = model.getMz();
203  if(alsMZ == alsMZ_cache && Mz == Mz_cache) {
204  if (mu == this-> mu && M == this->M && scheme == this->scheme)
205  return (*Evol(order));
206  }
207  alsMZ_cache = alsMZ;
208  Mz_cache = Mz;
209 
210  if (M < mu) {
211  std::stringstream out;
212  out << "M = " << M << " < mu = " << mu;
213  throw out.str();
214  }
215 
216  setScales(mu, M); // also assign evol to identity
217 
218  double m_down = mu;
219  double m_up = model.AboveTh(m_down);
220  double nf = model.Nf(m_down);
221  while (m_up < M) {
222  DC1EvolBuras(m_down, m_up, nf, scheme);
223  DC1PenguinThresholds(m_up, order);
224  m_down = m_up;
225  m_up = model.AboveTh(m_down);
226  nf += 1.;
227  }
228  DC1EvolBuras(m_down, M, nf, scheme);
229  return (*Evol(order));
230 
231  }
232 
233 void EvolDC1Buras::DC1EvolBuras(double mu, double M, double nf, schemes scheme)
234 {
235 
236  gslpp::matrix<double> resLO(dim, 0.), resNLO(dim, 0.), resNNLO(dim, 0.);
237 
238  int L = 6 - (int) nf;
239  double alsM = model.Als(M) / 4. / M_PI;
240  double alsmu = model.Als(mu) / 4. / M_PI;
241 
242  double eta = alsM / alsmu;
243 
244  for (unsigned int k = 0; k < dim; k++) {
245  double etap = pow(eta, a[L][k] / 2. / model.Beta0(nf));
246  for (unsigned int i = 0; i < dim; i++){
247  for (unsigned int j = 0; j < dim; j++) {
248  resNNLO(i, j) += 0.;
249 
250  if(fabs(e(i).real() - e(j).real() + 2. * model.Beta0(nf))>0.000000000001) {
251  resNLO(i, j) += c[L][i][j][k] * etap * alsmu;
252  resNLO(i, j) += d[L][i][j][k] * etap * alsM;
253  }
254  else{
255  resNLO(i, j) += - c[L][i][j][k] * etap * alsmu * log(eta);
256  }
257  resLO(i, j) += b[L][i][j][k] * etap;
258  }
259  }
260  }
261  switch(order) {
262  case NNLO:
263  *elem[NNLO] = 0.;
264  case NLO:
265  *elem[NLO] = (*elem[LO]) * resNLO + (*elem[NLO]) * resLO;
266  case LO:
267  *elem[LO] = (*elem[LO]) * resLO;
268  break;
269  case FULLNNLO:
270  case FULLNLO:
271  default:
272  throw std::runtime_error("Error in EvolDC1Buras::DC1EvolBuras()");
273  }
274 
275  }
276 
277 gslpp::matrix<double> EvolDC1Buras::StrongThresholds() const
278 {
279 
280 // entries of the threshold matrix for the evolution at the NLO
281 
282 gslpp::matrix<double> deltarsT(dim,0.);
283 
284 deltarsT(2,3) = 5./27.;
285 deltarsT(2,5) = 5./27.;
286 deltarsT(3,3) = -5./9.;
287 deltarsT(3,5) = -5./9.;
288 deltarsT(4,3) = 5./27.;
289 deltarsT(4,5) = 5./27.;
290 deltarsT(5,3) = -5./9.;
291 deltarsT(5,5) = -5./9.;
292 
293 return(deltarsT);
294 
295 }
296 
297 void EvolDC1Buras::DC1PenguinThresholds(double M, orders order)
298 {
299 
300  double alsM = model.Als(M) / 4. / M_PI;
301  gslpp::matrix<double> drsT(dim,0.);
302  drsT = alsM * StrongThresholds();
303  *elem[NLO] = (*elem[NLO]) + (*elem[LO]) * drsT;
304  }
EvolDC1Buras::d
double d[3][10][10][10]
Definition: EvolDC1Buras.h:74
EvolDC1Buras::e
gslpp::vector< gslpp::complex > e
Definition: EvolDC1Buras.h:91
EvolDC1Buras::s_s
gslpp::matrix< gslpp::complex > s_s
Definition: EvolDC1Buras.h:90
EvolDC1Buras.h
gslpp::matrix< double >
A class for constructing and defining operations on real matrices.
Definition: gslpp_matrix_double.h:48
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
QCD::Nf
double Nf(const double mu) const
The number of active flavour at scale .
Definition: QCD.cpp:438
StandardModel::getAlsMz
double getAlsMz() const
A get method to access the value of .
Definition: StandardModel.h:730
LO
Definition: OrderScheme.h:33
StandardModel.h
NDR
Definition: OrderScheme.h:21
gslpp::matrix< double >::transpose
matrix< double > transpose() const
Definition: gslpp_matrix_double.cpp:166
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
EvolDC1Buras::nu
int nu
Definition: EvolDC1Buras.h:67
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:477
EvolDC1Buras::jv
gslpp::matrix< gslpp::complex > jv
Definition: EvolDC1Buras.h:90
QCD::Beta0
double Beta0(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:466
EvolDC1Buras::jss
gslpp::matrix< gslpp::complex > jss
Definition: EvolDC1Buras.h:90
schemes
schemes
An enum type for regularization schemes.
Definition: OrderScheme.h:19
EvolDC1Buras::gg
gslpp::matrix< gslpp::complex > gg
Definition: EvolDC1Buras.h:90
EvolDC1Buras::js
gslpp::matrix< gslpp::complex > js
Definition: EvolDC1Buras.h:90
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
EvolDC1Buras::~EvolDC1Buras
virtual ~EvolDC1Buras()
EvolDC1Buras destructor.
Definition: EvolDC1Buras.cpp:78
EvolDC1Buras::Mz_cache
double Mz_cache
Definition: EvolDC1Buras.h:94
EvolDC1Buras::vi
gslpp::matrix< gslpp::complex > vi
Definition: EvolDC1Buras.h:90
NNLO
Definition: OrderScheme.h:35
EvolDC1Buras::alsMZ_cache
double alsMZ_cache
Definition: EvolDC1Buras.h:93
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
EvolDC1Buras::a
double a[3][10]
Definition: EvolDC1Buras.h:74
EvolDC1Buras::dim
unsigned int dim
Definition: EvolDC1Buras.h:92
EvolDC1Buras::EvolDC1Buras
EvolDC1Buras(unsigned int dim_i, schemes scheme, orders order, const StandardModel &model)
EvolDC1Buras constructor.
Definition: EvolDC1Buras.cpp:14
EvolDC1Buras::h
gslpp::matrix< gslpp::complex > h
Definition: EvolDC1Buras.h:90
EvolDC1Buras::c
double c[3][10][10][10]
Definition: EvolDC1Buras.h:74
orders
orders
An enum type for orders in QCD.
Definition: OrderScheme.h:31
LRI
Definition: OrderScheme.h:23
StandardModel::getMz
double getMz() const
A get method to access the mass of the boson .
Definition: StandardModel.h:721
EvolDC1Buras::StrongThresholds
gslpp::matrix< double > StrongThresholds() const
a method returning the matrix threshold for the QCD penguins at the NLO
Definition: EvolDC1Buras.cpp:277
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
EvolDC1Buras::AnomalousDimension_DC1_Buras
gslpp::matrix< double > AnomalousDimension_DC1_Buras(orders order, unsigned int n_u, unsigned int n_d) const
a method returning the anomalous dimension matrix given in the standard basis
Definition: EvolDC1Buras.cpp:81
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
EvolDC1Buras::b
double b[3][10][10][10]
Definition: EvolDC1Buras.h:74
EvolDC1Buras::DC1EvolBuras
gslpp::matrix< double > & DC1EvolBuras(double mu, double M, orders order, schemes scheme=NDR)
a method returning the evolutor related to the high scale and the low scale
FULLNNLO
Definition: OrderScheme.h:38
EvolDC1Buras::model
const StandardModel & model
Definition: EvolDC1Buras.h:75
gslpp::matrix< double >::eigensystem
void eigensystem(matrix< complex > &U, vector< complex > &S)
Definition: gslpp_matrix_double.cpp:280
EvolDC1Buras::vij
gslpp::matrix< gslpp::complex > vij
Definition: EvolDC1Buras.h:90
FULLNLO
Definition: OrderScheme.h:37
WilsonTemplate< gslpp::matrix< double > >::elem
gslpp::matrix< double > * elem[MAXORDER_QED+1]
Definition: WilsonTemplate.h:114
EvolDC1Buras::v
gslpp::matrix< gslpp::complex > v
Definition: EvolDC1Buras.h:90
EvolDC1Buras::nd
int nd
Definition: EvolDC1Buras.h:67
EvolDC1Buras::DC1PenguinThresholds
void DC1PenguinThresholds(double M, orders order)
a void type method for the implementation of the NLO threshold effects in the evolutor
Definition: EvolDC1Buras.cpp:297
EvolDC1Buras::jssv
gslpp::matrix< gslpp::complex > jssv
Definition: EvolDC1Buras.h:90