(Generated on Fri Feb 19 2016 13:25:08 by doxygen 1.8.9.1)
QCD.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2012 HEPfit Collaboration
3  * All rights reserved.
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 <gsl/gsl_sf.h>
18 #include <TF1.h>
19 #include <Math/WrappedTF1.h>
20 #include <Math/BrentRootFinder.h>
21 #include "QCD.h"
22 
23 const std::string QCD::QCDvars[NQCDvars] = {
24  "AlsM", "MAls",
25  "mup", "mdown", "mcharm", "mstrange", "mtop", "mbottom",
26  "muc", "mub", "mut",
27  "MK0", "MKp", "MD", "MBd", "MBp", "MBs", "MKstar", "Mphi",
28  "tKl", "tKp", "tBd", "tBp", "tBs", "tKstar", "tphi",
29  "DGs_Gs",
30  "FK", "FD", "FBs", "FBsoFBd", "FKstar", "FKstarp", "Fphi",
31  "BK1", "BK2", "BK3", "BK4", "BK5", "BKscale", "BKscheme",
32  "BD1", "BD2", "BD3", "BD4", "BD5", "BDscale", "BDscheme",
33  "BBsoBBd",
34  "BBs1", "BBs2", "BBs3", "BBs4", "BBs5", "BBsscale", "BBsscheme",
35  "BK(1/2)1", "BK(1/2)2", "BK(1/2)3", "BK(1/2)4", "BK(1/2)5",
36  "BK(1/2)6", "BK(1/2)7", "BK(1/2)8", "BK(1/2)9", "BK(1/2)10",
37  "BK(3/2)1", "BK(3/2)2", "BK(3/2)3", "BK(3/2)4", "BK(3/2)5",
38  "BK(3/2)6", "BK(3/2)7", "BK(3/2)8", "BK(3/2)9", "BK(3/2)10",
39  "BKd_scale", "BKd_scheme",
40  "ReA2_Kd", "ReA0_Kd", "Omega_eta_etap",
41  "Br_Kp_P0enu", "Br_Kp_munu", "Br_B_Xcenu", "DeltaP_cu", "IB_Kl", "IB_Kp",
42  "a_0V", "a_1V", "a_2V", "MRV", "a_0A0", "a_1A0", "a_2A0", "MRA0", "a_0A1", "a_1A1", "a_2A1", "MRA1", "a_0A12", "a_1A12", "a_2A12", "MRA12",
43  "a_0T1", "a_1T1", "a_2T1", "MRT1", "a_0T2", "a_1T2", "a_2T2", "MRT2", "a_0T23", "a_1T23", "a_2T23", "MRT23",
44  "a_0Vphi", "a_1Vphi", "a_2Vphi", "MRVphi", "a_0A0phi", "a_1A0phi", "a_2A0phi", "MRA0phi", "a_0A1phi", "a_1A1phi", "a_2A1phi", "MRA1phi", "a_0A12phi", "a_1A12phi", "a_2A12phi",
45  "MRA12phi", "a_0T1phi", "a_1T1phi", "a_2T1phi", "MRT1phi", "a_0T2phi", "a_1T2phi", "a_2T2phi", "MRT2phi", "a_0T23phi", "a_1T23phi", "a_2T23phi", "MRT23phi",
46  "absh_0", "absh_p", "absh_m", "argh_0", "argh_p", "argh_m",
47  "absh_0_1", "absh_p_1", "absh_m_1", "argh_0_1", "argh_p_1", "argh_m_1",
48  "absh_0_2", "absh_p_2", "absh_m_2", "argh_0_2", "argh_p_2", "argh_m_2",
49  "r_1_fplus", "r_2_fplus", "m_fit2_fplus", "r_1_fT", "r_2_fT", "m_fit2_fT", "r_2_f0", "m_fit2_f0",
50  "absh_0_MP", "argh_0_MP", "absh_0_1_MP", "argh_0_1_MP",
51  "bsgamma_E0", "BLNPcorr", "Gambino_mukin", "Gambino_BRsem", "Gambino_Mbkin", "Gambino_Mcatmuc", "Gambino_mupi2", "Gambino_rhoD3", "Gambino_muG2", "Gambino_rhoLS3",
52  "lambdaB", "alpha1kst", "alpha2kst"
53  //"r_2A0", "r_2T1", "r_2T2", "r_2A0phi", "r_2T1phi", "r_2T2phi"
54 };
55 
56 QCD::QCD()
57 : BBs(5), BBd(5), BD(5), BK(5), BKd1(10), BKd3(10)
58 {
59  Nc = 3.;
60  CF = Nc / 2. - 1. / (2. * Nc);
61  // Particle(std::string name, double mass, double mass_scale = 0., double width = 0., double charge = 0.,double isospin = 0.);
62  quarks[UP] = Particle("UP", 0., 2., 0., 2. / 3., .5);
63  quarks[CHARM] = Particle("CHARM", 0., 0., 0., 2. / 3., .5);
64  quarks[TOP] = Particle("TOP", 0., 0., 0., 2. / 3., .5);
65  quarks[DOWN] = Particle("DOWN", 0., 2., 0., -1. / 3., -.5);
66  quarks[STRANGE] = Particle("STRANGE", 0., 2., 0., -1. / 3., -.5);
67  quarks[BOTTOM] = Particle("BOTTOM", 0., 0., 0., -1. / 3., -.5);
68  zeta2 = gsl_sf_zeta_int(2);
69  zeta3 = gsl_sf_zeta_int(3);
70  for (int i = 0; i < CacheSize; i++) {
71  for (int j = 0; j < 8; j++)
72  als_cache[j][i] = 0.;
73  for (int j = 0; j < 4; j++)
74  logLambda5_cache[j][i] = 0.;
75  for (int j = 0; j < 10; j++)
76  mrun_cache[j][i] = 0.;
77  for (int j = 0; j < 5; j++)
78  mp2mbar_cache[j][i] = 0.;
79  }
80 
81  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("AlsM", boost::cref(AlsM)));
82  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MAls", boost::cref(MAls)));
83  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("mup", boost::cref(quarks[UP].getMass())));
84  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("mdown", boost::cref(quarks[DOWN].getMass())));
85  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("mcharm", boost::cref(quarks[CHARM].getMass())));
86  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("mstrange", boost::cref(quarks[STRANGE].getMass())));
87  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("mtop", boost::cref(mtpole)));
88  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("mbottom", boost::cref(quarks[BOTTOM].getMass())));
89  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("muc", boost::cref(muc)));
90  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("mub", boost::cref(mub)));
91  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("mut", boost::cref(mut)));
92  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MK0", boost::cref(mesons[K_0].getMass())));
93  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MKp", boost::cref(mesons[K_P].getMass())));
94  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MD", boost::cref(mesons[D_0].getMass())));
95  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MBd", boost::cref(mesons[B_D].getMass())));
96  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MBp", boost::cref(mesons[B_P].getMass())));
97  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MBs", boost::cref(mesons[B_S].getMass())));
98  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MKstar", boost::cref(mesons[K_star].getMass())));
99  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("Mphi", boost::cref(mesons[PHI].getMass())));
100  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("tKl", boost::cref(mesons[K_0].getWidth())));
101  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("tKp", boost::cref(mesons[K_P].getWidth())));
102  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("tBd", boost::cref(mesons[B_D].getWidth())));
103  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("tBp", boost::cref(mesons[B_P].getWidth())));
104  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("tBs", boost::cref(mesons[B_S].getWidth())));
105  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("DGs_Gs", boost::cref(mesons[B_S].getDgamma_gamma())));
106  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("tKstar", boost::cref(mesons[K_star].getWidth())));
107  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("tphi", boost::cref(mesons[PHI].getWidth())));
108  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("FK", boost::cref(mesons[K_0].getDecayconst())));
109  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("FD", boost::cref(mesons[D_0].getDecayconst())));
110  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("FBs", boost::cref(mesons[B_S].getDecayconst())));
111  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("FKstar", boost::cref(mesons[K_star].getDecayconst())));
112  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("FKstarp", boost::cref(FKstarp)));
113  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("Fphi", boost::cref(mesons[PHI].getDecayconst())));
114  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("FBsoFBd", boost::cref(FBsoFBd)));
115  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK1", boost::cref(BK.getBpars()(0))));
116  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK2", boost::cref(BK.getBpars()(1))));
117  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK3", boost::cref(BK.getBpars()(2))));
118  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK4", boost::cref(BK.getBpars()(3))));
119  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK5", boost::cref(BK.getBpars()(4))));
120  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BKscale", boost::cref(BK.getMu())));
121  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BD1", boost::cref(BD.getBpars()(0))));
122  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BD2", boost::cref(BD.getBpars()(1))));
123  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BD3", boost::cref(BD.getBpars()(2))));
124  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BD4", boost::cref(BD.getBpars()(3))));
125  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BD5", boost::cref(BD.getBpars()(4))));
126  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BDscale", boost::cref(BD.getMu())));
127  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BBsoBBd", boost::cref(BBsoBBd)));
128  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BBs1", boost::cref(BBs.getBpars()(0))));
129  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BBs2", boost::cref(BBs.getBpars()(1))));
130  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BBs3", boost::cref(BBs.getBpars()(2))));
131  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BBs4", boost::cref(BBs.getBpars()(3))));
132  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BBs5", boost::cref(BBs.getBpars()(4))));
133  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BBsscale", boost::cref(BBs.getMu())));
134  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(1/2)1", boost::cref(BKd1.getBpars()(0))));
135  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(1/2)2", boost::cref(BKd1.getBpars()(1))));
136  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(1/2)3", boost::cref(BKd1.getBpars()(2))));
137  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(1/2)4", boost::cref(BKd1.getBpars()(3))));
138  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(1/2)5", boost::cref(BKd1.getBpars()(4))));
139  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(1/2)6", boost::cref(BKd1.getBpars()(5))));
140  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(1/2)7", boost::cref(BKd1.getBpars()(6))));
141  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(1/2)8", boost::cref(BKd1.getBpars()(7))));
142  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(1/2)9", boost::cref(BKd1.getBpars()(8))));
143  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(1/2)10", boost::cref(BKd1.getBpars()(9))));
144  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(3/2)1", boost::cref(BKd3.getBpars()(0))));
145  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(3/2)2", boost::cref(BKd3.getBpars()(1))));
146  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(3/2)3", boost::cref(BKd3.getBpars()(2))));
147  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(3/2)4", boost::cref(BKd3.getBpars()(3))));
148  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(3/2)5", boost::cref(BKd3.getBpars()(4))));
149  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(3/2)6", boost::cref(BKd3.getBpars()(5))));
150  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(3/2)7", boost::cref(BKd3.getBpars()(6))));
151  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(3/2)8", boost::cref(BKd3.getBpars()(7))));
152  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(3/2)9", boost::cref(BKd3.getBpars()(8))));
153  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BK(3/2)10", boost::cref(BKd3.getBpars()(9))));
154  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BKd_scale", boost::cref(BKd1.getMu())));
155  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("ReA2_Kd", boost::cref(ReA2_Kd)));
156  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("ReA0_Kd", boost::cref(ReA0_Kd)));
157  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("Omega_eta_etap", boost::cref(Omega_eta_etap)));
158  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("Br_Kp_P0enu", boost::cref(Br_Kp_P0enu)));
159  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("Br_Kp_munu", boost::cref(Br_Kp_munu)));
160  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("Br_B_Xcenu", boost::cref(Br_B_Xcenu)));
161  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("DeltaP_cu", boost::cref(DeltaP_cu)));
162  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("IB_Kl", boost::cref(IB_Kl)));
163  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("IB_Kp", boost::cref(IB_Kp)));
164  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("absh_0", boost::cref(absh_0)));
165  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("absh_p", boost::cref(absh_p)));
166  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("absh_m", boost::cref(absh_m)));
167  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("argh_0", boost::cref(argh_0)));
168  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("argh_p", boost::cref(argh_p)));
169  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("argh_m", boost::cref(argh_m)));
170  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("absh_0_1", boost::cref(absh_0_1)));
171  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("absh_p_1", boost::cref(absh_p_1)));
172  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("absh_m_1", boost::cref(absh_m_1)));
173  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("argh_0_1", boost::cref(argh_0_1)));
174  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("argh_p_1", boost::cref(argh_p_1)));
175  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("argh_m_1", boost::cref(argh_m_1)));
176  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("absh_0_2", boost::cref(absh_0_2)));
177  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("absh_p_2", boost::cref(absh_p_2)));
178  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("absh_m_2", boost::cref(absh_m_2)));
179  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("argh_0_2", boost::cref(argh_0_2)));
180  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("argh_p_2", boost::cref(argh_p_2)));
181  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("argh_m_2", boost::cref(argh_m_2)));
182  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("absh_0_MP", boost::cref(absh_0_MP)));
183  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("argh_0_MP", boost::cref(argh_0_MP)));
184  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("absh_0_1_MP", boost::cref(absh_0_1_MP)));
185  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("argh_0_1_MP", boost::cref(argh_0_1_MP)));
186  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_0V", boost::cref(a_0V)));
187  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_1V", boost::cref(a_1V)));
188  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_2V", boost::cref(a_2V)));
189  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MRV", boost::cref(MRV)));
190  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_0A0", boost::cref(a_0A0)));
191  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_1A0", boost::cref(a_1A0)));
192  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_2A0", boost::cref(a_2A0)));
193  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MRA0", boost::cref(MRA0)));
194  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_0A1", boost::cref(a_0A1)));
195  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_1A1", boost::cref(a_1A1)));
196  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_2A1", boost::cref(a_2A1)));
197  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MRA1", boost::cref(MRA1)));
198  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_0A12", boost::cref(a_0A12)));
199  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_1A12", boost::cref(a_1A12)));
200  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_2A12", boost::cref(a_2A12)));
201  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MRA12", boost::cref(MRA12)));
202  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_0T1", boost::cref(a_0T1)));
203  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_1T1", boost::cref(a_1T1)));
204  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_2T1", boost::cref(a_2T1)));
205  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MRT1", boost::cref(MRT1)));
206  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_0T2", boost::cref(a_0T2)));
207  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_1T2", boost::cref(a_1T2)));
208  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_2T2", boost::cref(a_2T2)));
209  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MRT2", boost::cref(MRT2)));
210  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_0T23", boost::cref(a_0T23)));
211  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_1T23", boost::cref(a_1T23)));
212  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_2T23", boost::cref(a_2T23)));
213  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MRT23", boost::cref(MRT23)));
214  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_0Vphi", boost::cref(a_0Vphi)));
215  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_1Vphi", boost::cref(a_1Vphi)));
216  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_2Vphi", boost::cref(a_2Vphi)));
217  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MRVphi", boost::cref(MRVphi)));
218  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_0A0phi", boost::cref(a_0A0phi)));
219  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_1A0phi", boost::cref(a_1A0phi)));
220  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_2A0phi", boost::cref(a_2A0phi)));
221  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MRA0phi", boost::cref(MRA0phi)));
222  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_0A1phi", boost::cref(a_0A1phi)));
223  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_1A1phi", boost::cref(a_1A1phi)));
224  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_2A1phi", boost::cref(a_2A1phi)));
225  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MRA1phi", boost::cref(MRA1phi)));
226  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_0A12phi", boost::cref(a_0A12phi)));
227  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_1A12phi", boost::cref(a_1A12phi)));
228  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_2A12phi", boost::cref(a_2A12phi)));
229  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MRA12phi", boost::cref(MRA12phi)));
230  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_0T1phi", boost::cref(a_0T1phi)));
231  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_1T1phi", boost::cref(a_1T1phi)));
232  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_2T1phi", boost::cref(a_2T1phi)));
233  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MRT1phi", boost::cref(MRT1phi)));
234  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_0T2phi", boost::cref(a_0T2phi)));
235  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_1T2phi", boost::cref(a_1T2phi)));
236  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_2T2phi", boost::cref(a_2T2phi)));
237  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MRT2phi", boost::cref(MRT2phi)));
238  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_0T23phi", boost::cref(a_0T23phi)));
239  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_1T23phi", boost::cref(a_1T23phi)));
240  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("a_2T23phi", boost::cref(a_2T23phi)));
241  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("MRT23phi", boost::cref(MRT23phi)));
242  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("r_1_fplus", boost::cref(r_1_fplus)));
243  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("r_2_fplus", boost::cref(r_2_fplus)));
244  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("m_fit2_fplus", boost::cref(m_fit2_fplus)));
245  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("r_1_fT", boost::cref(r_1_fT)));
246  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("r_2_fT", boost::cref(r_2_fT)));
247  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("m_fit2_fT", boost::cref(m_fit2_fT)));
248  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("r_2_f0", boost::cref(r_2_f0)));
249  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("m_fit2_f0", boost::cref(m_fit2_f0)));
250  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("bsgamma_E0", boost::cref(bsgamma_E0)));
251  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("BLNPcorr", boost::cref(BLNPcorr)));
252  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("Gambino_mukin", boost::cref(Gambino_mukin)));
253  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("Gambino_BRsem", boost::cref(Gambino_BRsem)));
254  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("Gambino_Mbkin", boost::cref(Gambino_Mbkin)));
255  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("Gambino_Mcatmuc", boost::cref(Gambino_Mcatmuc)));
256  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("Gambino_mupi2", boost::cref(Gambino_mupi2)));
257  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("Gambino_rhoD3", boost::cref(Gambino_rhoD3)));
258  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("Gambino_muG2", boost::cref(Gambino_muG2)));
259  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("Gambino_rhoLS3", boost::cref(Gambino_rhoLS3)));
260  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("lambdaB", boost::cref(mesons[B_D].getLambdaM())));
261  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("alpha1kst", boost::cref(mesons[K_star].getGegenalpha(0))));
262  ModelParamMap.insert(std::pair<std::string, boost::reference_wrapper<const double> >("alpha2kst", boost::cref(mesons[K_star].getGegenalpha(1))));
263 
264  unknownParameterWarning = true;
265 }
266 
267 std::string QCD::orderToString(const orders order) const
268 {
269  switch (order) {
270  case LO:
271  return "LO";
272  case NLO:
273  return "NLO";
274  case FULLNLO:
275  return "FULLNLO";
276  case NNLO:
277  return "NNLO";
278  case FULLNNLO:
279  return "FULLNNLO";
280  default:
281  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::orderToString().");
282  }
283 }
284 
286 
287 bool QCD::Init(const std::map<std::string, double>& DPars)
288 {
289  bool check = CheckParameters(DPars);
290  if (!check) return (check);
291  check *= Update(DPars);
292  unknownParameterWarning = false;
293  return (check);
294 }
295 
296 bool QCD::PreUpdate()
297 {
298  requireYu = false;
299  requireYd = false;
300  computeBd = false;
301  computeFBd = false;
302  computemt = false;
303 
304  return (true);
305 }
306 
307 bool QCD::Update(const std::map<std::string, double>& DPars)
308 {
309  if (!PreUpdate()) return (false);
310 
311  UpdateError = false;
312 
313  for (std::map<std::string, double>::const_iterator it = DPars.begin(); it != DPars.end(); it++)
314  setParameter(it->first, it->second);
315 
316  if (UpdateError) return (false);
317 
318  if (!PostUpdate()) return (false);
319 
320  return (true);
321 }
322 
323 bool QCD::PostUpdate()
324 {
325  if (computeFBd)
326  mesons[B_D].setDecayconst(mesons[B_S].getDecayconst() / FBsoFBd);
327  if (computeBd)
328  BBd.setBpars(0, BBs.getBpars()(0) / BBsoBBd);
329  if (computemt) {
330  quarks[TOP].setMass(Mp2Mbar(mtpole, FULLNNLO));
331  quarks[TOP].setMass_scale(quarks[TOP].getMass());
332  }
333 
334  myh_0 = gslpp::complex(absh_0,argh_0,true);
335  myh_p = gslpp::complex(absh_p,argh_p,true);
336  myh_m = gslpp::complex(absh_m,argh_m,true);
337  myh_0_1 = gslpp::complex(absh_0_1,argh_0_1,true);
338  myh_p_1 = gslpp::complex(absh_p_1,argh_p_1,true);
339  myh_m_1 = gslpp::complex(absh_m_1,argh_m_1,true);
340  myh_0_2 = gslpp::complex(absh_0_2,argh_0_2,true);
341  myh_p_2 = gslpp::complex(absh_p_2,argh_p_2,true);
342  myh_m_2 = gslpp::complex(absh_m_2,argh_m_2,true);
343 
344  return (true);
345 }
346 
347 void QCD::setParameter(const std::string name, const double& value)
348 {
349  if (name.compare("AlsM") == 0) {
350  AlsM = value;
351  computemt = true;
352  requireYu = true;
353  requireYd = true;
354  } else if (name.compare("MAls") == 0) {
355  MAls = value;
356  computemt = true;
357  requireYu = true;
358  requireYd = true;
359  } else if (name.compare("mup") == 0) {
360  if (value < MEPS) UpdateError = true;
361  quarks[UP].setMass(value);
362  requireYu = true;
363  } else if (name.compare("mdown") == 0) {
364  if (value < MEPS) UpdateError = true;
365  quarks[DOWN].setMass(value);
366  requireYd = true;
367  } else if (name.compare("mcharm") == 0) {
368  quarks[CHARM].setMass(value);
369  quarks[CHARM].setMass_scale(value);
370  requireYu = true;
371  } else if (name.compare("mstrange") == 0) {
372  if (value < MEPS) UpdateError = true;
373  quarks[STRANGE].setMass(value);
374  requireYd = true;
375  } else if (name.compare("mtop") == 0) {
376  mtpole = value;
377  requireYu = true;
378  computemt = true;
379  } else if (name.compare("mbottom") == 0) {
380  quarks[BOTTOM].setMass(value);
381  quarks[BOTTOM].setMass_scale(value);
382  requireYd = true;
383  } else if (name.compare("muc") == 0)
384  muc = value;
385  else if (name.compare("mub") == 0)
386  mub = value;
387  else if (name.compare("mut") == 0)
388  mut = value;
389  else if (name.compare("MK0") == 0)
390  mesons[K_0].setMass(value);
391  else if (name.compare("MKp") == 0)
392  mesons[K_P].setMass(value);
393  else if (name.compare("MD") == 0)
394  mesons[D_0].setMass(value);
395  else if (name.compare("MBd") == 0)
396  mesons[B_D].setMass(value);
397  else if (name.compare("MBp") == 0)
398  mesons[B_P].setMass(value);
399  else if (name.compare("MBs") == 0)
400  mesons[B_S].setMass(value);
401  else if (name.compare("MKstar") == 0)
402  mesons[K_star].setMass(value);
403  else if (name.compare("Mphi") == 0)
404  mesons[PHI].setMass(value);
405  else if (name.compare("tKl") == 0)
406  mesons[K_0].setLifetime(value);
407  else if (name.compare("tKp") == 0)
408  mesons[K_P].setLifetime(value);
409  else if (name.compare("tBd") == 0)
410  mesons[B_D].setLifetime(value);
411  else if (name.compare("tBp") == 0)
412  mesons[B_P].setLifetime(value);
413  else if (name.compare("tBs") == 0)
414  mesons[B_S].setLifetime(value);
415  else if (name.compare("DGs_Gs") == 0)
416  mesons[B_S].setDgamma_gamma(value);
417  else if (name.compare("tKstar") == 0)
418  mesons[K_star].setLifetime(value);
419  else if (name.compare("tphi") == 0)
420  mesons[PHI].setLifetime(value);
421  //else if(name.compare("FP")==0) {
422  // mesons[P_0].setDecayconst(value);
423  // mesons[P_P].setDecayconst(value);
424  //}
425  else if (name.compare("FK") == 0)
426  mesons[K_0].setDecayconst(value);
427  else if (name.compare("FD") == 0)
428  mesons[D_0].setDecayconst(value);
429  else if (name.compare("FKstar") == 0)
430  mesons[K_star].setDecayconst(value);
431  else if (name.compare("FKstarp") == 0)
432  FKstarp = value;
433  else if (name.compare("Fphi") == 0)
434  mesons[PHI].setDecayconst(value);
435  else if (name.compare("FBs") == 0) {
436  mesons[B_S].setDecayconst(value);
437  computeFBd = true;
438  } else if (name.compare("FBsoFBd") == 0) {
439  FBsoFBd = value;
440  computeFBd = true;
441  } else if (name.compare("BK1") == 0) {
442  BK.setBpars(0, value);
443  BKB0 = value;
444  } else if (name.compare("BK2") == 0) {
445  BK.setBpars(1, value);
446  BKB1 = value;
447  } else if (name.compare("BK3") == 0) {
448  BK.setBpars(2, value);
449  BKB2 = value;
450  } else if (name.compare("BK4") == 0) {
451  BK.setBpars(3, value);
452  BKB3 = value;
453  } else if (name.compare("BK5") == 0) {
454  BK.setBpars(4, value);
455  BKB4 = value;
456  } else if (name.compare("BKscale") == 0) {
457  BK.setMu(value);
458  BKscale = value;
459  } else if (name.compare("BKscheme") == 0)
460  BK.setScheme((schemes) value);
461  else if (name.compare("BD1") == 0) {
462  BD.setBpars(0, value);
463  BDB0 = value;
464  } else if (name.compare("BD2") == 0) {
465  BD.setBpars(1, value);
466  BDB1 = value;
467  } else if (name.compare("BD3") == 0) {
468  BD.setBpars(2, value);
469  BDB2 = value;
470  } else if (name.compare("BD4") == 0) {
471  BD.setBpars(3, value);
472  BDB3 = value;
473  } else if (name.compare("BD5") == 0) {
474  BD.setBpars(4, value);
475  BDB4 = value;
476  } else if (name.compare("BDscale") == 0) {
477  BD.setMu(value);
478  BDscale = value;
479  } else if (name.compare("BDscheme") == 0)
480  BD.setScheme((schemes) value);
481  else if (name.compare("BBsoBBd") == 0) {
482  BBsoBBd = value;
483  computeBd = true;
484  } else if (name.compare("BBs1") == 0) {
485  BBs.setBpars(0, value);
486  BBsB0 = value;
487  computeBd = true;
488  } else if (name.compare("BBs2") == 0) {
489  BBd.setBpars(1, value);
490  BBs.setBpars(1, value);
491  BBsB1 = value;
492  //BBdB1 = value;
493  } else if (name.compare("BBs3") == 0) {
494  BBd.setBpars(2, value);
495  BBs.setBpars(2, value);
496  BBsB2 = value;
497  //BBdB2 = value;
498  } else if (name.compare("BBs4") == 0) {
499  BBd.setBpars(3, value);
500  BBs.setBpars(3, value);
501  BBsB3 = value;
502  //BBdB3 = value;
503  } else if (name.compare("BBs5") == 0) {
504  BBd.setBpars(4, value);
505  BBs.setBpars(4, value);
506  BBsB4 = value;
507  //BBdB4 = value;
508  } else if (name.compare("BBsscale") == 0) {
509  BBd.setMu(value);
510  BBs.setMu(value);
511  BBsscale = value;
512  //BBdscale = value;
513  } else if (name.compare("BBsscheme") == 0) {
514  BBd.setScheme((schemes) value);
515  BBs.setScheme((schemes) value);
516  } else if (name.compare("BK(1/2)1") == 0) {
517  BKd1.setBpars(0, value);
518  BKd1B0 = value;
519  } else if (name.compare("BK(1/2)2") == 0) {
520  BKd1.setBpars(1, value);
521  BKd1B1 = value;
522  } else if (name.compare("BK(1/2)3") == 0) {
523  BKd1.setBpars(2, value);
524  BKd1B2 = value;
525  } else if (name.compare("BK(1/2)4") == 0) {
526  BKd1.setBpars(3, value);
527  BKd1B3 = value;
528  } else if (name.compare("BK(1/2)5") == 0) {
529  BKd1.setBpars(4, value);
530  BKd1B4 = value;
531  } else if (name.compare("BK(1/2)6") == 0){
532  BKd1.setBpars(5, value);
533  BKd1B5 = value;
534  } else if (name.compare("BK(1/2)7") == 0){
535  BKd1.setBpars(6, value);
536  BKd1B6 = value;
537  } else if (name.compare("BK(1/2)8") == 0) {
538  BKd1.setBpars(7, value);
539  BKd1B7 = value;
540  } else if (name.compare("BK(1/2)9") == 0) {
541  BKd1.setBpars(8, value);
542  BKd1B8 = value;
543  } else if (name.compare("BK(1/2)10") == 0) {
544  BKd1.setBpars(9, value);
545  BKd1B9 = value;
546  } else if (name.compare("BK(3/2)1") == 0) {
547  BKd3.setBpars(0, value);
548  BKd3B0 = value;
549  } else if (name.compare("BK(3/2)2") == 0) {
550  BKd3.setBpars(1, value);
551  BKd3B1 = value;
552  } else if (name.compare("BK(3/2)3") == 0) {
553  BKd3.setBpars(2, value);
554  BKd3B2 = value;
555  } else if (name.compare("BK(3/2)4") == 0) {
556  BKd3.setBpars(3, value);
557  BKd3B3 = value;
558  } else if (name.compare("BK(3/2)5") == 0) {
559  BKd3.setBpars(4, value);
560  BKd3B4 = value;
561  } else if (name.compare("BK(3/2)6") == 0) {
562  BKd3.setBpars(5, value);
563  BKd3B5 = value;
564  } else if (name.compare("BK(3/2)7") == 0) {
565  BKd3.setBpars(6, value);
566  BKd3B6 = value;
567  } else if (name.compare("BK(3/2)8") == 0) {
568  BKd3.setBpars(7, value);
569  BKd3B7 = value;
570  } else if (name.compare("BK(3/2)9") == 0) {
571  BKd3.setBpars(8, value);
572  BKd3B8 = value;
573  } else if (name.compare("BK(3/2)10") == 0) {
574  BKd3.setBpars(9, value);
575  BKd3B9 = value;
576  } else if (name.compare("BKd_scale") == 0) {
577  BKd1.setMu(value);
578  BKd3.setMu(value);
579  BKd_scale = value;
580  } else if (name.compare("BKd_scheme") == 0) {
581  BKd1.setScheme((schemes) value);
582  BKd3.setScheme((schemes) value);
583  } else if (name.compare("ReA0_Kd") == 0)
584  ReA0_Kd = value;
585  else if (name.compare("ReA2_Kd") == 0)
586  ReA2_Kd = value;
587  else if (name.compare("Omega_eta_etap") == 0)
588  Omega_eta_etap = value;
589  else if (name.compare("Br_Kp_P0enu") == 0)
590  Br_Kp_P0enu = value;
591  else if (name.compare("Br_Kp_munu") == 0)
592  Br_Kp_munu = value;
593  else if (name.compare("Br_B_Xcenu") == 0)
594  Br_B_Xcenu = value;
595  else if (name.compare("DeltaP_cu") == 0)
596  DeltaP_cu = value;
597  else if (name.compare("IB_Kl") == 0)
598  IB_Kl = value;
599  else if (name.compare("IB_Kp") == 0)
600  IB_Kp = value;
601  else if (name.compare("absh_0") == 0)
602  absh_0 = value;
603  else if (name.compare("absh_p") == 0)
604  absh_p = value;
605  else if (name.compare("absh_m") == 0)
606  absh_m = value;
607  else if (name.compare("argh_0") == 0)
608  argh_0 = value;
609  else if (name.compare("argh_p") == 0)
610  argh_p = value;
611  else if (name.compare("argh_m") == 0)
612  argh_m = value;
613  else if (name.compare("absh_0_1") == 0)
614  absh_0_1 = value;
615  else if (name.compare("absh_p_1") == 0)
616  absh_p_1 = value;
617  else if (name.compare("absh_m_1") == 0)
618  absh_m_1 = value;
619  else if (name.compare("argh_0_1") == 0)
620  argh_0_1 = value;
621  else if (name.compare("argh_p_1") == 0)
622  argh_p_1 = value;
623  else if (name.compare("argh_m_1") == 0)
624  argh_m_1 = value;
625  else if (name.compare("absh_0_2") == 0)
626  absh_0_2 = value;
627  else if (name.compare("absh_p_2") == 0)
628  absh_p_2 = value;
629  else if (name.compare("absh_m_2") == 0)
630  absh_m_2 = value;
631  else if (name.compare("argh_0_2") == 0)
632  argh_0_2 = value;
633  else if (name.compare("argh_p_2") == 0)
634  argh_p_2 = value;
635  else if (name.compare("argh_m_2") == 0)
636  argh_m_2 = value;
637  else if (name.compare("absh_0_MP") == 0)
638  absh_0_MP = value;
639  else if (name.compare("argh_0_MP") == 0)
640  argh_0_MP = value;
641  else if (name.compare("absh_0_1_MP") == 0)
642  absh_0_1_MP = value;
643  else if (name.compare("argh_0_1_MP") == 0)
644  argh_0_1_MP = value;
645  else if (name.compare("a_0V") == 0)
646  a_0V = value;
647  else if (name.compare("a_1V") == 0)
648  a_1V = value;
649  else if (name.compare("a_2V") == 0)
650  a_2V = value;
651  else if (name.compare("MRV") == 0)
652  MRV = value;
653  else if (name.compare("a_0A0") == 0)
654  a_0A0 = value;
655  else if (name.compare("a_1A0") == 0)
656  a_1A0 = value;
657  else if (name.compare("a_2A0") == 0)
658  a_2A0 = value;
659  else if (name.compare("MRA0") == 0)
660  MRA0 = value;
661  else if (name.compare("a_0A1") == 0)
662  a_0A1 = value;
663  else if (name.compare("a_1A1") == 0)
664  a_1A1 = value;
665  else if (name.compare("a_2A1") == 0)
666  a_2A1 = value;
667  else if (name.compare("MRA1") == 0)
668  MRA1 = value;
669  else if (name.compare("a_0A12") == 0)
670  a_0A12 = value;
671  else if (name.compare("a_1A12") == 0)
672  a_1A12 = value;
673  else if (name.compare("a_2A12") == 0)
674  a_2A12 = value;
675  else if (name.compare("MRA12") == 0)
676  MRA12 = value;
677  else if (name.compare("a_0T1") == 0)
678  a_0T1 = value;
679  else if (name.compare("a_1T1") == 0)
680  a_1T1 = value;
681  else if (name.compare("a_2T1") == 0)
682  a_2T1 = value;
683  else if (name.compare("MRT1") == 0)
684  MRT1 = value;
685  else if (name.compare("a_0T2") == 0)
686  a_0T2 = value;
687  else if (name.compare("a_1T2") == 0)
688  a_1T2 = value;
689  else if (name.compare("a_2T2") == 0)
690  a_2T2 = value;
691  else if (name.compare("MRT2") == 0)
692  MRT2 = value;
693  else if (name.compare("a_0T23") == 0)
694  a_0T23 = value;
695  else if (name.compare("a_1T23") == 0)
696  a_1T23 = value;
697  else if (name.compare("a_2T23") == 0)
698  a_2T23 = value;
699  else if (name.compare("MRT23") == 0)
700  MRT23 = value;
701  else if (name.compare("a_0Vphi") == 0)
702  a_0Vphi = value;
703  else if (name.compare("a_1Vphi") == 0)
704  a_1Vphi = value;
705  else if (name.compare("a_2Vphi") == 0)
706  a_2Vphi = value;
707  else if (name.compare("MRVphi") == 0)
708  MRVphi = value;
709  else if (name.compare("a_0A0phi") == 0)
710  a_0A0phi = value;
711  else if (name.compare("a_1A0phi") == 0)
712  a_1A0phi = value;
713  else if (name.compare("a_2A0phi") == 0)
714  a_2A0phi = value;
715  else if (name.compare("MRA0phi") == 0)
716  MRA0phi = value;
717  else if (name.compare("a_0A1phi") == 0)
718  a_0A1phi = value;
719  else if (name.compare("a_1A1phi") == 0)
720  a_1A1phi = value;
721  else if (name.compare("a_2A1phi") == 0)
722  a_2A1phi = value;
723  else if (name.compare("MRA1phi") == 0)
724  MRA1phi = value;
725  else if (name.compare("a_0A12phi") == 0)
726  a_0A12phi = value;
727  else if (name.compare("a_1A12phi") == 0)
728  a_1A12phi = value;
729  else if (name.compare("a_2A12phi") == 0)
730  a_2A12phi = value;
731  else if (name.compare("MRA12phi") == 0)
732  MRA12phi = value;
733  else if (name.compare("a_0T1phi") == 0)
734  a_0T1phi = value;
735  else if (name.compare("a_1T1phi") == 0)
736  a_1T1phi = value;
737  else if (name.compare("a_2T1phi") == 0)
738  a_2T1phi = value;
739  else if (name.compare("MRT1phi") == 0)
740  MRT1phi = value;
741  else if (name.compare("a_0T2phi") == 0)
742  a_0T2phi = value;
743  else if (name.compare("a_1T2phi") == 0)
744  a_1T2phi = value;
745  else if (name.compare("a_2T2phi") == 0)
746  a_2T2phi = value;
747  else if (name.compare("MRT2phi") == 0)
748  MRT2phi = value;
749  else if (name.compare("a_0T23phi") == 0)
750  a_0T23phi = value;
751  else if (name.compare("a_1T23phi") == 0)
752  a_1T23phi = value;
753  else if (name.compare("a_2T23phi") == 0)
754  a_2T23phi = value;
755  else if (name.compare("MRT23phi") == 0)
756  MRT23phi = value;
757  else if (name.compare("r_1_fplus") == 0)
758  r_1_fplus = value;
759  else if (name.compare("r_2_fplus") == 0)
760  r_2_fplus = value;
761  else if (name.compare("m_fit2_fplus") == 0)
762  m_fit2_fplus = value;
763  else if (name.compare("r_1_fT") == 0)
764  r_1_fT = value;
765  else if (name.compare("r_2_fT") == 0)
766  r_2_fT = value;
767  else if (name.compare("m_fit2_fT") == 0)
768  m_fit2_fT = value;
769  else if (name.compare("r_2_f0") == 0)
770  r_2_f0 = value;
771  else if (name.compare("m_fit2_f0") == 0)
772  m_fit2_f0 = value;
773  else if (name.compare("bsgamma_E0") == 0)
774  bsgamma_E0 = value;
775  else if (name.compare("BLNPcorr") == 0)
776  BLNPcorr = value;
777  else if (name.compare("Gambino_mukin") == 0)
778  Gambino_mukin = value;
779  else if (name.compare("Gambino_BRsem") == 0)
780  Gambino_BRsem = value;
781  else if (name.compare("Gambino_Mbkin") == 0)
782  Gambino_Mbkin = value;
783  else if (name.compare("Gambino_Mcatmuc") == 0)
784  Gambino_Mcatmuc = value;
785  else if (name.compare("Gambino_mupi2") == 0)
786  Gambino_mupi2 = value;
787  else if (name.compare("Gambino_rhoD3") == 0)
788  Gambino_rhoD3 = value;
789  else if (name.compare("Gambino_muG2") == 0)
790  Gambino_muG2 = value;
791  else if (name.compare("Gambino_rhoLS3") == 0)
792  Gambino_rhoLS3 = value;
793  else if (name.compare("lambdaB") == 0)
794  mesons[B_D].setLambdaM(value);
795  else if (name.compare("alpha1kst") == 0)
796  mesons[K_star].setGegenalpha(0,value);
797  else if (name.compare("alpha2kst") == 0)
798  mesons[K_star].setGegenalpha(1,value);
799  else
800  if (unknownParameterWarning)
801  std::cout << "WARNING: unknown parameter " << name << " in model initialization" << std::endl;
802 }
803 
804 bool QCD::CheckParameters(const std::map<std::string, double>& DPars)
805 {
806  for (int i = 0; i < NQCDvars; i++)
807  if (DPars.find(QCDvars[i]) == DPars.end()) {
808  std::cout << "missing mandatory QCD parameter " << QCDvars[i] << std::endl;
809  return false;
810  }
811  return true;
812 }
813 
815 
816 bool QCD::setFlag(const std::string name, const bool value)
817 {
818  std::cout << "WARNING: wrong name or value for ModelFlag " << name << std::endl;
819  return (false);
820 }
821 
822 bool QCD::setFlagStr(const std::string name, const std::string value)
823 {
824  std::cout << "WARNING: wrong name or value for ModelFlag " << name << std::endl;
825  return (false);
826 }
827 
828 bool QCD::CheckFlags() const
829 {
830  return (true);
831 }
832 
834 
835 double QCD::Thresholds(const int i) const
836 {
837  if (!(mut > mub && mub > muc))
838  throw std::runtime_error("inverted thresholds in QCD::Thresholds()!");
839 
840  switch (i) {
841  case 0: return (1.0E10);
842  case 1: return (mut);
843  case 2: return (mub);
844  case 3: return (muc);
845  default: return (0.);
846  }
847 }
848 
849 double QCD::AboveTh(const double mu) const
850 {
851  int i;
852  for (i = 4; i >= 0; i--)
853  if (mu < Thresholds(i)) return (Thresholds(i));
854 
855  throw std::runtime_error("Error in QCD::AboveTh()");
856 }
857 
858 double QCD::BelowTh(const double mu) const
859 {
860  int i;
861  for (i = 0; i < 5; i++)
862  if (mu >= Thresholds(i)) return (Thresholds(i));
863 
864  throw std::runtime_error("Error in QCD::BelowTh()");
865 }
866 
867 double QCD::Nf(const double mu) const
868 {
869  int i;
870  for (i = 1; i < 5; i++)
871  if (mu >= Thresholds(i))
872  return (7. - (double) i);
873 
874  throw std::runtime_error("Error in QCD::Nf()");
875 }
876 
877 void QCD::CacheShift(double cache[][CacheSize], int n) const
878 {
879  int i, j;
880  for (i = CacheSize - 1; i > 0; i--)
881  for (j = 0; j < n; j++)
882  cache[j][i] = cache[j][i - 1];
883 }
884 
886 
887 double QCD::Beta0(const double nf) const
888 {
889  return ( (11. * Nc - 2. * nf) / 3.);
890 }
891 
892 double QCD::Beta1(const double nf) const
893 {
894  return ( 34. / 3. * Nc * Nc - 10. / 3. * Nc * nf - 2. * CF * nf);
895 }
896 
897 double QCD::Beta2(const double nf) const
898 {
899  return ( 2857. / 54. * Nc * Nc * Nc + CF * CF * nf - 205. / 18. * CF * Nc * nf
900  - 1415. / 54. * Nc * Nc * nf + 11. / 9. * CF * nf * nf + 79. / 54. * Nc * nf * nf);
901 }
902 
903 double QCD::AlsWithInit(const double mu, const double alsi, const double mu_i,
904  const orders order) const
905 {
906  double nf = Nf(mu);
907  if (nf != Nf(mu_i))
908  throw std::runtime_error("Error in QCD::AlsWithInit().");
909 
910  double v = 1. - Beta0(nf) * alsi / 2. / M_PI * log(mu_i / mu);
911  switch (order) {
912  case LO:
913  return (alsi / v);
914  case FULLNLO:
915  return (alsi / v * (1. - Beta1(nf) / Beta0(nf) * alsi / 4. / M_PI * log(v) / v));
916  case NLO:
917  return (alsi / v * (-Beta1(nf) / Beta0(nf) * alsi / 4. / M_PI * log(v) / v));
918  default:
919  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Als(mu,alsi,mi,order).");
920  }
921 }
922 
923 double QCD::Als4(const double mu) const
924 {
925  double v = 1. - Beta0(4.) * AlsM / 2. / M_PI * log(MAls / mu);
926  return (AlsM / v * (1. - Beta1(4.) / Beta0(4.) * AlsM / 4. / M_PI * log(v) / v));
927 }
928 
929 double QCD::Alstilde5(const double mu) const
930 {
931  double mu_0 = MAls;
932  double alphatilde_e = alphaMz()/4./M_PI;
933  double alphatilde_s = AlsM/4./M_PI;
934  unsigned int nf = 5;
935 
936  double B00S = Beta0(nf), B10S = Beta1(nf), B20S = Beta2(nf), B30S = gsl_sf_zeta_int(3) * 352864./81. - 598391./1458,
937  B01S = -22./9., B11S = -308./27., B02S = 4945./243.;
938 
939  double B00E = 80./9., B01E = 176./9., B10E = 464./27.;
940 
941  double B10soB00s = B10S / B00S;
942  double B01soB00e = B01S/B00E;
943 
944  double vs= 1. + 2. * B00S * alphatilde_s * log(mu/ mu_0);
945  double ve= 1. - 2. * B00E * alphatilde_e * log(mu/ mu_0);
946  double ps= B00S * alphatilde_s /(B00S * alphatilde_s + B00E * alphatilde_e);
947 
948  double logve = log(ve);
949  double logvs = log(vs);
950  double logeos = log(ve/vs);
951  double logsoe = log(vs/ve);
952  double asovs = alphatilde_s/vs;
953  double aeove = alphatilde_e/ve;
954 
955  double result = 0;
956 
957  result = asovs - pow(asovs, 2) * (logvs * B10soB00s - logve * B01soB00e)
958  + pow(asovs, 3) * ((1. - vs) * B20S / B00S + B10soB00s * B10soB00s * (logvs * logvs - logvs
959  + vs - 1.) + B01soB00e * B01soB00e * logve * logve + (-2. * logvs * logve
960  + ps * ve * logve) * B01S * B10S/(B00E * B00S))
961  + pow(asovs, 4) * (0.5 * B30S *(1. - vs * vs)/ B00S + ((2. * vs - 3.) * logvs + vs * vs
962  - vs) * B20S * B10soB00s /(B00S) + B10soB00s * B10soB00s * B10soB00s * (- pow(logvs,3)
963  + 5. * pow(logvs,2) / 2. + 2. * (1. - vs) * logvs - (vs - 1.) * (vs - 1.)* 0.5))
964  + pow(asovs, 2) * (aeove) * ((ve - 1.) * B02S / B00E
965  + ps * ve * logeos * B11S /B00S +(logve - ve + 1.) * B01soB00e * B10E/(B00S)
966  + logvs * ve * ps * B01S * B10soB00s/(B00S) +(logsoe * ve * ps - logvs) * B01soB00e * B01E/( B00S));
967  return (result);
968 }
969 
970 double QCD::AlsWithLambda(const double mu, const double logLambda,
971  const orders order) const
972 {
973  double nf = Nf(mu);
974  double L = 2. * (log(mu) - logLambda);
975 
976  // LO contribution
977  double b0 = Beta0(nf);
978  double b0L = b0*L;
979  double alsLO = 4. * M_PI / b0L;
980  if (order == LO) return alsLO;
981 
982  // NLO contribution
983  double b1 = Beta1(nf);
984  double log_L = log(L);
985  double alsNLO = 4. * M_PI / b0L * (-b1 * log_L / b0 / b0L);
986  if (order == NLO) return alsNLO;
987  if (order == FULLNLO) return (alsLO + alsNLO);
988 
989  // NNLO contribution
990  double b2 = Beta2(nf);
991  double alsNNLO = 4. * M_PI / b0L * (1. / b0L / b0L
992  * (b1 * b1 / b0 / b0 * (log_L * log_L - log_L - 1.) + b2 / b0));
993  if (order == NNLO) return alsNNLO;
994  if (order == FULLNNLO) return (alsLO + alsNLO + alsNNLO);
995 
996  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::AlsWithLambda().");
997 }
998 
999 double QCD::AlsWithLambda(const double mu, const orders order) const
1000 {
1001  return AlsWithLambda(mu, logLambda(Nf(mu), order), order);
1002 }
1003 
1004 double QCD::Als(const double mu, const orders order) const
1005 {
1006  int i;
1007  for (i = 0; i < CacheSize; ++i)
1008  if ((mu == als_cache[0][i]) && ((double) order == als_cache[1][i]) &&
1009  (AlsM == als_cache[2][i]) && (MAls == als_cache[3][i]) &&
1010  (mut == als_cache[4][i]) && (mub == als_cache[5][i]) &&
1011  (muc == als_cache[6][i]))
1012  return als_cache[7][i];
1013 
1014  double nfmu = Nf(mu), nfz = Nf(MAls), mu_thre1, mu_thre2, Als_tmp, mf;
1015  double als;
1016 
1017  switch (order) {
1018  case LO:
1019  case FULLNLO:
1020  case NLO:
1021  if (nfmu == nfz)
1022  als = AlsWithInit(mu, AlsM, MAls, order);
1023  else if (nfmu > nfz) {
1024  if (order == NLO)
1025  throw std::runtime_error("NLO is not implemented in QCD::Als(mu,order).");
1026  if (nfmu == nfz + 1.) {
1027  mu_thre1 = AboveTh(MAls); // mut
1028  Als_tmp = AlsWithInit(mu_thre1 - MEPS, AlsM, MAls, order);
1029  if (order == FULLNLO) {
1030  mf = mtpole;
1031  Als_tmp = (1. + Als_tmp / M_PI * log(mu_thre1 / mf) / 3.) * Als_tmp;
1032  }
1033  als = AlsWithInit(mu, Als_tmp, mu_thre1 + MEPS, order);
1034  } else
1035  throw std::runtime_error("Error in QCD::Als(mu,order)");
1036  } else {
1037  if (order == NLO)
1038  throw std::runtime_error("NLO is not implemented in QCD::Als(mu,order).");
1039  if (nfmu == nfz - 1.) {
1040  mu_thre1 = BelowTh(MAls); // mub
1041  Als_tmp = AlsWithInit(mu_thre1 + MEPS, AlsM, MAls, order);
1042  if (order == FULLNLO) {
1043  mf = getQuarks(BOTTOM).getMass();
1044  Als_tmp = (1. - Als_tmp / M_PI * log(mu_thre1 / mf) / 3.) * Als_tmp;
1045  }
1046  als = AlsWithInit(mu, Als_tmp, mu_thre1 - MEPS, order);
1047  } else if (nfmu == nfz - 2.) {
1048  mu_thre1 = BelowTh(MAls); // mub
1049  mu_thre2 = AboveTh(mu); // muc
1050  Als_tmp = Als(mu_thre1 + MEPS, order);
1051  if (order == FULLNLO) {
1052  mf = getQuarks(BOTTOM).getMass();
1053  Als_tmp = (1. - Als_tmp / M_PI * log(mu_thre1 / mf) / 3.) * Als_tmp;
1054  }
1055  Als_tmp = AlsWithInit(mu_thre2, Als_tmp, mu_thre1 - MEPS, order);
1056  if (order == FULLNLO) {
1057  mf = getQuarks(CHARM).getMass();
1058  Als_tmp = (1. - Als_tmp / M_PI * log(mu_thre2 / mf) / 3.) * Als_tmp;
1059  }
1060  als = AlsWithInit(mu, Als_tmp, mu_thre2 - MEPS, order);
1061  } else
1062  throw std::runtime_error("Error in QCD::Als(mu,order)");
1063  }
1064  break;
1065  case FULLNNLO:
1066  case NNLO:
1067  /* alpha_s(mu) computed with Lambda_QCD for Nf=nfmu */
1068  als = AlsWithLambda(mu, order);
1069  break;
1070  default:
1071  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Als(mu,order).");
1072  }
1073 
1074  CacheShift(als_cache, 8);
1075  als_cache[0][0] = mu;
1076  als_cache[1][0] = (double) order;
1077  als_cache[2][0] = AlsM;
1078  als_cache[3][0] = MAls;
1079  als_cache[4][0] = mut;
1080  als_cache[5][0] = mub;
1081  als_cache[6][0] = muc;
1082  als_cache[7][0] = als;
1083 
1084  return als;
1085 }
1086 
1087 double QCD::ZeroNf6NLO(double *logLambda6, double *logLambda5_in) const
1088 {
1089  return ( AlsWithLambda(mut + 1.e-10, *logLambda6, FULLNLO)
1090  - AlsWithLambda(mut - 1.e-10, *logLambda5_in, FULLNLO));
1091 }
1092 
1093 double QCD::ZeroNf5(double *logLambda5, double *order) const
1094 {
1095  return ( AlsWithLambda(MAls, *logLambda5, (orders) * order) - AlsM);
1096 }
1097 
1098 double QCD::ZeroNf4NLO(double *logLambda4, double *logLambda5_in) const
1099 {
1100  return ( AlsWithLambda(mub - 1.e-10, *logLambda4, FULLNLO)
1101  - AlsWithLambda(mub + 1.e-10, *logLambda5_in, FULLNLO));
1102 }
1103 
1104 double QCD::ZeroNf3NLO(double *logLambda3, double *logLambda4_in) const
1105 {
1106  return ( AlsWithLambda(muc - 1.e-10, *logLambda3, FULLNLO)
1107  - AlsWithLambda(muc + 1.e-10, *logLambda4_in, FULLNLO));
1108 }
1109 
1110 double QCD::logLambda5(orders order) const
1111 {
1112  if (order == NLO) order = FULLNLO;
1113  if (order == NNLO) order = FULLNNLO;
1114 
1115  for (int i = 0; i < CacheSize; ++i)
1116  if ((AlsM == logLambda5_cache[0][i])
1117  && (MAls == logLambda5_cache[1][i])
1118  && ((double) order == logLambda5_cache[2][i]))
1119  return logLambda5_cache[3][i];
1120 
1121  CacheShift(logLambda5_cache, 4);
1122  logLambda5_cache[0][0] = AlsM;
1123  logLambda5_cache[1][0] = MAls;
1124  logLambda5_cache[2][0] = (double) order;
1125 
1126  if (order == LO)
1127  logLambda5_cache[3][0] = log(MAls) - 2. * M_PI / Beta0(5.) / AlsM;
1128  else {
1129  double xmin = -4., xmax = -0.2;
1130  TF1 f = TF1("f", this, &QCD::ZeroNf5, xmin, xmax, 1, "QCD", "zeroNf5");
1131 
1132  ROOT::Math::WrappedTF1 wf1(f);
1133  double ledouble = (double) order;
1134  wf1.SetParameters(&ledouble);
1135 
1136  ROOT::Math::BrentRootFinder brf;
1137  brf.SetFunction(wf1, xmin, xmax);
1138 
1139  if (brf.Solve()) logLambda5_cache[3][0] = brf.Root();
1140  else
1141  throw std::runtime_error("Error in QCD::logLambda5()");
1142  }
1143  return ( logLambda5_cache[3][0]);
1144 }
1145 
1146 double QCD::logLambdaNLO(const double nfNEW, const double nfORG,
1147  const double logLambdaORG) const
1148 {
1149  for (int i = 0; i < CacheSize; ++i)
1150  if ((AlsM == logLambdaNLO_cache[0][i])
1151  && (MAls == logLambdaNLO_cache[1][i])
1152  && (mut == logLambdaNLO_cache[2][i])
1153  && (mub == logLambdaNLO_cache[3][i])
1154  && (muc == logLambdaNLO_cache[4][i])
1155  && (nfNEW == logLambdaNLO_cache[5][i])
1156  && (nfORG == logLambdaNLO_cache[6][i])
1157  && (logLambdaORG == logLambdaNLO_cache[7][i]))
1158  return logLambdaNLO_cache[8][i];
1159 
1160  CacheShift(logLambdaNLO_cache, 9);
1161  logLambdaNLO_cache[0][0] = AlsM;
1162  logLambdaNLO_cache[1][0] = MAls;
1163  logLambdaNLO_cache[2][0] = mut;
1164  logLambdaNLO_cache[3][0] = mub;
1165  logLambdaNLO_cache[4][0] = muc;
1166  logLambdaNLO_cache[5][0] = nfNEW;
1167  logLambdaNLO_cache[6][0] = nfORG;
1168  logLambdaNLO_cache[7][0] = logLambdaORG;
1169 
1170  double xmin = -4., xmax = -0.2;
1171 
1172  TF1 f;
1173  if (nfNEW == 6. && nfORG == 5.) {
1174  f = TF1("f", this, &QCD::ZeroNf6NLO, xmin, xmax, 1, "QCD", "zeroNf6NLO");
1175  } else if (nfNEW == 4. && nfORG == 5.) {
1176  f = TF1("f", this, &QCD::ZeroNf4NLO, xmin, xmax, 1, "QCD", "zeroNf4NLO");
1177  } else if (nfNEW == 3. && nfORG == 4.) {
1178  f = TF1("f", this, &QCD::ZeroNf3NLO, xmin, xmax, 1, "QCD", "zeroNf3NLO");
1179  } else
1180  throw std::runtime_error("Error in QCD::logLambdaNLO()");
1181 
1182  ROOT::Math::WrappedTF1 wf1(f);
1183  wf1.SetParameters(&logLambdaORG);
1184 
1185  ROOT::Math::BrentRootFinder brf;
1186  brf.SetFunction(wf1, xmin, xmax);
1187 
1188  if (brf.Solve()) logLambdaNLO_cache[8][0] = brf.Root();
1189  else
1190  throw std::runtime_error("Error in QCD::logLambdaNLO()");
1191 
1192  return ( logLambdaNLO_cache[8][0]);
1193 }
1194 
1195 double QCD::logLambda(const double muMatching, const double mf,
1196  const double nfNEW, const double nfORG,
1197  const double logLambdaORG, orders order) const
1198 {
1199  if (fabs(nfNEW - nfORG) != 1.)
1200  throw std::runtime_error("Error in QCD::logLambda()");
1201  if (order == NLO) order = FULLNLO;
1202  if (order == NNLO) order = FULLNNLO;
1203 
1204  /* We do not use the codes below for FULLNLO, since threshold corrections
1205  * can be regarded as an NNLO effect as long as setting the matching scale
1206  * to be close to the mass scale of the decoupling quark. In order to use
1207  * the relation als^{nf+1} = als^{nf} exactly, we use logLambdaNLO method.
1208  */
1209  if (order == FULLNLO)
1210  return logLambdaNLO(nfNEW, nfORG, logLambdaORG);
1211 
1212  double logMuMatching = log(muMatching);
1213  double L = 2. * (logMuMatching - logLambdaORG);
1214  double rNEW = 0.0, rORG = 0.0, log_mu2_mf2 = 0.0, log_L = 0.0;
1215  double C1 = 0.0, C2 = 0.0; // threshold corrections
1216  double logLambdaNEW;
1217 
1218  // LO contribution
1219  logLambdaNEW = 1. / 2. / Beta0(nfNEW)
1220  *(Beta0(nfNEW) - Beta0(nfORG)) * L + logLambdaORG;
1221 
1222  // NLO contribution
1223  if (order == FULLNLO || order == FULLNNLO) {
1224  rNEW = Beta1(nfNEW) / Beta0(nfNEW);
1225  rORG = Beta1(nfORG) / Beta0(nfORG);
1226  log_mu2_mf2 = 2. * (logMuMatching - log(mf));
1227  log_L = log(L);
1228  if (nfNEW < nfORG)
1229  C1 = 2. / 3. * log_mu2_mf2;
1230  else
1231  C1 = -2. / 3. * log_mu2_mf2;
1232  logLambdaNEW += 1. / 2. / Beta0(nfNEW)
1233  *((rNEW - rORG) * log_L
1234  - rNEW * log(Beta0(nfNEW) / Beta0(nfORG)) - C1);
1235  }
1236 
1237  // NNLO contribution
1238  if (order == FULLNNLO) {
1239  if (nfNEW == 5. && nfORG == 6.)
1240  C2 = -16. * (log_mu2_mf2 * log_mu2_mf2 / 36. - 19. / 24. * log_mu2_mf2 - 7. / 24.);
1241  else if (nfNEW == 6. && nfORG == 5.)
1242  C2 = -16. * (log_mu2_mf2 * log_mu2_mf2 / 36. + 19. / 24. * log_mu2_mf2 + 7. / 24.);
1243  else {
1244  if (nfNEW < nfORG)
1245  C2 = -16. * (log_mu2_mf2 * log_mu2_mf2 / 36. - 19. / 24. * log_mu2_mf2 + 11. / 72.);
1246  else
1247  C2 = -16. * (log_mu2_mf2 * log_mu2_mf2 / 36. + 19. / 24. * log_mu2_mf2 - 11. / 72.);
1248  }
1249  logLambdaNEW += 1. / 2. / Beta0(nfNEW) / Beta0(nfORG) / L
1250  * (rORG * (rNEW - rORG) * log_L + rNEW * rNEW - rORG * rORG
1251  - Beta2(nfNEW) / Beta0(nfNEW) + Beta2(nfORG) / Beta0(nfORG)
1252  + rNEW * C1 - C1 * C1 - C2);
1253  }
1254 
1255  return logLambdaNEW;
1256 }
1257 
1258 double QCD::logLambda(const double nf, orders order) const
1259 {
1260  if (order == NLO) order = FULLNLO;
1261  if (order == NNLO) order = FULLNNLO;
1262 
1263  double muMatching, mf, logLambdaORG, logLambdaNEW;
1264  if (nf == 5.)
1265  return logLambda5(order);
1266  else if (nf == 6.) {
1267  muMatching = Thresholds(1); // mut
1268  /* matching condition from Nf=5 to Nf=6 is given in terms of the top pole mass. */
1269  mf = mtpole; // top pole mass
1270  return logLambda(muMatching, mf, 6., 5., logLambda5(order), order);
1271  } else if (nf == 4. || nf == 3.) {
1272  muMatching = Thresholds(2); // mub
1273  mf = getQuarks(BOTTOM).getMass(); // m_b(m_b)
1274  logLambdaORG = logLambda5(order);
1275  logLambdaNEW = logLambda(muMatching, mf, 4., 5., logLambdaORG, order);
1276  if (nf == 3.) {
1277  muMatching = Thresholds(3); // muc
1278  mf = getQuarks(CHARM).getMass(); // m_c(m_c)
1279  logLambdaORG = logLambdaNEW;
1280  logLambdaNEW = logLambda(muMatching, mf, 3., 4., logLambdaORG, order);
1281  }
1282  return logLambdaNEW;
1283  } else
1284  throw std::runtime_error("Error in QCD::logLambda()");
1285 }
1286 
1288 
1289 double QCD::Gamma0(const double nf) const
1290 {
1291  return ( 6. * CF);
1292 }
1293 
1294 double QCD::Gamma1(const double nf) const
1295 {
1296  return ( CF * (3. * CF + 97. / 3. * Nc - 10. / 3. * nf));
1297 }
1298 
1299 double QCD::Gamma2(const double nf) const
1300 {
1301  return ( 129. * CF * CF * CF - 129. / 2. * CF * CF * Nc + 11413. / 54. * CF * Nc * Nc
1302  + CF * CF * nf * (-46. + 48. * zeta3) + CF * Nc * nf * (-556. / 27. - 48. * zeta3)
1303  - 70. / 27. * CF * nf * nf);
1304 }
1305 
1306 double QCD::threCorrForMass(const double nf_f, const double nf_i) const
1307 {
1308  if (fabs(nf_f - nf_i) != 1.)
1309  throw std::runtime_error("Error in QCD::threCorrForMass()");
1310 
1311  double mu_threshold, mf, log_mu2_mf2;
1312  if (nf_f > nf_i) {
1313  if (nf_f == 6.) {
1314  mu_threshold = mut;
1315  mf = quarks[TOP].getMass(); // m_t(m_t)
1316  } else if (nf_f == 5.) {
1317  mu_threshold = mub;
1318  mf = quarks[BOTTOM].getMass(); // m_b(m_b)
1319  } else if (nf_f == 4.) {
1320  mu_threshold = muc;
1321  mf = quarks[CHARM].getMass(); // m_c(m_c)
1322  } else
1323  throw std::runtime_error("Error in QCD::threCorrForMass()");
1324  log_mu2_mf2 = 2. * log(mu_threshold / mf);
1325  return (1. + pow(Als(mu_threshold - MEPS, FULLNNLO) / M_PI, 2.)
1326  *(-log_mu2_mf2 * log_mu2_mf2 / 12. + 5. / 36. * log_mu2_mf2 - 89. / 432.));
1327  } else {
1328  if (nf_i == 6.) {
1329  mu_threshold = mut;
1330  mf = quarks[TOP].getMass(); // m_t(m_t)
1331  } else if (nf_i == 5.) {
1332  mu_threshold = mub;
1333  mf = quarks[BOTTOM].getMass(); // m_b(m_b)
1334  } else if (nf_i == 4.) {
1335  mu_threshold = muc;
1336  mf = quarks[CHARM].getMass(); // m_c(m_c)
1337  } else
1338  throw std::runtime_error("Error in QCD::threCorrForMass()");
1339  log_mu2_mf2 = 2. * log(mu_threshold / mf);
1340  return (1. + pow(Als(mu_threshold + MEPS, FULLNNLO) / M_PI, 2.)
1341  *(log_mu2_mf2 * log_mu2_mf2 / 12. - 5. / 36. * log_mu2_mf2 + 89. / 432.));
1342  }
1343 }
1344 
1345 double QCD::Mrun(const double mu, const double m, const orders order) const
1346 {
1347  return Mrun(mu, m, m, order);
1348 }
1349 
1350 double QCD::Mrun(const double mu_f, const double mu_i, const double m,
1351  const orders order) const
1352 {
1353  // Note: When the scale evolves across a flavour threshold, the definitions
1354  // of the outputs for "NLO" and "NNLO" become complicated.
1355 
1356  int i;
1357  for (i = 0; i < CacheSize; ++i) {
1358  if ((mu_f == mrun_cache[0][i]) && (mu_i == mrun_cache[1][i]) &&
1359  (m == mrun_cache[2][i]) && ((double) order == mrun_cache[3][i]) &&
1360  (AlsM == mrun_cache[4][i]) && (MAls == mrun_cache[5][i]) &&
1361  (mut == mrun_cache[6][i]) && (mub == mrun_cache[7][i]) &&
1362  (muc == mrun_cache[8][i]))
1363  return mrun_cache[9][i];
1364  }
1365 
1366  double nfi = Nf(mu_i), nff = Nf(mu_f);
1367  double mu_threshold, mu_threshold2, mu_threshold3, mrun;
1368  if (nff == nfi)
1369  mrun = MrunTMP(mu_f, mu_i, m, order);
1370  else if (nff > nfi) {
1371  if (order == NLO || order == NNLO)
1372  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Mrun(mu_f,mu_i,m,order)");
1373  mu_threshold = AboveTh(mu_i);
1374  mrun = MrunTMP(mu_threshold - MEPS, mu_i, m, order);
1375  if (order == FULLNNLO)
1376  mrun *= threCorrForMass(nfi + 1., nfi); // threshold corrections
1377  if (nff == nfi + 1.) {
1378  mrun = MrunTMP(mu_f, mu_threshold + MEPS, mrun, order);
1379  } else if (nff == nfi + 2.) {
1380  mu_threshold2 = BelowTh(mu_f);
1381  mrun = MrunTMP(mu_threshold2 - MEPS, mu_threshold + MEPS, mrun, order);
1382  if (order == FULLNNLO)
1383  mrun *= threCorrForMass(nff, nfi + 1.); // threshold corrections
1384  mrun = MrunTMP(mu_f, mu_threshold2 + MEPS, mrun, order);
1385  } else if (nff == nfi + 3.) {
1386  mu_threshold2 = AboveTh(mu_threshold);
1387  mrun = MrunTMP(mu_threshold2 - MEPS, mu_threshold + MEPS, mrun, order);
1388  if (order == FULLNNLO)
1389  mrun *= threCorrForMass(nfi + 2., nfi + 1.); // threshold corrections
1390  mu_threshold3 = BelowTh(mu_f);
1391  mrun = MrunTMP(mu_threshold3 - MEPS, mu_threshold2 + MEPS, mrun, order);
1392  if (order == FULLNNLO)
1393  mrun *= threCorrForMass(nff, nfi + 2.); // threshold corrections
1394  mrun = MrunTMP(mu_f, mu_threshold3 + MEPS, mrun, order);
1395  } else
1396  throw std::runtime_error("Error in QCD::Mrun(mu_f,mu_i,m,order)");
1397  } else {
1398  if (order == NLO || order == NNLO)
1399  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Mrun(mu_f,mu_i,m,order)");
1400  mu_threshold = BelowTh(mu_i);
1401  mrun = MrunTMP(mu_threshold + MEPS, mu_i, m, order);
1402  if (order == FULLNNLO)
1403  mrun *= threCorrForMass(nfi - 1., nfi); // threshold corrections
1404  if (nff == nfi - 1.)
1405  mrun = MrunTMP(mu_f, mu_threshold - MEPS, mrun, order);
1406  else if (nff == nfi - 2.) {
1407  mu_threshold2 = AboveTh(mu_f);
1408  mrun = MrunTMP(mu_threshold2 + MEPS, mu_threshold - MEPS, mrun, order);
1409  if (order == FULLNNLO)
1410  mrun *= threCorrForMass(nff, nfi - 1.); // threshold corrections
1411  mrun = MrunTMP(mu_f, mu_threshold2 - MEPS, mrun, order);
1412  } else
1413  throw std::runtime_error("Error in QCD::Mrun(mu_f,mu_i,m,order)");
1414  }
1415 
1416  if (mrun < 0.0) {
1417  std::stringstream out;
1418  out << "QCD::Mrun(): A quark mass becomes tachyonic in QCD::Mrun("
1419  << mu_f << ", " << mu_i << ", " << m << ", " << orderToString(order) << ")"
1420  << std::endl
1421  << " Als(" << mu_i << ", " << orderToString(order) << ")/(4pi)="
1422  << Als(mu_i, order) / (4. * M_PI) << std::endl
1423  << " Als(" << mu_f << ", " << orderToString(order) << ")/(4pi)="
1424  << Als(mu_f, order) / (4. * M_PI);
1425  throw std::runtime_error(out.str());
1426  }
1427 
1428  CacheShift(mrun_cache, 10);
1429  mrun_cache[0][0] = mu_f;
1430  mrun_cache[1][0] = mu_i;
1431  mrun_cache[2][0] = m;
1432  mrun_cache[3][0] = (double) order;
1433  mrun_cache[4][0] = AlsM;
1434  mrun_cache[5][0] = MAls;
1435  mrun_cache[6][0] = mut;
1436  mrun_cache[7][0] = mub;
1437  mrun_cache[8][0] = muc;
1438  mrun_cache[9][0] = mrun;
1439 
1440  return mrun;
1441 }
1442 
1443 double QCD::MrunTMP(const double mu_f, const double mu_i, const double m,
1444  const orders order) const
1445 {
1446  double nf = Nf(mu_f);
1447  if (nf != Nf(mu_i))
1448  throw std::runtime_error("Error in QCD::MrunTMP().");
1449 
1450  // alpha_s/(4pi)
1451  orders orderForAls;
1452  if (order == LO) orderForAls = LO;
1453  if (order == NLO || order == FULLNLO) orderForAls = FULLNLO;
1454  if (order == NNLO || order == FULLNNLO) orderForAls = FULLNNLO;
1455  double ai = Als(mu_i, orderForAls) / (4. * M_PI);
1456  double af = Als(mu_f, orderForAls) / (4. * M_PI);
1457 
1458  // LO contribution
1459  double b0 = Beta0(nf), g0 = Gamma0(nf);
1460  double mLO = m * pow(af / ai, g0 / (2. * b0));
1461  if (order == LO) return mLO;
1462 
1463  // NLO contribution
1464  double b1 = Beta1(nf), g1 = Gamma1(nf);
1465  double A1 = g1 / (2. * b0) - b1 * g0 / (2. * b0 * b0);
1466  double mNLO = mLO * A1 * (af - ai);
1467  if (order == NLO) return mNLO;
1468  if (order == FULLNLO) return (mLO + mNLO);
1469 
1470  // NNLO contribution
1471  double b2 = Beta2(nf), g2 = Gamma2(nf);
1472  double A2 = b1 * b1 * g0 / (2. * b0 * b0 * b0) - b2 * g0 / (2. * b0 * b0) - b1 * g1 / (2. * b0 * b0) + g2 / (2. * b0);
1473  double mNNLO = mLO * (A1 * A1 / 2. * (af - ai)*(af - ai) + A2 / 2. * (af * af - ai * ai));
1474  if (order == NNLO) return mNNLO;
1475  if (order == FULLNNLO) return (mLO + mNLO + mNNLO);
1476 
1477  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::MrunTMP()");
1478 }
1479 
1480 double QCD::Mrun4(const double mu_f, const double mu_i, const double m) const
1481 {
1482  double nf = 4.;
1483 
1484  // alpha_s/(4pi)
1485  double ai = Als4(mu_i) / (4. * M_PI);
1486  double af = Als4(mu_f) / (4. * M_PI);
1487 
1488  // LO contribution
1489  double b0 = Beta0(nf), g0 = Gamma0(nf);
1490  double mLO = m * pow(af / ai, g0 / (2. * b0));
1491 
1492  // NLO contribution
1493  double b1 = Beta1(nf), g1 = Gamma1(nf);
1494  double A1 = g1 / (2. * b0) - b1 * g0 / (2. * b0 * b0);
1495  double mNLO = mLO * A1 * (af - ai);
1496  return (mLO + mNLO);
1497 
1498 }
1499 
1501 
1502 double QCD::Mbar2Mp(const double mbar, const orders order) const
1503 {
1504  // LO contribution
1505  double MpLO = mbar;
1506  if (order == LO) return MpLO;
1507 
1508  // alpha_s(mbar)/pi
1509  orders orderForAls;
1510  if (order == NLO || order == FULLNLO) orderForAls = FULLNLO;
1511  if (order == NNLO || order == FULLNNLO) orderForAls = FULLNNLO;
1512  double a = Als(mbar + MEPS, orderForAls) / M_PI;
1513 
1514  // NLO contribution
1515  double MpNLO = mbar * 4. / 3. * a;
1516  if (order == NLO) return MpNLO;
1517  if (order == FULLNLO) return (MpLO + MpNLO);
1518 
1519  // NNLO contribution
1520  double nl, x;
1521  if (mbar < 3.)
1522  throw std::runtime_error("QCD::Mbar2Mp() can convert only top and bottom masses");
1523  else if (mbar < 50.) {
1524  // for the b quark
1525  nl = 4.;
1526  /* simply adding m_s(2 GeV) and m_c(m_c) */
1527  x = (quarks[STRANGE].getMass() + quarks[CHARM].getMass()) / mbar;
1528  } else {
1529  // for the top quark
1530  nl = 5.;
1531  /* simply adding m_s(2 GeV), m_c(m_c) and m_b(m_b) */
1532  x = (quarks[STRANGE].getMass() + quarks[CHARM].getMass()
1533  + quarks[BOTTOM].getMass()) / mbar;
1534  }
1535  double Delta = M_PI * M_PI / 8. * x - 0.597 * x * x + 0.230 * x * x*x;
1536  double MpNNLO = mbar * (307. / 32. + 2. * zeta2 + 2. / 3. * zeta2 * log(2.0) - zeta3 / 6.
1537  - nl / 3. * (zeta2 + 71. / 48.) + 4. / 3. * Delta) * a*a;
1538  if (order == NNLO) return MpNNLO;
1539  if (order == FULLNNLO) return (MpLO + MpNLO + MpNNLO);
1540 
1541  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Mbar2Mp().");
1542 }
1543 
1544 double QCD::Mp2MbarTMP(double *mu, double *params) const
1545 {
1546  double mp = params[0];
1547  orders order = (orders) params[1];
1548  return (mp - Mbar2Mp(*mu, order));
1549 }
1550 
1551 double QCD::Mp2Mbar(const double mp, const orders order) const
1552 {
1553  if (order == NLO || order == NNLO)
1554  throw std::runtime_error(orderToString(order) + " is not implemented in QCD::Mp2Mbar().");
1555 
1556  int i;
1557  double ms = quarks[STRANGE].getMass(), mc = quarks[CHARM].getMass();
1558  double alsmp = Als(mp, order);
1559  for (i = 0; i < CacheSize; ++i)
1560  if (alsmp == mp2mbar_cache[0][i] && ms == mp2mbar_cache[1][i] &&
1561  mc == mp2mbar_cache[2][i] && (double) order == mp2mbar_cache[3][i])
1562  return mp2mbar_cache[4][i];
1563 
1564  CacheShift(mp2mbar_cache, 5);
1565  mp2mbar_cache[0][0] = alsmp;
1566  mp2mbar_cache[1][0] = ms;
1567  mp2mbar_cache[2][0] = mc;
1568  mp2mbar_cache[3][0] = (double) order;
1569 
1570  TF1 f("f", this, &QCD::Mp2MbarTMP, mp / 2., 2. * mp, 2, "QCD", "mp2mbara");
1571 
1572  ROOT::Math::WrappedTF1 wf1(f);
1573  double params[2];
1574  params[0] = mp;
1575  params[1] = (double) order;
1576  wf1.SetParameters(params);
1577 
1578  ROOT::Math::BrentRootFinder brf;
1579 
1580  brf.SetFunction(wf1, .7 * mp, 1.3 * mp);
1581  if (brf.Solve())
1582  mp2mbar_cache[4][0] = brf.Root();
1583  else
1584  throw std::runtime_error("error in QCD::mp2mbar");
1585 
1586  return (mp2mbar_cache[4][0]);
1587 }
1588 
1589 double QCD::MS2DRqmass(const double MSbar) const
1590 {
1591  return (MSbar / (1. + Als(MSbar, FULLNLO) / 4. / M_PI * CF));
1592 }
1593 
1594 double QCD::MS2DRqmass(const double MSscale, const double MSbar) const
1595 {
1596  return (MSbar / (1. + Als(MSscale, FULLNLO) / 4. / M_PI * CF));
1597 }
QCD::logLambda
double logLambda(const double nf, orders order) const
Computes with nf flavours in GeV.
Definition: QCD.cpp:1195
QCD::K_P
Definition: QCD.h:717
Particle
A class for particles.
Definition: Particle.h:26
QCD::NQCDvars
static const int NQCDvars
The number of model parameters in QCD.
Definition: QCD.h:739
QCD::Beta1
double Beta1(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:892
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:970
orders
orders
An enum type for orders in QCD.
Definition: OrderScheme.h:31
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
QCD::PostUpdate
virtual bool PostUpdate()
The post-update method for QCD.
Definition: QCD.cpp:323
QCD::Update
virtual bool Update(const std::map< std::string, double > &DPars)
The update method for QCD.
Definition: QCD.cpp:307
QCD::setFlagStr
virtual bool setFlagStr(const std::string name, const std::string value)
A method to set a flag of QCD.
Definition: QCD.cpp:822
QCD::Als
double Als(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:1004
QCD::K_0
Definition: QCD.h:716
QCD::D_0
Definition: QCD.h:718
Model::UpdateError
bool UpdateError
A boolean set to false if update is successful.
Definition: Model.h:192
QCD::UP
Definition: QCD.h:731
QCD::setFlag
virtual bool setFlag(const std::string name, const bool value)
A method to set a flag of QCD.
Definition: QCD.cpp:816
NLO
Definition: OrderScheme.h:34
QCD::CheckFlags
virtual bool CheckFlags() const
A method to check the sanity of the set of model flags.
Definition: QCD.cpp:828
schemes
schemes
An enum type for regularization schemes.
Definition: OrderScheme.h:19
FULLNNLO
Definition: OrderScheme.h:37
QCD::TOP
Definition: QCD.h:735
NNLO
Definition: OrderScheme.h:35
QCD::BelowTh
double BelowTh(const double mu) const
The active flavour threshold below the scale as defined in QCD::Thresholds().
Definition: QCD.cpp:858
QCD::AboveTh
double AboveTh(const double mu) const
The active flavour threshold above the scale as defined in QCD::Thresholds().
Definition: QCD.cpp:849
QCD::PreUpdate
virtual bool PreUpdate()
The pre-update method for QCD.
Definition: QCD.cpp:296
QCD::QCD
QCD()
Constructor.
Definition: QCD.cpp:56
QCD.h
LO
Definition: OrderScheme.h:33
QCD::CheckParameters
virtual bool CheckParameters(const std::map< std::string, double > &DPars)
A method to check if all the mandatory parameters for StandardModel have been provided in model initi...
Definition: QCD.cpp:804
QCD::B_P
Definition: QCD.h:720
QCD::DOWN
Definition: QCD.h:732
QCD::alphaMz
virtual double alphaMz() const =0
QCD::B_S
Definition: QCD.h:721
FULLNLO
Definition: OrderScheme.h:36
QCD::getQuarks
Particle getQuarks(const quark q) const
A get method to access a quark as an object of the type Particle.
Definition: QCD.h:869
MEPS
#define MEPS
Definition: QCD.h:14
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:903
QCD::PHI
Definition: QCD.h:722
Particle::getMass
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
QCD::B_D
Definition: QCD.h:719
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
QCD::Init
virtual bool Init(const std::map< std::string, double > &DPars)
Initializes the QCD parameters found in the argument.
Definition: QCD.cpp:287
QCD::Nf
double Nf(const double mu) const
The number of active flavour at scale .
Definition: QCD.cpp:867
QCD::orderToString
std::string orderToString(const orders order) const
Converts an object of the enum type "orders" to the corresponding string.
Definition: QCD.cpp:267
gslpp::complex
A class for defining operations on and functions of complex numbers.
Definition: gslpp_complex.h:35
QCD::Beta0
double Beta0(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:887
QCD::QCDvars
static const std::string QCDvars[NQCDvars]
An array containing the labels under which all QCD parameters are stored in a vector of ModelParamete...
Definition: QCD.h:745
QCD::Als4
double Als4(const double mu) const
The value of at any scale with the number of flavours .
Definition: QCD.cpp:923
Model::ModelParamMap
std::map< std::string, boost::reference_wrapper< const double > > ModelParamMap
Definition: Model.h:200
QCD::Beta2
double Beta2(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:897
QCD::Thresholds
double Thresholds(const int i) const
For accessing the active flavour threshold scales.
Definition: QCD.cpp:835
Model::setParameter
virtual void setParameter(const std::string name, const double &value)=0
A method to set the value of a parameter of the model.