LHAPDF is hosted by Hepforge, IPPP Durham
LHAPDF  6.2.0
/tests/testalphas.cc
// Example program to test alpha_s calculation
#include "LHAPDF/LHAPDF.h"
#include "boost/assign.hpp"
#include <iostream>
#include <fstream>
#include <limits>
using namespace LHAPDF;
using namespace std;
using namespace boost::assign;
int main() {
// Set up three standalone AlphaS solvers:
// Can set order of QCD (up to 4)
// 0 returns a constant value -- needs to be set by as_ana.setAlphaSMZ(double value);
as_ana.setOrderQCD(4);
// Set quark masses for evolution (both transition points and gradients by default)
as_ana.setQuarkMass(1, 0.0017);
as_ana.setQuarkMass(2, 0.0041);
as_ana.setQuarkMass(3, 0.1);
as_ana.setQuarkMass(4, 1.29);
as_ana.setQuarkMass(5, 4.1);
as_ana.setQuarkMass(6, 172.5);
as_ana.setLambda(3, 0.339);
as_ana.setLambda(4, 0.296);
as_ana.setLambda(5, 0.213);
// Can override quark masses for Nf transitions thresholds by setting them explicitly
// You can't mix the two: if you set one flavor threshold explicitly you need to set all of them.
//as_ode.setQuarkThreshold(6, 650);
//as_ode.setQuarkThreshold(5, 10);
//as_ode.setQuarkThreshold(4, 2);
//as_ode.setQuarkThreshold(3, 0.3);
//as_ode.setQuarkThreshold(2, 0.1);
//as_ode.setQuarkThreshold(1, 0.08);
// Uncomment to use fixed flavor scheme for analytic solver
// as_ana.setFlavorScheme(AlphaS::FIXED, 5);
AlphaS_ODE as_ode;
// As above: order = 0 returns
// constant value set by
// as_ode.setAlphaSMZ(double value);
as_ode.setOrderQCD(5);
as_ode.setMZ(91);
as_ode.setAlphaSMZ(0.118);
// as_ode.setMassReference(4.1);
// as_ode.setAlphaSReference(0.21);
as_ode.setQuarkMass(1, 0.0017);
as_ode.setQuarkMass(2, 0.0041);
as_ode.setQuarkMass(3, 0.1);
as_ode.setQuarkMass(4, 1.29);
as_ode.setQuarkMass(5, 4.1);
as_ode.setQuarkMass(6, 172.5);
// Can override quark masses for Nf transitions thresholds by setting them explicitly
// You can't mix the two: if you set one flavor threshold explicitly you need to set all of them.
// as_ode.setQuarkThreshold(6, 650);
// as_ode.setQuarkThreshold(5, 10);
// as_ode.setQuarkThreshold(4, 2);
// as_ode.setQuarkThreshold(3, 0.3);
// as_ode.setQuarkThreshold(2, 0.1);
// as_ode.setQuarkThreshold(1, 0.08);
// Uncomment to use fixed flavor scheme for ODE solver
// as_ode.setFlavorScheme(AlphaS::FIXED, 4);
AlphaS_Ipol as_ipol;
vector<double> qs; qs += 1.300000e+00, 1.300000e+00, 1.560453e+00, 1.873087e+00, 2.248357e+00, 2.698811e+00, 3.239513e+00, 3.888544e+00, 4.667607e+00, 5.602754e+00, 6.725257e+00, 8.072650e+00, 9.689992e+00, 1.163137e+01, 1.396169e+01, 1.675889e+01, 2.011651e+01, 2.414681e+01, 2.898459e+01, 3.479160e+01, 4.176203e+01, 5.012899e+01, 6.017224e+01, 7.222765e+01, 8.669834e+01, 1.040682e+02, 1.249181e+02, 1.499452e+02, 1.799865e+02, 2.160465e+02, 2.593310e+02, 3.112875e+02, 3.736534e+02, 4.485143e+02, 5.383733e+02, 6.462355e+02, 7.757077e+02, 9.311194e+02, 1.117668e+03, 1.341590e+03, 1.610376e+03, 1.933012e+03, 2.320287e+03, 2.785153e+03, 3.343154e+03, 4.012949e+03, 4.816936e+03, 4.816936e+03;
vector<double> alphas; alphas += 4.189466e-01, 4.189466e-01, 3.803532e-01, 3.482705e-01, 3.211791e-01, 2.979983e-01, 2.779383e-01, 2.604087e-01, 2.451922e-01, 2.324903e-01, 2.210397e-01, 2.106640e-01, 2.012187e-01, 1.925841e-01, 1.846600e-01, 1.773623e-01, 1.706194e-01, 1.643704e-01, 1.585630e-01, 1.531520e-01, 1.480981e-01, 1.433671e-01, 1.389290e-01, 1.347574e-01, 1.308290e-01, 1.271232e-01, 1.236215e-01, 1.203076e-01, 1.171667e-01, 1.141857e-01, 1.113525e-01, 1.086566e-01, 1.060881e-01, 1.036382e-01, 1.012990e-01, 9.906295e-02, 9.692353e-02, 9.487457e-02, 9.291044e-02, 9.102599e-02, 8.921646e-02, 8.747747e-02, 8.580498e-02, 8.419524e-02, 8.264479e-02, 8.115040e-02, 7.970910e-02, 7.970910e-02;
as_ipol.setQValues(qs);
as_ipol.setAlphaSValues(alphas);
// Can interpolate ODE with given knots in Q
// as_ode.setQValues(qs);
// Test these solvers and the CT10nlo PDF's default behaviours:
PDF* pdf = mkPDF("CT10", 0);
const double inf = numeric_limits<double>::infinity();
ofstream fa("alphas_ana.dat"), fo("alphas_ode.dat"), fi("alphas_ipol.dat"), fc("alphas_ct10nlo.dat");
cout << endl;
for (double log10q = -0.5; log10q < 3; log10q += 0.05) {
const double q = pow(10, log10q);
const double as_ana_q = as_ana.alphasQ(q);
cout << fixed;
cout << "Q = " << setprecision(3) << q << " GeV" << endl;
cout << "Analytical solution: " << setprecision(3) << setw(6) << ( (as_ana_q > 2) ? inf : as_ana_q )
<< " num flavs = " << as_ana.numFlavorsQ(q) << endl;
fa << q << " " << as_ana_q << endl;
const double as_ode_q = as_ode.alphasQ(q);
cout << "ODE solution: " << setprecision(3) << setw(6) << ( (as_ode_q > 2) ? inf : as_ode_q )
<< " num flavs = " << as_ode.numFlavorsQ(q) << endl;
fo << q << " " << as_ode_q << endl;
/*
const double as_ipol_q = as_ipol.alphasQ(q);
cout << "Interpolated solution: " << setprecision(3) << setw(6) << ( (as_ipol_q > 2) ? inf : as_ipol_q ) << endl;
fi << q << " " << as_ipol_q << endl;
const double as_ct10_q = pdf->alphasQ(q);
cout << "CT10 solution: " << setprecision(3) << setw(6) << as_ct10_q << endl;
fc << q << " " << as_ct10_q << endl;
const double as_ct10_q_2 = pdf->alphaS().alphasQ(q);
cout << "CT10 AlphaS solution: " << setprecision(3) << setw(6) << as_ct10_q_2
<< " agrees = " << boolalpha << (as_ct10_q == as_ct10_q_2) << endl;
cout << endl;
*/
}
fa.close(); fo.close(); fi.close(); fc.close();
delete pdf;
return 0;
}