a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
QCD.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 <iostream>
9 #include <stdlib.h>
10 #include <sstream>
11 #include <math.h>
12 #include <map>
13 #include <stdexcept>
14 #include <gsl/gsl_errno.h>
15 #include <gsl/gsl_math.h>
16 #include <gsl/gsl_roots.h>
17 #include <TF1.h>
18 #include <Math/WrappedTF1.h>
19 #include <Math/BrentRootFinder.h>
20 #include <algorithm>
21 #include "QCD.h"
22 #include "gslpp_special_functions.h"
23 
24 std::string QCD::QCDvars[NQCDvars] = {
25  "AlsM", "MAls",
26  "mup", "mdown", "mcharm", "mstrange", "mtop", "mbottom",
27  "muc", "mub", "mut"
28 };
29 
30 QCD::QCD()
31 {
32  FlagCsi = true;
33  computeFBd = false;
34  computeFBp = false;
35  computeBd = false;
36  computeBs = false;
37  Nc = 3.;
38  TF = 0.5;
39  CF = Nc / 2. - 1. / (2. * Nc);
40  CA = Nc;
41  dFdA_NA = Nc*(Nc*Nc+6.)/48.;
42  dAdA_NA = Nc*Nc*(Nc*Nc+36.)/24.;
43  dFdF_NA = (Nc*Nc-6.+18./Nc/Nc)/96.;
44  NA = Nc*Nc-1.;
45 
46  // Particle(std::string name, double mass, double mass_scale = 0., double width = 0., double charge = 0.,double isospin = 0.);
47  quarks[UP] = Particle("UP", 0., 2., 0., 2. / 3., .5);
48  quarks[CHARM] = Particle("CHARM", 0., 0., 0., 2. / 3., .5);
49  quarks[TOP] = Particle("TOP", 0., 0., 0., 2. / 3., .5);
50  quarks[DOWN] = Particle("DOWN", 0., 2., 0., -1. / 3., -.5);
51  quarks[STRANGE] = Particle("STRANGE", 0., 2., 0., -1. / 3., -.5);
52  quarks[BOTTOM] = Particle("BOTTOM", 0., 0., 0., -1. / 3., -.5);
53 
54  zeta2 = gslpp_special_functions::zeta(2);
55  zeta3 = gslpp_special_functions::zeta(3);
56  for (int i = 0; i < CacheSize; i++) {
57  for (int j = 0; j < 8; j++)
58  als_cache[j][i] = 0.;
59  for (int j = 0; j < 4; j++)
60  logLambda5_cache[j][i] = 0.;
61  for (int j = 0; j < 10; j++)
62  mrun_cache[j][i] = 0.;
63  for (int j = 0; j < 5; j++)
64  mp2mbar_cache[j][i] = 0.;
65  }
66 
67  ModelParamMap.insert(std::make_pair("AlsM", std::cref(AlsM)));
68  ModelParamMap.insert(std::make_pair("MAls", std::cref(MAls)));
69  ModelParamMap.insert(std::make_pair("mup", std::cref(quarks[UP].getMass())));
70  ModelParamMap.insert(std::make_pair("mdown", std::cref(quarks[DOWN].getMass())));
71  ModelParamMap.insert(std::make_pair("mcharm", std::cref(quarks[CHARM].getMass())));
72  ModelParamMap.insert(std::make_pair("mstrange", std::cref(quarks[STRANGE].getMass())));
73  ModelParamMap.insert(std::make_pair("mtop", std::cref(mtpole)));
74  ModelParamMap.insert(std::make_pair("mbottom", std::cref(quarks[BOTTOM].getMass())));
75  ModelParamMap.insert(std::make_pair("muc", std::cref(muc)));
76  ModelParamMap.insert(std::make_pair("mub", std::cref(mub)));
77  ModelParamMap.insert(std::make_pair("mut", std::cref(mut)));
78 
79  unknownParameterWarning = true;
80  realorder = LO;
81 }
82 
83 std::string QCD::orderToString(const orders order) const
84 {
85  switch (order) {
86  case LO:
87  return "LO";
88  case NLO:
89  return "NLO";
90  case FULLNLO:
91  return "FULLNLO";
92  case NNLO:
93  return "NNLO";
94  case FULLNNLO:
95  return "FULLNNLO";
96  case NNNLO:
97  return "NNNLO";
98  case FULLNNNLO:
99  return "FULLNNNLO";
100  default:
101  throw std::runtime_error("QCD::orderToString(): order not implemented.");
102  }
103 }
104 
106 
107 bool QCD::Init(const std::map<std::string, double>& DPars)
108 {
109  bool check = CheckParameters(DPars);
110  if (!check) return (check);
111  check *= Update(DPars);
112  unknownParameterWarning = false;
113  return (check);
114 }
115 
116 bool QCD::PreUpdate()
117 {
118  requireYu = false;
119  requireYd = false;
120  computemt = false;
121 
122  return (true);
123 }
124 
125 bool QCD::Update(const std::map<std::string, double>& DPars)
126 {
127  if (!PreUpdate()) return (false);
128 
129  UpdateError = false;
130 
131  for (std::map<std::string, double>::const_iterator it = DPars.begin(); it != DPars.end(); it++)
132  {
133  setParameter(it->first, it->second);
134  }
135 
136  if (UpdateError) return (false);
137 
138  if (!PostUpdate()) return (false);
139 
140  return (true);
141 }
142 
143 bool QCD::PostUpdate()
144 {
145  if (computeFBd) mesonsMap.at(QCD::B_D).setDecayconst(mesonsMap.at(QCD::B_S).getDecayconst() / mesonsMap.at(QCD::B_D).getFBsoFBd());
146  if (computeFBp) mesonsMap.at(QCD::B_P).setDecayconst(mesonsMap.at(QCD::B_S).getDecayconst() / mesonsMap.at(QCD::B_P).getFBsoFBd()); /**** FOR NOW FB+ = FBd ****/
147  if (computeBs && FlagCsi)
148  {
149  BParameterMap.at("BBs").setBpars(0, BParameterMap.at("BBs").getFBsSqrtBBs1() * BParameterMap.at("BBs").getFBsSqrtBBs1() / mesonsMap.at(QCD::B_S).getDecayconst() / mesonsMap.at(QCD::B_S).getDecayconst());
150  BParameterMap.at("BBs").setBpars(1, BParameterMap.at("BBs").getFBsSqrtBBs2() * BParameterMap.at("BBs").getFBsSqrtBBs2() / mesonsMap.at(QCD::B_S).getDecayconst() / mesonsMap.at(QCD::B_S).getDecayconst());
151  BParameterMap.at("BBs").setBpars(2, BParameterMap.at("BBs").getFBsSqrtBBs3() * BParameterMap.at("BBs").getFBsSqrtBBs3() / mesonsMap.at(QCD::B_S).getDecayconst() / mesonsMap.at(QCD::B_S).getDecayconst());
152  BParameterMap.at("BBs").setBpars(3, BParameterMap.at("BBs").getFBsSqrtBBs4() * BParameterMap.at("BBs").getFBsSqrtBBs4() / mesonsMap.at(QCD::B_S).getDecayconst() / mesonsMap.at(QCD::B_S).getDecayconst());
153  BParameterMap.at("BBs").setBpars(4, BParameterMap.at("BBs").getFBsSqrtBBs5() * BParameterMap.at("BBs").getFBsSqrtBBs5() / mesonsMap.at(QCD::B_S).getDecayconst() / mesonsMap.at(QCD::B_S).getDecayconst());
154  }
155  if (computeBd) {
156  if (FlagCsi) {
157  BParameterMap.at("BBd").setBpars(0, mesonsMap.at(QCD::B_D).getFBsoFBd() * mesonsMap.at(QCD::B_D).getFBsoFBd() * BParameterMap.at("BBs").getBpars()(0) / BParameterMap.at("BBd").getcsi() / BParameterMap.at("BBd").getcsi());
158  BParameterMap.at("BBd").setBBsoBBd(BParameterMap.at("BBs").getBpars()(0) / BParameterMap.at("BBd").getBpars()(0));
159  BParameterMap.at("BBd").setBpars(1, BParameterMap.at("BBd").getFBdSqrtBBd2() * BParameterMap.at("BBd").getFBdSqrtBBd2() / mesonsMap.at(QCD::B_D).getDecayconst() / mesonsMap.at(QCD::B_D).getDecayconst());
160  BParameterMap.at("BBd").setBpars(2, BParameterMap.at("BBd").getFBdSqrtBBd3() * BParameterMap.at("BBd").getFBdSqrtBBd3() / mesonsMap.at(QCD::B_D).getDecayconst() / mesonsMap.at(QCD::B_D).getDecayconst());
161  BParameterMap.at("BBd").setBpars(3, BParameterMap.at("BBd").getFBdSqrtBBd4() * BParameterMap.at("BBd").getFBdSqrtBBd4() / mesonsMap.at(QCD::B_D).getDecayconst() / mesonsMap.at(QCD::B_D).getDecayconst());
162  BParameterMap.at("BBd").setBpars(4, BParameterMap.at("BBd").getFBdSqrtBBd5() * BParameterMap.at("BBd").getFBdSqrtBBd5() / mesonsMap.at(QCD::B_D).getDecayconst() / mesonsMap.at(QCD::B_D).getDecayconst());
163  } else
164  BParameterMap.at("BBd").setBpars(0, BParameterMap.at("BBs").getBpars()(0) / BParameterMap.at("BBd").getBBsoBBd());
165  }
166  if (computemt) {
167  quarks[TOP].setMass(Mp2Mbar(mtpole, FULLNNLO));
168 #if SUSYFIT_DEBUG & 2
169 // std::cout << "WARNING: using NLO mt for debugging purpose!"<< std::endl;
170 // quarks[TOP].setMass(Mp2Mbar(mtpole, FULLNLO));
171  mut = quarks[TOP].getMass();
172  std::cout << "postupdate: " << mtpole << std::endl;
173  std::cout << "postupdate: " << Mp2Mbar(mtpole, FULLNNLO) << std::endl;
174  std::cout << "mut set to " << mut << std::endl;
175 #endif
176  quarks[TOP].setMass_scale(quarks[TOP].getMass());
177  }
178 
179  return (true);
180 }
181 
182 void QCD::addParameters(std::vector<std::string> params_i)
183 {
184  for (std::vector<std::string>::iterator it = params_i.begin(); it < params_i.end(); it++) {
185  if (optionalParameters.find(*it) == optionalParameters.end()){
186  optionalParameters[*it] = 0.;
187  ModelParamMap.insert(std::make_pair(*it, std::cref(optionalParameters[*it])));
188  }
189  }
190 }
191 
192 void QCD::initializeBParameter(std::string name_i) const
193 {
194  if (BParameterMap.find(name_i) != BParameterMap.end()) return;
195 
196  if (name_i.compare("BBs") == 0 || name_i.compare("BBd") == 0) {
197  BParameterMap.insert(std::pair<std::string, BParameter >("BBs", BParameter(5, "BBs")));
198  BParameterMap.at("BBs").setFlagCsi(FlagCsi);
199  BParameterMap.at("BBs").ModelParameterMapInsert(ModelParamMap);
200  computeBs = true;
201  initializeMeson(QCD::B_S);
202 
203  BParameterMap.insert(std::pair<std::string, BParameter >("BBd", BParameter(5, "BBd")));
204  BParameterMap.at("BBd").setFlagCsi(FlagCsi);
205  BParameterMap.at("BBd").ModelParameterMapInsert(ModelParamMap);
206  computeBd = true;
207  initializeMeson(QCD::B_D);
208  return;
209  }
210  if (name_i.compare("BD") == 0) {
211  BParameterMap.insert(std::pair<std::string, BParameter >(name_i, BParameter(5, name_i)));
212  BParameterMap.at(name_i).ModelParameterMapInsert(ModelParamMap);
213  initializeMeson(QCD::D_0);
214  return;
215  }
216  if (name_i.compare("BK") == 0) {
217  BParameterMap.insert(std::pair<std::string, BParameter >(name_i, BParameter(5, name_i)));
218  BParameterMap.at(name_i).ModelParameterMapInsert(ModelParamMap);
219  initializeMeson(QCD::K_0);
220  return;
221  }
222  if (name_i.compare("BKd1") == 0) {
223  BParameterMap.insert(std::pair<std::string, BParameter >(name_i, BParameter(10, name_i)));
224  BParameterMap.at(name_i).ModelParameterMapInsert(ModelParamMap);
225  initializeMeson(QCD::K_0);
226  return;
227  }
228  if (name_i.compare("BKd3") == 0) {
229  BParameterMap.insert(std::pair<std::string, BParameter >(name_i, BParameter(10, name_i)));
230  BParameterMap.at(name_i).ModelParameterMapInsert(ModelParamMap);
231  initializeMeson(QCD::K_0);
232  return;
233  }
234 }
235 
236 void QCD::initializeMeson(const QCD::meson meson_i) const
237 {
238  if (mesonsMap.find(meson_i) != mesonsMap.end()) return;
239 
240  mesonsMap.insert(std::pair<const QCD::meson, Meson>(meson_i, Meson()));
241 
242  if (meson_i == QCD::P_0) mesonsMap.at(meson_i).setName("P_0");
243  else if (meson_i == QCD::P_P) mesonsMap.at(meson_i).setName("P_P");
244  else if (meson_i == QCD::K_0) mesonsMap.at(meson_i).setName("K_0");
245  else if (meson_i == QCD::K_P) mesonsMap.at(meson_i).setName("K_P");
246  else if (meson_i == QCD::D_0) mesonsMap.at(meson_i).setName("D_0");
247  else if (meson_i == QCD::D_P) mesonsMap.at(meson_i).setName("D_P");
248  else if (meson_i == QCD::B_D) mesonsMap.at(meson_i).setName("B_D");
249  else if (meson_i == QCD::B_P) mesonsMap.at(meson_i).setName("B_P");
250  else if (meson_i == QCD::B_S) mesonsMap.at(meson_i).setName("B_S");
251  else if (meson_i == QCD::B_C) mesonsMap.at(meson_i).setName("B_C");
252  else if (meson_i == QCD::PHI) mesonsMap.at(meson_i).setName("PHI");
253  else if (meson_i == QCD::K_star) mesonsMap.at(meson_i).setName("K_star");
254  else if (meson_i == QCD::K_star_P) mesonsMap.at(meson_i).setName("K_star_P");
255  else if (meson_i == QCD::D_star_P) mesonsMap.at(meson_i).setName("D_star_P");
256  else if (meson_i == QCD::RHO) mesonsMap.at(meson_i).setName("RHO");
257  else if (meson_i == QCD::RHO_P) mesonsMap.at(meson_i).setName("RHO_P");
258  else if (meson_i == QCD::OMEGA) mesonsMap.at(meson_i).setName("OMEGA");
259  else {
260  std::stringstream out;
261  out << meson_i;
262  throw std::runtime_error("QCD::initializeMeson() meson " + out.str() + " not implemented");
263  }
264 
265  if (meson_i == QCD::B_D) computeFBd = true;
266  if (meson_i == QCD::B_P) computeFBp = true;
267 
268  if ((computeFBd || computeFBp) && (mesonsMap.find(QCD::B_S) == mesonsMap.end())) initializeMeson(QCD::B_S);
269 
270  mesonsMap.at(meson_i).ModelParameterMapInsert(ModelParamMap);
271 }
272 
273 void QCD::setParameter(const std::string name, const double& value)
274 {
275  int notABparameter = 0;
276  int notAMesonParameter = 0;
277 
278  if (name.compare("AlsM") == 0) {
279  AlsM = value;
280  computemt = true;
281  requireYu = true;
282  requireYd = true;
283  } else if (name.compare("MAls") == 0) {
284  MAls = value;
285  computemt = true;
286  requireYu = true;
287  requireYd = true;
288  } else if (name.compare("mup") == 0) {
289  if (value < MEPS) UpdateError = true;
290  quarks[UP].setMass(value);
291  requireYu = true;
292  } else if (name.compare("mdown") == 0) {
293  if (value < MEPS) UpdateError = true;
294  quarks[DOWN].setMass(value);
295  requireYd = true;
296  } else if (name.compare("mcharm") == 0) {
297  quarks[CHARM].setMass(value);
298  quarks[CHARM].setMass_scale(value);
299  requireYu = true;
300  } else if (name.compare("mstrange") == 0) {
301  if (value < MEPS) UpdateError = true;
302  quarks[STRANGE].setMass(value);
303  requireYd = true;
304  } else if (name.compare("mtop") == 0) {
305  mtpole = value;
306  requireYu = true;
307  computemt = true;
308  } else if (name.compare("mbottom") == 0) {
309  quarks[BOTTOM].setMass(value);
310  quarks[BOTTOM].setMass_scale(value);
311  requireYd = true;
312  } else if (name.compare("muc") == 0)
313  muc = value;
314  else if (name.compare("mub") == 0)
315  mub = value;
316  else if (name.compare("mut") == 0)
317  mut = value;
318  else if (optionalParameters.find(name) != optionalParameters.end())
319  setOptionalParameter(name, value);
320  else {
321  if (!BParameterMap.empty())
322  for (std::map<std::string, BParameter>::iterator it = BParameterMap.begin(); it != BParameterMap.end(); it++)
323  if(it->second.setParameter(name, value))
324  notABparameter += 1;
325  if (!mesonsMap.empty())
326  for (std::map<const QCD::meson, Meson>::iterator it = mesonsMap.begin(); it != mesonsMap.end(); it++)
327  if(it->second.setParameter(name, value))
328  notAMesonParameter += 1;
329  if (unknownParameterWarning && !isSliced && notABparameter == 0 && notAMesonParameter == 0)
330  if (std::find(unknownParameters.begin(), unknownParameters.end(), name) == unknownParameters.end()) unknownParameters.push_back(name);
331  }
332 
333 }
334 
335 bool QCD::CheckParameters(const std::map<std::string, double>& DPars)
336 {
337  for (int i = 0; i < NQCDvars; i++)
338  if (DPars.find(QCDvars[i]) == DPars.end()) {
339  std::cout << "ERROR: missing mandatory QCD parameter " << QCDvars[i] << std::endl;
340  raiseMissingModelParameterCount();
341  addMissingModelParameter(QCDvars[i]);
342  }
343  for (std::map<std::string, double>::iterator it = optionalParameters.begin(); it != optionalParameters.end(); it++) {
344  if (DPars.find(it->first) == DPars.end()) {
345  std::cout << "ERROR: missing optional parameter " << it->first << std::endl;
346  raiseMissingModelParameterCount();
347  addMissingModelParameter(it->first);
348  }
349  }
350  if (!BParameterMap.empty()) {
351  for (std::map<std::string, BParameter>::iterator it = BParameterMap.begin(); it != BParameterMap.end(); it++) {
352  std::vector<std::string> parameters = it->second.parameterList(it->first);
353  for (std::vector<std::string>::iterator it1 = parameters.begin(); it1 != parameters.end(); it1++) {
354  if (DPars.find(*it1) == DPars.end()) {
355  std::cout << "ERROR: missing parameter for " << it->first << ": " << *it1 << std::endl;
356  raiseMissingModelParameterCount();
357  addMissingModelParameter(*it1);
358  }
359  }
360  }
361  }
362  if (!mesonsMap.empty()) {
363  for (std::map<const QCD::meson, Meson>::iterator it = mesonsMap.begin(); it != mesonsMap.end(); it++) {
364  std::vector<std::string> parameters = it->second.parameterList(it->second.getName());
365  for (std::vector<std::string>::iterator it1 = parameters.begin(); it1 != parameters.end(); it1++) {
366  if (DPars.find(*it1) == DPars.end()) {
367  std::cout << "ERROR: missing parameter for " << mesonsMap.at(it->first).getName() << ": " << *it1 << std::endl;
368  raiseMissingModelParameterCount();
369  addMissingModelParameter(*it1);
370  }
371  }
372  }
373  }
374  if (getMissingModelParametersCount() > 0) return false;
375  else return true;
376 }
377 
379 
380 bool QCD::setFlag(const std::string name, const bool value)
381 {
382  bool res = false;
383  if (name.compare("FlagCsi") == 0) {
384  FlagCsi = value;
385  if (computeBs) BParameterMap.at("BBs").setFlagCsi(FlagCsi);
386  if (computeBd) BParameterMap.at("BBd").setFlagCsi(FlagCsi);
387  res = true;
388  }
389 
390  return res;
391 }
392 
393 bool QCD::setFlagStr(const std::string name, const std::string value)
394 {
395  std::cout << "WARNING: wrong name or value for ModelFlag " << name << std::endl;
396  return (false);
397 }
398 
399 bool QCD::CheckFlags() const
400 {
401  return (true);
402 }
403 
405 
406 double QCD::Thresholds(const int i) const
407 {
408  if (!(mut > mub && mub > muc))
409  throw std::runtime_error("inverted thresholds in QCD::Thresholds()!");
410 
411  switch (i) {
412  case 0: return (1.0E10);
413  case 1: return (mut);
414  case 2: return (mub);
415  case 3: return (muc);
416  default: return (0.);
417  }
418 }
419 
420 double QCD::AboveTh(const double mu) const
421 {
422  int i;
423  for (i = 4; i >= 0; i--)
424  if (mu < Thresholds(i)) return (Thresholds(i));
425 
426  throw std::runtime_error("Error in QCD::AboveTh()");
427 }
428 
429 double QCD::BelowTh(const double mu) const
430 {
431  int i;
432  for (i = 0; i < 5; i++)
433  if (mu >= Thresholds(i)) return (Thresholds(i));
434 
435  throw std::runtime_error("Error in QCD::BelowTh()");
436 }
437 
438 double QCD::Nf(const double mu) const
439 {
440  int i;
441  for (i = 1; i < 5; i++)
442  if (mu >= Thresholds(i))
443  return (7. - (double) i);
444 
445  throw std::runtime_error("Error in QCD::Nf()");
446 }
447 
448 void QCD::CacheShift(double cache[][CacheSize], int n) const
449 {
450  int i, j;
451  for (i = CacheSize - 1; i > 0; i--)
452  for (j = 0; j < n; j++)
453  cache[j][i] = cache[j][i - 1];
454 }
455 
456 void QCD::CacheShift(int cache[][CacheSize], int n) const
457 {
458  int i, j;
459  for (i = CacheSize - 1; i > 0; i--)
460  for (j = 0; j < n; j++)
461  cache[j][i] = cache[j][i - 1];
462 }
463 
465 
466 double QCD::Beta0(const double nf) const
467 {
468  return ( (11. * CA - 4. * TF * nf) / 3. );
469 }
470 
471 double QCD::Beta1(const double nf) const
472 {
473  return ( 34./3. * CA * CA - ( 20./3. * CA + 4. * CF )* TF * nf);
474 }
475 
476 double QCD::Beta2(const double nf) const
477 {
478  return ( 2857./54. * CA * CA * CA - (1415./27. * CA * CA + 205./9. * CF * CA -
479  2. * CF * CF) * TF * nf +
480  (158./27. * CA + 44./9. * CF ) * TF * TF * nf * nf );
481 }
482 
483 // Czakon, hep-ph/0411261, eq. (14)
484 double QCD::Beta3(const double nf) const
485 {
486  return ( CA * CF * TF * TF * nf * nf * (17152./243. + 448./9. * zeta3) +
487  CA * CF * CF * TF * nf * (-4204./27. + 352./9. * zeta3) +
488  424./243. * CA * TF * TF * TF * nf * nf * nf + CA * CA * CF * TF * nf *
489  (7073./243. - 656./9. * zeta3) + CA * CA * TF * TF * nf * nf * (7930./81.+
490  224./9. * zeta3) + 1232./243. * CF * TF * TF * TF * nf * nf * nf +
491  CA * CA * CA * TF * nf * (-39143./81. + 136./3. * zeta3) + CA * CA * CA *
492  CA * (150653./486. - 44./9. * zeta3) + CF * CF * TF * TF * nf * nf *
493  (1352./27. - 704./9. * zeta3) + 46. * CF * CF * CF * TF * nf +
494  nf * dFdA_NA * (512./9. - 1664./3. * zeta3) + nf * nf * dFdF_NA * (
495  -704./9. + 512./3. * zeta3) + dAdA_NA * (-80./9. + 704./3. * zeta3) );
496 }
497 
498 double QCD::AlsWithInit(const double mu, const double alsi, const double mu_i,
499  const orders order) const
500 {
501  double nf = Nf(mu_i);
502 // double nf = Nf(mu);
503 // if (nf != Nf(mu_i))
504 // throw std::runtime_error("Error in QCD::AlsWithInit().");
505 
506  double b1_b0 = Beta1(nf)/Beta0(nf);
507  double v = 1. - Beta0(nf) * alsi / 2. / M_PI * log(mu_i / mu);
508  double logv = log(v);
509 
510  switch (order) {
511  case LO:
512  return (alsi / v);
513  case NLO:
514  return (- alsi * alsi / 4. / M_PI / v / v * b1_b0 * logv );
515  case NNLO:
516  return (alsi * alsi * alsi / 4. / 4. / M_PI /M_PI / v / v / v * (
517  Beta2(nf) / Beta0(nf) * (1. - v) + b1_b0 * b1_b0 * (logv * logv -
518  logv + v - 1.)));
519  case NNNLO:
520  return (alsi * alsi * alsi * alsi / 4. / 4. / 4. / M_PI /M_PI / M_PI /
521  v / v / v / v * ( Beta3(nf) / Beta0(nf) * (1. - v * v) / 2. +
522  b1_b0 * Beta2(nf) / Beta0(nf) * ((2. * v - 3.) * logv + v * v -
523  v) + b1_b0 * b1_b0 * b1_b0 * (-logv * logv * logv + 2.5 *
524  logv * logv + 2. * (1. - v) * logv - 0.5 * (v - 1.) * (v - 1.))));
525  case FULLNLO:
526  return (AlsWithInit(mu,alsi,mu_i,LO) + AlsWithInit(mu,alsi,mu_i,NLO));
527  case FULLNNLO:
528  return(AlsWithInit(mu,alsi,mu_i,LO)+AlsWithInit(mu,alsi,mu_i,NLO)+AlsWithInit(mu,alsi,mu_i,NNLO));
529  case FULLNNNLO:
530  return(AlsWithInit(mu,alsi,mu_i,LO)+AlsWithInit(mu,alsi,mu_i,NLO)+AlsWithInit(mu,alsi,mu_i,NNLO)+AlsWithInit(mu,alsi,mu_i,NNNLO));
531  default:
532  throw std::runtime_error("QCD::AlsWithInit(): " + orderToString(order) + " is not implemented.");
533  }
534 }
535 
536 double QCD::Als4(const double mu) const
537 {
538  double v = 1. - Beta0(4.) * AlsM / 2. / M_PI * log(MAls / mu);
539  return (AlsM / v * (1. - Beta1(4.) / Beta0(4.) * AlsM / 4. / M_PI * log(v) / v));
540 }
541 
542 double QCD::AlsWithLambda(const double mu, const double logLambda,
543  const orders order) const
544 {
545  double nf = Nf(mu);
546  double L = 2. * (log(mu) - logLambda);
547 
548  // LO contribution
549  double b0 = Beta0(nf);
550  double b0L = b0*L;
551  double alsLO = 4. * M_PI / b0L;
552  if (order == LO) return alsLO;
553 
554  // NLO contribution
555  double b1 = Beta1(nf);
556  double log_L = log(L);
557  double alsNLO = 4. * M_PI / b0L * (-b1 * log_L / b0 / b0L);
558  if (order == NLO) return alsNLO;
559  if (order == FULLNLO) return (alsLO + alsNLO);
560 
561  // NNLO contribution
562  double b2 = Beta2(nf);
563  double alsNNLO = 4. * M_PI / b0L * (1. / b0L / b0L
564  * (b1 * b1 / b0 / b0 * (log_L * log_L - log_L - 1.) + b2 / b0));
565  if (order == NNLO) return alsNNLO;
566  if (order == FULLNNLO) return (alsLO + alsNLO + alsNNLO);
567 
568  // NNNLO contribution
569  double b3 = Beta3(nf);
570  double alsNNNLO = 4. * M_PI / b0L * (-1. / b0L / b0L / b0L
571  * (b1 * b1 * b1 / b0 / b0 / b0 * (log_L * log_L * log_L - 5./2. * log_L * log_L -2. * log_L - 0.5) + 3. * b1 * b2 * log_L / b0 / b0 - 0.5 * b3 / b0));
572  if (order == NNNLO) return alsNNNLO;
573  if (order == FULLNNNLO) return (alsLO + alsNLO + alsNNLO + alsNNNLO);
574 
575  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::AlsWithLambda().");
576 }
577 
578 double QCD::AlsWithLambda(const double mu, const orders order) const
579 {
580  return AlsWithLambda(mu, logLambda(Nf(mu), order), order);
581 }
582 
583 double QCD::NfThresholdCorrections(double mu, double M, double als, int nf, orders order) const
584 {
585  double lmM = 2.*log(mu/M), res = 0.;
586  switch(order)
587  {
588  case FULLNNNLO:
589  res += als*als*als/M_PI/M_PI/M_PI*(-564731./124416. + 82043./27648.*gslpp_special_functions::zeta(3) +
590  2191./576.*lmM + 511./576.*lmM*lmM + lmM*lmM*lmM/216. + ((double) nf - 1.) * (2633./31104. - 67./576.*lmM + lmM*lmM/36.));
591  case FULLNNLO:
592  res += als*als/M_PI/M_PI*(-11./72. + 19./24.*lmM + lmM*lmM/36.);
593  case FULLNLO:
594  res += als/M_PI*lmM/6.;
595  case LO:
596  break;
597  default:
598  throw std::runtime_error("QCD::ThresholdCorrections(): order not implemented.");
599  }
600  return(res);
601 }
602 
603 orders QCD::FullOrder(orders order) const
604 {
605  switch(order)
606  {
607  case LO:
608  return(LO);
609  case NLO:
610  return(FULLNLO);
611  case NNLO:
612  return(FULLNNLO);
613  case NNNLO:
614  return(FULLNNNLO);
615  default:
616  throw std::runtime_error("QCD::FullOrder(): " + orderToString(order) + " is not implemented.");
617  }
618 }
619 
620 double QCD::MassOfNf(int nf) const
621 {
622  switch(nf)
623  {
624  case 6:
625  return(quarks[TOP].getMass());
626  case 5:
627  return(quarks[BOTTOM].getMass());
628  case 4:
629  return(quarks[CHARM].getMass());
630  case 3:
631  return(quarks[STRANGE].getMass());
632  default:
633  throw std::runtime_error("QCD::MassOfNf(): no running masses for light quarks");
634  }
635 }
636 
637 double QCD::Als(const double mu, const orders order, bool Nf_thr) const
638 {
639  switch (order)
640  {
641  case LO:
642  realorder = order;
643  return AlsByOrder(mu, LO, Nf_thr);
644  case FULLNLO:
645  realorder = order;
646  return (AlsByOrder(mu, LO, Nf_thr) + AlsByOrder(mu, NLO, Nf_thr));
647  case FULLNNLO:
648  realorder = order;
649  return (AlsByOrder(mu, LO, Nf_thr) + AlsByOrder(mu, NLO, Nf_thr) + AlsByOrder(mu, NNLO, Nf_thr));
650  case FULLNNNLO:
651  realorder = order;
652  return (AlsByOrder(mu, LO, Nf_thr) + AlsByOrder(mu, NLO, Nf_thr) + AlsByOrder(mu, NNLO, Nf_thr) + AlsByOrder(mu, NNNLO, Nf_thr));
653  default:
654  throw std::runtime_error("QCD::Als(): " + orderToString(order) + " is not implemented.");
655  }
656 }
657 
658 double QCD::AlsByOrder(const double mu, const orders order, bool Nf_thr) const
659 {
660  int i, nfAls = (int) Nf(MAls), nfmu = Nf_thr ? (int) Nf(mu) : nfAls;
661  double als, alstmp, mutmp;
662  orders fullord;
663 
664  for (i = 0; i < CacheSize; ++i)
665  if ((mu == als_cache[0][i]) && ((double) order == als_cache[1][i]) &&
666  (AlsM == als_cache[2][i]) && (MAls == als_cache[3][i]) &&
667  (mut == als_cache[4][i]) && (mub == als_cache[5][i]) &&
668  (muc == als_cache[6][i]) && (Nf_thr == (bool) als_cache[7][i]))
669  return als_cache[8][i];
670 
671  switch (order)
672  {
673  case LO:
674  case NLO:
675  case NNLO:
676  case NNNLO:
677  if (nfAls == nfmu)
678  als = AlsWithInit(mu, AlsM, MAls, order);
679  fullord = FullOrder(order);
680  if (nfAls > nfmu) {
681  mutmp = BelowTh(MAls);
682  alstmp = AlsWithInit(mutmp, AlsM, MAls, realorder);
683  alstmp *= (1. - NfThresholdCorrections(mutmp, MassOfNf(nfAls), alstmp, nfAls, fullord));
684  for (i = nfAls - 1; i > nfmu; i--) {
685  mutmp = BelowTh(mutmp - MEPS);
686  alstmp = AlsWithInit(mutmp, alstmp, AboveTh(mutmp) - MEPS, realorder);
687  alstmp *= (1. - NfThresholdCorrections(mutmp, MassOfNf(i), alstmp, i, fullord));
688  }
689  als = AlsWithInit(mu, alstmp, AboveTh(mu) - MEPS, order);
690  }
691 
692  if (nfAls < nfmu) {
693  mutmp = AboveTh(MAls) - MEPS;
694  alstmp = AlsWithInit(mutmp, AlsM, MAls, realorder);
695  alstmp *= (1. + NfThresholdCorrections(mutmp, MassOfNf(nfAls + 1), alstmp, nfAls + 1, fullord));
696  for (i = nfAls + 1; i < nfmu; i++) {
697  mutmp = AboveTh(mutmp) - MEPS;
698  alstmp = AlsWithInit(mutmp, alstmp, BelowTh(mutmp) + MEPS, realorder);
699  alstmp *= (1. + NfThresholdCorrections(mutmp, MassOfNf(i + 1), alstmp, i + 1, fullord));
700  }
701  als = AlsWithInit(mu, alstmp, BelowTh(mu) + MEPS, order);
702  }
703 
704  CacheShift(als_cache, 9);
705  als_cache[0][0] = mu;
706  als_cache[1][0] = (double) order;
707  als_cache[2][0] = AlsM;
708  als_cache[3][0] = MAls;
709  als_cache[4][0] = mut;
710  als_cache[5][0] = mub;
711  als_cache[6][0] = muc;
712  als_cache[7][0] = (int) Nf_thr;
713  als_cache[8][0] = als;
714 
715  return als;
716  default:
717  throw std::runtime_error("QCD::Als(): " + orderToString(order) + " is not implemented.");
718  }
719 }
720 
721 double QCD::AlsOLD(const double mu, const orders order) const
722 {
723  int i;
724  for (i = 0; i < CacheSize; ++i)
725  if ((mu == als_cache[0][i]) && ((double) order == als_cache[1][i]) &&
726  (AlsM == als_cache[2][i]) && (MAls == als_cache[3][i]) &&
727  (mut == als_cache[4][i]) && (mub == als_cache[5][i]) &&
728  (muc == als_cache[6][i]))
729  return als_cache[7][i];
730 
731  double nfmu = Nf(mu), nfz = Nf(MAls), mu_thre1, mu_thre2, Als_tmp, mf;
732  double als;
733 
734  switch (order) {
735  case LO:
736  case FULLNLO:
737  case NLO:
738  if (nfmu == nfz)
739  als = AlsWithInit(mu, AlsM, MAls, order);
740  else if (nfmu > nfz) {
741  if (order == NLO)
742  throw std::runtime_error("NLO is not implemented in QCD::Als(mu,order).");
743  if (nfmu == nfz + 1.) {
744  mu_thre1 = AboveTh(MAls); // mut
745  Als_tmp = AlsWithInit(mu_thre1 - MEPS, AlsM, MAls, order);
746  if (order == FULLNLO) {
747  mf = getQuarks(TOP).getMass(); //mf = mtpole;
748  Als_tmp = (1. + Als_tmp / M_PI * log(mu_thre1 / mf) / 3.) * Als_tmp;
749  }
750  als = AlsWithInit(mu, Als_tmp, mu_thre1 + MEPS, order);
751  } else
752  throw std::runtime_error("Error in QCD::Als(mu,order)");
753  } else {
754  if (order == NLO)
755  throw std::runtime_error("NLO is not implemented in QCD::Als(mu,order).");
756  if (nfmu == nfz - 1.) {
757  mu_thre1 = BelowTh(MAls); // mub
758  Als_tmp = AlsWithInit(mu_thre1 + MEPS, AlsM, MAls, order);
759  if (order == FULLNLO) {
760  mf = getQuarks(BOTTOM).getMass();
761  Als_tmp = (1. - Als_tmp / M_PI * log(mu_thre1 / mf) / 3.) * Als_tmp;
762  }
763  als = AlsWithInit(mu, Als_tmp, mu_thre1 - MEPS, order);
764  } else if (nfmu == nfz - 2.) {
765  mu_thre1 = BelowTh(MAls); // mub
766  mu_thre2 = AboveTh(mu); // muc
767  Als_tmp = Als(mu_thre1 + MEPS, order);
768  if (order == FULLNLO) {
769  mf = getQuarks(BOTTOM).getMass();
770  Als_tmp = (1. - Als_tmp / M_PI * log(mu_thre1 / mf) / 3.) * Als_tmp;
771  }
772  Als_tmp = AlsWithInit(mu_thre2, Als_tmp, mu_thre1 - MEPS, order);
773  if (order == FULLNLO) {
774  mf = getQuarks(CHARM).getMass();
775  Als_tmp = (1. - Als_tmp / M_PI * log(mu_thre2 / mf) / 3.) * Als_tmp;
776  }
777  als = AlsWithInit(mu, Als_tmp, mu_thre2 - MEPS, order);
778  } else
779  throw std::runtime_error("Error in QCD::Als(mu,order)");
780  }
781  break;
782  case NNLO:
783 // case NNNLO:
784  case FULLNNLO:
785 // case FULLNNNLO:
786  /* alpha_s(mu) computed with Lambda_QCD for Nf=nfmu */
787  als = AlsWithLambda(mu, order);
788  break;
789  default:
790  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Als(mu,order).");
791  }
792 
793  CacheShift(als_cache, 8);
794  als_cache[0][0] = mu;
795  als_cache[1][0] = (double) order;
796  als_cache[2][0] = AlsM;
797  als_cache[3][0] = MAls;
798  als_cache[4][0] = mut;
799  als_cache[5][0] = mub;
800  als_cache[6][0] = muc;
801  als_cache[7][0] = als;
802 
803  return als;
804 }
805 
806 double QCD::ZeroNf6NLO(double *logLambda6, double *logLambda5_in) const
807 {
808  return ( AlsWithLambda(mut + 1.e-10, *logLambda6, FULLNLO)
809  - AlsWithLambda(mut - 1.e-10, *logLambda5_in, FULLNLO));
810 }
811 
812 double QCD::ZeroNf5(double *logLambda5, double *order) const
813 {
814  return ( AlsWithLambda(MAls, *logLambda5, (orders) * order) - AlsM);
815 }
816 
817 double QCD::ZeroNf4NLO(double *logLambda4, double *logLambda5_in) const
818 {
819  return ( AlsWithLambda(mub - 1.e-10, *logLambda4, FULLNLO)
820  - AlsWithLambda(mub + 1.e-10, *logLambda5_in, FULLNLO));
821 }
822 
823 double QCD::ZeroNf3NLO(double *logLambda3, double *logLambda4_in) const
824 {
825  return ( AlsWithLambda(muc - 1.e-10, *logLambda3, FULLNLO)
826  - AlsWithLambda(muc + 1.e-10, *logLambda4_in, FULLNLO));
827 }
828 
829 double QCD::logLambda5(orders order) const
830 {
831  if (order == NLO) order = FULLNLO;
832  if (order == NNLO) order = FULLNNLO;
833 
834  for (int i = 0; i < CacheSize; ++i)
835  if ((AlsM == logLambda5_cache[0][i])
836  && (MAls == logLambda5_cache[1][i])
837  && ((double) order == logLambda5_cache[2][i]))
838  return logLambda5_cache[3][i];
839 
840  CacheShift(logLambda5_cache, 4);
841  logLambda5_cache[0][0] = AlsM;
842  logLambda5_cache[1][0] = MAls;
843  logLambda5_cache[2][0] = (double) order;
844 
845  if (order == LO)
846  logLambda5_cache[3][0] = log(MAls) - 2. * M_PI / Beta0(5.) / AlsM;
847  else {
848  double xmin = -4., xmax = -0.2;
849  TF1 f = TF1("f", this, &QCD::ZeroNf5, xmin, xmax, 1, "QCD", "zeroNf5");
850 
851  ROOT::Math::WrappedTF1 wf1(f);
852  double ledouble = (double) order;
853  wf1.SetParameters(&ledouble);
854 
855  ROOT::Math::BrentRootFinder brf;
856  brf.SetFunction(wf1, xmin, xmax);
857 
858  if (brf.Solve()) logLambda5_cache[3][0] = brf.Root();
859  else
860  throw std::runtime_error("Error in QCD::logLambda5()");
861  }
862  return ( logLambda5_cache[3][0]);
863 }
864 
865 double QCD::logLambdaNLO(const double nfNEW, const double nfORG,
866  const double logLambdaORG) const
867 {
868  for (int i = 0; i < CacheSize; ++i)
869  if ((AlsM == logLambdaNLO_cache[0][i])
870  && (MAls == logLambdaNLO_cache[1][i])
871  && (mut == logLambdaNLO_cache[2][i])
872  && (mub == logLambdaNLO_cache[3][i])
873  && (muc == logLambdaNLO_cache[4][i])
874  && (nfNEW == logLambdaNLO_cache[5][i])
875  && (nfORG == logLambdaNLO_cache[6][i])
876  && (logLambdaORG == logLambdaNLO_cache[7][i]))
877  return logLambdaNLO_cache[8][i];
878 
879  CacheShift(logLambdaNLO_cache, 9);
880  logLambdaNLO_cache[0][0] = AlsM;
881  logLambdaNLO_cache[1][0] = MAls;
882  logLambdaNLO_cache[2][0] = mut;
883  logLambdaNLO_cache[3][0] = mub;
884  logLambdaNLO_cache[4][0] = muc;
885  logLambdaNLO_cache[5][0] = nfNEW;
886  logLambdaNLO_cache[6][0] = nfORG;
887  logLambdaNLO_cache[7][0] = logLambdaORG;
888 
889  double xmin = -4., xmax = -0.2;
890 
891  TF1 f;
892  if (nfNEW == 6. && nfORG == 5.) {
893  f = TF1("f", this, &QCD::ZeroNf6NLO, xmin, xmax, 1, "QCD", "zeroNf6NLO");
894  } else if (nfNEW == 4. && nfORG == 5.) {
895  f = TF1("f", this, &QCD::ZeroNf4NLO, xmin, xmax, 1, "QCD", "zeroNf4NLO");
896  } else if (nfNEW == 3. && nfORG == 4.) {
897  f = TF1("f", this, &QCD::ZeroNf3NLO, xmin, xmax, 1, "QCD", "zeroNf3NLO");
898  } else
899  throw std::runtime_error("Error in QCD::logLambdaNLO()");
900 
901  ROOT::Math::WrappedTF1 wf1(f);
902  wf1.SetParameters(&logLambdaORG);
903 
904  ROOT::Math::BrentRootFinder brf;
905  brf.SetFunction(wf1, xmin, xmax);
906 
907  if (brf.Solve()) logLambdaNLO_cache[8][0] = brf.Root();
908  else
909  throw std::runtime_error("Error in QCD::logLambdaNLO()");
910 
911  return ( logLambdaNLO_cache[8][0]);
912 }
913 
914 double QCD::logLambda(const double muMatching, const double mf,
915  const double nfNEW, const double nfORG,
916  const double logLambdaORG, orders order) const
917 {
918  if (fabs(nfNEW - nfORG) != 1.)
919  throw std::runtime_error("Error in QCD::logLambda()");
920  if (order == NLO) order = FULLNLO;
921  if (order == NNLO) order = FULLNNLO;
922 
923  /* We do not use the codes below for FULLNLO, since threshold corrections
924  * can be regarded as an NNLO effect as long as setting the matching scale
925  * to be close to the mass scale of the decoupling quark. In order to use
926  * the relation als^{nf+1} = als^{nf} exactly, we use logLambdaNLO method.
927  */
928  if (order == FULLNLO)
929  return logLambdaNLO(nfNEW, nfORG, logLambdaORG);
930 
931  double logMuMatching = log(muMatching);
932  double L = 2. * (logMuMatching - logLambdaORG);
933  double rNEW = 0.0, rORG = 0.0, log_mu2_mf2 = 0.0, log_L = 0.0;
934  double C1 = 0.0, C2 = 0.0; // threshold corrections
935  double logLambdaNEW;
936 
937  // LO contribution
938  logLambdaNEW = 1. / 2. / Beta0(nfNEW)
939  *(Beta0(nfNEW) - Beta0(nfORG)) * L + logLambdaORG;
940 
941  // NLO contribution
942  if (order == FULLNLO || order == FULLNNLO) {
943  rNEW = Beta1(nfNEW) / Beta0(nfNEW);
944  rORG = Beta1(nfORG) / Beta0(nfORG);
945  log_mu2_mf2 = 2. * (logMuMatching - log(mf));
946  log_L = log(L);
947  if (nfNEW < nfORG)
948  C1 = 2. / 3. * log_mu2_mf2;
949  else
950  C1 = -2. / 3. * log_mu2_mf2;
951  logLambdaNEW += 1. / 2. / Beta0(nfNEW)
952  *((rNEW - rORG) * log_L
953  - rNEW * log(Beta0(nfNEW) / Beta0(nfORG)) - C1);
954  }
955 
956  // NNLO contribution
957  if (order == FULLNNLO) {
958  if (nfNEW == 5. && nfORG == 6.)
959  C2 = -16. * (log_mu2_mf2 * log_mu2_mf2 / 36. - 19. / 24. * log_mu2_mf2 - 7. / 24.);
960  else if (nfNEW == 6. && nfORG == 5.)
961  C2 = -16. * (log_mu2_mf2 * log_mu2_mf2 / 36. + 19. / 24. * log_mu2_mf2 + 7. / 24.);
962  else {
963  if (nfNEW < nfORG)
964  C2 = -16. * (log_mu2_mf2 * log_mu2_mf2 / 36. - 19. / 24. * log_mu2_mf2 + 11. / 72.);
965  else
966  C2 = -16. * (log_mu2_mf2 * log_mu2_mf2 / 36. + 19. / 24. * log_mu2_mf2 - 11. / 72.);
967  }
968  logLambdaNEW += 1. / 2. / Beta0(nfNEW) / Beta0(nfORG) / L
969  * (rORG * (rNEW - rORG) * log_L + rNEW * rNEW - rORG * rORG
970  - Beta2(nfNEW) / Beta0(nfNEW) + Beta2(nfORG) / Beta0(nfORG)
971  + rNEW * C1 - C1 * C1 - C2);
972  }
973 
974  return logLambdaNEW;
975 }
976 
977 double QCD::logLambda(const double nf, orders order) const
978 {
979  if (order == NLO) order = FULLNLO;
980  if (order == NNLO) order = FULLNNLO;
981 
982  double muMatching, mf, logLambdaORG, logLambdaNEW;
983  if (nf == 5.)
984  return logLambda5(order);
985  else if (nf == 6.) {
986  muMatching = Thresholds(1); // mut
987  /* matching condition from Nf=5 to Nf=6 is given in terms of the top pole mass. */
988  mf = mtpole; // top pole mass
989  return logLambda(muMatching, mf, 6., 5., logLambda5(order), order);
990  } else if (nf == 4. || nf == 3.) {
991  muMatching = Thresholds(2); // mub
992  mf = getQuarks(BOTTOM).getMass(); // m_b(m_b)
993  logLambdaORG = logLambda5(order);
994  logLambdaNEW = logLambda(muMatching, mf, 4., 5., logLambdaORG, order);
995  if (nf == 3.) {
996  muMatching = Thresholds(3); // muc
997  mf = getQuarks(CHARM).getMass(); // m_c(m_c)
998  logLambdaORG = logLambdaNEW;
999  logLambdaNEW = logLambda(muMatching, mf, 3., 4., logLambdaORG, order);
1000  }
1001  return logLambdaNEW;
1002  } else
1003  throw std::runtime_error("Error in QCD::logLambda()");
1004 }
1005 
1007 
1008 double QCD::Gamma0(const double nf) const
1009 {
1010  return ( 6. * CF);
1011 }
1012 
1013 double QCD::Gamma1(const double nf) const
1014 {
1015  return ( CF * (3. * CF + 97. / 3. * Nc - 10. / 3. * nf));
1016 }
1017 
1018 double QCD::Gamma2(const double nf) const
1019 {
1020  return ( 129. * CF * CF * CF - 129. / 2. * CF * CF * Nc + 11413. / 54. * CF * Nc * Nc
1021  + CF * CF * nf * (-46. + 48. * zeta3) + CF * Nc * nf * (-556. / 27. - 48. * zeta3)
1022  - 70. / 27. * CF * nf * nf);
1023 }
1024 
1025 double QCD::threCorrForMass(const double nf_f, const double nf_i) const
1026 {
1027  if (fabs(nf_f - nf_i) != 1.)
1028  throw std::runtime_error("Error in QCD::threCorrForMass()");
1029 
1030  double mu_threshold, mf, log_mu2_mf2;
1031  if (nf_f > nf_i) {
1032  if (nf_f == 6.) {
1033  mu_threshold = mut;
1034  mf = quarks[TOP].getMass(); // m_t(m_t)
1035  } else if (nf_f == 5.) {
1036  mu_threshold = mub;
1037  mf = quarks[BOTTOM].getMass(); // m_b(m_b)
1038  } else if (nf_f == 4.) {
1039  mu_threshold = muc;
1040  mf = quarks[CHARM].getMass(); // m_c(m_c)
1041  } else
1042  throw std::runtime_error("Error in QCD::threCorrForMass()");
1043  log_mu2_mf2 = 2. * log(mu_threshold / mf);
1044  return (1. + pow(Als(mu_threshold - MEPS, FULLNNLO) / M_PI, 2.)
1045  *(-log_mu2_mf2 * log_mu2_mf2 / 12. + 5. / 36. * log_mu2_mf2 - 89. / 432.));
1046  } else {
1047  if (nf_i == 6.) {
1048  mu_threshold = mut;
1049  mf = quarks[TOP].getMass(); // m_t(m_t)
1050  } else if (nf_i == 5.) {
1051  mu_threshold = mub;
1052  mf = quarks[BOTTOM].getMass(); // m_b(m_b)
1053  } else if (nf_i == 4.) {
1054  mu_threshold = muc;
1055  mf = quarks[CHARM].getMass(); // m_c(m_c)
1056  } else
1057  throw std::runtime_error("Error in QCD::threCorrForMass()");
1058  log_mu2_mf2 = 2. * log(mu_threshold / mf);
1059  return (1. + pow(Als(mu_threshold + MEPS, FULLNNLO) / M_PI, 2.)
1060  *(log_mu2_mf2 * log_mu2_mf2 / 12. - 5. / 36. * log_mu2_mf2 + 89. / 432.));
1061  }
1062 }
1063 
1064 double QCD::Mrun(const double mu, const double m, const orders order) const
1065 {
1066  return Mrun(mu, m, m, order);
1067 }
1068 
1069 double QCD::Mrun(const double mu_f, const double mu_i, const double m,
1070  const orders order) const
1071 {
1072  // Note: When the scale evolves across a flavour threshold, the definitions
1073  // of the outputs for "NLO" and "NNLO" become complicated.
1074 
1075  int i;
1076  for (i = 0; i < CacheSize; ++i) {
1077  if ((mu_f == mrun_cache[0][i]) && (mu_i == mrun_cache[1][i]) &&
1078  (m == mrun_cache[2][i]) && ((double) order == mrun_cache[3][i]) &&
1079  (AlsM == mrun_cache[4][i]) && (MAls == mrun_cache[5][i]) &&
1080  (mut == mrun_cache[6][i]) && (mub == mrun_cache[7][i]) &&
1081  (muc == mrun_cache[8][i]))
1082  return mrun_cache[9][i];
1083  }
1084 
1085  double nfi = Nf(mu_i), nff = Nf(mu_f);
1086  double mu_threshold, mu_threshold2, mu_threshold3, mrun;
1087  if (nff == nfi)
1088  mrun = MrunTMP(mu_f, mu_i, m, order);
1089  else if (nff > nfi) {
1090  if (order == NLO || order == NNLO)
1091  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Mrun(mu_f,mu_i,m,order)");
1092  mu_threshold = AboveTh(mu_i);
1093  mrun = MrunTMP(mu_threshold - MEPS, mu_i, m, order);
1094  if (order == FULLNNLO)
1095  mrun *= threCorrForMass(nfi + 1., nfi); // threshold corrections
1096  if (nff == nfi + 1.) {
1097  mrun = MrunTMP(mu_f, mu_threshold + MEPS, mrun, order);
1098  } else if (nff == nfi + 2.) {
1099  mu_threshold2 = BelowTh(mu_f);
1100  mrun = MrunTMP(mu_threshold2 - MEPS, mu_threshold + MEPS, mrun, order);
1101  if (order == FULLNNLO)
1102  mrun *= threCorrForMass(nff, nfi + 1.); // threshold corrections
1103  mrun = MrunTMP(mu_f, mu_threshold2 + MEPS, mrun, order);
1104  } else if (nff == nfi + 3.) {
1105  mu_threshold2 = AboveTh(mu_threshold);
1106  mrun = MrunTMP(mu_threshold2 - MEPS, mu_threshold + MEPS, mrun, order);
1107  if (order == FULLNNLO)
1108  mrun *= threCorrForMass(nfi + 2., nfi + 1.); // threshold corrections
1109  mu_threshold3 = BelowTh(mu_f);
1110  mrun = MrunTMP(mu_threshold3 - MEPS, mu_threshold2 + MEPS, mrun, order);
1111  if (order == FULLNNLO)
1112  mrun *= threCorrForMass(nff, nfi + 2.); // threshold corrections
1113  mrun = MrunTMP(mu_f, mu_threshold3 + MEPS, mrun, order);
1114  } else
1115  throw std::runtime_error("Error in QCD::Mrun(mu_f,mu_i,m,order)");
1116  } else {
1117  if (order == NLO || order == NNLO)
1118  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Mrun(mu_f,mu_i,m,order)");
1119  mu_threshold = BelowTh(mu_i);
1120  mrun = MrunTMP(mu_threshold + MEPS, mu_i, m, order);
1121  if (order == FULLNNLO)
1122  mrun *= threCorrForMass(nfi - 1., nfi); // threshold corrections
1123  if (nff == nfi - 1.)
1124  mrun = MrunTMP(mu_f, mu_threshold - MEPS, mrun, order);
1125  else if (nff == nfi - 2.) {
1126  mu_threshold2 = AboveTh(mu_f);
1127  mrun = MrunTMP(mu_threshold2 + MEPS, mu_threshold - MEPS, mrun, order);
1128  if (order == FULLNNLO)
1129  mrun *= threCorrForMass(nff, nfi - 1.); // threshold corrections
1130  mrun = MrunTMP(mu_f, mu_threshold2 - MEPS, mrun, order);
1131  } else
1132  throw std::runtime_error("Error in QCD::Mrun(mu_f,mu_i,m,order)");
1133  }
1134 
1135  if (mrun < 0.0) {
1136  std::stringstream out;
1137  out << "QCD::Mrun(): A quark mass becomes tachyonic in QCD::Mrun("
1138  << mu_f << ", " << mu_i << ", " << m << ", " << orderToString(order) << ")"
1139  << std::endl
1140  << " Als(" << mu_i << ", " << orderToString(order) << ")/(4pi)="
1141  << Als(mu_i, order) / (4. * M_PI) << std::endl
1142  << " Als(" << mu_f << ", " << orderToString(order) << ")/(4pi)="
1143  << Als(mu_f, order) / (4. * M_PI);
1144  throw std::runtime_error(out.str());
1145  }
1146 
1147  CacheShift(mrun_cache, 10);
1148  mrun_cache[0][0] = mu_f;
1149  mrun_cache[1][0] = mu_i;
1150  mrun_cache[2][0] = m;
1151  mrun_cache[3][0] = (double) order;
1152  mrun_cache[4][0] = AlsM;
1153  mrun_cache[5][0] = MAls;
1154  mrun_cache[6][0] = mut;
1155  mrun_cache[7][0] = mub;
1156  mrun_cache[8][0] = muc;
1157  mrun_cache[9][0] = mrun;
1158 
1159  return mrun;
1160 }
1161 
1162 double QCD::MrunTMP(const double mu_f, const double mu_i, const double m,
1163  const orders order) const
1164 {
1165  double nf = Nf(mu_f);
1166  if (nf != Nf(mu_i))
1167  throw std::runtime_error("Error in QCD::MrunTMP().");
1168 
1169  // alpha_s/(4pi)
1170  orders orderForAls;
1171  if (order == LO) orderForAls = LO;
1172  if (order == NLO || order == FULLNLO) orderForAls = FULLNLO;
1173  if (order == NNLO || order == FULLNNLO) orderForAls = FULLNNLO;
1174  double ai = Als(mu_i, orderForAls) / (4. * M_PI);
1175  double af = Als(mu_f, orderForAls) / (4. * M_PI);
1176 
1177  // LO contribution
1178  double b0 = Beta0(nf), g0 = Gamma0(nf);
1179  double mLO = m * pow(af / ai, g0 / (2. * b0));
1180  if (order == LO) return mLO;
1181 
1182  // NLO contribution
1183  double b1 = Beta1(nf), g1 = Gamma1(nf);
1184  double A1 = g1 / (2. * b0) - b1 * g0 / (2. * b0 * b0);
1185  double mNLO = mLO * A1 * (af - ai);
1186  if (order == NLO) return mNLO;
1187  if (order == FULLNLO) return (mLO + mNLO);
1188 
1189  // NNLO contribution
1190  double b2 = Beta2(nf), g2 = Gamma2(nf);
1191  double A2 = b1 * b1 * g0 / (2. * b0 * b0 * b0) - b2 * g0 / (2. * b0 * b0) - b1 * g1 / (2. * b0 * b0) + g2 / (2. * b0);
1192  double mNNLO = mLO * (A1 * A1 / 2. * (af - ai)*(af - ai) + A2 / 2. * (af * af - ai * ai));
1193  if (order == NNLO) return mNNLO;
1194  if (order == FULLNNLO) return (mLO + mNLO + mNNLO);
1195 
1196  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::MrunTMP()");
1197 }
1198 
1199 double QCD::Mrun4(const double mu_f, const double mu_i, const double m) const
1200 {
1201  double nf = 4.;
1202 
1203  // alpha_s/(4pi)
1204  double ai = Als4(mu_i) / (4. * M_PI);
1205  double af = Als4(mu_f) / (4. * M_PI);
1206 
1207  // LO contribution
1208  double b0 = Beta0(nf), g0 = Gamma0(nf);
1209  double mLO = m * pow(af / ai, g0 / (2. * b0));
1210 
1211  // NLO contribution
1212  double b1 = Beta1(nf), g1 = Gamma1(nf);
1213  double A1 = g1 / (2. * b0) - b1 * g0 / (2. * b0 * b0);
1214  double mNLO = mLO * A1 * (af - ai);
1215  return (mLO + mNLO);
1216 
1217 }
1218 
1220 
1221 double QCD::Mbar2Mp(const double mbar, const orders order) const
1222 {
1223  // LO contribution
1224  double MpLO = mbar;
1225  if (order == LO) return MpLO;
1226 
1227  // alpha_s(mbar)/pi
1228  orders orderForAls;
1229  if (order == NLO || order == FULLNLO) orderForAls = FULLNLO;
1230  if (order == NNLO || order == FULLNNLO) orderForAls = FULLNNLO;
1231  double a = Als(mbar + MEPS, orderForAls) / M_PI;
1232 
1233  // NLO contribution
1234  double MpNLO = mbar * 4. / 3. * a;
1235  if (order == NLO) return MpNLO;
1236  if (order == FULLNLO) return (MpLO + MpNLO);
1237 
1238  // NNLO contribution
1239  double nl, x;
1240  if (mbar < 3.)
1241  throw std::runtime_error("QCD::Mbar2Mp() can convert only top and bottom masses");
1242  else if (mbar < 50.) {
1243  // for the b quark
1244  nl = 4.;
1245  /* simply adding m_s(2 GeV) and m_c(m_c) */
1246  x = (quarks[STRANGE].getMass() + quarks[CHARM].getMass()) / mbar;
1247  } else {
1248  // for the top quark
1249  nl = 5.;
1250  /* simply adding m_s(2 GeV), m_c(m_c) and m_b(m_b) */
1251  x = (quarks[STRANGE].getMass() + quarks[CHARM].getMass()
1252  + quarks[BOTTOM].getMass()) / mbar;
1253  }
1254  double Delta = M_PI * M_PI / 8. * x - 0.597 * x * x + 0.230 * x * x*x;
1255  double MpNNLO = mbar * (307. / 32. + 2. * zeta2 + 2. / 3. * zeta2 * log(2.0) - zeta3 / 6.
1256  - nl / 3. * (zeta2 + 71. / 48.) + 4. / 3. * Delta) * a*a;
1257  if (order == NNLO) return MpNNLO;
1258  if (order == FULLNNLO) return (MpLO + MpNLO + MpNNLO);
1259 
1260  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Mbar2Mp().");
1261 }
1262 
1263 double QCD::Mp2MbarTMP(double *mu, double *params) const
1264 {
1265  double mp = params[0];
1266  orders order = (orders) params[1];
1267  return (mp - Mbar2Mp(*mu, order));
1268 }
1269 
1270 double QCD::Mp2Mbar(const double mp, const orders order) const
1271 {
1272  if (order == NLO || order == NNLO)
1273  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Mp2Mbar().");
1274 
1275  int i;
1276  double ms = quarks[STRANGE].getMass(), mc = quarks[CHARM].getMass();
1277  double alsmp = Als(mp, order);
1278  for (i = 0; i < CacheSize; ++i)
1279  if (alsmp == mp2mbar_cache[0][i] && ms == mp2mbar_cache[1][i] &&
1280  mc == mp2mbar_cache[2][i] && (double) order == mp2mbar_cache[3][i])
1281  return mp2mbar_cache[4][i];
1282 
1283  CacheShift(mp2mbar_cache, 5);
1284  mp2mbar_cache[0][0] = alsmp;
1285  mp2mbar_cache[1][0] = ms;
1286  mp2mbar_cache[2][0] = mc;
1287  mp2mbar_cache[3][0] = (double) order;
1288 
1289  TF1 f("f", this, &QCD::Mp2MbarTMP, mp / 2., 2. * mp, 2, "QCD", "mp2mbara");
1290 
1291  ROOT::Math::WrappedTF1 wf1(f);
1292  double params[2];
1293  params[0] = mp;
1294  params[1] = (double) order;
1295  wf1.SetParameters(params);
1296 
1297  ROOT::Math::BrentRootFinder brf;
1298 
1299  brf.SetFunction(wf1, .7 * mp, 1.3 * mp);
1300  if (brf.Solve())
1301  mp2mbar_cache[4][0] = brf.Root();
1302  else
1303  throw std::runtime_error("error in QCD::mp2mbar");
1304 
1305  return (mp2mbar_cache[4][0]);
1306 }
1307 
1308 double QCD::MS2DRqmass(const double MSbar) const
1309 {
1310  return (MSbar / (1. + Als(MSbar, FULLNLO) / 4. / M_PI * CF));
1311 }
1312 
1313 double QCD::MS2DRqmass(const double MSscale, const double MSbar) const
1314 {
1315  return (MSbar / (1. + Als(MSscale, FULLNLO) / 4. / M_PI * CF));
1316 }
gslpp_special_functions::zeta
double zeta(int i)
Definition: gslpp_special_functions.cpp:20
BParameter
A class for the bag parameters.
Definition: BParameter.h:151
QCD::ZeroNf6NLO
double ZeroNf6NLO(double *logLambda6, double *logLambda5_in) const
A member for calculating the difference in across the six-five flavour threshold using AlsWithLambda...
Definition: QCD.cpp:806
QCD::computeFBd
bool computeFBd
Switch for computing from .
Definition: QCD.h:941
QCD::logLambdaNLO
double logLambdaNLO(const double nfNEW, const double nfORG, const double logLambdaORG) const
used for computation of at FULLNLO.
Definition: QCD.cpp:865
QCD::BOTTOM
Definition: QCD.h:329
QCD::ZeroNf5
double ZeroNf5(double *logLambda5, double *order) const
A member for calculating the difference in using AlsWithLambda() and the input value of given as a ...
Definition: QCD.cpp:812
Meson
A class for mesons.
Definition: Meson.h:315
QCD::CacheSize
static const int CacheSize
Defines the depth of the cache.
Definition: QCD.h:945
QCD::Als
virtual double Als(const double mu, const orders order=FULLNLO, bool Nf_thr=true) const
Definition: QCD.cpp:637
QCD::mub
double mub
The threshold between five- and four-flavour theory in GeV.
Definition: QCD.h:929
QCD::MS2DRqmass
double MS2DRqmass(const double MSscale, const double MSbar) const
Converts a quark mass from the scheme to the scheme.
Definition: QCD.cpp:1313
QCD::AlsByOrder
virtual double AlsByOrder(const double mu, const orders order=FULLNLO, bool Nf_thr=true) const
Definition: QCD.cpp:658
Particle
A class for particles.
Definition: Particle.h:26
QCD::unknownParameterWarning
bool unknownParameterWarning
A flag to stop the unknown parameter warning after the first time.
Definition: QCD.h:951
QCD::RHO_P
Definition: QCD.h:352
QCD::MAls
double MAls
The mass scale in GeV at which the strong coupling measurement is provided.
Definition: QCD.h:926
QCD::setOptionalParameter
void setOptionalParameter(std::string name, double value)
A method to set the parameter value for the parameters that are specific to only one set of observabl...
Definition: QCD.h:458
QCD::B_S
Definition: QCD.h:345
QCD::Beta1
double Beta1(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:471
QCD::initializeBParameter
void initializeBParameter(std::string name_i) const
A method to initialize B Parameter and the corresponding meson.
Definition: QCD.cpp:192
QCD::mesonsMap
std::map< const QCD::meson, Meson > mesonsMap
The map of defined mesons.
Definition: QCD.h:954
QCD::CheckFlags
virtual bool CheckFlags() const
A method to check the sanity of the set of model flags.
Definition: QCD.cpp:399
QCD::Nf
double Nf(const double mu) const
The number of active flavour at scale .
Definition: QCD.cpp:438
FULLNNNLO
Definition: OrderScheme.h:39
QCD::initializeMeson
void initializeMeson(QCD::meson meson_i) const
A method to initialize a meson.
Definition: QCD.cpp:236
NNNLO
Definition: OrderScheme.h:36
QCD::computeFBp
bool computeFBp
Switch for computing from .
Definition: QCD.h:942
QCD::computemt
bool computemt
Switch for computing the mass of the top quark.
Definition: QCD.h:920
QCD::dFdA_NA
double dFdA_NA
Definition: QCD.h:933
LO
Definition: OrderScheme.h:33
QCD::UP
Definition: QCD.h:324
QCD::D_P
Definition: QCD.h:342
QCD::Mp2Mbar
double Mp2Mbar(const double mp, const orders order=FULLNNLO) const
Converts a quark pole mass to the corresponding mass .
Definition: QCD.cpp:1270
Model::addMissingModelParameter
void addMissingModelParameter(const std::string &missingParameterName)
Definition: Model.h:232
QCD::addParameters
void addParameters(std::vector< std::string > params_i)
A method to add parameters that are specific to only one set of observables.
Definition: QCD.cpp:182
QCD::CHARM
Definition: QCD.h:326
QCD::TF
double TF
Definition: QCD.h:933
QCD::als_cache
double als_cache[9][CacheSize]
Cache for .
Definition: QCD.h:946
QCD::dFdF_NA
double dFdF_NA
Definition: QCD.h:933
QCD::logLambda
double logLambda(const double nf, orders order) const
Computes with nf flavours in GeV.
Definition: QCD.cpp:977
QCD::FlagCsi
bool FlagCsi
A flag to determine whether and or (false) and (default, true) are used as inputs.
Definition: QCD.h:955
QCD::setParameter
virtual void setParameter(const std::string name, const double &value)
A method to set the value of a parameter of QCD.
Definition: QCD.cpp:273
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
QCD::P_P
Definition: QCD.h:338
QCD::orderToString
std::string orderToString(const orders order) const
Converts an object of the enum type "orders" to the corresponding string.
Definition: QCD.cpp:83
Model::UpdateError
bool UpdateError
A boolean set to false if update is successful.
Definition: Model.h:254
QCD::AlsOLD
double AlsOLD(const double mu, const orders order=FULLNLO) const
Computes the running strong coupling in the scheme. In the cases of LO, NLO and FULLNNLO,...
Definition: QCD.cpp:721
QCD::MassOfNf
double MassOfNf(int nf) const
The Mbar mass of the heaviest quark in the theory with Nf active flavour.
Definition: QCD.cpp:620
QCD::mp2mbar_cache
double mp2mbar_cache[5][CacheSize]
Cache for pole mass to msbar mass conversion.
Definition: QCD.h:950
QCD::mtpole
double mtpole
The pole mass of the top quark.
Definition: QCD.h:927
QCD::Mbar2Mp
double Mbar2Mp(const double mbar, const orders order=FULLNNLO) const
Converts the mass to the pole mass.
Definition: QCD.cpp:1221
QCD::mrun_cache
double mrun_cache[10][CacheSize]
Cache for running quark mass.
Definition: QCD.h:949
Model::ModelParamMap
std::map< std::string, std::reference_wrapper< const double > > ModelParamMap
Definition: Model.h:262
QCD::CheckParameters
virtual bool CheckParameters(const std::map< std::string, double > &DPars)
A method to check if all the mandatory parameters for QCD have been provided in model initialization.
Definition: QCD.cpp:335
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
QCD::PostUpdate
virtual bool PostUpdate()
The post-update method for QCD.
Definition: QCD.cpp:143
QCD::Init
virtual bool Init(const std::map< std::string, double > &DPars)
Initializes the QCD parameters found in the argument.
Definition: QCD.cpp:107
QCD::muc
double muc
The threshold between four- and three-flavour theory in GeV.
Definition: QCD.h:930
QCD::Als4
double Als4(const double mu) const
The value of at any scale with the number of flavours .
Definition: QCD.cpp:536
Particle::getMass
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
QCD::BParameterMap
std::map< std::string, BParameter > BParameterMap
Definition: QCD.h:937
QCD::Update
virtual bool Update(const std::map< std::string, double > &DPars)
The update method for QCD.
Definition: QCD.cpp:125
QCD::RHO
Definition: QCD.h:351
QCD::dAdA_NA
double dAdA_NA
Definition: QCD.h:933
QCD::K_star
Definition: QCD.h:348
QCD::meson
meson
An enum type for mesons.
Definition: QCD.h:336
QCD::B_P
Definition: QCD.h:344
QCD::TOP
Definition: QCD.h:328
QCD::OMEGA
Definition: QCD.h:353
Particle::setMass_scale
void setMass_scale(double mass_scale)
A set method to fix the scale at which the particle mass is defined.
Definition: Particle.h:142
QCD::zeta3
double zeta3
computed with the GSL.
Definition: QCD.h:940
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
QCD::Gamma1
double Gamma1(const double nf) const
The coefficient used to compute the running of a mass.
Definition: QCD.cpp:1013
Model::getMissingModelParametersCount
unsigned int getMissingModelParametersCount()
Definition: Model.h:247
Model::raiseMissingModelParameterCount
void raiseMissingModelParameterCount()
Definition: Model.h:242
Model::isSliced
bool isSliced
A boolean set to true if the current istance is a slice of an extended object.
Definition: Model.h:264
QCD::Beta3
double Beta3(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:484
QCD::Gamma0
double Gamma0(const double nf) const
The coefficient used to compute the running of a mass.
Definition: QCD.cpp:1008
QCD::zeta2
double zeta2
computed with the GSL.
Definition: QCD.h:939
QCD::threCorrForMass
double threCorrForMass(const double nf_f, const double nf_i) const
The threshold correction for running of a mass when crossing a flavour threshold.
Definition: QCD.cpp:1025
NNLO
Definition: OrderScheme.h:35
QCD::B_C
Definition: QCD.h:346
gslpp_special_functions.h
QCD::CacheShift
void CacheShift(double cache[][5], int n) const
A member used to manage the caching for this class.
QCD::getQuarks
Particle getQuarks(const QCD::quark q) const
A get method to access a quark as an object of the type Particle.
Definition: QCD.h:534
QCD::optionalParameters
std::map< std::string, double > optionalParameters
A map for containing the list and values of the parameters that are used only by a specific set of ob...
Definition: QCD.h:952
QCD::D_0
Definition: QCD.h:341
QCD::K_P
Definition: QCD.h:340
QCD::BelowTh
double BelowTh(const double mu) const
The active flavour threshold below the scale as defined in QCD::Thresholds().
Definition: QCD.cpp:429
QCD::requireYd
bool requireYd
Switch for generating the Yukawa couplings to the down-type quarks.
Definition: QCD.h:922
QCD.h
QCD::QCD
QCD()
Constructor.
Definition: QCD.cpp:30
QCD::Mp2MbarTMP
double Mp2MbarTMP(double *mu, double *params) const
The member used for finding the numerical solution to the pole mass from the mass.
Definition: QCD.cpp:1263
Particle::setMass
void setMass(double mass)
A set method to fix the particle mass.
Definition: Particle.h:70
QCD::D_star_P
Definition: QCD.h:350
QCD::QCDvars
static std::string QCDvars[NQCDvars]
An array containing the labels under which all QCD parameters are stored in a vector of ModelParamete...
Definition: QCD.h:363
orders
orders
An enum type for orders in QCD.
Definition: OrderScheme.h:31
QCD::computeBs
bool computeBs
Switch for computing from .
Definition: QCD.h:944
QCD::PreUpdate
virtual bool PreUpdate()
The pre-update method for QCD.
Definition: QCD.cpp:116
QCD::Thresholds
double Thresholds(const int i) const
For accessing the active flavour threshold scales.
Definition: QCD.cpp:406
QCD::logLambdaNLO_cache
double logLambdaNLO_cache[9][CacheSize]
Definition: QCD.h:948
QCD::Gamma2
double Gamma2(const double nf) const
The coefficient used to compute the running of a mass.
Definition: QCD.cpp:1018
QCD::P_0
Definition: QCD.h:337
QCD::K_star_P
Definition: QCD.h:349
QCD::NA
double NA
Definition: QCD.h:933
QCD::K_0
Definition: QCD.h:339
QCD::ZeroNf4NLO
double ZeroNf4NLO(double *logLambda4, double *logLambda5_in) const
A member for calculating the difference in across the four-five flavour threshold using AlsWithLambd...
Definition: QCD.cpp:817
QCD::Mrun
double Mrun(const double mu, const double m, const orders order=FULLNNLO) const
Computes a running quark mass from .
Definition: QCD.cpp:1064
QCD::STRANGE
Definition: QCD.h:327
QCD::requireYu
bool requireYu
Switch for generating the Yukawa couplings to the up-type quarks.
Definition: QCD.h:921
QCD::FullOrder
orders FullOrder(orders order) const
Return the FULLORDER enum corresponding to order.
Definition: QCD.cpp:603
QCD::PHI
Definition: QCD.h:347
QCD::computeBd
bool computeBd
Switch for computing from .
Definition: QCD.h:943
QCD::setFlagStr
virtual bool setFlagStr(const std::string name, const std::string value)
A method to set a flag of QCD.
Definition: QCD.cpp:393
QCD::B_D
Definition: QCD.h:343
Model::name
std::string name
The name of the model.
Definition: Model.h:267
QCD::AboveTh
double AboveTh(const double mu) const
The active flavour threshold above the scale as defined in QCD::Thresholds().
Definition: QCD.cpp:420
QCD::NfThresholdCorrections
double NfThresholdCorrections(double mu, double M, double als, int nf, orders order) const
Threshold corrections in matching with from eq. (34) of hep-ph/0512060.
Definition: QCD.cpp:583
QCD::Nc
double Nc
The number of colours.
Definition: QCD.h:932
QCD::ZeroNf3NLO
double ZeroNf3NLO(double *logLambda3, double *logLambda4_in) const
A member for calculating the difference in across the three-four flavour threshold using AlsWithLamb...
Definition: QCD.cpp:823
NLO
Definition: OrderScheme.h:34
QCD::AlsWithLambda
double AlsWithLambda(const double mu, const orders order) const
Computes the running strong coupling in the scheme with the use of .
Definition: QCD.cpp:578
QCD::realorder
orders realorder
Definition: QCD.h:956
QCD::CA
double CA
Definition: QCD.h:933
QCD::unknownParameters
std::vector< std::string > unknownParameters
A vector for containing the names of the parameters that are not being used but specified in the conf...
Definition: QCD.h:953
QCD::AlsWithInit
double AlsWithInit(const double mu, const double alsi, const double mu_i, const orders order) const
Computes the running strong coupling from in the scheme, where it is forbidden to across a flavour...
Definition: QCD.cpp:498
QCD::CF
double CF
Definition: QCD.h:933
FULLNNLO
Definition: OrderScheme.h:38
QCD::Mrun4
double Mrun4(const double mu_f, const double mu_i, const double m) const
The running of a mass with the number of flavours .
Definition: QCD.cpp:1199
QCD::DOWN
Definition: QCD.h:325
QCD::logLambda5_cache
double logLambda5_cache[4][CacheSize]
Definition: QCD.h:947
QCD::setFlag
virtual bool setFlag(const std::string name, const bool value)
A method to set a flag of QCD.
Definition: QCD.cpp:380
FULLNLO
Definition: OrderScheme.h:37
QCD::MrunTMP
double MrunTMP(const double mu_f, const double mu_i, const double m, const orders order) const
A function to calculate the running of the mass between flavour thresholds.
Definition: QCD.cpp:1162
QCD::NQCDvars
static const int NQCDvars
The number of model parameters in QCD.
Definition: QCD.h:357
QCD::quarks
Particle quarks[6]
The vector of all SM quarks.
Definition: QCD.h:934
QCD::logLambda5
double logLambda5(orders order) const
for .
Definition: QCD.cpp:829
QCD::mut
double mut
The threshold between six- and five-flavour theory in GeV.
Definition: QCD.h:928
QCD::AlsM
double AlsM
The strong coupling constant at the mass scale MAls, .
Definition: QCD.h:925