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
24std::string QCD::QCDvars[NQCDvars] = {
25 "AlsM", "MAls",
26 "mup", "mdown", "mcharm", "mstrange", "mtop", "mbottom",
27 "muc", "mub", "mut"};
28
29QCD::QCD()
30{
31 setModelName("QCD");
32 FlagCsi = true;
33 std::cout << "Warning: use of pole top mass as input is now deprecated. By default, the MSbar mass is used." << std::endl;
34 std::cout << "To use the pole mass, set FlagMtPole to true. (to be removed in future releases)" << std::endl;
35 FlagMtPole = false;
36 FlagMpole2MbarNumeric = true;
37 computeFBd = false;
38 computeFBp = false;
39 computeBd = false;
40 computeBs = false;
41 computemt = false;
42 requireYu = false;
43 requireYd = false;
44 Nc = 3.;
45 TF = 0.5;
46 CF = Nc / 2. - 1. / (2. * Nc);
47 CA = Nc;
48 dFdA_NA = Nc * (Nc * Nc + 6.) / 48.;
49 dAdA_NA = Nc * Nc * (Nc * Nc + 36.) / 24.;
50 dFdF_NA = (Nc * Nc - 6. + 18. / Nc / Nc) / 96.;
51 NA = Nc * Nc - 1.;
52
53 // Particle(std::string name, double mass, double mass_scale = 0., double width = 0., double charge = 0.,double isospin = 0.);
54 quarks[UP] = Particle("UP", 0., 2., 0., 2. / 3., .5);
55 quarks[CHARM] = Particle("CHARM", 0., 0., 0., 2. / 3., .5);
56 quarks[TOP] = Particle("TOP", 0., 0., 0., 2. / 3., .5);
57 quarks[DOWN] = Particle("DOWN", 0., 2., 0., -1. / 3., -.5);
58 quarks[STRANGE] = Particle("STRANGE", 0., 2., 0., -1. / 3., -.5);
59 quarks[BOTTOM] = Particle("BOTTOM", 0., 0., 0., -1. / 3., -.5);
60
61 zeta2 = gslpp_special_functions::zeta(2);
62 zeta3 = gslpp_special_functions::zeta(3);
63 for (int i = 0; i < CacheSize; i++)
64 {
65 for (int j = 0; j < 9; j++)
66 als_cache[j][i] = 0.;
67 for (int j = 0; j < 4; j++)
68 logLambda5_cache[j][i] = 0.;
69 for (int j = 0; j < 10; j++)
70 mrun_cache[j][i] = 0.;
71 for (int j = 0; j < 6; j++)
72 mp2mbar_cache[j][i] = 0.;
73 }
74
75 ModelParamMap.insert(std::make_pair("AlsM", std::cref(AlsM)));
76 ModelParamMap.insert(std::make_pair("MAls", std::cref(MAls)));
77 ModelParamMap.insert(std::make_pair("mup", std::cref(quarks[UP].getMass())));
78 ModelParamMap.insert(std::make_pair("mdown", std::cref(quarks[DOWN].getMass())));
79 ModelParamMap.insert(std::make_pair("mcharm", std::cref(quarks[CHARM].getMass())));
80 ModelParamMap.insert(std::make_pair("mstrange", std::cref(quarks[STRANGE].getMass())));
81 if (FlagMtPole)
82 ModelParamMap.insert(std::make_pair("mtop", std::cref(mtpole)));
83 else
84 ModelParamMap.insert(std::make_pair("mtop", std::cref(quarks[TOP].getMass())));
85
86 ModelParamMap.insert(std::make_pair("mbottom", std::cref(quarks[BOTTOM].getMass())));
87 ModelParamMap.insert(std::make_pair("muc", std::cref(muc)));
88 ModelParamMap.insert(std::make_pair("mub", std::cref(mub)));
89 ModelParamMap.insert(std::make_pair("mut", std::cref(mut)));
90
91 unknownParameterWarning = true;
92 realorder = LO;
93}
94
95const std::string QCD::orderToString(const orders order) const
96{
97 switch (order)
98 {
99 case LO:
100 return "LO";
101 case NLO:
102 return "NLO";
103 case FULLNLO:
104 return "FULLNLO";
105 case NNLO:
106 return "NNLO";
107 case FULLNNLO:
108 return "FULLNNLO";
109 case NNNLO:
110 return "NNNLO";
111 case FULLNNNLO:
112 return "FULLNNNLO";
113 default:
114 throw std::runtime_error("QCD::orderToString(): order not implemented.");
115 }
116}
117
119
120bool QCD::Init(const std::map<std::string, double> &DPars)
121{
122 bool check = CheckParameters(DPars);
123 if (!check)
124 return (check);
125 check *= Update(DPars) * QCDsuccess;
126 unknownParameterWarning = false;
127 return (check);
128}
129
130bool QCD::PreUpdate()
131{
132 computemt = false;
133 return (true);
134}
135
136bool QCD::Update(const std::map<std::string, double> &DPars)
137{
138 QCDsuccess = true;
139 if (!PreUpdate())
140 return (false);
141
142 UpdateError = false;
143
144 for (std::map<std::string, double>::const_iterator it = DPars.begin(); it != DPars.end(); it++)
145 {
146 setParameter(it->first, it->second);
147 }
148
149 if (UpdateError || !QCDsuccess)
150 return (false);
151
152 if (!PostUpdate() || !QCDsuccess)
153 return (false);
154
155 return (true);
156}
157
158bool QCD::PostUpdate()
159{
160 if (computeFBd)
161 mesonsMap.at(QCD::B_D).setDecayconst(mesonsMap.at(QCD::B_S).getDecayconst() / mesonsMap.at(QCD::B_D).getFBsoFBd());
162 if (computeFBp)
163 mesonsMap.at(QCD::B_P).setDecayconst(mesonsMap.at(QCD::B_S).getDecayconst() / mesonsMap.at(QCD::B_P).getFBsoFBd()); /**** FOR NOW FB+ = FBd ****/
164 if (computeBs && FlagCsi)
165 {
166 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());
167 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());
168 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());
169 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());
170 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());
171 }
172 if (computeBd)
173 {
174 if (FlagCsi)
175 {
176 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());
177 BParameterMap.at("BBd").setBBsoBBd(BParameterMap.at("BBs").getBpars()(0) / BParameterMap.at("BBd").getBpars()(0));
178 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());
179 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());
180 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());
181 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());
182 }
183 else
184 BParameterMap.at("BBd").setBpars(0, BParameterMap.at("BBs").getBpars()(0) / BParameterMap.at("BBd").getBBsoBBd());
185 }
186 if (computemt)
187 {
188 if(FlagMtPole){
189 quarks[TOP].setMass(163.); //Cannot use MSbar mass for Alphas thresholds since it is not yet updated
190 quarks[TOP].setMass(Mp2Mbar(mtpole, TOP, FULLNNLO));
191 quarks[TOP].setMass_scale(quarks[TOP].getMass());
192 }
193 else
194 mtpole = Mbar2Mp(quarks[TOP].getMass(), TOP, FULLNNLO);
195 }
196 return (QCDsuccess);
197}
198
199void QCD::addParameters(std::vector<std::string> params_i)
200{
201 for (std::vector<std::string>::iterator it = params_i.begin(); it < params_i.end(); it++)
202 {
203 if (optionalParameters.find(*it) == optionalParameters.end())
204 {
205 optionalParameters[*it] = 0.;
206 ModelParamMap.insert(std::make_pair(*it, std::cref(optionalParameters[*it])));
207 }
208 }
209}
210
211void QCD::initializeBParameter(std::string name_i) const
212{
213 if (BParameterMap.find(name_i) != BParameterMap.end())
214 return;
215
216 if (name_i.compare("BBs") == 0 || name_i.compare("BBd") == 0)
217 {
218 BParameterMap.insert(std::pair<std::string, BParameter>("BBs", BParameter(5, "BBs")));
219 BParameterMap.at("BBs").setFlagCsi(FlagCsi);
220 BParameterMap.at("BBs").ModelParameterMapInsert(ModelParamMap);
221 computeBs = true;
222 initializeMeson(QCD::B_S);
223
224 BParameterMap.insert(std::pair<std::string, BParameter>("BBd", BParameter(5, "BBd")));
225 BParameterMap.at("BBd").setFlagCsi(FlagCsi);
226 BParameterMap.at("BBd").ModelParameterMapInsert(ModelParamMap);
227 computeBd = true;
228 initializeMeson(QCD::B_D);
229 return;
230 }
231 if (name_i.compare("BBs_subleading") == 0)
232 {
233 BParameterMap.insert(std::pair<std::string, BParameter>(name_i, BParameter(2, name_i)));
234 BParameterMap.at(name_i).ModelParameterMapInsert(ModelParamMap);
235 return;
236 }
237 if (name_i.compare("BBd_subleading") == 0)
238 {
239 BParameterMap.insert(std::pair<std::string, BParameter>(name_i, BParameter(2, name_i)));
240 BParameterMap.at(name_i).ModelParameterMapInsert(ModelParamMap);
241 return;
242 }
243 if (name_i.compare("BD") == 0)
244 {
245 BParameterMap.insert(std::pair<std::string, BParameter>(name_i, BParameter(5, name_i)));
246 BParameterMap.at(name_i).ModelParameterMapInsert(ModelParamMap);
247 initializeMeson(QCD::D_0);
248 return;
249 }
250 if (name_i.compare("BK") == 0)
251 {
252 BParameterMap.insert(std::pair<std::string, BParameter>(name_i, BParameter(5, name_i)));
253 BParameterMap.at(name_i).ModelParameterMapInsert(ModelParamMap);
254 initializeMeson(QCD::K_0);
255 return;
256 }
257 if (name_i.compare("BKd1") == 0)
258 {
259 BParameterMap.insert(std::pair<std::string, BParameter>(name_i, BParameter(10, name_i)));
260 BParameterMap.at(name_i).ModelParameterMapInsert(ModelParamMap);
261 initializeMeson(QCD::K_0);
262 initializeMeson(QCD::K_P);
263 initializeMeson(QCD::K_S);
264 initializeMeson(QCD::P_0);
265 initializeMeson(QCD::P_P);
266 return;
267 }
268 if (name_i.compare("BKd3") == 0)
269 {
270 BParameterMap.insert(std::pair<std::string, BParameter>(name_i, BParameter(10, name_i)));
271 BParameterMap.at(name_i).ModelParameterMapInsert(ModelParamMap);
272 initializeMeson(QCD::K_0);
273 initializeMeson(QCD::K_P);
274 initializeMeson(QCD::P_0);
275 initializeMeson(QCD::P_P);
276 return;
277 }
278}
279
280void QCD::initializeMeson(const QCD::meson meson_i) const
281{
282 if (mesonsMap.find(meson_i) != mesonsMap.end())
283 return;
284
285 mesonsMap.insert(std::pair<const QCD::meson, Meson>(meson_i, Meson()));
286
287 if (meson_i == QCD::P_0)
288 mesonsMap.at(meson_i).setName("P_0");
289 else if (meson_i == QCD::P_P)
290 mesonsMap.at(meson_i).setName("P_P");
291 else if (meson_i == QCD::K_0)
292 mesonsMap.at(meson_i).setName("K_0");
293 else if (meson_i == QCD::K_P)
294 mesonsMap.at(meson_i).setName("K_P");
295 else if (meson_i == QCD::K_S)
296 mesonsMap.at(meson_i).setName("K_S");
297 else if (meson_i == QCD::D_0)
298 mesonsMap.at(meson_i).setName("D_0");
299 else if (meson_i == QCD::D_P)
300 mesonsMap.at(meson_i).setName("D_P");
301 else if (meson_i == QCD::D_S)
302 mesonsMap.at(meson_i).setName("D_S");
303 else if (meson_i == QCD::B_D)
304 mesonsMap.at(meson_i).setName("B_D");
305 else if (meson_i == QCD::B_P)
306 mesonsMap.at(meson_i).setName("B_P");
307 else if (meson_i == QCD::B_S)
308 mesonsMap.at(meson_i).setName("B_S");
309 else if (meson_i == QCD::B_C)
310 mesonsMap.at(meson_i).setName("B_C");
311 else if (meson_i == QCD::PHI)
312 mesonsMap.at(meson_i).setName("PHI");
313 else if (meson_i == QCD::K_star)
314 mesonsMap.at(meson_i).setName("K_star");
315 else if (meson_i == QCD::K_star_P)
316 mesonsMap.at(meson_i).setName("K_star_P");
317 else if (meson_i == QCD::D_star_P)
318 mesonsMap.at(meson_i).setName("D_star_P");
319 else if (meson_i == QCD::RHO)
320 mesonsMap.at(meson_i).setName("RHO");
321 else if (meson_i == QCD::RHO_P)
322 mesonsMap.at(meson_i).setName("RHO_P");
323 else if (meson_i == QCD::OMEGA)
324 mesonsMap.at(meson_i).setName("OMEGA");
325 else
326 {
327 std::stringstream out;
328 out << meson_i;
329 throw std::runtime_error("QCD::initializeMeson() meson " + out.str() + " not implemented");
330 }
331
332 if (meson_i == QCD::B_D)
333 computeFBd = true;
334 if (meson_i == QCD::B_P)
335 computeFBp = true;
336
337 if ((computeFBd || computeFBp) && (mesonsMap.find(QCD::B_S) == mesonsMap.end()))
338 initializeMeson(QCD::B_S);
339
340 mesonsMap.at(meson_i).ModelParameterMapInsert(ModelParamMap);
341}
342
343void QCD::setParameter(const std::string name, const double &value)
344{
345 int notABparameter = 0;
346 int notAMesonParameter = 0;
347
348 if (name.compare("AlsM") == 0)
349 {
350 AlsM = value;
351 computemt = true;
352 }
353 else if (name.compare("MAls") == 0)
354 {
355 MAls = value;
356 computemt = true;
357 }
358 else if (name.compare("mup") == 0)
359 {
360 if (value < MEPS)
361 UpdateError = true;
362 quarks[UP].setMass(value);
363 }
364 else if (name.compare("mdown") == 0)
365 {
366 if (value < MEPS)
367 UpdateError = true;
368 quarks[DOWN].setMass(value);
369 }
370 else if (name.compare("mcharm") == 0)
371 {
372 quarks[CHARM].setMass(value);
373 quarks[CHARM].setMass_scale(value);
374 }
375 else if (name.compare("mstrange") == 0)
376 {
377 if (value < MEPS)
378 UpdateError = true;
379 quarks[STRANGE].setMass(value);
380 }
381 else if (name.compare("mtop") == 0)
382 {
383 if (FlagMtPole)
384 mtpole = value;
385 else
386 {
387 quarks[TOP].setMass(value);
388 quarks[TOP].setMass_scale(value);
389 }
390 computemt = true;
391 }
392 else if (name.compare("mbottom") == 0)
393 {
394 quarks[BOTTOM].setMass(value);
395 quarks[BOTTOM].setMass_scale(value);
396 }
397 else if (name.compare("muc") == 0)
398 muc = value;
399 else if (name.compare("mub") == 0)
400 mub = value;
401 else if (name.compare("mut") == 0)
402 mut = value;
403 else if (optionalParameters.find(name) != optionalParameters.end())
404 setOptionalParameter(name, value);
405 else
406 {
407 if (!BParameterMap.empty())
408 for (std::map<std::string, BParameter>::iterator it = BParameterMap.begin(); it != BParameterMap.end(); it++)
409 if (it->second.setParameter(name, value))
410 notABparameter += 1;
411 if (!mesonsMap.empty())
412 for (std::map<const QCD::meson, Meson>::iterator it = mesonsMap.begin(); it != mesonsMap.end(); it++)
413 if (it->second.setParameter(name, value))
414 notAMesonParameter += 1;
415 if (unknownParameterWarning && !isSliced && notABparameter == 0 && notAMesonParameter == 0)
416 if (std::find(unknownParameters.begin(), unknownParameters.end(), name) == unknownParameters.end())
417 unknownParameters.push_back(name);
418 }
419}
420
421bool QCD::CheckParameters(const std::map<std::string, double> &DPars)
422{
423 for (int i = 0; i < NQCDvars; i++)
424 if (DPars.find(QCDvars[i]) == DPars.end())
425 {
426 std::cout << "ERROR: missing mandatory QCD parameter " << QCDvars[i] << std::endl;
427 raiseMissingModelParameterCount();
428 addMissingModelParameter(QCDvars[i]);
429 }
430 for (std::map<std::string, double>::iterator it = optionalParameters.begin(); it != optionalParameters.end(); it++)
431 {
432 if (DPars.find(it->first) == DPars.end())
433 {
434 std::cout << "ERROR: missing optional parameter " << it->first << std::endl;
435 raiseMissingModelParameterCount();
436 addMissingModelParameter(it->first);
437 }
438 }
439 if (!BParameterMap.empty())
440 {
441 for (std::map<std::string, BParameter>::iterator it = BParameterMap.begin(); it != BParameterMap.end(); it++)
442 {
443 std::vector<std::string> parameters = it->second.parameterList(it->first);
444 for (std::vector<std::string>::iterator it1 = parameters.begin(); it1 != parameters.end(); it1++)
445 {
446 if (DPars.find(*it1) == DPars.end())
447 {
448 std::cout << "ERROR: missing parameter for " << it->first << ": " << *it1 << std::endl;
449 raiseMissingModelParameterCount();
450 addMissingModelParameter(*it1);
451 }
452 }
453 }
454 }
455 if (!mesonsMap.empty())
456 {
457 for (std::map<const QCD::meson, Meson>::iterator it = mesonsMap.begin(); it != mesonsMap.end(); it++)
458 {
459 std::vector<std::string> parameters = it->second.parameterList(it->second.getName());
460 for (std::vector<std::string>::iterator it1 = parameters.begin(); it1 != parameters.end(); it1++)
461 {
462 if (DPars.find(*it1) == DPars.end())
463 {
464 std::cout << "ERROR: missing parameter for " << mesonsMap.at(it->first).getName() << ": " << *it1 << std::endl;
465 raiseMissingModelParameterCount();
466 addMissingModelParameter(*it1);
467 }
468 }
469 }
470 }
471 if (getMissingModelParametersCount() > 0)
472 return false;
473 else
474 return true;
475}
476
478
479bool QCD::setFlag(const std::string name, const bool value)
480{
481 bool res = false;
482 if (name.compare("FlagCsi") == 0)
483 {
484 FlagCsi = value;
485 if (computeBs)
486 BParameterMap.at("BBs").setFlagCsi(FlagCsi);
487 if (computeBd)
488 BParameterMap.at("BBd").setFlagCsi(FlagCsi);
489 res = true;
490 }
491 else if (name.compare("FlagMtPole") == 0)
492 {
493 FlagMtPole = value;
494 if (FlagMtPole) {
495 ModelParamMap.at("mtop")=std::cref(mtpole);
496 std::cout << "WARNING: use of pole top mass as input is now deprecated. This feature will be removed in future releases. " << std::endl;
497 std::cout << "Using anyway the pole mass as requested." << std::endl;
498 }
499 else
500 ModelParamMap.at("mtop")=std::cref(quarks[TOP].getMass());
501 res = true;
502 }
503 else if (name.compare("FlagMpole2MbarNumeric") == 0)
504 {
505 FlagMpole2MbarNumeric = value;
506 res = true;
507 }
508 return res;
509}
510
511bool QCD::setFlagStr(const std::string name, const std::string value)
512{
513 std::cout << "WARNING: wrong name or value for ModelFlag " << name << std::endl;
514 return (false);
515}
516
517bool QCD::CheckFlags() const
518{
519 return (true);
520}
521
523
524const double QCD::Thresholds(const int i) const
525{
526 if (!(mut > mub && mub > muc))
527 {
528 QCDsuccess = false;
529 std::cout << "inverted thresholds in QCD::Thresholds()! QCDsuccess set to false" << std::endl;
530 }
531
532 switch (i)
533 {
534 case 0:
535 return (1.0E10);
536 case 1:
537 return (mut);
538 case 2:
539 return (mub);
540 case 3:
541 return (muc);
542 default:
543 return (0.);
544 }
545}
546
547const double QCD::AboveTh(const double mu) const
548{
549 int i;
550 for (i = 4; i >= 0; i--)
551 if (mu < Thresholds(i))
552 return (Thresholds(i));
553
554 QCDsuccess = false;
555 std::cout << "Error in QCD::AboveTh()! QCDsuccess set to false" << std::endl;
556 return (0.);
557}
558
559const double QCD::BelowTh(const double mu) const
560{
561 int i;
562 for (i = 0; i < 5; i++)
563 if (mu >= Thresholds(i))
564 return (Thresholds(i));
565
566 QCDsuccess = false;
567 std::cout << "Error in QCD::BelowTh()! QCDsuccess set to false" << std::endl;
568 return (0.);
569}
570
571const double QCD::Nf(const double mu) const
572{
573 int i;
574 for (i = 1; i < 5; i++)
575 if (mu >= Thresholds(i))
576 return (7. - (double)i);
577
578 QCDsuccess = false;
579 std::cout << "Error in QCD::Nf(" << mu <<")! QCDsuccess set to false" << std::endl;
580 return (0.);
581}
582
583void QCD::CacheShift(double cache[][CacheSize], int n) const
584{
585 int i, j;
586 for (i = CacheSize - 1; i > 0; i--)
587 for (j = 0; j < n; j++)
588 cache[j][i] = cache[j][i - 1];
589}
590
591void QCD::CacheShift(int cache[][CacheSize], int n) const
592{
593 int i, j;
594 for (i = CacheSize - 1; i > 0; i--)
595 for (j = 0; j < n; j++)
596 cache[j][i] = cache[j][i - 1];
597}
598
600
601const double QCD::Beta0(const double nf) const
602{
603 return ((11. * CA - 4. * TF * nf) / 3.);
604}
605
606const double QCD::Beta1(const double nf) const
607{
608 return (34. / 3. * CA * CA - (20. / 3. * CA + 4. * CF) * TF * nf);
609}
610
611const double QCD::Beta2(const double nf) const
612{
613 return (2857. / 54. * CA * CA * CA - (1415. / 27. * CA * CA + 205. / 9. * CF * CA - 2. * CF * CF) * TF * nf +
614 (158. / 27. * CA + 44. / 9. * CF) * TF * TF * nf * nf);
615}
616
617// Czakon, hep-ph/0411261, eq. (14)
618const double QCD::Beta3(const double nf) const
619{
620 return (CA * CF * TF * TF * nf * nf * (17152. / 243. + 448. / 9. * zeta3) +
621 CA * CF * CF * TF * nf * (-4204. / 27. + 352. / 9. * zeta3) +
622 424. / 243. * CA * TF * TF * TF * nf * nf * nf + CA * CA * CF * TF * nf * (7073. / 243. - 656. / 9. * zeta3) + CA * CA * TF * TF * nf * nf * (7930. / 81. + 224. / 9. * zeta3) + 1232. / 243. * CF * TF * TF * TF * nf * nf * nf +
623 CA * CA * CA * TF * nf * (-39143. / 81. + 136. / 3. * zeta3) + CA * CA * CA * CA * (150653. / 486. - 44. / 9. * zeta3) + CF * CF * TF * TF * nf * nf * (1352. / 27. - 704. / 9. * zeta3) + 46. * CF * CF * CF * TF * nf +
624 nf * dFdA_NA * (512. / 9. - 1664. / 3. * zeta3) + nf * nf * dFdF_NA * (-704. / 9. + 512. / 3. * zeta3) + dAdA_NA * (-80. / 9. + 704. / 3. * zeta3));
625}
626
627const double QCD::AlsWithInit(const double mu, const double alsi, const double mu_i, const int nf,
628 const orders order) const
629{
630
631 double b1_b0 = Beta1((double)nf) / Beta0((double)nf);
632 double v = 1. - Beta0((double)nf) * alsi / 2. / M_PI * log(mu_i / mu);
633 double logv = log(v);
634
635 switch (order)
636 {
637 case LO:
638 return (alsi / v);
639 case NLO:
640 return (-alsi * alsi / 4. / M_PI / v / v * b1_b0 * logv);
641 case NNLO:
642 return (alsi * alsi * alsi / 4. / 4. / M_PI / M_PI / v / v / v * (Beta2((double)nf) / Beta0((double)nf) * (1. - v) + b1_b0 * b1_b0 * (logv * logv - logv + v - 1.)));
643 case NNNLO:
644 return (alsi * alsi * alsi * alsi / 4. / 4. / 4. / M_PI / M_PI / M_PI /
645 v / v / v / v * (Beta3((double)nf) / Beta0((double)nf) * (1. - v * v) / 2. + b1_b0 * Beta2((double)nf) / Beta0((double)nf) * ((2. * v - 3.) * logv + v * v - v) + b1_b0 * b1_b0 * b1_b0 * (-logv * logv * logv + 2.5 * logv * logv + 2. * (1. - v) * logv - 0.5 * (v - 1.) * (v - 1.))));
646 case FULLNLO:
647 return (AlsWithInit(mu, alsi, mu_i, nf, LO) + AlsWithInit(mu, alsi, mu_i, nf, NLO));
648 case FULLNNLO:
649 return (AlsWithInit(mu, alsi, mu_i, nf, LO) + AlsWithInit(mu, alsi, mu_i, nf, NLO) + AlsWithInit(mu, alsi, mu_i, nf, NNLO));
650 case FULLNNNLO:
651 return (AlsWithInit(mu, alsi, mu_i, nf, LO) + AlsWithInit(mu, alsi, mu_i, nf, NLO) + AlsWithInit(mu, alsi, mu_i, nf, NNLO) + AlsWithInit(mu, alsi, mu_i, nf, NNNLO));
652 default:
653 throw std::runtime_error("QCD::AlsWithInit(): " + orderToString(order) + " is not implemented.");
654 }
655}
656
657const double QCD::Als4(const double mu) const
658{
659 double v = 1. - Beta0(4.) * AlsM / 2. / M_PI * log(MAls / mu);
660 return (AlsM / v * (1. - Beta1(4.) / Beta0(4.) * AlsM / 4. / M_PI * log(v) / v));
661}
662
663const double QCD::AlsWithLambda(const double mu, const double logLambda,
664 const orders order) const
665{
666 double nf = Nf(mu);
667 double L = 2. * (log(mu) - logLambda);
668
669 // LO contribution
670 double b0 = Beta0(nf);
671 double b0L = b0 * L;
672 double alsLO = 4. * M_PI / b0L;
673 if (order == LO)
674 return alsLO;
675
676 // NLO contribution
677 double b1 = Beta1(nf);
678 double log_L = log(L);
679 double alsNLO = 4. * M_PI / b0L * (-b1 * log_L / b0 / b0L);
680 if (order == NLO)
681 return alsNLO;
682 if (order == FULLNLO)
683 return (alsLO + alsNLO);
684
685 // NNLO contribution
686 double b2 = Beta2(nf);
687 double alsNNLO = 4. * M_PI / b0L * (1. / b0L / b0L * (b1 * b1 / b0 / b0 * (log_L * log_L - log_L - 1.) + b2 / b0));
688 if (order == NNLO)
689 return alsNNLO;
690 if (order == FULLNNLO)
691 return (alsLO + alsNLO + alsNNLO);
692
693 // NNNLO contribution
694 double b3 = Beta3(nf);
695 double alsNNNLO = 4. * M_PI / b0L * (-1. / b0L / b0L / b0L * (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));
696 if (order == NNNLO)
697 return alsNNNLO;
698 if (order == FULLNNNLO)
699 return (alsLO + alsNLO + alsNNLO + alsNNNLO);
700
701 throw std::runtime_error(orderToString(order) + " is not implemented in QCD::AlsWithLambda().");
702}
703
704const double QCD::AlsWithLambda(const double mu, const orders order) const
705{
706 return AlsWithLambda(mu, logLambda(Nf(mu), order), order);
707}
708
709const double QCD::NfThresholdCorrections(double mu, double M, double als, int nf, orders order) const
710{
711 double lmM = 2. * log(mu / M), res = 0.;
712 switch (order)
713 {
714 case FULLNNNLO:
715 res += als * als * als / M_PI / M_PI / M_PI * (-564731. / 124416. + 82043. / 27648. * gslpp_special_functions::zeta(3) + 2191. / 576. * lmM + 511. / 576. * lmM * lmM + lmM * lmM * lmM / 216. + ((double)nf - 1.) * (2633. / 31104. - 281. / 1728. * lmM));
716 case FULLNNLO:
717 res += als * als / M_PI / M_PI * (-11. / 72. + 19. / 24. * lmM + lmM * lmM / 36.);
718 case FULLNLO:
719 res += als / M_PI * lmM / 6.;
720 case LO:
721 break;
722 default:
723 throw std::runtime_error("QCD::ThresholdCorrections(): order not implemented.");
724 }
725 return (res);
726}
727
728const orders QCD::FullOrder(orders order) const
729{
730 switch (order)
731 {
732 case LO:
733 return (LO);
734 case NLO:
735 return (FULLNLO);
736 case NNLO:
737 return (FULLNNLO);
738 case NNNLO:
739 return (FULLNNNLO);
740 default:
741 throw std::runtime_error("QCD::FullOrder(): " + orderToString(order) + " is not implemented.");
742 }
743}
744
745const double QCD::MassOfNf(int nf) const
746{
747 switch (nf)
748 {
749 case 6:
750 return (quarks[TOP].getMass());
751 case 5:
752 return (quarks[BOTTOM].getMass());
753 case 4:
754 return (quarks[CHARM].getMass());
755 case 3:
756 return (quarks[STRANGE].getMass());
757 default:
758 throw std::runtime_error("QCD::MassOfNf(): no running masses for light quarks");
759 }
760}
761
762const double QCD::Als(const double mu, const orders order, const bool Nf_thr) const
763{
764 switch (order)
765 {
766 case LO:
767 realorder = order;
768 return AlsByOrder(mu, LO, Nf_thr);
769 case FULLNLO:
770 realorder = order;
771 return (AlsByOrder(mu, LO, Nf_thr) + AlsByOrder(mu, NLO, Nf_thr));
772 case FULLNNLO:
773 realorder = order;
774 return (AlsByOrder(mu, LO, Nf_thr) + AlsByOrder(mu, NLO, Nf_thr) + AlsByOrder(mu, NNLO, Nf_thr));
775 case FULLNNNLO:
776 realorder = order;
777 return (AlsByOrder(mu, LO, Nf_thr) + AlsByOrder(mu, NLO, Nf_thr) + AlsByOrder(mu, NNLO, Nf_thr) + AlsByOrder(mu, NNNLO, Nf_thr));
778 default:
779 throw std::runtime_error("QCD::Als(): " + orderToString(order) + " is not implemented.");
780 }
781}
782
783const double QCD::Als(const double mu, const int Nf, const orders order) const
784{
785 switch (order)
786 {
787 case LO:
788 realorder = order;
789 return AlsByOrder(mu, Nf, LO);
790 case FULLNLO:
791 realorder = order;
792 return (AlsByOrder(mu, Nf, LO) + AlsByOrder(mu, Nf, NLO));
793 case FULLNNLO:
794 realorder = order;
795 return (AlsByOrder(mu, Nf, LO) + AlsByOrder(mu, Nf, NLO) + AlsByOrder(mu, Nf, NNLO));
796 case FULLNNNLO:
797 realorder = order;
798 return (AlsByOrder(mu, Nf, LO) + AlsByOrder(mu, Nf, NLO) + AlsByOrder(mu, Nf, NNLO) + AlsByOrder(mu, Nf, NNNLO));
799 default:
800 throw std::runtime_error("QCD::Als(): " + orderToString(order) + " is not implemented.");
801 }
802}
803
804const double QCD::AlsByOrder(const double mu, const orders order, const bool Nf_thr) const
805{
806 int i, nfAls = (int)Nf(MAls), nfmu = Nf_thr ? (int)Nf(mu) : nfAls;
807 double als, alstmp, mutmp;
808 orders fullord;
809
810 for (i = 0; i < CacheSize; ++i)
811 if ((mu == als_cache[0][i]) && ((double)order == als_cache[1][i]) &&
812 (AlsM == als_cache[2][i]) && (MAls == als_cache[3][i]) &&
813 (mut == als_cache[4][i]) && (mub == als_cache[5][i]) &&
814 (muc == als_cache[6][i]) && (Nf_thr == (bool)als_cache[7][i]))
815 return als_cache[8][i];
816
817 switch (order)
818 {
819 case LO:
820 case NLO:
821 case NNLO:
822 case NNNLO:
823 if (nfAls == nfmu)
824 als = AlsWithInit(mu, AlsM, MAls, nfmu, order);
825 fullord = FullOrder(order);
826 if (nfAls > nfmu)
827 {
828 mutmp = BelowTh(MAls);
829 alstmp = AlsWithInit(mutmp, AlsM, MAls, nfAls, realorder);
830 alstmp *= (1. - NfThresholdCorrections(mutmp, MassOfNf(nfAls), alstmp, nfAls, fullord));
831 for (i = nfAls - 1; i > nfmu; i--)
832 {
833 mutmp = BelowTh(mutmp - MEPS);
834 alstmp = AlsWithInit(mutmp, alstmp, AboveTh(mutmp) - MEPS, i, realorder);
835 alstmp *= (1. - NfThresholdCorrections(mutmp, MassOfNf(i), alstmp, i, fullord));
836 }
837 als = AlsWithInit(mu, alstmp, AboveTh(mu) - MEPS, nfmu, order);
838 }
839
840 if (nfAls < nfmu)
841 {
842 mutmp = AboveTh(MAls) - MEPS;
843 alstmp = AlsWithInit(mutmp, AlsM, MAls, nfAls, realorder);
844 alstmp *= (1. + NfThresholdCorrections(mutmp, MassOfNf(nfAls + 1), alstmp, nfAls + 1, fullord));
845 for (i = nfAls + 1; i < nfmu; i++)
846 {
847 mutmp = AboveTh(mutmp) - MEPS;
848 alstmp = AlsWithInit(mutmp, alstmp, BelowTh(mutmp) + MEPS, i, realorder);
849 alstmp *= (1. + NfThresholdCorrections(mutmp, MassOfNf(i + 1), alstmp, i + 1, fullord));
850 }
851 als = AlsWithInit(mu, alstmp, BelowTh(mu) + MEPS, nfmu, order);
852 }
853
854 CacheShift(als_cache, 9);
855 als_cache[0][0] = mu;
856 als_cache[1][0] = (double)order;
857 als_cache[2][0] = AlsM;
858 als_cache[3][0] = MAls;
859 als_cache[4][0] = mut;
860 als_cache[5][0] = mub;
861 als_cache[6][0] = muc;
862 als_cache[7][0] = (int)Nf_thr;
863 als_cache[8][0] = als;
864
865 return als;
866 default:
867 throw std::runtime_error("QCD::Als(): " + orderToString(order) + " is not implemented.");
868 }
869}
870
871const double QCD::AlsOLD(const double mu, const orders order) const
872{
873 int i;
874 for (i = 0; i < CacheSize; ++i)
875 if ((mu == als_cache[0][i]) && ((double)order == als_cache[1][i]) &&
876 (AlsM == als_cache[2][i]) && (MAls == als_cache[3][i]) &&
877 (mut == als_cache[4][i]) && (mub == als_cache[5][i]) &&
878 (muc == als_cache[6][i]) && -1 == (int)als_cache[7][i])
879 return als_cache[8][i];
880
881 double nfmu = Nf(mu), nfz = Nf(MAls), mu_thre1, mu_thre2, Als_tmp, mf;
882 double als;
883
884 switch (order)
885 {
886 case LO:
887 case FULLNLO:
888 case NLO:
889 if (nfmu == nfz)
890 als = AlsWithInit(mu, AlsM, MAls, nfmu, order);
891 else if (nfmu > nfz)
892 {
893 if (order == NLO)
894 throw std::runtime_error("NLO is not implemented in QCD::Als(mu,order).");
895 if (nfmu == nfz + 1.)
896 {
897 mu_thre1 = AboveTh(MAls); // mut
898 Als_tmp = AlsWithInit(mu_thre1 - MEPS, AlsM, MAls, nfz, order);
899 if (order == FULLNLO)
900 {
901 mf = getQuarks(TOP).getMass(); // mf = mtpole;
902 Als_tmp = (1. + Als_tmp / M_PI * log(mu_thre1 / mf) / 3.) * Als_tmp;
903 }
904 als = AlsWithInit(mu, Als_tmp, mu_thre1 + MEPS, nfmu, order);
905 }
906 else
907 {
908 QCDsuccess = false;
909 std::cout << "Error in QCD::Als(mu,order)! QCDsuccess set to false" << std::endl;
910 als = 0.;
911 }
912 }
913 else
914 {
915 if (order == NLO)
916 throw std::runtime_error("NLO is not implemented in QCD::Als(mu,order).");
917 if (nfmu == nfz - 1.)
918 {
919 mu_thre1 = BelowTh(MAls); // mub
920 Als_tmp = AlsWithInit(mu_thre1 + MEPS, AlsM, MAls, nfz, order);
921 if (order == FULLNLO)
922 {
923 mf = getQuarks(BOTTOM).getMass();
924 Als_tmp = (1. - Als_tmp / M_PI * log(mu_thre1 / mf) / 3.) * Als_tmp;
925 }
926 als = AlsWithInit(mu, Als_tmp, mu_thre1 - MEPS, nfmu, order);
927 }
928 else if (nfmu == nfz - 2.)
929 {
930 mu_thre1 = BelowTh(MAls); // mub
931 mu_thre2 = AboveTh(mu); // muc
932 Als_tmp = Als(mu_thre1 + MEPS, order);
933 if (order == FULLNLO)
934 {
935 mf = getQuarks(BOTTOM).getMass();
936 Als_tmp = (1. - Als_tmp / M_PI * log(mu_thre1 / mf) / 3.) * Als_tmp;
937 }
938 Als_tmp = AlsWithInit(mu_thre2, Als_tmp, mu_thre1 - MEPS, nfz-1, order);
939 if (order == FULLNLO)
940 {
941 mf = getQuarks(CHARM).getMass();
942 Als_tmp = (1. - Als_tmp / M_PI * log(mu_thre2 / mf) / 3.) * Als_tmp;
943 }
944 als = AlsWithInit(mu, Als_tmp, mu_thre2 - MEPS, nfmu, order);
945 }
946 else
947 {
948 QCDsuccess = false;
949 std::cout << "Error in QCD::Als(mu,order)! QCDsuccess set to false" << std::endl;
950 als = 0.;
951 }
952 }
953 break;
954 case NNLO:
955 // case NNNLO:
956 case FULLNNLO:
957 // case FULLNNNLO:
958 /* alpha_s(mu) computed with Lambda_QCD for Nf=nfmu */
959 als = AlsWithLambda(mu, order);
960 break;
961 default:
962 throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Als(mu,order).");
963 }
964
965 CacheShift(als_cache, 9);
966 als_cache[0][0] = mu;
967 als_cache[1][0] = (double)order;
968 als_cache[2][0] = AlsM;
969 als_cache[3][0] = MAls;
970 als_cache[4][0] = mut;
971 als_cache[5][0] = mub;
972 als_cache[6][0] = muc;
973 als_cache[7][0] = -1;
974 als_cache[8][0] = als;
975
976 return als;
977}
978
979const double QCD::AlsByOrder(const double mu, const int Nf_in, const orders order) const
980{
981 int i, nfAls = (int)Nf(MAls), nfmu = Nf_in;
982 double als, alstmp, mutmp;
983 orders fullord;
984
985 for (i = 0; i < CacheSize; ++i)
986 if ((mu == als_cache[0][i]) && ((double)order == als_cache[1][i]) &&
987 (AlsM == als_cache[2][i]) && (MAls == als_cache[3][i]) &&
988 (mut == als_cache[4][i]) && (mub == als_cache[5][i]) &&
989 (muc == als_cache[6][i]) && (Nf_in == (int)als_cache[7][i]))
990 return als_cache[8][i];
991
992 switch (order)
993 {
994 case LO:
995 case NLO:
996 case NNLO:
997 case NNNLO:
998 if (nfAls == nfmu)
999 als = AlsWithInit(mu, AlsM, MAls, nfmu, order);
1000 fullord = FullOrder(order);
1001 if (nfAls > nfmu)
1002 {
1003 mutmp = BelowTh(MAls);
1004 alstmp = AlsWithInit(mutmp, AlsM, MAls, nfAls, realorder);
1005 alstmp *= (1. - NfThresholdCorrections(mutmp, MassOfNf(nfAls), alstmp, nfAls, fullord));
1006 for (i = nfAls - 1; i > nfmu; i--)
1007 {
1008 mutmp = BelowTh(mutmp - MEPS);
1009 alstmp = AlsWithInit(mutmp, alstmp, AboveTh(mutmp) - MEPS, i, realorder);
1010 alstmp *= (1. - NfThresholdCorrections(mutmp, MassOfNf(i), alstmp, i, fullord));
1011 }
1012 als = AlsWithInit(mu, alstmp, mutmp, nfmu, order);
1013 }
1014
1015 if (nfAls < nfmu)
1016 {
1017 mutmp = AboveTh(MAls) - MEPS;
1018 alstmp = AlsWithInit(mutmp, AlsM, MAls, nfAls, realorder);
1019 alstmp *= (1. + NfThresholdCorrections(mutmp, MassOfNf(nfAls + 1), alstmp, nfAls + 1, fullord));
1020 for (i = nfAls + 1; i < nfmu; i++)
1021 {
1022 double oldmu = mutmp;
1023 mutmp = AboveTh(mutmp + 2. * MEPS) - MEPS;
1024 alstmp = AlsWithInit(mutmp, alstmp, oldmu, i, realorder);
1025 alstmp *= (1. + NfThresholdCorrections(mutmp, MassOfNf(i + 1), alstmp, i + 1, fullord));
1026 }
1027 als = AlsWithInit(mu, alstmp, mutmp + 2. * MEPS, nfmu, order);
1028 }
1029
1030 CacheShift(als_cache, 9);
1031 als_cache[0][0] = mu;
1032 als_cache[1][0] = (double)order;
1033 als_cache[2][0] = AlsM;
1034 als_cache[3][0] = MAls;
1035 als_cache[4][0] = mut;
1036 als_cache[5][0] = mub;
1037 als_cache[6][0] = muc;
1038 als_cache[7][0] = Nf_in;
1039 als_cache[8][0] = als;
1040
1041 return als;
1042 default:
1043 throw std::runtime_error("QCD::Als(): " + orderToString(order) + " is not implemented.");
1044 }
1045}
1046
1047const double QCD::ZeroNf6NLO(double *logLambda6, double *logLambda5_in) const
1048{
1049 return (AlsWithLambda(mut + 1.e-10, *logLambda6, FULLNLO) - AlsWithLambda(mut - 1.e-10, *logLambda5_in, FULLNLO));
1050}
1051
1052const double QCD::ZeroNf5(double *logLambda5, double *order) const
1053{
1054 return (AlsWithLambda(MAls, *logLambda5, (orders)*order) - AlsM);
1055}
1056
1057const double QCD::ZeroNf4NLO(double *logLambda4, double *logLambda5_in) const
1058{
1059 return (AlsWithLambda(mub - 1.e-10, *logLambda4, FULLNLO) - AlsWithLambda(mub + 1.e-10, *logLambda5_in, FULLNLO));
1060}
1061
1062const double QCD::ZeroNf3NLO(double *logLambda3, double *logLambda4_in) const
1063{
1064 return (AlsWithLambda(muc - 1.e-10, *logLambda3, FULLNLO) - AlsWithLambda(muc + 1.e-10, *logLambda4_in, FULLNLO));
1065}
1066
1067const double QCD::logLambda5(orders order) const
1068{
1069 if (order == NLO)
1070 order = FULLNLO;
1071 if (order == NNLO)
1072 order = FULLNNLO;
1073
1074 for (int i = 0; i < CacheSize; ++i)
1075 if ((AlsM == logLambda5_cache[0][i]) && (MAls == logLambda5_cache[1][i]) && ((double)order == logLambda5_cache[2][i]))
1076 return logLambda5_cache[3][i];
1077
1078 CacheShift(logLambda5_cache, 4);
1079 logLambda5_cache[0][0] = AlsM;
1080 logLambda5_cache[1][0] = MAls;
1081 logLambda5_cache[2][0] = (double)order;
1082
1083 if (order == LO)
1084 logLambda5_cache[3][0] = log(MAls) - 2. * M_PI / Beta0(5.) / AlsM;
1085 else
1086 {
1087 double xmin = -4., xmax = -0.2;
1088 TF1 f = TF1("f", this, &QCD::ZeroNf5, xmin, xmax, 1, "QCD", "zeroNf5");
1089
1090 ROOT::Math::WrappedTF1 wf1(f);
1091 double ledouble = (double)order;
1092 wf1.SetParameters(&ledouble);
1093
1094 ROOT::Math::BrentRootFinder brf;
1095 brf.SetFunction(wf1, xmin, xmax);
1096
1097 if (brf.Solve())
1098 logLambda5_cache[3][0] = brf.Root();
1099 else
1100 {
1101 QCDsuccess = false;
1102 std::cout << "Error in QCD::logLambda5()! QCDsuccess set to false" << std::endl;
1103 }
1104 }
1105 return (logLambda5_cache[3][0]);
1106}
1107
1108const double QCD::logLambdaNLO(const double nfNEW, const double nfORG,
1109 const double logLambdaORG) const
1110{
1111 for (int i = 0; i < CacheSize; ++i)
1112 if ((AlsM == logLambdaNLO_cache[0][i]) && (MAls == logLambdaNLO_cache[1][i]) && (mut == logLambdaNLO_cache[2][i]) && (mub == logLambdaNLO_cache[3][i]) && (muc == logLambdaNLO_cache[4][i]) && (nfNEW == logLambdaNLO_cache[5][i]) && (nfORG == logLambdaNLO_cache[6][i]) && (logLambdaORG == logLambdaNLO_cache[7][i]))
1113 return logLambdaNLO_cache[8][i];
1114
1115 CacheShift(logLambdaNLO_cache, 9);
1116 logLambdaNLO_cache[0][0] = AlsM;
1117 logLambdaNLO_cache[1][0] = MAls;
1118 logLambdaNLO_cache[2][0] = mut;
1119 logLambdaNLO_cache[3][0] = mub;
1120 logLambdaNLO_cache[4][0] = muc;
1121 logLambdaNLO_cache[5][0] = nfNEW;
1122 logLambdaNLO_cache[6][0] = nfORG;
1123 logLambdaNLO_cache[7][0] = logLambdaORG;
1124
1125 double xmin = -4., xmax = -0.2;
1126
1127 TF1 f;
1128 if (nfNEW == 6. && nfORG == 5.)
1129 {
1130 f = TF1("f", this, &QCD::ZeroNf6NLO, xmin, xmax, 1, "QCD", "zeroNf6NLO");
1131 }
1132 else if (nfNEW == 4. && nfORG == 5.)
1133 {
1134 f = TF1("f", this, &QCD::ZeroNf4NLO, xmin, xmax, 1, "QCD", "zeroNf4NLO");
1135 }
1136 else if (nfNEW == 3. && nfORG == 4.)
1137 {
1138 f = TF1("f", this, &QCD::ZeroNf3NLO, xmin, xmax, 1, "QCD", "zeroNf3NLO");
1139 }
1140 else
1141 {
1142 QCDsuccess = false;
1143 std::cout << "Error in QCD::logLambdaNLO()! QCDsuccess set to false" << std::endl;
1144 }
1145
1146 ROOT::Math::WrappedTF1 wf1(f);
1147 wf1.SetParameters(&logLambdaORG);
1148
1149 ROOT::Math::BrentRootFinder brf;
1150 brf.SetFunction(wf1, xmin, xmax);
1151
1152 if (brf.Solve())
1153 logLambdaNLO_cache[8][0] = brf.Root();
1154 else
1155 {
1156 QCDsuccess = false;
1157 std::cout << "Error in QCD::logLambdaNLO()! QCDsuccess set to false" << std::endl;
1158 }
1159
1160 return (logLambdaNLO_cache[8][0]);
1161}
1162
1163const double QCD::logLambda(const double muMatching, const double mf,
1164 const double nfNEW, const double nfORG,
1165 const double logLambdaORG, orders order) const
1166{
1167 if (fabs(nfNEW - nfORG) != 1.)
1168 {
1169 QCDsuccess = false;
1170 std::cout << "Error in QCD::logLambda()! QCDsuccess set to false" << std::endl;
1171 }
1172
1173 if (order == NLO)
1174 order = FULLNLO;
1175 if (order == NNLO)
1176 order = FULLNNLO;
1177
1178 /* We do not use the codes below for FULLNLO, since threshold corrections
1179 * can be regarded as an NNLO effect as long as setting the matching scale
1180 * to be close to the mass scale of the decoupling quark. In order to use
1181 * the relation als^{nf+1} = als^{nf} exactly, we use logLambdaNLO method.
1182 */
1183 if (order == FULLNLO)
1184 return logLambdaNLO(nfNEW, nfORG, logLambdaORG);
1185
1186 double logMuMatching = log(muMatching);
1187 double L = 2. * (logMuMatching - logLambdaORG);
1188 double rNEW = 0.0, rORG = 0.0, log_mu2_mf2 = 0.0, log_L = 0.0;
1189 double C1 = 0.0, C2 = 0.0; // threshold corrections
1190 double logLambdaNEW;
1191
1192 // LO contribution
1193 logLambdaNEW = 1. / 2. / Beta0(nfNEW) * (Beta0(nfNEW) - Beta0(nfORG)) * L + logLambdaORG;
1194
1195 // NLO contribution
1196 if (order == FULLNLO || order == FULLNNLO)
1197 {
1198 rNEW = Beta1(nfNEW) / Beta0(nfNEW);
1199 rORG = Beta1(nfORG) / Beta0(nfORG);
1200 log_mu2_mf2 = 2. * (logMuMatching - log(mf));
1201 log_L = log(L);
1202 if (nfNEW < nfORG)
1203 C1 = 2. / 3. * log_mu2_mf2;
1204 else
1205 C1 = -2. / 3. * log_mu2_mf2;
1206 logLambdaNEW += 1. / 2. / Beta0(nfNEW) * ((rNEW - rORG) * log_L - rNEW * log(Beta0(nfNEW) / Beta0(nfORG)) - C1);
1207 }
1208
1209 // NNLO contribution
1210 if (order == FULLNNLO)
1211 {
1212 if (nfNEW == 5. && nfORG == 6.)
1213 C2 = -16. * (log_mu2_mf2 * log_mu2_mf2 / 36. - 19. / 24. * log_mu2_mf2 - 7. / 24.);
1214 else if (nfNEW == 6. && nfORG == 5.)
1215 C2 = -16. * (log_mu2_mf2 * log_mu2_mf2 / 36. + 19. / 24. * log_mu2_mf2 + 7. / 24.);
1216 else
1217 {
1218 if (nfNEW < nfORG)
1219 C2 = -16. * (log_mu2_mf2 * log_mu2_mf2 / 36. - 19. / 24. * log_mu2_mf2 + 11. / 72.);
1220 else
1221 C2 = -16. * (log_mu2_mf2 * log_mu2_mf2 / 36. + 19. / 24. * log_mu2_mf2 - 11. / 72.);
1222 }
1223 logLambdaNEW += 1. / 2. / Beta0(nfNEW) / Beta0(nfORG) / L * (rORG * (rNEW - rORG) * log_L + rNEW * rNEW - rORG * rORG - Beta2(nfNEW) / Beta0(nfNEW) + Beta2(nfORG) / Beta0(nfORG) + rNEW * C1 - C1 * C1 - C2);
1224 }
1225
1226 return logLambdaNEW;
1227}
1228
1229const double QCD::logLambda(const double nf, orders order) const
1230{
1231 if (order == NLO)
1232 order = FULLNLO;
1233 if (order == NNLO)
1234 order = FULLNNLO;
1235
1236 double muMatching, mf, logLambdaORG, logLambdaNEW;
1237 if (nf == 5.)
1238 return logLambda5(order);
1239 else if (nf == 6.)
1240 {
1241 muMatching = Thresholds(1); // mut
1242 /* matching condition from Nf=5 to Nf=6 is given in terms of the top pole mass. */
1243 mf = mtpole; // top pole mass
1244 return logLambda(muMatching, mf, 6., 5., logLambda5(order), order);
1245 }
1246 else if (nf == 4. || nf == 3.)
1247 {
1248 muMatching = Thresholds(2); // mub
1249 mf = getQuarks(BOTTOM).getMass(); // m_b(m_b)
1250 logLambdaORG = logLambda5(order);
1251 logLambdaNEW = logLambda(muMatching, mf, 4., 5., logLambdaORG, order);
1252 if (nf == 3.)
1253 {
1254 muMatching = Thresholds(3); // muc
1255 mf = getQuarks(CHARM).getMass(); // m_c(m_c)
1256 logLambdaORG = logLambdaNEW;
1257 logLambdaNEW = logLambda(muMatching, mf, 3., 4., logLambdaORG, order);
1258 }
1259 return logLambdaNEW;
1260 }
1261 else
1262 {
1263 QCDsuccess = false;
1264 std::cout << "Error in QCD::logLambda()! QCDsuccess set to false" << std::endl;
1265 }
1266 return 0.;
1267}
1268
1270
1271const double QCD::Gamma0(const double nf) const
1272{
1273 return (6. * CF);
1274}
1275
1276const double QCD::Gamma1(const double nf) const
1277{
1278 return (CF * (3. * CF + 97. / 3. * Nc - 10. / 3. * nf));
1279}
1280
1281const double QCD::Gamma2(const double nf) const
1282{
1283 return (129. * CF * CF * CF - 129. / 2. * CF * CF * Nc + 11413. / 54. * CF * Nc * Nc + CF * CF * nf * (-46. + 48. * zeta3) + CF * Nc * nf * (-556. / 27. - 48. * zeta3) - 70. / 27. * CF * nf * nf);
1284}
1285
1286const double QCD::threCorrForMass(const double nf_f, const double nf_i) const
1287{
1288 if (fabs(nf_f - nf_i) != 1.)
1289 {
1290 QCDsuccess = false;
1291 std::cout << "Error in QCD::threCorrForMass()! QCDsuccess set to false" << std::endl;
1292 return 0.;
1293 }
1294
1295 // hep-ph/9712201v2 eq. (48)
1296 double mu_threshold, mf, log_mu2_mf2, AlsoverPi;
1297 if (nf_f > nf_i)
1298 {
1299 if (nf_f == 6.)
1300 {
1301 mu_threshold = mut;
1302 mf = quarks[TOP].getMass(); // m_t(m_t)
1303 }
1304 else if (nf_f == 5.)
1305 {
1306 mu_threshold = mub;
1307 mf = quarks[BOTTOM].getMass(); // m_b(m_b)
1308 }
1309 else if (nf_f == 4.)
1310 {
1311 mu_threshold = muc;
1312 mf = quarks[CHARM].getMass(); // m_c(m_c)
1313 }
1314 else
1315 {
1316 QCDsuccess = false;
1317 std::cout << "Error in QCD::threCorrForMass()! QCDsuccess set to false" << std::endl;
1318 return 0.;
1319 }
1320 log_mu2_mf2 = 2. * log(mu_threshold / mf);
1321 AlsoverPi = Als(mu_threshold - MEPS, FULLNNLO) / M_PI;
1322 return (1. + AlsoverPi * AlsoverPi * (-log_mu2_mf2 * log_mu2_mf2 / 12. + 5. / 36. * log_mu2_mf2 - 89. / 432.));
1323 }
1324 else
1325 {
1326 if (nf_i == 6.)
1327 {
1328 mu_threshold = mut;
1329 mf = quarks[TOP].getMass(); // m_t(m_t)
1330 }
1331 else if (nf_i == 5.)
1332 {
1333 mu_threshold = mub;
1334 mf = quarks[BOTTOM].getMass(); // m_b(m_b)
1335 }
1336 else if (nf_i == 4.)
1337 {
1338 mu_threshold = muc;
1339 mf = quarks[CHARM].getMass(); // m_c(m_c)
1340 }
1341 else
1342 {
1343 QCDsuccess = false;
1344 std::cout << "Error in QCD::threCorrForMass()! QCDsuccess set to false" << std::endl;
1345 return 0.;
1346 }
1347 log_mu2_mf2 = 2. * log(mu_threshold / mf);
1348 AlsoverPi = Als(mu_threshold + MEPS, FULLNNLO) / M_PI;
1349 return (1. + AlsoverPi * AlsoverPi * (log_mu2_mf2 * log_mu2_mf2 / 12. - 5. / 36. * log_mu2_mf2 + 89. / 432.));
1350 }
1351}
1352
1353const double QCD::Mrun(const double mu, const double m, const quark q, const orders order) const
1354{
1355 return Mrun(mu, m, m, q, order);
1356}
1357
1358const double QCD::Mrun(const double mu_f, const double mu_i, const double m, const quark q,
1359 const orders order) const
1360{
1361 // Note: When the scale evolves across a flavour threshold, the definitions
1362 // of the outputs for "NLO" and "NNLO" become complicated.
1363
1364 int i;
1365 for (i = 0; i < CacheSize; ++i)
1366 {
1367 if ((mu_f == mrun_cache[0][i]) && (mu_i == mrun_cache[1][i]) &&
1368 (m == mrun_cache[2][i]) && ((double)order == mrun_cache[3][i]) &&
1369 ((double)q == mrun_cache[4][i]) &&
1370 (AlsM == mrun_cache[5][i]) && (MAls == mrun_cache[6][i]) &&
1371 (mut == mrun_cache[7][i]) && (mub == mrun_cache[8][i]) &&
1372 (muc == mrun_cache[9][i]))
1373 return mrun_cache[10][i];
1374 }
1375
1376 double nfi = Nf(mu_i), nff = Nf(mu_f);
1377 double nfq;
1378 if (q == TOP)
1379 nfq = 6.;
1380 else if (q == BOTTOM)
1381 nfq = 5.;
1382 else if (q == CHARM)
1383 nfq = 4.;
1384 else
1385 nfq = 3.;
1386
1387 nfi = (nfi > nfq) ? nfi : nfq;
1388 nff = (nff > nfq) ? nff : nfq;
1389
1390 double mu_threshold, mu_threshold2, mu_threshold3, mrun;
1391 if (nff == nfi)
1392 mrun = MrunTMP(mu_f, mu_i, m, lround(nfi), order);
1393 else if (nff > nfi)
1394 {
1395 if (order == NLO || order == NNLO)
1396 throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Mrun(mu_f,mu_i,m,q,order)");
1397 mu_threshold = Nf(mu_i) < nfi ? Thresholds(6-lround(nfi)): AboveTh(mu_i) ;
1398 mrun = MrunTMP(mu_threshold - MEPS, mu_i, m, lround(nfi), order);
1399 if (order == FULLNNLO)
1400 mrun *= threCorrForMass(nfi + 1., nfi); // threshold corrections
1401 if (nff == nfi + 1.)
1402 {
1403 mrun = MrunTMP(mu_f, mu_threshold + MEPS, mrun, lround(nff), order);
1404 }
1405 else if (nff == nfi + 2.)
1406 {
1407 mu_threshold2 = BelowTh(mu_f);
1408 mrun = MrunTMP(mu_threshold2 - MEPS, mu_threshold + MEPS, mrun, lround(nfi) + 1, order);
1409 if (order == FULLNNLO)
1410 mrun *= threCorrForMass(nff, nfi + 1.); // threshold corrections
1411 mrun = MrunTMP(mu_f, mu_threshold2 + MEPS, mrun, lround(nff), order);
1412 }
1413 else if (nff == nfi + 3.)
1414 {
1415 mu_threshold2 = AboveTh(mu_threshold);
1416 mrun = MrunTMP(mu_threshold2 - MEPS, mu_threshold + MEPS, mrun, lround(nfi) + 1, order);
1417 if (order == FULLNNLO)
1418 mrun *= threCorrForMass(nfi + 2., nfi + 1.); // threshold corrections
1419 mu_threshold3 = BelowTh(mu_f);
1420 mrun = MrunTMP(mu_threshold3 - MEPS, mu_threshold2 + MEPS, mrun, lround(nff) - 1, order);
1421 if (order == FULLNNLO)
1422 mrun *= threCorrForMass(nff, nfi + 2.); // threshold corrections
1423 mrun = MrunTMP(mu_f, mu_threshold3 + MEPS, mrun, lround(nff), order);
1424 }
1425 else
1426 {
1427 QCDsuccess = false;
1428 std::cout << "Error in QCD::Mrun(mu_f,mu_i,m,order)! QCDsuccess set to false" << std::endl;
1429 mrun = 0.;
1430 }
1431 }
1432 else
1433 {
1434 if (order == NLO || order == NNLO)
1435 throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Mrun(mu_f,mu_i,m,order)");
1436 mu_threshold = BelowTh(mu_i);
1437 mrun = MrunTMP(mu_threshold + MEPS, mu_i, m, lround(nfi), order);
1438 if (order == FULLNNLO)
1439 mrun *= threCorrForMass(nfi - 1., nfi); // threshold corrections
1440 if (nff == nfi - 1.)
1441 mrun = MrunTMP(mu_f, mu_threshold - MEPS, mrun, lround(nff), order);
1442 else if (nff == nfi - 2.)
1443 {
1444 mu_threshold2 = Nf(mu_f) < nff ? Thresholds(6-lround(nff)) : AboveTh(mu_f);
1445 mrun = MrunTMP(mu_threshold2 + MEPS, mu_threshold - MEPS, mrun, lround(nfi) - 1, order);
1446 if (order == FULLNNLO)
1447 mrun *= threCorrForMass(nff, nfi - 1.); // threshold corrections
1448 mrun = MrunTMP(mu_f, mu_threshold2 - MEPS, mrun, lround(nff), order);
1449 }
1450 else
1451 {
1452 QCDsuccess = false;
1453 std::cout << "Error in QCD::Mrun(mu_f,mu_i,m,order)! QCDsuccess set to false" << std::endl;
1454 mrun = 0.;
1455 }
1456 }
1457
1458 if (mrun < 0.0)
1459 {
1460 std::stringstream out;
1461 out << "QCD::Mrun(): A quark mass becomes tachyonic in QCD::Mrun("
1462 << mu_f << ", " << mu_i << ", " << m << ", " << orderToString(order) << ")"
1463 << std::endl
1464 << " Als(" << mu_i << ", " << orderToString(order) << ")/(4pi)="
1465 << Als(mu_i, order) / (4. * M_PI) << std::endl
1466 << " Als(" << mu_f << ", " << orderToString(order) << ")/(4pi)="
1467 << Als(mu_f, order) / (4. * M_PI);
1468 QCDsuccess = false;
1469 std::cout << "Error in QCD::Mrun(mu_f,mu_i,m,order)! QCDsuccess set to false" << std::endl;
1470 mrun = 0.;
1471 }
1472
1473 CacheShift(mrun_cache, 11);
1474 mrun_cache[0][0] = mu_f;
1475 mrun_cache[1][0] = mu_i;
1476 mrun_cache[2][0] = m;
1477 mrun_cache[3][0] = (double)order;
1478 mrun_cache[4][0] = (double)q;
1479 mrun_cache[5][0] = AlsM;
1480 mrun_cache[6][0] = MAls;
1481 mrun_cache[7][0] = mut;
1482 mrun_cache[8][0] = mub;
1483 mrun_cache[9][0] = muc;
1484 mrun_cache[10][0] = mrun;
1485
1486 return mrun;
1487}
1488
1489const double QCD::MrunTMP(const double mu_f, const double mu_i, const double m, const int nf,
1490 const orders order) const
1491{
1492
1493 // alpha_s/(4pi)
1494 orders orderForAls;
1495 if (order == LO)
1496 orderForAls = LO;
1497 if (order == NLO || order == FULLNLO)
1498 orderForAls = FULLNLO;
1499 if (order == NNLO || order == FULLNNLO)
1500 orderForAls = FULLNNLO;
1501 double ai = Als(mu_i, nf, orderForAls) / (4. * M_PI);
1502 double af = Als(mu_f, nf, orderForAls) / (4. * M_PI);
1503
1504 // LO contribution
1505 double b0 = Beta0((double) nf), g0 = Gamma0((double) nf);
1506 double mLO = m * pow(af / ai, g0 / (2. * b0));
1507 if (order == LO)
1508 return mLO;
1509
1510 // NLO contribution
1511 double b1 = Beta1((double) nf), g1 = Gamma1((double) nf);
1512 double A1 = g1 / (2. * b0) - b1 * g0 / (2. * b0 * b0);
1513 double mNLO = mLO * A1 * (af - ai);
1514 if (order == NLO)
1515 return mNLO;
1516 if (order == FULLNLO)
1517 return (mLO + mNLO);
1518
1519 // NNLO contribution: Chetyrkin, hep-ph/9703278v1 eq. (15)
1520 double b2 = Beta2((double) nf), g2 = Gamma2((double) nf);
1521 double A2 = b1 * b1 * g0 / (2. * b0 * b0 * b0) - b2 * g0 / (2. * b0 * b0) - b1 * g1 / (2. * b0 * b0) + g2 / (2. * b0);
1522 double mNNLO = mLO * (A1 * A1 / 2. * (af - ai) * (af - ai) + A2 / 2. * (af * af - ai * ai));
1523 if (order == NNLO)
1524 return mNNLO;
1525 if (order == FULLNNLO)
1526 return (mLO + mNLO + mNNLO);
1527
1528 throw std::runtime_error(orderToString(order) + " is not implemented in QCD::MrunTMP()");
1529}
1530
1531const double QCD::Mrun4(const double mu_f, const double mu_i, const double m) const
1532{
1533 double nf = 4.;
1534
1535 // alpha_s/(4pi)
1536 double ai = Als4(mu_i) / (4. * M_PI);
1537 double af = Als4(mu_f) / (4. * M_PI);
1538
1539 // LO contribution
1540 double b0 = Beta0(nf), g0 = Gamma0(nf);
1541 double mLO = m * pow(af / ai, g0 / (2. * b0));
1542
1543 // NLO contribution
1544 double b1 = Beta1(nf), g1 = Gamma1(nf);
1545 double A1 = g1 / (2. * b0) - b1 * g0 / (2. * b0 * b0);
1546 double mNLO = mLO * A1 * (af - ai);
1547 return (mLO + mNLO);
1548}
1549
1551
1552const double QCD::Mbar2Mp(const double mbar, const quark q, const orders order) const
1553{
1554 int Nf_m;
1555 double nl;
1556 std::vector<double> x = {0.0};
1557
1558 if (q == CHARM)
1559 Nf_m = 4;
1560 else if (q == BOTTOM)
1561 Nf_m = 5;
1562 else if (q == TOP)
1563 Nf_m = 6;
1564 else
1565 throw std::runtime_error("QCD::Mbar2Mp() can convert only top, bottom and charm masses");
1566
1567 nl = Nf_m-1;
1568
1569 // LO contribution
1570 double MpLO = mbar;
1571 if (order == LO)
1572 return MpLO;
1573
1574 // alpha_s(mbar)/pi
1575 orders orderForAls;
1576 if (order == NLO || order == FULLNLO)
1577 orderForAls = FULLNLO;
1578 if (order == NNLO || order == FULLNNLO)
1579 orderForAls = FULLNNLO;
1580 double a = Als(mbar, Nf_m, orderForAls) / M_PI;
1581
1582 // NLO contribution
1583 double MpNLO = mbar * 4. / 3. * a;
1584 if (order == NLO)
1585 return MpNLO;
1586 if (order == FULLNLO)
1587 return (MpLO + MpNLO);
1588
1589 // NNLO contribution
1590
1591 x.push_back(quarks[STRANGE].getMass() / mbar);
1592 if (nl >3)
1593 x.push_back(quarks[CHARM].getMass() / mbar);
1594 if (nl > 4)
1595 x.push_back(quarks[BOTTOM].getMass() / mbar);
1596
1597 double Delta = 0.;
1598 for (int i = 0; i < x.size(); i++)
1599 Delta += DeltaMass(x[i]);
1600
1601 double MpNNLO = mbar * (307. / 32. + 2. * zeta2 + 2. / 3. * zeta2 * log(2.0) - zeta3 / 6. - nl / 3. * (zeta2 + 71. / 48.) + 4. / 3. * Delta) * a * a;
1602
1603 if (order == NNLO)
1604 return MpNNLO;
1605 if (order == FULLNNLO)
1606 return (MpLO + MpNLO + MpNNLO);
1607
1608 throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Mbar2Mp().");
1609}
1610
1611double QCD::DeltaMass(double x) const
1612{
1613 return M_PI * M_PI / 8. * x - 0.597 * x * x + 0.230 * x * x * x;
1614}
1615
1616
1617const double QCD::Mp2MbarTMP(double *mu, double *params) const
1618{
1619 double mp = params[0];
1620 quark q = (quark)params[1];
1621 orders order = (orders)params[2];
1622 return (mp - Mbar2Mp(*mu, q, order));
1623}
1624
1625const double QCD::Mp2Mbar_bar(const double mp, const quark q, const orders order) const
1626{
1627 if (order == NLO || order == NNLO)
1628 throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Mp2Mbar().");
1629
1630 int i;
1631 double ms = quarks[STRANGE].getMass(), mc = quarks[CHARM].getMass(),
1632 mb = quarks[BOTTOM].getMass();
1633 double alsmp = Als(mp, order);
1634 for (i = 0; i < CacheSize; ++i)
1635 if (alsmp == mp2mbar_cache[0][i] && ms == mp2mbar_cache[1][i] &&
1636 mc == mp2mbar_cache[2][i] && mb == mp2mbar_cache[3][i]
1637 && (double)order == mp2mbar_cache[4][i])
1638 return mp2mbar_cache[5][i];
1639
1640 TF1 f("f", this, &QCD::Mp2MbarTMP, mp / 2., 2. * mp, 3, "QCD", "Mp2MbarTMP");
1641
1642 ROOT::Math::WrappedTF1 wf1(f);
1643 double params[3];
1644 params[0] = mp;
1645 params[1] = (double)q;
1646 params[2] = (double)order;
1647 wf1.SetParameters(params);
1648
1649 ROOT::Math::BrentRootFinder brf;
1650
1651 brf.SetNpx(5);
1652 brf.SetFunction(wf1, .5 * mp, 2. * mp);
1653 if (brf.Solve())
1654 {
1655 CacheShift(mp2mbar_cache, 6);
1656 mp2mbar_cache[0][0] = alsmp;
1657 mp2mbar_cache[1][0] = ms;
1658 mp2mbar_cache[2][0] = mc;
1659 mp2mbar_cache[3][0] = mb;
1660 mp2mbar_cache[4][0] = (double)order;
1661 mp2mbar_cache[5][0] = brf.Root();
1662 }
1663 else
1664 {
1665 QCDsuccess = false;
1666 std::cout << "Error in QCD::Mp2Mbar()! QCDsuccess set to false" << std::endl;
1667 }
1668
1669 return (mp2mbar_cache[5][0]);
1670}
1671
1672const double QCD::Mp2Mbar(const double mp, const quark q, const orders order) const
1673{
1674 if (q == TOP)
1675// return Mp2Mbar_pole(mp, order);
1676 return Mp2Mbar_bar(mp, q, order);
1677 else if (q == BOTTOM || q == CHARM)
1678 return Mp2Mbar_bar(mp, q, order);
1679 else
1680 throw std::runtime_error("QCD::Mp2Mbar() can convert only top, bottom and charm masses");
1681}
1682
1683const double QCD::Mp2Mbar_pole(const double mp, const orders order) const
1684{
1685 int Nfmp = 6;
1686 // LO contribution
1687 double mbarLO = mp;
1688 if (order == LO)
1689 return mbarLO;
1690
1691 // alpha_s(mp)/pi
1692 orders orderForAls;
1693 if (order == NLO || order == FULLNLO)
1694 orderForAls = FULLNLO;
1695 if (order == NNLO || order == FULLNNLO)
1696 orderForAls = FULLNNLO;
1697 double a = Als(mp, Nfmp, orderForAls) / M_PI;
1698
1699 // NLO contribution
1700 double mbarNLO = mp * (-4. / 3.) * a;
1701 if (order == NLO)
1702 return mbarNLO;
1703 if (order == FULLNLO)
1704 return (mbarLO + mbarNLO);
1705
1706 double mbarNNLO = mp * ((-1111. * CA * CF) / 384. + (199. * CF * CF) / 128. + (143. * CF * TF) / 96. +
1707 (71. * CF * ((double)Nfmp - 1.) * TF) / 96. + (CA * CF * zeta2) / 2. - (15. * CF * CF * zeta2) / 8. -
1708 (3. * CA * CF * log(2.) * zeta2) / 2. + 3. * CF * CF * log(2.) * zeta2 - CF * TF * zeta2 +
1709 (CF * ((double)Nfmp - 1.) * TF * zeta2) / 2. + (3. * CA * CF * zeta3) / 8. - (3. * CF * CF * zeta3) / 4.) * a * a;
1710
1711 //for the top, the correction from light quark masses is negligible
1712 if (order == NNLO)
1713 return mbarNNLO;
1714 if (order == FULLNNLO)
1715 return (mbarLO + mbarNLO + mbarNNLO);
1716
1717 throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Mp2Mbar_new().");
1718}
1719
1720const double QCD::MS2DRqmass(const double MSbar) const
1721{
1722 return (MSbar / (1. + Als(MSbar, FULLNLO) / 4. / M_PI * CF));
1723}
1724
1725const double QCD::MS2DRqmass(const double MSscale, const double MSbar) const
1726{
1727 return (MSbar / (1. + Als(MSscale, FULLNLO) / 4. / M_PI * CF));
1728}
1729
1730const double QCD::Mofmu2MbarTMP(double *mu, double *params) const
1731{
1732 double mofmu = params[0];
1733 double muI = params[1];
1734 quark q = (quark)params[2];
1735 return (*mu - Mrun(*mu, muI, mofmu, q));
1736}
1737
1738const double QCD::Mofmu2Mbar(const double m, const double mu, const quark q) const
1739{
1740
1741 // First move to the right region
1742
1743 double mutmp=0.;
1744 switch (q)
1745 {
1746 case TOP:
1747 mutmp = mut;
1748 break;
1749 case BOTTOM:
1750 mutmp = mub;
1751 break;
1752 case CHARM:
1753 mutmp = muc;
1754 break;
1755 default:
1756 QCDsuccess = false;
1757 std::cout << "Error in QCD::Mofmu2Mbar()! QCDsuccess set to false" << std::endl;
1758 break;
1759 }
1760 double mlow = Mrun(mutmp, mu, m, q);
1761 TF1 f("f", this, &QCD::Mofmu2MbarTMP, mlow / 2., 2. * mlow, 3, "QCD", "mofmu2mbara");
1762
1763 ROOT::Math::WrappedTF1 wf1(f);
1764
1765 double params[3];
1766 params[0] = mlow;
1767 params[1] = mutmp;
1768 params[2] = (double)q;
1769 wf1.SetParameters(params);
1770
1771 ROOT::Math::BrentRootFinder brf;
1772 brf.SetNpx(5);
1773 brf.SetFunction(wf1, .5 * mlow, 2. * mlow);
1774 if (brf.Solve())
1775 mlow = brf.Root();
1776 else
1777 {
1778 QCDsuccess = false;
1779 std::cout << "Error in QCD::Mofmu2Mbar()! QCDsuccess set to false" << std::endl;
1780 mlow = 0.;
1781 }
1782
1783 return (mlow);
1784}
DPars
std::map< std::string, double > DPars
Definition: Minimal.cpp:11
FULLNNNLO
@ FULLNNNLO
Definition: OrderScheme.h:40
NNLO
@ NNLO
Definition: OrderScheme.h:36
LO
@ LO
Definition: OrderScheme.h:34
NNNLO
@ NNNLO
Definition: OrderScheme.h:37
NLO
@ NLO
Definition: OrderScheme.h:35
FULLNNLO
@ FULLNNLO
Definition: OrderScheme.h:39
FULLNLO
@ FULLNLO
Definition: OrderScheme.h:38
BParameter
A class for the bag parameters.
Definition: BParameter.h:151
Meson
A class for mesons.
Definition: Meson.h:315
Model::addMissingModelParameter
void addMissingModelParameter(const std::string &missingParameterName)
Definition: Model.h:250
Model::ModelParamMap
std::map< std::string, std::reference_wrapper< const double > > ModelParamMap
Definition: Model.h:280
Model::name
std::string name
The name of the model.
Definition: Model.h:285
Model::setModelName
void setModelName(const std::string name)
A method to set the name of the model.
Definition: Model.h:50
Model::isSliced
bool isSliced
A boolean set to true if the current istance is a slice of an extended object.
Definition: Model.h:282
Model::UpdateError
bool UpdateError
A boolean set to false if update is successful.
Definition: Model.h:272
Model::raiseMissingModelParameterCount
void raiseMissingModelParameterCount()
Definition: Model.h:260
Model::getMissingModelParametersCount
unsigned int getMissingModelParametersCount()
Definition: Model.h:265
Particle
A class for particles.
Definition: Particle.h:26
Particle::getMass
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
Particle::setMass
void setMass(double mass)
A set method to fix the particle mass.
Definition: Particle.h:70
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::mp2mbar_cache
double mp2mbar_cache[6][CacheSize]
Cache for pole mass to msbar mass conversion.
Definition: QCD.h:1043
QCD::CacheSize
static const int CacheSize
Defines the depth of the cache.
Definition: QCD.h:1038
QCD::Als4
const double Als4(const double mu) const
The value of at any scale with the number of flavours .
Definition: QCD.cpp:657
QCD::logLambdaNLO
const double logLambdaNLO(const double nfNEW, const double nfORG, const double logLambdaORG) const
used for computation of at FULLNLO.
Definition: QCD.cpp:1108
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:421
QCD::meson
meson
An enum type for mesons.
Definition: QCD.h:336
QCD::OMEGA
@ OMEGA
Definition: QCD.h:355
QCD::K_P
@ K_P
Definition: QCD.h:340
QCD::PHI
@ PHI
Definition: QCD.h:348
QCD::K_star
@ K_star
Definition: QCD.h:349
QCD::D_P
@ D_P
Definition: QCD.h:342
QCD::B_C
@ B_C
Definition: QCD.h:347
QCD::B_P
@ B_P
Definition: QCD.h:345
QCD::K_star_P
@ K_star_P
Definition: QCD.h:350
QCD::P_0
@ P_0
Definition: QCD.h:337
QCD::RHO_P
@ RHO_P
Definition: QCD.h:354
QCD::D_S
@ D_S
Definition: QCD.h:343
QCD::D_star_P
@ D_star_P
Definition: QCD.h:352
QCD::K_0
@ K_0
Definition: QCD.h:339
QCD::K_S
@ K_S
Definition: QCD.h:351
QCD::B_D
@ B_D
Definition: QCD.h:344
QCD::RHO
@ RHO
Definition: QCD.h:353
QCD::D_0
@ D_0
Definition: QCD.h:341
QCD::P_P
@ P_P
Definition: QCD.h:338
QCD::B_S
@ B_S
Definition: QCD.h:346
QCD::requireYu
bool requireYu
Switch for generating the Yukawa couplings to the up-type quarks.
Definition: QCD.h:1012
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:199
QCD::ZeroNf3NLO
const 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:1062
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:460
QCD::Mp2Mbar_pole
const double Mp2Mbar_pole(const double mp, const orders order=FULLNNLO) const
Converts the top pole mass to the corresponding mass using eq. (16) of hep-ph/0004189.
Definition: QCD.cpp:1683
QCD::mut
double mut
The threshold between six- and five-flavour theory in GeV.
Definition: QCD.h:1021
QCD::MAls
double MAls
The mass scale in GeV at which the strong coupling measurement is provided.
Definition: QCD.h:1019
QCD::Mrun4
const 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:1531
QCD::CF
double CF
Definition: QCD.h:1026
QCD::MrunTMP
const double MrunTMP(const double mu_f, const double mu_i, const double m, const int nf, const orders order) const
A function to calculate the running of the mass between flavour thresholds.
Definition: QCD.cpp:1489
QCD::PostUpdate
virtual bool PostUpdate()
The post-update method for QCD.
Definition: QCD.cpp:158
QCD::BParameterMap
std::map< std::string, BParameter > BParameterMap
Definition: QCD.h:1030
QCD::Beta2
const double Beta2(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:611
QCD::computeBd
bool computeBd
Switch for computing from .
Definition: QCD.h:1036
QCD::NQCDvars
static const int NQCDvars
The number of model parameters in QCD.
Definition: QCD.h:359
QCD::logLambdaNLO_cache
double logLambdaNLO_cache[9][CacheSize]
Definition: QCD.h:1041
QCD::FlagCsi
bool FlagCsi
A flag to determine whether and or (false) and (default, true) are used as inputs.
Definition: QCD.h:1048
QCD::unknownParameterWarning
bool unknownParameterWarning
A flag to stop the unknown parameter warning after the first time.
Definition: QCD.h:1044
QCD::FlagMpole2MbarNumeric
bool FlagMpole2MbarNumeric
A flag to determine whether the pole mass to mass conversion is done numerically.
Definition: QCD.h:1015
QCD::zeta2
double zeta2
computed with the GSL.
Definition: QCD.h:1032
QCD::Nc
double Nc
The number of colours.
Definition: QCD.h:1025
QCD::zeta3
double zeta3
computed with the GSL.
Definition: QCD.h:1033
QCD::muc
double muc
The threshold between four- and three-flavour theory in GeV.
Definition: QCD.h:1023
QCD::ZeroNf6NLO
const 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:1047
QCD::Gamma0
const double Gamma0(const double nf) const
The coefficient used to compute the running of a mass.
Definition: QCD.cpp:1271
QCD::realorder
orders realorder
Definition: QCD.h:1049
QCD::AlsWithLambda
const double AlsWithLambda(const double mu, const orders order) const
Computes the running strong coupling in the scheme with the use of .
Definition: QCD.cpp:704
QCD::CA
double CA
Definition: QCD.h:1026
QCD::FlagMtPole
bool FlagMtPole
A flag to determine whether the pole mass of the top quark is used as input.
Definition: QCD.h:1014
QCD::PreUpdate
virtual bool PreUpdate()
The pre-update method for QCD.
Definition: QCD.cpp:130
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:365
QCD::Mp2MbarTMP
const 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:1617
QCD::MS2DRqmass
const double MS2DRqmass(const double MSscale, const double MSbar) const
Converts a quark mass from the scheme to the scheme.
Definition: QCD.cpp:1725
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:1045
QCD::mesonsMap
std::map< const QCD::meson, Meson > mesonsMap
The map of defined mesons.
Definition: QCD.h:1047
QCD::mrun_cache
double mrun_cache[11][CacheSize]
Cache for running quark mass.
Definition: QCD.h:1042
QCD::setFlag
virtual bool setFlag(const std::string name, const bool value)
A method to set a flag of QCD.
Definition: QCD.cpp:479
QCD::logLambda
const double logLambda(const double nf, orders order) const
Computes with nf flavours in GeV.
Definition: QCD.cpp:1229
QCD::Beta1
const double Beta1(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:606
QCD::Gamma1
const double Gamma1(const double nf) const
The coefficient used to compute the running of a mass.
Definition: QCD.cpp:1276
QCD::Mrun
const double Mrun(const double mu, const double m, const quark q, const orders order=FULLNNLO) const
Computes a running quark mass from .
Definition: QCD.cpp:1353
QCD::QCD
QCD()
Constructor.
Definition: QCD.cpp:29
QCD::quark
quark
An enum type for quarks.
Definition: QCD.h:323
QCD::UP
@ UP
Definition: QCD.h:324
QCD::BOTTOM
@ BOTTOM
Definition: QCD.h:329
QCD::TOP
@ TOP
Definition: QCD.h:328
QCD::DOWN
@ DOWN
Definition: QCD.h:325
QCD::STRANGE
@ STRANGE
Definition: QCD.h:327
QCD::CHARM
@ CHARM
Definition: QCD.h:326
QCD::ZeroNf5
const 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:1052
QCD::computeFBd
bool computeFBd
Switch for computing from .
Definition: QCD.h:1034
QCD::setFlagStr
virtual bool setFlagStr(const std::string name, const std::string value)
A method to set a flag of QCD.
Definition: QCD.cpp:511
QCD::CheckFlags
virtual bool CheckFlags() const
A method to check the sanity of the set of model flags.
Definition: QCD.cpp:517
QCD::Init
virtual bool Init(const std::map< std::string, double > &DPars)
Initializes the QCD parameters found in the argument.
Definition: QCD.cpp:120
QCD::Mofmu2Mbar
const double Mofmu2Mbar(const double m, const double mu, const quark q) const
Converts a quark running mass at an arbitrary scale to the corresponding mass .
Definition: QCD.cpp:1738
QCD::Mp2Mbar
const double Mp2Mbar(const double mp, const quark q, orders order=FULLNNLO) const
Converts a quark pole mass to the corresponding mass .
Definition: QCD.cpp:1672
QCD::BelowTh
const double BelowTh(const double mu) const
The active flavour threshold below the scale as defined in QCD::Thresholds().
Definition: QCD.cpp:559
QCD::logLambda5
const double logLambda5(orders order) const
for .
Definition: QCD.cpp:1067
QCD::TF
double TF
Definition: QCD.h:1026
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:343
QCD::Beta0
const double Beta0(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:601
QCD::als_cache
double als_cache[9][CacheSize]
Cache for .
Definition: QCD.h:1039
QCD::orderToString
const std::string orderToString(const orders order) const
Converts an object of the enum type "orders" to the corresponding string.
Definition: QCD.cpp:95
QCD::Mofmu2MbarTMP
const double Mofmu2MbarTMP(double *mu, double *params) const
The member used for finding the numerical solution to the mass from mass.
Definition: QCD.cpp:1730
QCD::Beta3
const double Beta3(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:618
QCD::computeFBp
bool computeFBp
Switch for computing from .
Definition: QCD.h:1035
QCD::dFdF_NA
double dFdF_NA
Definition: QCD.h:1026
QCD::computemt
bool computemt
Switch for computing the mass of the top quark.
Definition: QCD.h:1011
QCD::requireYd
bool requireYd
Switch for generating the Yukawa couplings to the down-type quarks.
Definition: QCD.h:1013
QCD::Thresholds
const double Thresholds(const int i) const
For accessing the active flavour threshold scales.
Definition: QCD.cpp:524
QCD::dAdA_NA
double dAdA_NA
Definition: QCD.h:1026
QCD::AboveTh
const double AboveTh(const double mu) const
The active flavour threshold above the scale as defined in QCD::Thresholds().
Definition: QCD.cpp:547
QCD::initializeBParameter
void initializeBParameter(std::string name_i) const
A method to initialize B Parameter and the corresponding meson.
Definition: QCD.cpp:211
QCD::QCDsuccess
bool QCDsuccess
Definition: QCD.h:995
QCD::Mp2Mbar_bar
const double Mp2Mbar_bar(const double mp, const quark q, const orders order=FULLNNLO) const
Converts the bottom pole mass to the corresponding mass by iteration.
Definition: QCD.cpp:1625
QCD::Nf
const double Nf(const double mu) const
The number of active flavour at scale .
Definition: QCD.cpp:571
QCD::AlsWithInit
const double AlsWithInit(const double mu, const double alsi, const double mu_i, const int nf, const orders order) const
Computes the running strong coupling from in the scheme, where it is forbidden to across a flavour...
Definition: QCD.cpp:627
QCD::NfThresholdCorrections
const 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:709
QCD::FullOrder
const orders FullOrder(orders order) const
Return the FULLORDER enum corresponding to order.
Definition: QCD.cpp:728
QCD::logLambda5_cache
double logLambda5_cache[4][CacheSize]
Definition: QCD.h:1040
QCD::getQuarks
const Particle & getQuarks(const QCD::quark q) const
A get method to access a quark as an object of the type Particle.
Definition: QCD.h:536
QCD::dFdA_NA
double dFdA_NA
Definition: QCD.h:1026
QCD::AlsByOrder
const double AlsByOrder(const double mu, const orders order=FULLNLO, bool Nf_thr=true) const
Definition: QCD.cpp:804
QCD::DeltaMass
double DeltaMass(double x) const
A method to compute the correction due to light quarks to the conversion between pole and mass.
Definition: QCD.cpp:1611
QCD::quarks
Particle quarks[6]
The vector of all SM quarks.
Definition: QCD.h:1027
QCD::Update
virtual bool Update(const std::map< std::string, double > &DPars)
The update method for QCD.
Definition: QCD.cpp:136
QCD::MassOfNf
const double MassOfNf(int nf) const
The Mbar mass of the heaviest quark in the theory with Nf active flavour.
Definition: QCD.cpp:745
QCD::Als
const double Als(const double mu, const orders order=FULLNLO, const bool Nf_thr=true) const
Definition: QCD.cpp:762
QCD::mtpole
double mtpole
The pole mass of the top quark.
Definition: QCD.h:1020
QCD::threCorrForMass
const 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:1286
QCD::Gamma2
const double Gamma2(const double nf) const
The coefficient used to compute the running of a mass.
Definition: QCD.cpp:1281
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:1046
QCD::initializeMeson
void initializeMeson(QCD::meson meson_i) const
A method to initialize a meson.
Definition: QCD.cpp:280
QCD::AlsM
double AlsM
The strong coupling constant at the mass scale MAls, .
Definition: QCD.h:1018
QCD::ZeroNf4NLO
const 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:1057
QCD::AlsOLD
const 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:871
QCD::computeBs
bool computeBs
Switch for computing from .
Definition: QCD.h:1037
QCD::mub
double mub
The threshold between five- and four-flavour theory in GeV.
Definition: QCD.h:1022
QCD::Mbar2Mp
const double Mbar2Mp(const double mbar, const quark q, const orders order=FULLNNLO) const
Converts the mass to the pole mass.
Definition: QCD.cpp:1552
QCD::NA
double NA
Definition: QCD.h:1026
QCD::CacheShift
void CacheShift(double cache[][5], int n) const
A member used to manage the caching for this class.
mtpole
A class for , the pole mass of the top quark.
Definition: masses.h:164
orders
orders
An enum type for orders in QCD.
Definition: OrderScheme.h:33