Page MenuHomeHEPForge
Files F8310207

No OneTemporary
Actions

View Options

diff --git a/src/AlphaS.cc b/src/AlphaS.cc
--- a/src/AlphaS.cc
+++ b/src/AlphaS.cc
@@ -1,102 +1,102 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/AlphaS.h"
#include "LHAPDF/Utils.h"
namespace LHAPDF {
// Base class constructor for default param setup
AlphaS::AlphaS() {
_qcdorder = 4;
_mz = 91.1876;
_alphas_mz = 0.118;
_flavorscheme = VARIABLE;
_fixflav = -1;
_customref = false;
}
// Calculate the number of active quark flavours at energy scale Q2
int AlphaS::numFlavorsQ2(double q2) const {
if ( _flavorscheme == FIXED ) return _fixflav;
int nf = 0;
/// Use quark masses if flavour threshold not set explicitly
if ( _flavorthresholds.empty() ) {
for (int it = 1; it <= 6; ++it) {
std::map<int, double>::const_iterator element = _quarkmasses.find(it);
if ( element == _quarkmasses.end() ) continue;
if ( sqr(element->second) < q2 ) nf = it;
}
} else {
for (int it = 1; it <= 6; ++it) {
std::map<int, double>::const_iterator element = _flavorthresholds.find(it);
if ( element == _flavorthresholds.end() ) continue;
if ( sqr(element->second) < q2 ) nf = it;
}
}
if ( _fixflav != -1 && nf > _fixflav ) nf = _fixflav;
return nf;
}
// Calculate a beta function given the number of active flavours
double AlphaS::_beta(int i, int nf) const {
if (i == 0) return (double) 0.875352187 - 0.053051647*nf; //(33 - 2*nf)/(12*M_PI)
if (i == 1) return (double) 0.6459225457 - 0.0802126037*nf; //(153 - 19*nf)/(24*sqr(M_PI))
if (i == 2) return (double) 0.719864327 - 0.140904490*nf + 0.00303291339*nf*nf; //(2857 - (5033 / 9.0)*nf + (325 / 27.0)*sqr(nf))/(128*sqr(M_PI)*M_PI)
if (i == 3) return (double) 1.172686 - 0.2785458*nf + 0.01624467*nf*nf + 0.0000601247*nf*nf*nf;
// ( (149753/6.) + 3564*ZETA_3 - ((1078361/162.) + (6502/27.)*ZETA_3)*nf +
// ((50065/162.) + (6472/81.)*ZETA_3)*sqr(nf) + (1093/729.)*sqr(nf)*nf)/(256*sqr(M_PI)*sqr(M_PI))
throw Exception("Invalid index " + to_str(i) + " for requested beta function");
}
// Calculate beta functions given the number of active flavours
vector<double> AlphaS::_betas(int nf) const {
vector<double> rtn; rtn.reserve(4);
for (int i = 0; i < 4; ++i) rtn.push_back(_beta(i, nf));
return rtn;
}
// Set a quark mass, explicitly giving its ID
void AlphaS::setQuarkMass(int id, double value) {
if (abs(id) > 6 || id == 0)
throw Exception("Invalid ID " + to_str(id) + " for quark given (should be 1-6).");
_quarkmasses[abs(id)] = value;
}
// Set a flavour threshold, explicitly giving its ID
void AlphaS::setQuarkThreshold(int id, double value) {
if (abs(id) > 6 || id == 0)
throw Exception("Invalid ID " + to_str(id) + " for flavour threshold given (should be 1-6).");
_flavorthresholds[abs(id)] = value;
}
// Get a quark mass by ID
double AlphaS::quarkMass(int id) const {
std::map<int, double>::const_iterator quark = _quarkmasses.find(abs(id));
if ( quark == _quarkmasses.end() )
throw Exception("Quark mass " + to_str(id) + " not set!");
return quark->second;
}
// Get a quark mass by ID
double AlphaS::quarkThreshold(int id) const {
std::map<int, double>::const_iterator threshold = _flavorthresholds.find(abs(id));
if ( threshold == _flavorthresholds.end() )
throw Exception("Flavour threshold " + to_str(id) + " not set!");
return threshold->second;
}
void AlphaS::setFlavorScheme(FlavorScheme scheme, int nf) {
if( scheme == FIXED && nf == -1 ) throw Exception("You need to define the number of flavors when using a fixed scheme!");
_flavorscheme = scheme;
_fixflav = nf;
}
}
diff --git a/src/AlphaS_Analytic.cc b/src/AlphaS_Analytic.cc
--- a/src/AlphaS_Analytic.cc
+++ b/src/AlphaS_Analytic.cc
@@ -1,134 +1,134 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/AlphaS.h"
#include "LHAPDF/Utils.h"
namespace LHAPDF {
// Calculate the number of active quark flavours at energy scale Q2.
// Respects min/max nf
/// Currently returns the active number of flavors,
/// not the number actually used for lambdaQCD
/// (in case only lambda3 and lambda5 are defined and
/// we are in the 4 flavour range, we use lambda3 but this returns 4)
/// @todo Is this the "correct" behaviour?
int AlphaS_Analytic::numFlavorsQ2(double q2) const {
if ( _flavorscheme == FIXED ) return _fixflav;
int nf = _nfmin;
/// Use quark masses if flavour threshold not set explicitly
if ( _flavorthresholds.empty() ) {
for ( int it = _nfmin; it <= _nfmax; ++it ) {
std::map<int, double>::const_iterator element = _quarkmasses.find(it);
if ( element == _quarkmasses.end() ) continue;
if ( sqr(element->second) < q2 ) nf = it;
}
} else {
for ( int it = _nfmin; it <= _nfmax; ++it ) {
std::map<int, double>::const_iterator element = _flavorthresholds.find(it);
if ( element == _flavorthresholds.end() ) continue;
if ( sqr(element->second) < q2 ) nf = it;
}
}
if ( _fixflav != -1 && nf > _fixflav ) nf = _fixflav;
return nf;
}
// Set lambda_i && recalculate nfmax and nfmin
void AlphaS_Analytic::setLambda(unsigned int i, double lambda) {
_lambdas[i] = lambda;
_setFlavors();
}
// Recalculate nfmax and nfmin after a new lambda has been set
void AlphaS_Analytic::_setFlavors() {
for (int it = 0; it <= 6; ++it) {
std::map<int, double>::iterator element = _lambdas.find(it);
if ( element == _lambdas.end() ) continue;
_nfmin = it;
break;
}
for (int it = 6; it >= 0; --it) {
std::map<int, double>::iterator element = _lambdas.find(it);
if ( element == _lambdas.end() ) continue;
_nfmax = it;
break;
}
}
// Return the correct lambda for a given number of active flavours
// Uses recursion to find the closest defined-but-lower lambda for the given
// number of active flavours
// If a fixed flavor scheme is used, require the correct lambda to be set
double AlphaS_Analytic::_lambdaQCD(int nf) const {
if ( _flavorscheme == FIXED ) {
std::map<int, double>::const_iterator lambda = _lambdas.find(_fixflav);
if ( lambda == _lambdas.end() ) throw Exception("Set lambda(" + to_str(_fixflav) + ") when using a fixed " + to_str(_fixflav) + " flavor scheme.");
return lambda->second;
} else {
if ( nf < 0 ) throw Exception("Requested lambdaQCD for " + to_str(nf) + " number of flavours.");
std::map<int, double>::const_iterator lambda = _lambdas.find(nf);
if ( lambda == _lambdas.end() ) return _lambdaQCD(nf-1);
return lambda->second;
}
}
// Calculate alpha_s(Q2) by an analytic approximation
double AlphaS_Analytic::alphasQ2(double q2) const {
/// Get the number of active flavours and corresponding LambdaQCD
/// Should support any number of active flavors as long as the
/// corresponding lambas are set
if ( _lambdas.empty() ) throw Exception("You need to set at least one lambda value to calculate alpha_s by analytic means!");
const int nf = this->numFlavorsQ2(q2);
const double lambdaQCD = _lambdaQCD(nf);
if (q2 <= lambdaQCD * lambdaQCD) return std::numeric_limits<double>::max();
// Get beta coeffs for the number of active (above threshold) quark flavours at energy Q
const std::vector<double> beta = _betas(nf);
const double beta02 = sqr(beta[0]);
const double beta12 = sqr(beta[1]);
// Pre-calculate ln(Q2/lambdaQCD) and expansion term y = 1/ln(Q2/lambdaQCD)
const double x = q2 / (lambdaQCD*lambdaQCD);
const double lnx = log(x);
const double lnlnx = log(lnx);
const double lnlnx2 = lnlnx * lnlnx;
const double lnlnx3 = lnlnx * lnlnx * lnlnx;
const double y = 1 / lnx;
// Calculate terms up to qcdorder = 4
// A bit messy because the actual expressions are
// quite messy...
/// @todo Is it okay to use _alphas_mz as the constant value?
if(_qcdorder == 0) return _alphas_mz;
const double A = 1 / beta[0];
const double a_0 = 1;
double tmp = a_0;
if (_qcdorder > 1) {
const double a_1 = beta[1] * lnlnx / beta02;
tmp -= a_1 * y;
}
if (_qcdorder > 2) {
const double B = beta12 / (beta02 * beta02);
const double a_20 = lnlnx2 - lnlnx;
const double a_21 = beta[2] * beta[0] / beta12;
const double a_22 = 1;
tmp += B * y*y * (a_20 + a_21 - a_22);
}
if (_qcdorder > 3) {
const double C = 1. / (beta02 * beta02 * beta02);
const double a_30 = (beta12 * beta[1]) * (lnlnx3 - (5/2.) * lnlnx2 - 2 * lnlnx + 0.5);
const double a_31 = 3 * beta[0] * beta[1] * beta[2] * lnlnx;
const double a_32 = 0.5 * beta02 * beta[3];
tmp -= C * y*y*y * (a_30 + a_31 - a_32);
}
const double alphaS = A * y * tmp;
return alphaS;
}
}
diff --git a/src/AlphaS_Ipol.cc b/src/AlphaS_Ipol.cc
--- a/src/AlphaS_Ipol.cc
+++ b/src/AlphaS_Ipol.cc
@@ -1,130 +1,130 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/AlphaS.h"
#include "LHAPDF/Utils.h"
namespace LHAPDF {
void AlphaS_Ipol::setQ2Values(const std::vector<double>& q2s) {
_q2s = q2s;
}
void AlphaS_Ipol::setAlphaSValues(const std::vector<double>& as) {
_as = as;
}
/// @note This is const so it can be called silently from a const method
void AlphaS_Ipol::_setup_grids() const {
if (!_knotarrays.empty())
throw LogicError("AlphaS interpolation subgrids being initialized a second time!");
if (_q2s.size() != _as.size())
throw MetadataError("AlphaS value and Q interpolation arrays are differently sized");
// Walk along the Q2 vector, making subgrids at each boundary
double prevQ2 = _q2s.front();
vector<double> q2s, as;
size_t combined_lenq2s = 0; //< For consistency checking
for (size_t i = 0; i <= _q2s.size(); ++i) { //< The iteration to len+1 is intentional
// Get current Q2 and alpha_s points, faked if i > vector to force syncing
const double currQ2 = (i != _q2s.size()) ? _q2s[i] : _q2s.back();
const double currAS = (i != _q2s.size()) ? _as[i] : -1;
// If the Q2 value is repeated, sync the current subgrid and start a new one.
// Note special treatment for the first and last points in q2s.
if (abs(currQ2 - prevQ2) < numeric_limits<double>::epsilon()) {
// Sync current subgrid as as AlphaSArray
if (i != 0) {
_knotarrays[q2s.front()] = AlphaSArray(q2s, as);
combined_lenq2s += q2s.size();
}
// Reset temporary vectors
q2s.clear(); q2s.reserve(_q2s.size() - i);
as.clear(); as.reserve(_q2s.size() - i);
}
// Append current value to temporary vectors
q2s.push_back(currQ2);
as.push_back(currAS);
prevQ2 = currQ2;
}
if (combined_lenq2s != _q2s.size())
throw AlphaSError("Sum of alpha_s subgrid sizes does not match input knot array ("
+ to_str(combined_lenq2s) + " vs. " + to_str(_q2s.size()) + ")");
}
double AlphaS_Ipol::_interpolateCubic(double T, double VL, double VDL, double VH, double VDH) const {
// Pre-calculate powers of T
const double t2 = T*T;
const double t3 = t2*T;
// Calculate left point
const double p0 = (2*t3 - 3*t2 + 1)*VL;
const double m0 = (t3 - 2*t2 + T)*VDL;
// Calculate right point
const double p1 = (-2*t3 + 3*t2)*VH;
const double m1 = (t3 - t2)*VDH;
return abs(p0 + m0 + p1 + m1) < 2. ? p0 + m0 + p1 + m1 : std::numeric_limits<double>::max();
}
// Interpolate alpha_s from tabulated points in Q2 via metadata
double AlphaS_Ipol::alphasQ2(double q2) const {
assert(q2 >= 0);
// Using base 10 for logs to get constant gradient extrapolation in
// a log 10 - log 10 plot
if (q2 < _q2s.front()) {
// Remember to take situations where the first knot also is a
// flavor threshold into account
double dlogq2, dlogas;
unsigned int next_point = 1;
while ( _q2s[0] == _q2s[next_point] ) next_point++;
dlogq2 = log10( _q2s[next_point] / _q2s[0] );
dlogas = log10( _as[next_point] / _as[0] );
const double loggrad = dlogas / dlogq2;
return _as[0] * pow( q2/_q2s[0] , loggrad );
}
if (q2 > _q2s.back()) return _as.back();
// If this is the first valid query, set up the ipol grids
if (_knotarrays.empty()) _setup_grids();
// Retrieve the appropriate subgrid
map<double, AlphaSArray>::const_iterator it = --(_knotarrays.upper_bound(q2));
const AlphaSArray& arr = it->second;
// Get the Q/alpha_s index on this array which is *below* this Q point
const size_t i = arr.iq2below(q2);
// Calculate derivatives
double didlogq2, di1dlogq2;
if ( i == 0 ) {
didlogq2 = arr.ddlogq_forward(i);
di1dlogq2 = arr.ddlogq_central(i+1);
} else if ( i == arr.logq2s().size()-2 ) {
didlogq2 = arr.ddlogq_central(i);
di1dlogq2 = arr.ddlogq_backward(i+1);
} else {
didlogq2 = arr.ddlogq_central(i);
di1dlogq2 = arr.ddlogq_central(i+1);
}
// Calculate alpha_s
const double dlogq2 = arr.logq2s()[i+1] - arr.logq2s()[i];
const double tlogq2 = (log(q2) - arr.logq2s()[i]) / dlogq2;
return _interpolateCubic( tlogq2,
arr.alphas()[i], didlogq2*dlogq2,
arr.alphas()[i+1], di1dlogq2*dlogq2 );
}
}
diff --git a/src/AlphaS_ODE.cc b/src/AlphaS_ODE.cc
--- a/src/AlphaS_ODE.cc
+++ b/src/AlphaS_ODE.cc
@@ -1,313 +1,313 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/AlphaS.h"
#include "LHAPDF/Utils.h"
#include "boost/bind.hpp"
namespace LHAPDF {
// Calculate first order derivative, dy/dt, as it appears in the differential equation
double AlphaS_ODE::_derivative(double t, double y, const vector<double>& beta) const {
if ( _qcdorder == 0 ) return 0;
const double b0 = beta[0];
double d = (b0*y*y);
if ( _qcdorder == 1 ) return - d / t;
const double b1 = beta[1];
d += (b1*y*y*y);
if ( _qcdorder == 2 ) return - d / t;
const double b2 = beta[2];
d += (b2*y*y*y*y);
if ( _qcdorder == 3 ) return - d / t;
const double b3 = beta[3];
d += (b3*y*y*y*y*y);
return - d / t;
}
// Calculate decoupling for transition from num. flavour = ni -> nf
double AlphaS_ODE::_decouple(double y, double t, unsigned int ni, unsigned int nf) const {
if ( ni == nf || _qcdorder == 0 ) return 1.;
double as = y / M_PI;
unsigned int heavyQuark = nf > ni ? nf : ni;
std::map<int, double>::const_iterator quark = _quarkmasses.find(heavyQuark);
if ( quark == _quarkmasses.end() ) throw AlphaSError("Quark masses are not set, required for using the ODE solver with a variable flavor scheme.");
const double qmass = quark->second;
double lnmm = log(t/sqr(qmass));
double as1 = 0, as2 = 0, as3 = 0, as4 = 0;
if ( ni > nf ) {
as1 = - 0.166666*lnmm*as;
as2 = (0.152778 - 0.458333*lnmm + 0.0277778*lnmm*lnmm)*as*as;
as3 = (0.972057 - 0.0846515*nf + (- 1.65799 + 0.116319*nf)*lnmm +
(0.0920139 - 0.0277778*nf)*lnmm*lnmm - 0.00462963*lnmm*lnmm*lnmm)*as*as*as;
as4 = (5.17035 - 1.00993*nf - 0.0219784*nf*nf + (- 8.42914 + 1.30983*nf + 0.0367852*nf*nf)*lnmm +
(0.629919 - 0.143036*nf + 0.00371335*nf*nf)*lnmm*lnmm + (-0.181617 - 0.0244985*nf + 0.00308642*nf*nf)*lnmm*lnmm*lnmm +
0.000771605*lnmm*lnmm*lnmm*lnmm)*as*as*as*as;
} else {
as1 = 0.166667*lnmm*as;
as2 = (- 0.152778 + 0.458333*lnmm + 0.0277778*lnmm*lnmm)*as*as;
as3 = (- 0.972057 + 0.0846515*ni + (1.53067 - 0.116319*ni)*lnmm +
(0.289931 + 0.0277778*ni)*lnmm*lnmm + 0.00462963*lnmm*lnmm*lnmm)*as*as*as;
as4 = (- 5.10032 + 1.00993*ni + 0.0219784*ni*ni + (7.03696 - 1.22518*ni - 0.0367852*ni*ni)*lnmm +
(1.59462 + 0.0267168*ni + 0.00371335*ni*ni)*lnmm*lnmm + (0.280575 + 0.0522762*ni - 0.00308642*ni*ni)*lnmm*lnmm*lnmm +
0.000771605*lnmm*lnmm*lnmm*lnmm)*as*as*as*as;
}
double decoupling = 1.;
decoupling += as1;
if ( _qcdorder == 1 ) return decoupling;
decoupling += as2;
if ( _qcdorder == 2 ) return decoupling;
decoupling += as3;
if ( _qcdorder == 3 ) return decoupling;
decoupling += as4;
return decoupling;
}
// Calculate the next step, using recursion to achieve
// adaptive step size. Passing by reference explained
// below.
void AlphaS_ODE::_rk4(double& t, double& y, double h,
const double allowed_change, const vector<double>& bs) const {
// Determine increments in y based on the slopes of the function at the
// beginning, midpoint, and end of the interval
const double k1 = h * _derivative(t, y, bs);
const double k2 = h * _derivative(t + h/2.0, y + k1/2.0, bs);
const double k3 = h * _derivative(t + h/2.0, y + k2/2.0, bs);
const double k4 = h * _derivative(t + h, y + k3, bs);
const double change = (k1 + 2*k2 + 2*k3 + k4) / 6.0;
// Only call recursively if Q2 > 1 GeV (constant step
// size after this)
if (t > 1. && fabs(change) > allowed_change) {
_rk4(t, y, h/2., allowed_change, bs);
} else {
y += change;
t += h;
}
}
// Solve for alpha_s(q2) using alpha_s = y and Q2 = t as starting points
// Pass y and t by reference and change them to avoid having to
// return them in a container -- bit confusing but should be more
// efficient
void AlphaS_ODE::_solve(double q2, double& t, double& y,
const double& allowed_relative, double h, double accuracy) const {
if ( q2 == t ) return;
while (fabs(q2 - t) > accuracy) {
/// Can make the allowed change smaller as the q2 scale gets greater, does not seem necessary
const double allowed_change = allowed_relative;
/// Mechanism to shrink the steps if accuracy < stepsize and we are close to Q2
if (fabs(h) > accuracy && fabs(q2 - t)/h < 10 && t > 1.) h = accuracy/2.1;
/// Take constant steps for Q2 < 1 GeV
if (fabs(h) > 0.01 && t < 1.) { accuracy = 0.0051; h = 0.01; }
// Check if we are going to run forward or backwards in energy scale towards target Q2.
/// @todo C++11's std::copysign would be perfect here
if ((q2 < t && sgn(h) > 0) || (q2 > t && sgn(h) < 0)) h *= -1;
// Get beta coefficients
const vector<double> bs = _betas(numFlavorsQ2(t));
// Calculate next step
_rk4(t, y, h, allowed_change, bs);
if (y > 2.) { y = std::numeric_limits<double>::max(); break; }
}
}
/// Interpolate to get Alpha_S if the ODE has been solved,
/// otherwise solve ODE from scratch
double AlphaS_ODE::alphasQ2(double q2) const {
// Tabulate ODE solutions for interpolation and return interpolated value
_interpolate();
return _ipol.alphasQ2(q2);
// // Or directly return the ODE result (for testing)
// double h = 2;
// const double allowed_relative = 0.01;
// const double faccuracy = 0.01;
// double accuracy = faccuracy;
// double t = sqr(_mz); // starting point
// double y = _alphas_mz; // starting value
// _solve(q2, t, y, allowed_relative, h, accuracy); // solve ODE
// return y;
}
// Solve the differential equation in alphaS using an implementation of RK4
void AlphaS_ODE::_interpolate() const {
if ( _calculated ) return;
// Initial step size
double h = 2.0;
/// This the the relative error allowed for the adaptive step size.
const double allowed_relative = 0.01;
/// Accuracy of Q2 (error in Q2 within this / 2)
double accuracy = 0.001;
// Run in Q2 using RK4 algorithm until we are within our defined accuracy
double t;
double y;
if (_customref) {
t = sqr(_mreference); // starting point
y = _alphas_reference; // starting value
} else {
t = sqr(_mz); // starting point
y = _alphas_mz; // starting value
}
// If a vector of knots in q2 has been given, solve for those.
// This creates a default grid which should be overkill for most
// purposes
if ( _q2s.empty() ) {
for (int q = 1; (q/10.) < 1; ++q) {
_q2s.push_back(sqr(q/10.));
}
for (int q = 4; (q/4.) < _mz; ++q) {
_q2s.push_back(sqr(q/4.));
}
for (int q = ceil(_mz/4); (4*q) < 1000; ++q) {
_q2s.push_back(sqr(4*q));
}
for (int q = (1000/50); (50*q) < 2000; ++q) {
_q2s.push_back(sqr(50*q));
}
if ( _flavorthresholds.empty() ) {
for (int it = 4; it <= 6; ++it) {
std::map<int, double>::const_iterator element = _quarkmasses.find(it);
if ( element == _quarkmasses.end() ) continue;
_q2s.push_back(sqr(element->second)); _q2s.push_back(sqr(element->second));
}
} else {
for (int it = 4; it <= 6; ++it) {
std::map<int, double>::const_iterator element = _flavorthresholds.find(it);
if ( element == _flavorthresholds.end() ) continue;
_q2s.push_back(sqr(element->second)); _q2s.push_back(sqr(element->second));
}
}
sort(_q2s.begin(),_q2s.end());
}
// If for some reason the highest q2 knot is below m_{Z},
// force a knot there anyway (since we know it, might as well
// use it)
if ( _q2s[_q2s.size()-1] < sqr(_mz) ) _q2s.push_back(sqr(_mz));
// Find the index of the knot right below m_{Z}
unsigned int index_of_mz_lower = 0;
while ( _q2s[index_of_mz_lower + 1] < sqr(_mz) ) {
if ( index_of_mz_lower == _q2s.size() -1 ) break;
index_of_mz_lower++;
}
vector<pair<int, double> > grid; // for storing in correct order
grid.reserve(_q2s.size());
double low_lim = 0;
double last_val = -1;
bool threshold = false;
// We do this by starting from m_{Z}, going down to the lowest q2,
// and then jumping back up to m_{Z} to avoid calculating things twice
for ( int ind = index_of_mz_lower; ind >= 0; --ind) {
const double q2 = _q2s[ind];
// Deal with cases with two identical adjacent points (thresholds) by decreasing step size,
// allowed errors, and accuracy.
if ( ind != 0 && ind != 1 ) {
if ( q2 == _q2s[ind-1] ) {
last_val = q2;
threshold = true;
_solve(q2, t, y, allowed_relative/5, h/5, accuracy/5);
grid.push_back(make_pair(ind, y));
y = y * _decouple(y, t, numFlavorsQ2(_q2s[ind+1]), numFlavorsQ2(_q2s[ind-2]));
// Define divergence after y > 2. -- we have no accuracy after that any way
if ( y > 2. ) { low_lim = q2; }
continue;
}
}
// If q2 is lower than a value that already diverged, it will also diverge
if ( q2 < low_lim ) {
grid.push_back(make_pair(ind, std::numeric_limits<double>::max()));
continue;
// If last point was the same we don't need to recalculate
} else if ( q2 == last_val ) {
grid.push_back(make_pair(ind, y));
continue;
// Else calculate
} else {
last_val = q2;
if ( threshold ) { _solve(q2, t, y, allowed_relative/5, h/5, accuracy/5); threshold = false; }
else { _solve(q2, t, y, allowed_relative, h, accuracy); }
grid.push_back(make_pair(ind, y));
// Define divergence after y > 2. -- we have no accuracy after that any way
if ( y > 2. ) { low_lim = q2; }
}
}
if (_customref) {
t = sqr(_mreference); // starting point
y = _alphas_reference; // starting value
} else {
t = sqr(_mz); // starting point
y = _alphas_mz; // starting value
}
for ( size_t ind = index_of_mz_lower + 1; ind < _q2s.size(); ++ind) {
double q2 = _q2s[ind];
// Deal with cases with two identical adjacent points (thresholds) by decreasing step size,
// allowed errors, and accuracy.
if ( ind != _q2s.size() - 1 && ind != _q2s.size() - 2 ) {
if ( q2 == _q2s[ind+1] ) {
last_val = q2;
_solve(q2, t, y, allowed_relative/5, h/5, accuracy/5);
grid.push_back(make_pair(ind, y));
y = y * _decouple(y, t, numFlavorsQ2(_q2s[ind-1]), numFlavorsQ2(_q2s[ind+2]));
// Define divergence after y > 2. -- we have no accuracy after that any way
if ( y > 2. ) { low_lim = q2; }
continue;
}
}
// If q2 is lower than a value that already diverged, it will also diverge
if ( q2 < low_lim ) {
grid.push_back(make_pair(ind, std::numeric_limits<double>::max()));
continue;
// If last point was the same we don't need to recalculate
} else if ( q2 == last_val ) {
grid.push_back(make_pair(ind, y));
continue;
// Else calculate
} else {
last_val = q2;
_solve(q2, t, y, allowed_relative, h, accuracy);
grid.push_back(make_pair(ind, y));
// Define divergence after y > 2. -- we have no accuracy after that any way
if ( y > 2. ) { low_lim = q2; }
}
}
std::sort(grid.begin(), grid.end(),
boost::bind(&std::pair<int, double>::first, _1) < boost::bind(&std::pair<int, double>::first, _2));
vector<double> alphas;
alphas.reserve(_q2s.size());
for ( size_t x = 0; x < grid.size(); ++x ) {
// cout << sqrt(_q2s.at(x)) << " " << grid.at(x).second << endl;
alphas.push_back(grid.at(x).second);
}
_ipol.setQ2Values(_q2s);
_ipol.setAlphaSValues(alphas);
_calculated = true;
}
}
diff --git a/src/BicubicInterpolator.cc b/src/BicubicInterpolator.cc
--- a/src/BicubicInterpolator.cc
+++ b/src/BicubicInterpolator.cc
@@ -1,126 +1,126 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/BicubicInterpolator.h"
#include <iostream>
namespace LHAPDF {
namespace { // Unnamed namespace
// One-dimensional linear interpolation for y(x)
inline double _interpolateLinear(double x, double xl, double xh, double yl, double yh) {
assert(x >= xl);
assert(xh >= x);
return yl + (x - xl) / (xh - xl) * (yh - yl);
}
// One-dimensional cubic interpolation
inline double _interpolateCubic(double T, double VL, double VDL, double VH, double VDH) {
// Pre-calculate powers of T
const double t2 = T*T;
const double t3 = t2*T;
// Calculate left point
const double p0 = (2*t3 - 3*t2 + 1)*VL;
const double m0 = (t3 - 2*t2 + T)*VDL;
// Calculate right point
const double p1 = (-2*t3 + 3*t2)*VH;
const double m1 = (t3 - t2)*VDH;
return p0 + m0 + p1 + m1;
}
// Provides d/dx at all grid locations
double _ddx(const KnotArray1F& subgrid, size_t ix, size_t iq2) {
/// @todo Re-order this if so that branch prediction will favour the "normal" central case
if (ix == 0) { //< If at leftmost edge, use forward difference
return (subgrid.xf(ix+1, iq2) - subgrid.xf(ix, iq2)) / (subgrid.xs()[ix+1] - subgrid.xs()[ix]);
} else if (ix == subgrid.xs().size() - 1) { //< If at rightmost edge, use backward difference
return (subgrid.xf(ix, iq2) - subgrid.xf(ix-1, iq2)) / (subgrid.xs()[ix] - subgrid.xs()[ix-1]);
} else { //< If central, use the central difference
const double lddx = (subgrid.xf(ix, iq2) - subgrid.xf(ix-1, iq2)) / (subgrid.xs()[ix] - subgrid.xs()[ix-1]);
const double rddx = (subgrid.xf(ix+1, iq2) - subgrid.xf(ix, iq2)) / (subgrid.xs()[ix+1] - subgrid.xs()[ix]);
return (lddx + rddx) / 2.0;
}
}
}
double BicubicInterpolator::_interpolateXQ2(const KnotArray1F& subgrid, double x, size_t ix, double q2, size_t iq2) const {
if (subgrid.logxs().size() < 4)
throw GridError("PDF subgrids are required to have at least 4 x-knots for use with BicubicInterpolator");
if (subgrid.logq2s().size() < 4) {
if (subgrid.logq2s().size() > 1) {
// Fallback to BilinearInterpolator if either 2 or 3 Q2-knots
// First interpolate in x
const double f_ql = _interpolateLinear(x, subgrid.xs()[ix], subgrid.xs()[ix+1], subgrid.xf(ix, iq2), subgrid.xf(ix+1, iq2));
const double f_qh = _interpolateLinear(x, subgrid.xs()[ix], subgrid.xs()[ix+1], subgrid.xf(ix, iq2+1), subgrid.xf(ix+1, iq2+1));
// Then interpolate in Q2, using the x-ipol results as anchor points
return _interpolateLinear(q2, subgrid.q2s()[iq2], subgrid.q2s()[iq2+1], f_ql, f_qh);
} else throw GridError("PDF subgrids are required to have at least 2 Q2-knots for use with BicubicInterpolator");
}
/// @todo Allow interpolation right up to the borders of the grid in Q2 and x... the last inter-knot range is currently broken
/// @todo Also treat the x top/bottom edges carefully, cf. the Q2 ones
// Distance parameters
const double dx = subgrid.xs()[ix+1] - subgrid.xs()[ix];
const double tx = (x - subgrid.xs()[ix]) / dx;
/// @todo Only compute these if the +1 and +2 indices are guaranteed to be valid
const double dq_0 = subgrid.q2s()[iq2] - subgrid.q2s()[iq2-1];
const double dq_1 = subgrid.q2s()[iq2+1] - subgrid.q2s()[iq2];
const double dq_2 = subgrid.q2s()[iq2+2] - subgrid.q2s()[iq2+1];
const double dq = dq_1;
const double tq = (q2 - subgrid.q2s()[iq2]) / dq;
// Points in Q2
double vl = _interpolateCubic(tx, subgrid.xf(ix, iq2), _ddx(subgrid, ix, iq2) * dx,
subgrid.xf(ix+1, iq2), _ddx(subgrid, ix+1, iq2) * dx);
double vh = _interpolateCubic(tx, subgrid.xf(ix, iq2+1), _ddx(subgrid, ix, iq2+1) * dx,
subgrid.xf(ix+1, iq2+1), _ddx(subgrid, ix+1, iq2+1) * dx);
// Derivatives in Q2
double vdl, vdh;
if (iq2 == 0) {
// Forward difference for lower q
vdl = (vh - vl) / dq_1;
// Central difference for higher q
double vhh = _interpolateCubic(tx, subgrid.xf(ix, iq2+2), _ddx(subgrid, ix, iq2+2) * dx,
subgrid.xf(ix+1, iq2+2), _ddx(subgrid, ix+1, iq2+2) * dx);
vdh = (vdl + (vhh - vh)/dq_2) / 2.0;
}
else if (iq2+1 == subgrid.q2s().size()-1) {
// Backward difference for higher q
vdh = (vh - vl) / dq_1;
// Central difference for lower q
double vll = _interpolateCubic(tx, subgrid.xf(ix, iq2-1), _ddx(subgrid, ix, iq2-1) * dx,
subgrid.xf(ix+1, iq2-1), _ddx(subgrid, ix+1, iq2-1) * dx);
vdl = (vdh + (vl - vll)/dq_0) / 2.0;
}
else {
// Central difference for both q
double vll = _interpolateCubic(tx, subgrid.xf(ix, iq2-1), _ddx(subgrid, ix, iq2-1) * dx,
subgrid.xf(ix+1, iq2-1), _ddx(subgrid, ix+1, iq2-1) * dx);
vdl = ( (vh - vl)/dq_1 + (vl - vll)/dq_0 ) / 2.0;
double vhh = _interpolateCubic(tx, subgrid.xf(ix, iq2+2), _ddx(subgrid, ix, iq2+2) * dx,
subgrid.xf(ix+1, iq2+2), _ddx(subgrid, ix+1, iq2+2) * dx);
vdh = ( (vh - vl)/dq_1 + (vhh - vh)/dq_2 ) / 2.0;
}
vdl *= dq;
vdh *= dq;
return _interpolateCubic(tq, vl, vdl, vh, vdh);
}
}
diff --git a/src/BilinearInterpolator.cc b/src/BilinearInterpolator.cc
--- a/src/BilinearInterpolator.cc
+++ b/src/BilinearInterpolator.cc
@@ -1,36 +1,36 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/BilinearInterpolator.h"
namespace LHAPDF {
namespace { // Unnamed namespace
// One-dimensional linear interpolation for y(x)
inline double _interpolateLinear(double x, double xl, double xh, double yl, double yh) {
assert(x >= xl);
assert(xh >= x);
return yl + (x - xl) / (xh - xl) * (yh - yl);
}
}
double BilinearInterpolator::_interpolateXQ2(const KnotArray1F& subgrid, double x, size_t ix, double q2, size_t iq2) const {
if (subgrid.logxs().size() < 2)
throw GridError("PDF subgrids are required to have at least 2 x-knots for use with BilinearInterpolator");
if (subgrid.logq2s().size() < 2)
throw GridError("PDF subgrids are required to have at least 2 Q2-knots for use with BilinearInterpolator");
// First interpolate in x
const double f_ql = _interpolateLinear(x, subgrid.xs()[ix], subgrid.xs()[ix+1], subgrid.xf(ix, iq2), subgrid.xf(ix+1, iq2));
const double f_qh = _interpolateLinear(x, subgrid.xs()[ix], subgrid.xs()[ix+1], subgrid.xf(ix, iq2+1), subgrid.xf(ix+1, iq2+1));
// Then interpolate in Q2, using the x-ipol results as anchor points
return _interpolateLinear(q2, subgrid.q2s()[iq2], subgrid.q2s()[iq2+1], f_ql, f_qh);
}
}
diff --git a/src/CompositePDF.cc b/src/CompositePDF.cc
--- a/src/CompositePDF.cc
+++ b/src/CompositePDF.cc
@@ -1,32 +1,32 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/CompositePDF.h"
#include "LHAPDF/Factories.h"
using namespace std;
namespace LHAPDF {
typedef pair<string, int> SetNameMem;
CompositePDF::CompositePDF(const vector<SetNameMem>& setnames_members) {
for (const SetNameMem& setname_member : setnames_members) {
PDF* pdf = mkPDF(setname_member.first, setname_member.second);
addConstituentPDF(pdf);
}
}
CompositePDF::CompositePDF(const vector<int>& lhaids) {
for (int lhapdfid : lhaids) {
PDF* pdf = mkPDF(lhapdfid);
addConstituentPDF(pdf);
}
}
}
diff --git a/src/Config.cc b/src/Config.cc
--- a/src/Config.cc
+++ b/src/Config.cc
@@ -1,25 +1,24 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
-
#include "LHAPDF/Config.h"
#include "LHAPDF/Version.h"
using namespace std;
namespace LHAPDF {
Config::~Config() {
// Emit citation information at the end of the job, via the Config destructor
// std::cout << "CONFIG DESTRUCTION" << std::endl;
if (verbosity() > 0) {
cout << "Thanks for using LHAPDF " << version() << ". Please make sure to cite the paper:\n";
cout << " Eur.Phys.J. C75 (2015) 3, 132 (http://arxiv.org/abs/1412.7420)" << endl;
}
}
}
diff --git a/src/ContinuationExtrapolator.cc b/src/ContinuationExtrapolator.cc
--- a/src/ContinuationExtrapolator.cc
+++ b/src/ContinuationExtrapolator.cc
@@ -1,124 +1,124 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/ContinuationExtrapolator.h"
#include "LHAPDF/GridPDF.h"
namespace LHAPDF {
namespace { // Unnamed namespace
// One-dimensional linear extrapolation for y(x).
// Extrapolate in log(x) rather than just in x.
inline double _extrapolateLinear(double x, double xl, double xh, double yl, double yh) {
if (yl > 1e-3 && yh > 1e-3) {
// If yl and yh are sufficiently positive, keep y positive by extrapolating log(y).
return exp(log(yl) + (log(x) - log(xl)) / (log(xh) - log(xl)) * (log(yh) - log(yl)));
} else {
// Otherwise just extrapolate y itself.
return yl + (log(x) - log(xl)) / (log(xh) - log(xl)) * (yh - yl);
}
}
}
double ContinuationExtrapolator::extrapolateXQ2(int id, double x, double q2) const {
// The ContinuationExtrapolator provides an implementation of the extrapolation used in
// the MSTW standalone code (and LHAPDFv5 when using MSTW sets), G. Watt, October 2014.
const size_t nxknots = pdf().xKnots().size(); // total number of x knots (all subgrids)
const size_t nq2knots = pdf().q2Knots().size(); // total number of q2 knots (all subgrids)
const double xMin = pdf().xKnots()[0]; // first x knot
const double xMin1 = pdf().xKnots()[1]; // second x knot
const double xMax = pdf().xKnots()[nxknots-1]; // last x knot
const double q2Min = pdf().q2Knots()[0]; // first q2 knot
const double q2Max1 = pdf().q2Knots()[nq2knots-2]; // second-last q2 knot
const double q2Max = pdf().q2Knots()[nq2knots-1]; // last q2 knot
double fxMin, fxMin1, fq2Max, fq2Max1, fq2Min, fq2Min1, xpdf, anom;
if (x < xMin && (q2 >= q2Min && q2 <= q2Max)) {
// Extrapolation in small x only.
fxMin = pdf().interpolator().interpolateXQ2(id, xMin, q2); // PDF at (xMin,q2)
fxMin1 = pdf().interpolator().interpolateXQ2(id, xMin1, q2); // PDF at (xMin1,q2)
xpdf = _extrapolateLinear(x, xMin, xMin1, fxMin, fxMin1); // PDF at (x,q2)
} else if ((x >= xMin && x <= xMax) && q2 > q2Max) {
// Extrapolation in large q2 only.
fq2Max = pdf().interpolator().interpolateXQ2(id, x, q2Max); // PDF at (x,q2Max)
fq2Max1 = pdf().interpolator().interpolateXQ2(id, x, q2Max1); // PDF at (x,q2Max1)
xpdf = _extrapolateLinear(q2, q2Max, q2Max1, fq2Max, fq2Max1); // PDF at (x,q2)
} else if (x < xMin && q2 > q2Max) {
// Extrapolation in large q2 AND small x.
fq2Max = pdf().interpolator().interpolateXQ2(id, xMin, q2Max); // PDF at (xMin,q2Max)
fq2Max1 = pdf().interpolator().interpolateXQ2(id, xMin, q2Max1); // PDF at (xMin,q2Max1)
fxMin = _extrapolateLinear(q2, q2Max, q2Max1, fq2Max, fq2Max1); // PDF at (xMin,q2)
fq2Max = pdf().interpolator().interpolateXQ2(id, xMin1, q2Max); // PDF at (xMin1,q2Max)
fq2Max1 = pdf().interpolator().interpolateXQ2(id, xMin1, q2Max1); // PDF at (xMin1,q2Max1)
fxMin1 = _extrapolateLinear(q2, q2Max, q2Max1, fq2Max, fq2Max1); // PDF at (xMin1,q2)
xpdf = _extrapolateLinear(x, xMin, xMin1, fxMin, fxMin1); // PDF at (x,q2)
} else if (q2 < q2Min && x <= xMax) {
// Extrapolation in small q2.
if (x < xMin) {
// Extrapolation also in small x.
fxMin = pdf().interpolator().interpolateXQ2(id, xMin, q2Min); // PDF at (xMin,q2Min)
fxMin1 = pdf().interpolator().interpolateXQ2(id, xMin1, q2Min); // PDF at (xMin1,q2Min)
fq2Min = _extrapolateLinear(x, xMin, xMin1, fxMin, fxMin1); // PDF at (x,q2Min)
fxMin = pdf().interpolator().interpolateXQ2(id, xMin, 1.01*q2Min); // PDF at (xMin,1.01*q2Min)
fxMin1 = pdf().interpolator().interpolateXQ2(id, xMin1, 1.01*q2Min); // PDF at (xMin1,1.01*q2Min)
fq2Min1 = _extrapolateLinear(x, xMin, xMin1, fxMin, fxMin1); // PDF at (x,1.01*q2Min)
} else {
// Usual interpolation in x.
fq2Min = pdf().interpolator().interpolateXQ2(id, x, q2Min); // PDF at (x,q2Min)
fq2Min1 = pdf().interpolator().interpolateXQ2(id, x, 1.01*q2Min); // PDF at (x,1.01*q2Min)
}
// Calculate the anomalous dimension, dlog(f)/dlog(q2),
// evaluated at q2Min. Then extrapolate the PDFs to low
// q2 < q2Min by interpolating the anomalous dimension between
// the value at q2Min and a value of 1 for q2 << q2Min.
// If value of PDF at q2Min is very small, just set
// anomalous dimension to 1 to prevent rounding errors.
// Impose minimum anomalous dimension of -2.5.
if (abs(fq2Min) >= 1e-5) {
// anom = dlog(f)/dlog(q2) = q2/f * df/dq2 evaluated at q2 = q2Min,
// where derivative df/dq2 = ( f(1.01*q2Min) - f(q2Min) ) / (0.01*q2Min).
anom = max( -2.5, (fq2Min1 - fq2Min) / fq2Min / 0.01 );
} else anom = 1.0;
// Interpolates between f(q2Min)*(q2/q2Min)^anom for q2 ~ q2Min and
// f(q2Min)*(q2/q2Min) for q2 << q2Min, i.e. PDFs vanish as q2 --> 0.
xpdf = fq2Min * pow( q2/q2Min, anom*q2/q2Min + 1.0 - q2/q2Min );
}
else throw LogicError("We shouldn't be able to get here!");
return xpdf;
}
}
diff --git a/src/ErrExtrapolator.cc b/src/ErrExtrapolator.cc
--- a/src/ErrExtrapolator.cc
+++ b/src/ErrExtrapolator.cc
@@ -1,18 +1,18 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/ErrExtrapolator.h"
#include "LHAPDF/Exceptions.h"
namespace LHAPDF {
double ErrExtrapolator::extrapolateXQ2(int id, double x, double q2) const {
throw RangeError("Point x=" + to_str(x) + ", Q2=" + to_str(q2) +
" is outside the PDF grid boundaries");
}
}
diff --git a/src/Factories.cc b/src/Factories.cc
--- a/src/Factories.cc
+++ b/src/Factories.cc
@@ -1,340 +1,340 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/Info.h"
#include "LHAPDF/Config.h"
#include "LHAPDF/PDFSet.h"
#include "LHAPDF/PDFInfo.h"
#include "LHAPDF/PDF.h"
#include "LHAPDF/GridPDF.h"
#include "LHAPDF/BilinearInterpolator.h"
#include "LHAPDF/BicubicInterpolator.h"
#include "LHAPDF/LogBilinearInterpolator.h"
#include "LHAPDF/LogBicubicInterpolator.h"
#include "LHAPDF/ErrExtrapolator.h"
#include "LHAPDF/NearestPointExtrapolator.h"
#include "LHAPDF/ContinuationExtrapolator.h"
#include "LHAPDF/AlphaS.h"
namespace LHAPDF {
Info& getConfig() {
return Config::get();
}
PDFSet& getPDFSet(const string& setname) {
static map<string, PDFSet> _sets;
map<string, PDFSet>::iterator it = _sets.find(setname);
if (it != _sets.end()) return it->second;
_sets[setname] = PDFSet(setname);
return _sets[setname];
}
PDFInfo* mkPDFInfo(const std::string& setname, int member) {
//return new Info(findpdfmempath(setname, member));
return new PDFInfo(setname, member);
}
PDFInfo* mkPDFInfo(int lhaid) {
const pair<string,int> setname_memid = lookupPDF(lhaid);
return mkPDFInfo(setname_memid.first, setname_memid.second);
}
PDF* mkPDF(const string& setname, int member) {
// Find the member data file and ensure that it exists
const string searchpath = findpdfmempath(setname, member);
if (searchpath.empty()) {
const int setsize = getPDFSet(setname).size();
if (member > setsize-1)
throw UserError("PDF " + setname + "/" + to_str(member) + " is out of the member range of set " + setname);
throw UserError("Can't find a valid PDF " + setname + "/" + to_str(member));
}
// First create an Info object to work out what format of PDF this is:
Info info(searchpath);
const string fmt = info.get_entry("Format");
// Then use the format information to call the appropriate concrete PDF constructor:
if (fmt == "lhagrid1") return new GridPDF(setname, member);
/// @todo Throw a deprecation error if format version is too old or new
throw FactoryError("No LHAPDF factory defined for format type '" + fmt + "'");
}
namespace { // Keep this in an unnamed namespace for now, until stable/final for API
// /// @brief Expand a string describing a chain of PDFs with wildcards into a list of definite chain strings
// ///
// /// Operation:
// /// SETA/1 * SETB/* or SETA/1 * SETB
// /// ->
// /// [SETA/1, SETB/*]
// /// ->
// /// [SETA/1 * SETB/0, SETA/1 * SETB/1, ...] (via PDFSet)
// ///
// /// @todo Handle more complex wildcards such as member ID ranges
// ///
// vector<string> expandPDFsStr(const string& pdfsstr) {
// // Split into multiplicative parts
// vector<string> parts;
// iter_split(parts, pdfsstr, boost::algorithm::first_finder(" * "));
// for (string& part : parts) part = trim(part);
// // Expand each part into a single PDF spec and append it
// vector<string> rtn;
// for (const string& part : parts) {
// // No member number corresponds to whole-set wildcard expansion in this context
// if (!contains(part, "/")) part += "/*";
// // Make a list of strings corresponding to the expansion of the current part
// vector<string> expandedpart;
// if (endswith(part, "/*")) {
// const string setname = part.substr(0, part.find("/*"));
// const PDFSet& set = getPDFSet(setname);
// for (size_t i = 0; i < set.size(); ++i)
// expandedpart.push_back(setname + "/" + to_str(i));
// } else {
// expandedpart.push_back(part);
// }
// // Append to or outer-product the return vector as appropriate
// if (expandedpart.size() == 1) {
// // If there were no wildcards to expand, no point in playing outer-product games
// for (const string& basestr : rtn) {
// const string newstr = (!basestr.empty() ? basestr + " * " : "") + expandedpart[0];
// }
// } else {
// // There was a wildcard to expand, so we work harder...
// vector<string> tmp;
// tmp.reserve(expandedpart.size() * rtn.size());
// for (const string& basestr : rtn) {
// for (const string& pdfstr : expandedpart) {
// const string newstr = (!basestr.empty() ? basestr + " * " : "") + pdfstr;
// tmp.push_back(newstr);
// }
// }
// rtn = tmp; //< Update rtn with the next level of Cartesian product from wildcard expansion
// }
// }
// return rtn;
// }
/// @brief Decode a single PDF member ID string into a setname,memid pair
///
/// SETA/1 SETA
/// -> ->
/// <SETA,1> <SETA,0>
pair<string,int> decodePDFStr(const string& pdfstr) {
int nmem = 0;
const size_t slashpos = pdfstr.find("/");
const string setname = trim(pdfstr.substr(0, slashpos));
try {
if (slashpos != string::npos) {
const string smem = pdfstr.substr(slashpos+1);
nmem = lexical_cast<int>(smem);
}
} catch (...) {
throw UserError("Could not parse PDF identity string " + pdfstr);
}
return make_pair(setname, nmem);
}
// /// @brief Decode a multiple PDF member ID string into a list of setname,memid pairs
// ///
// /// SETA/1 * SETB/0
// /// ->
// /// [<SETA,1>, <SETB,0>]
// vector< pair<string,int> > decodePDFsStr(const string& pdfsstr) {
// // Split into multiplicative parts
// vector<string> parts;
// iter_split(parts, pdfsstr, boost::algorithm::first_finder(" * "));
// // Create the list of set/id pairs
// vector< pair<string,int> > rtn;
// for (const string& part : parts) {
// rtn.push_back(decodePDFStr(part));
// }
// return rtn;
// }
}
PDF* mkPDF(const string& setname_nmem) {
const pair<string,int> idpair = decodePDFStr(setname_nmem);
return mkPDF(idpair.first, idpair.second);
}
PDF* mkPDF(int lhaid) {
const pair<string,int> setname_nmem = lookupPDF(lhaid);
return mkPDF(setname_nmem.first, setname_nmem.second);
}
void mkPDFs(const string& setname, vector<PDF*>& pdfs) {
getPDFSet(setname).mkPDFs(pdfs);
}
vector<PDF*> mkPDFs(const string& setname) {
return getPDFSet(setname).mkPDFs();
}
Interpolator* mkInterpolator(const string& name) {
// Convert name to lower case for comparisons
const string iname = to_lower(name);
if (iname == "linear")
return new BilinearInterpolator();
else if (iname == "cubic")
return new BicubicInterpolator();
else if (iname == "log")
return new LogBilinearInterpolator();
else if (iname == "logcubic")
return new LogBicubicInterpolator();
else
throw FactoryError("Undeclared interpolator requested: " + name);
}
Extrapolator* mkExtrapolator(const string& name) {
// Convert name to lower case for comparisons
const string iname = to_lower(name);
if (iname == "nearest")
return new NearestPointExtrapolator();
else if (iname == "error")
return new ErrExtrapolator();
else if (iname == "continuation")
return new ContinuationExtrapolator();
else
throw FactoryError("Undeclared extrapolator requested: " + name);
}
AlphaS* mkAlphaS(const Info& info) {
AlphaS* as = 0;
const string itype = to_lower(info.get_entry("AlphaS_Type"));
if (itype == "analytic") as = new AlphaS_Analytic();
else if (itype == "ode") as = new AlphaS_ODE();
else if (itype == "ipol") as = new AlphaS_Ipol();
else throw FactoryError("Undeclared AlphaS requested: " + itype);
// Configure the QCD params on this AlphaS
if (info.has_key("AlphaS_OrderQCD")) as->setOrderQCD(info.get_entry_as<int>("AlphaS_OrderQCD"));
/// @todo Fall back to generic OrderQCD? (Default value is 4 at the moment --KN)
/// Thresholds are where the flavor transition happens in the ODE/analytic solvers, AlphaS_MQ
/// should be given as \bar{MQ} (\bar{MQ}) if it is to be used for thresholds away from the mass.
/// Since you need the heavy quark mass to calculate the decoupling in the ODE solver
/// when threshold =/= mass, I've implemented this as AlphaS_Threshold* -> Threshold*, AlphaS_M* -> M*
if (info.has_key("AlphaS_ThresholdDown") && info.has_key("AlphaS_ThresholdUp") && info.has_key("AlphaS_ThresholdStrange")
&& info.has_key("AlphaS_ThresholdCharm") && info.has_key("AlphaS_ThresholdBottom") && info.has_key("AlphaS_ThresholdTop")) {
as->setQuarkThreshold(1, info.get_entry_as<double>("AlphaS_ThresholdDown"));
as->setQuarkThreshold(2, info.get_entry_as<double>("AlphaS_ThresholdUp"));
as->setQuarkThreshold(3, info.get_entry_as<double>("AlphaS_ThresholdStrange"));
as->setQuarkThreshold(4, info.get_entry_as<double>("AlphaS_ThresholdCharm"));
as->setQuarkThreshold(5, info.get_entry_as<double>("AlphaS_ThresholdBottom"));
as->setQuarkThreshold(6, info.get_entry_as<double>("AlphaS_ThresholdTop"));
}
else if (info.has_key("ThresholdDown") && info.has_key("ThresholdUp") && info.has_key("ThresholdStrange")
&& info.has_key("ThresholdCharm") && info.has_key("ThresholdBottom") && info.has_key("ThresholdTop")) {
as->setQuarkThreshold(1, info.get_entry_as<double>("ThresholdDown"));
as->setQuarkThreshold(2, info.get_entry_as<double>("ThresholdUp"));
as->setQuarkThreshold(3, info.get_entry_as<double>("ThresholdStrange"));
as->setQuarkThreshold(4, info.get_entry_as<double>("ThresholdCharm"));
as->setQuarkThreshold(5, info.get_entry_as<double>("ThresholdBottom"));
as->setQuarkThreshold(6, info.get_entry_as<double>("ThresholdTop"));
}
if (info.has_key("AlphaS_MDown") && info.has_key("AlphaS_MUp") && info.has_key("AlphaS_MStrange")
&& info.has_key("AlphaS_MCharm") && info.has_key("AlphaS_MBottom") && info.has_key("AlphaS_MTop")) {
as->setQuarkMass(1, info.get_entry_as<double>("AlphaS_MDown"));
as->setQuarkMass(2, info.get_entry_as<double>("AlphaS_MUp"));
as->setQuarkMass(3, info.get_entry_as<double>("AlphaS_MStrange"));
as->setQuarkMass(4, info.get_entry_as<double>("AlphaS_MCharm"));
as->setQuarkMass(5, info.get_entry_as<double>("AlphaS_MBottom"));
as->setQuarkMass(6, info.get_entry_as<double>("AlphaS_MTop"));
}
// This falls back to lhapdf.conf so should in theory never throw the MetadataError
else if (info.has_key("MDown") && info.has_key("MUp") && info.has_key("MStrange")
&& info.has_key("MCharm") && info.has_key("MBottom") && info.has_key("MTop")) {
as->setQuarkMass(1, info.get_entry_as<double>("MDown"));
as->setQuarkMass(2, info.get_entry_as<double>("MUp"));
as->setQuarkMass(3, info.get_entry_as<double>("MStrange"));
as->setQuarkMass(4, info.get_entry_as<double>("MCharm"));
as->setQuarkMass(5, info.get_entry_as<double>("MBottom"));
as->setQuarkMass(6, info.get_entry_as<double>("MTop"));
} else {
throw MetadataError("All quark masses required (either as AlphaS_MQ or MQ) for AlphaS.");
}
const string fscheme = to_lower(info.get_entry("AlphaS_FlavorScheme", info.get_entry("FlavorScheme", "variable"))); // default is VFNS
const int nflavs = info.get_entry_as<int>("AlphaS_NumFlavors", info.get_entry_as<int>("NumFlavors", 5)); // default is 5 flavour evolution
if (fscheme == "fixed") as->setFlavorScheme(AlphaS::FIXED, nflavs);
else if (fscheme == "variable") as->setFlavorScheme(AlphaS::VARIABLE, nflavs);
else as->setFlavorScheme(AlphaS::VARIABLE, 5); // default fallback mode
// Required parameter settings for each calculation mode
if (as->type() == "ode") {
/// @todo Handle FFNS / VFNS
if ( (!info.has_key("AlphaS_MZ") || !info.has_key("MZ")) || (!info.has_key("AlphaS_MassReference") || !info.has_key("AlphaS_Reference")) )
throw MetadataError("Requested ODE AlphaS but there is no reference point given: define either AlphaS_MZ and MZ, or AlphaS_MassReference and AlphaS_Reference. The latter is given preference if both are defined.");
if (info.has_key("AlphaS_MZ")) as->setAlphaSMZ(info.get_entry_as<double>("AlphaS_MZ"));
if (info.has_key("MZ"))as->setMZ(info.get_entry_as<double>("MZ"));
if (info.has_key("AlphaS_Reference"))as->setAlphaSReference(info.get_entry_as<double>("AlphaS_Reference"));
if (info.has_key("AlphaS_MassReference"))as->setMassReference(info.get_entry_as<double>("AlphaS_MassReference"));
if (info.has_key("AlphaS_Qs")) {
AlphaS_ODE* as_o = dynamic_cast<AlphaS_ODE*>(as);
if (info.has_key("AlphaS_Qs")) as_o->setQValues( info.get_entry_as< vector<double> >("AlphaS_Qs"));
}
}
else if (as->type() == "analytic") {
/// @todo Handle FFNS / VFNS
if (!info.has_key("AlphaS_Lambda5") && !info.has_key("AlphaS_Lambda4") && !info.has_key("AlphaS_Lambda3") )
throw MetadataError("Requested analytic AlphaS but the required parameters are not defined.");
if (info.has_key("AlphaS_Lambda3")) as->setLambda(3, info.get_entry_as<double>("AlphaS_Lambda3"));
if (info.has_key("AlphaS_Lambda4")) as->setLambda(4, info.get_entry_as<double>("AlphaS_Lambda4"));
if (info.has_key("AlphaS_Lambda5")) as->setLambda(5, info.get_entry_as<double>("AlphaS_Lambda5"));
}
else if (as->type() == "ipol") {
if (!info.has_key("AlphaS_Qs") || !info.has_key("AlphaS_Vals") )
throw MetadataError("Requested ipol AlphaS but the required parameters are not defined.");
AlphaS_Ipol* as_i = dynamic_cast<AlphaS_Ipol*>(as);
if (info.has_key("AlphaS_Qs")) as_i->setQValues( info.get_entry_as< vector<double> >("AlphaS_Qs"));
if (info.has_key("AlphaS_Vals")) as_i->setAlphaSValues( info.get_entry_as< vector<double> >("AlphaS_Vals"));
}
return as;
}
AlphaS* mkAlphaS(const std::string& setname) {
return mkAlphaS(getPDFSet(setname));
}
AlphaS* mkAlphaS(const std::string& setname, int member) {
unique_ptr<Info> info( mkPDFInfo(setname, member) );
return mkAlphaS(*info);
}
AlphaS* mkAlphaS(int lhaid) {
unique_ptr<Info> info( mkPDFInfo(lhaid) );
return mkAlphaS(*info);
}
}
diff --git a/src/GridPDF.cc b/src/GridPDF.cc
--- a/src/GridPDF.cc
+++ b/src/GridPDF.cc
@@ -1,198 +1,198 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/GridPDF.h"
#include "LHAPDF/Interpolator.h"
#include "LHAPDF/Factories.h"
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <stdexcept>
#include <cstring>
using namespace std;
namespace LHAPDF {
double GridPDF::_xfxQ2(int id, double x, double q2) const {
/// Decide whether to use interpolation or extrapolation... the sanity checks
/// are done in the public PDF::xfxQ2 function.
// cout << "From GridPDF[0]: x = " << x << ", Q2 = " << q2 << endl;
double xfx = 0;
if (inRangeXQ2(x, q2)) {
// cout << "From GridPDF[ipol]: x = " << x << ", Q2 = " << q2 << endl;
// cout << "Num subgrids = " << _knotarrays.size() << endl;
// int i = 0;
// for (std::map<double, KnotArrayNF>::const_iterator it = _knotarrays.begin(); it != _knotarrays.end(); ++it)
// cout << "#" << i++ << " from Q = " << sqrt(it->first) << endl;
xfx = interpolator().interpolateXQ2(id, x, q2);
} else {
// cout << "From GridPDF[xpol]: x = " << x << ", Q2 = " << q2 << endl;
xfx = extrapolator().extrapolateXQ2(id, x, q2);
}
return xfx;
}
namespace {
// A wrapper for std::strtod and std::strtol, for fast tokenizing when all
// input is guaranteed to be numeric (as in this data block). Based very
// closely on FastIStringStream by Gavin Salam.
class NumParser {
public:
// Constructor from char*
NumParser(const char* line=0) { reset(line); }
// Constructor from std::string
NumParser(const string& line) { reset(line); }
// Re-init to new line as char*
void reset(const char* line=0) {
_next = const_cast<char*>(line);
_new_next = _next;
_error = false;
}
// Re-init to new line as std::string
void reset(const string& line) { reset(line.c_str()); }
// Tokenizing stream operator (forwards to double and int specialisations)
template<class T> NumParser& operator>>(T& value) {
_get(value);
if (_new_next == _next) _error = true; // handy error condition behaviour!
_next = _new_next;
return *this;
}
// Allow use of operator>> in a while loop
operator bool() const { return !_error; }
private:
void _get(double& x) { x = std::strtod(_next, &_new_next); }
void _get(float& x) { x = std::strtof(_next, &_new_next); }
void _get(int& i) { i = std::strtol(_next, &_new_next, 10); } // force base 10!
char *_next, *_new_next;
bool _error;
};
}
void GridPDF::_loadData(const std::string& mempath) {
string line, prevline;
int iblock(0), iblockline(0), iline(0);
vector<double> xs, q2s;
vector<int> pids;
vector< vector<double> > ipid_xfs;
try {
ifstream file(mempath.c_str());
NumParser nparser; double ftoken; int itoken;
while (getline(file, line)) {
// Trim the current line to ensure that there is no effect of leading spaces, etc.
line = trim(line);
prevline = line; // used to test the last line after the while loop fails
// If the line is commented out, increment the line number but not the block line
iline += 1;
if (line.find("#") == 0) continue;
iblockline += 1;
if (line != "---") { // if we are not on a block separator line...
// Block 0 is the metadata, which we ignore here
if (iblock == 0) continue;
// Debug printout
// cout << iline << " = block line #" << iblockline << " => " << line << endl;
// Parse the data lines
nparser.reset(line);
if (iblockline == 1) { // x knots line
while (nparser >> ftoken) xs.push_back(ftoken);
if (xs.empty())
throw ReadError("Empty x knot array on line " + to_str(iline));
} else if (iblockline == 2) { // Q knots line
while (nparser >> ftoken) q2s.push_back(ftoken*ftoken); // note Q -> Q2
if (q2s.empty())
throw ReadError("Empty Q knot array on line " + to_str(iline));
} else if (iblockline == 3) { // internal flavor IDs ordering line
while (nparser >> itoken) pids.push_back(itoken);
// Check that each line has many tokens as there should be flavours
if (pids.size() != flavors().size())
throw ReadError("PDF grid data error on line " + to_str(iline) + ": " + to_str(pids.size()) +
" parton flavors declared but " + to_str(flavors().size()) + " expected from Flavors metadata");
/// @todo Handle sea/valence representations via internal pseudo-PIDs
} else {
if (iblockline == 4) { // on the first line of the xf block, resize the arrays
ipid_xfs.resize(pids.size());
const size_t subgridsize = xs.size() * q2s.size();
for (size_t ipid = 0; ipid < pids.size(); ++ipid) {
ipid_xfs[ipid].reserve(subgridsize);
}
}
size_t ipid = 0;
while (nparser >> ftoken) {
ipid_xfs[ipid].push_back(ftoken);
ipid += 1;
}
// Check that each line has many tokens as there should be flavours
if (ipid != pids.size())
throw ReadError("PDF grid data error on line " + to_str(iline) + ": " + to_str(ipid) +
" flavor entries seen but " + to_str(pids.size()) + " expected");
}
} else { // we *are* on a block separator line
// Check that the expected number of data lines were seen in the last block
if (iblock > 0 && iblockline - 1 != int(xs.size()*q2s.size()) + 3)
throw ReadError("PDF grid data error on line " + to_str(iline) + ": " +
to_str(iblockline-1) + " data lines were seen in block " + to_str(iblock-1) +
" but " + to_str(xs.size()*q2s.size() + 3) + " expected");
// Ignore block registration if we've just finished reading the 0th (metadata) block
if (iblock > 0) {
// Throw if the last subgrid block was of zero size
if (ipid_xfs.empty())
throw ReadError("Empty xf values array in data block " + to_str(iblock) + ", ending on line " + to_str(iline));
// Register data from the block into the GridPDF data structure
KnotArrayNF& arraynf = _knotarrays[q2s.front()]; //< Reference to newly created subgrid object
for (size_t ipid = 0; ipid < pids.size(); ++ipid) {
const int pid = pids[ipid];
// Create the 2D array with the x and Q2 knot positions
arraynf[pid] = KnotArray1F(xs, q2s);
// Populate the xf data array
arraynf[pid].xfs().assign(ipid_xfs[ipid].begin(), ipid_xfs[ipid].end());
}
}
// Increment/reset the block and line counters, subgrid arrays, etc.
iblock += 1;
iblockline = 0;
xs.clear(); q2s.clear();
for (size_t ipid = 0; ipid < pids.size(); ++ipid)
ipid_xfs[ipid].clear();
pids.clear();
}
}
// File reading finished: complain if it was not properly terminated
if (prevline != "---")
throw ReadError("Grid file " + mempath + " is not properly terminated: .dat files MUST end with a --- separator line");
// Error handling
} catch (Exception& e) {
throw;
} catch (std::exception& e) {
throw ReadError("Read error while parsing " + mempath + " as a GridPDF data file");
}
}
}
diff --git a/src/Info.cc b/src/Info.cc
--- a/src/Info.cc
+++ b/src/Info.cc
@@ -1,106 +1,106 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/Info.h"
#include "LHAPDF/PDFIndex.h"
#include "yaml-cpp/yaml.h"
#ifdef YAML_NAMESPACE
#define YAML YAML_NAMESPACE
#endif
namespace LHAPDF {
void Info::load(const string& filepath) {
// Complain if the path is empty
if (filepath.empty()) throw ReadError("Empty PDF file name given to Info::load");
// But complain if a non-empty path is provided, but it's invalid
if (!file_exists(filepath)) throw ReadError("PDF data file '" + filepath + "' not found");
// Read the YAML part of the file into the metadata map
try {
// Do the parsing "manually" up to the first doc delimiter
std::ifstream file(filepath.c_str());
YAML::Node doc;
#if YAMLCPP_API == 3
YAML::Parser parser(file);
parser.GetNextDocument(doc);
for (YAML::Iterator it = doc.begin(); it != doc.end(); ++it) {
string key, val;
it.first() >> key;
try {
// Assume the value is a scalar type -- it'll throw an exception if not
it.second() >> val;
} catch (const YAML::InvalidScalar& ex) {
// It's a list: process the entries individually into a comma-separated string
string subval;
for (size_t i = 0; i < it.second().size(); ++i) {
it.second()[i] >> subval;
val += subval + ((i < it.second().size()-1) ? "," : "");
}
}
//cout << key << ": " << val << endl;
_metadict[key] = val;
}
#elif YAMLCPP_API == 5
string docstr, line;
while (getline(file, line)) {
//cout << "@ " << line << endl;
if (line == "---") break;
docstr += line + "\n";
}
doc = YAML::Load(docstr);
for (YAML::const_iterator it = doc.begin(); it != doc.end(); ++it) {
const string key = it->first.as<string>();
const YAML::Node& val = it->second;
if (val.IsScalar()) {
// Scalar value
_metadict[key] = val.as<string>();
} else {
// Process the sequence entries into a comma-separated string
string seqstr = "";
for (size_t i = 0; i < val.size(); ++i)
seqstr += val[i].as<string>() + ((i < val.size()-1) ? "," : "");
_metadict[key] = seqstr;
}
}
#endif
} catch (const YAML::ParserException& ex) {
throw ReadError("YAML parse error in " + filepath + " :" + ex.what());
} catch (const LHAPDF::Exception& ex) {
throw;
} catch (const std::exception& ex) {
throw ReadError("Trouble when reading " + filepath + " :" + ex.what());
}
}
// /// @todo Only support loading via PDF set name and member ID, not explicit paths
// /// @todo Replace the loading of the set metadata into the member info with set-level Info singletons
// void Info::loadFull(const string& mempath) { //< @todo Need a better method name!
// // Extract the set name from the member data file path
// const string memberdata = findFile(mempath);
// if (memberdata.empty() || !file_exists(memberdata)) throw ReadError("Could not find PDF data file '" + mempath + "'");
// const string memname = basename(memberdata); //< Can use this to alternatively work out the set name...
// const string setdir = dirname(memberdata);
// const string setname = basename(setdir);
// path setinfo = findpdfsetinfopath(setname);
// // Load the set info
// if (file_exists(setinfo)) load(setinfo);
// // Load the member info (possibly overriding the set-level metadata)
// load(memberdata);
// }
}
diff --git a/src/Interpolator.cc b/src/Interpolator.cc
--- a/src/Interpolator.cc
+++ b/src/Interpolator.cc
@@ -1,29 +1,29 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/Interpolator.h"
#include "LHAPDF/GridPDF.h"
namespace LHAPDF {
double Interpolator::interpolateXQ2(int id, double x, double q2) const {
// Subgrid lookup
/// @todo Do this in two stages to cache the KnotArrayNF?
/// @todo Add flavour error checking
const KnotArray1F& subgrid = pdf().subgrid(id, q2);
// Index look-up
/// @todo Cache this index lookup for performance?
// cout << "From Ipol: x = " << x << ", Q2 = " << q2 << endl;
const size_t ix = subgrid.ixbelow(x);
const size_t iq2 = subgrid.iq2below(q2);
// cout << "ix = " << ix << ", iq2 = " << iq2 << ", xf[ix, iq2] = " << subgrid.xf(ix, iq2) << endl;
/// Call the overloaded interpolation routine on this subgrid
return _interpolateXQ2(subgrid, x, ix, q2, iq2);
}
}
diff --git a/src/LHAGlue.cc b/src/LHAGlue.cc
--- a/src/LHAGlue.cc
+++ b/src/LHAGlue.cc
@@ -1,1149 +1,1149 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/PDF.h"
#include "LHAPDF/PDFSet.h"
#include "LHAPDF/PDFIndex.h"
#include "LHAPDF/Factories.h"
#include "LHAPDF/Utils.h"
#include "LHAPDF/Paths.h"
#include "LHAPDF/Version.h"
#include "LHAPDF/LHAGlue.h"
#include <cstring>
using namespace std;
// We have to create and initialise some common blocks here for backwards compatibility
struct w50512 {
double qcdl4, qcdl5;
};
w50512 w50512_;
struct w50513 {
double xmin, xmax, q2min, q2max;
};
w50513 w50513_;
struct lhapdfr {
double qcdlha4, qcdlha5;
int nfllha;
};
lhapdfr lhapdfr_;
namespace { //< Unnamed namespace to restrict visibility to this file
/// @brief PDF object storage here is a smart pointer to ensure deletion of created PDFs
typedef std::shared_ptr<LHAPDF::PDF> PDFPtr;
/// @brief A struct for handling the active PDFs for the Fortran interface.
///
/// We operate in a string-based way, since maybe there will be sets with names, but no
/// index entry in pdfsets.index.
///
/// @todo Can we avoid the strings and just work via the LHAPDF ID and factory construction?
///
/// Smart pointers are used in the native map used for PDF member storage so
/// that they auto-delete if the PDFSetHandler that holds them goes out of
/// scope (i.e. is overwritten).
struct PDFSetHandler {
/// Default constructor
PDFSetHandler() : currentmem(0)
{ } //< It'll be stored in a map so we need one of these...
/// Constructor from a PDF set name
PDFSetHandler(const string& name)
: setname(name)
{
loadMember(0);
}
/// Constructor from a PDF set's LHAPDF ID code
PDFSetHandler(int lhaid) {
pair<string,int> set_mem = LHAPDF::lookupPDF(lhaid);
// First check that the lookup was successful, i.e. it was a valid ID for the LHAPDF6 set collection
if (set_mem.first.empty() || set_mem.second < 0)
throw LHAPDF::UserError("Could not find a valid PDF with LHAPDF ID = " + LHAPDF::to_str(lhaid));
// Try to load this PDF (checking that the member number is in the set's range is done in mkPDF, called by loadMember)
setname = set_mem.first;
loadMember(set_mem.second);
}
/// @brief Load a new PDF member
///
/// If it's already loaded, the existing object will not be reloaded.
void loadMember(int mem) {
if (mem < 0)
throw LHAPDF::UserError("Tried to load a negative PDF member ID: " + LHAPDF::to_str(mem) + " in set " + setname);
if (members.find(mem) == members.end())
members[mem] = PDFPtr(LHAPDF::mkPDF(setname, mem));
currentmem = mem;
}
/// Actively delete a PDF member to save memory
void unloadMember(int mem) {
members.erase(mem);
const int nextmem = (!members.empty()) ? members.begin()->first : 0;
loadMember(nextmem);
}
/// @brief Get a PDF member
///
/// Non-const because it can secretly load the member. Not that constness
/// matters in a Fortran interface utility function!
const PDFPtr member(int mem) {
loadMember(mem);
return members.find(mem)->second;
}
/// Get the currently active PDF member
///
/// Non-const because it can secretly load the member. Not that constness
/// matters in a Fortran interface utility function!
const PDFPtr activemember() {
return member(currentmem);
}
/// The currently active member in this set
int currentmem;
/// Name of this set
string setname;
/// Map of pointers to selected member PDFs
///
// /// It's mutable so that a "const" member-getting operation can implicitly
// /// load a new PDF object. Good idea / bad idea? Disabled for now.
// mutable map<int, PDFPtr> members;
map<int, PDFPtr> members;
};
/// Collection of active sets
static map<int, PDFSetHandler> ACTIVESETS;
/// The currently active set
int CURRENTSET = 0;
/// Useful C-string -> Fortran-string converter
// Credit: https://stackoverflow.com/questions/10163485/passing-char-arrays-from-c-to-fortran
void cstr_to_fstr(const char* cstring, char* fstring, std::size_t fstring_len) {
std::size_t inlen = std::strlen(cstring);
std::size_t cpylen = std::min(inlen, fstring_len);
// TODO: truncation error or warning
//if (inlen > fstring_len) FOOOOO();
std::copy(cstring, cstring+cpylen, fstring);
std::fill(fstring+cpylen, fstring+fstring_len, ' ');
}
}
string lhaglue_get_current_pdf(int nset) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
return "NONE";
CURRENTSET = nset;
return ACTIVESETS[nset].activemember()->set().name() + " (" +
LHAPDF::to_str(ACTIVESETS[nset].activemember()->lhapdfID()) + ")";
}
extern "C" {
// NEW FORTRAN INTERFACE FUNCTIONS
/// LHAPDF library version
void lhapdf_getversion_(char* s, size_t len) {
// strncpy(s, LHAPDF_VERSION, len);
cstr_to_fstr(LHAPDF_VERSION, s, len);
}
/// List of available PDF sets, returned as a space-separated string
void lhapdf_getpdfsetlist_(char* s, size_t len) {
string liststr;
for (const string& setname : LHAPDF::availablePDFSets()) {
if (!liststr.empty()) liststr += " ";
liststr += setname;
}
// strncpy(s, liststr.c_str(), len);
cstr_to_fstr(liststr.c_str(), s, len);
}
//////////////////
// LHAPDF5 / PDFLIB COMPATIBILITY INTERFACE FUNCTIONS
// System-level info
/// LHAPDF library version
void getlhapdfversion_(char* s, size_t len) {
// strncpy(s, LHAPDF_VERSION, len);
cstr_to_fstr(LHAPDF_VERSION, s, len);
}
/// Does nothing, only provided for backward compatibility
void lhaprint_(int& a) { }
/// Set LHAPDF parameters -- does nothing in LHAPDF6!
void setlhaparm_(const char* par, int parlength) {
/// @todo Can any Fortran LHAPDF params be usefully mapped?
}
/// Get LHAPDF parameters -- does nothing in LHAPDF6!
void getlhaparm_(int dummy, char* par, int parlength) {
/// @todo Can any Fortran LHAPDF params be usefully mapped?
cstr_to_fstr("", par, parlength);
}
/// Return a dummy max number of sets (there is no limitation now)
void getmaxnumsets_(int& nmax) {
nmax = 1000;
}
/// Set PDF data path
void setpdfpath_(const char* s, size_t len) {
/// @todo Works? Need to check C-string copying, null termination
char s2[1024];
s2[len] = '\0';
strncpy(s2, s, len);
LHAPDF::pathsPrepend(s2);
}
/// Get PDF data path (colon-separated if there is more than one element)
void getdatapath_(char* s, size_t len) {
/// @todo Works? Need to check Fortran string return, string macro treatment, etc.
string pathstr;
for (const string& path : LHAPDF::paths()) {
if (!pathstr.empty()) pathstr += ":";
pathstr += path;
}
// strncpy(s, pathstr.c_str(), len);
cstr_to_fstr(pathstr.c_str(), s, len);
}
// PDF initialisation and focus-switching
/// Load a PDF set
///
/// @todo Does this version actually take a *path*? What to do?
void initpdfsetm_(const int& nset, const char* setpath, int setpathlength) {
// Strip file extension for backward compatibility
string fullp = string(setpath, setpathlength);
// Remove trailing whitespace
fullp.erase( std::remove_if( fullp.begin(), fullp.end(), ::isspace ), fullp.end() );
// Use only items after the last /
const string pap = LHAPDF::dirname(fullp);
const string p = LHAPDF::basename(fullp);
// Prepend path to search area
LHAPDF::pathsPrepend(pap);
// Handle extensions
string path = LHAPDF::file_extn(p).empty() ? p : LHAPDF::file_stem(p);
/// @note We correct the misnamed CTEQ6L1/CTEQ6ll set name as a backward compatibility special case.
if (LHAPDF::to_lower(path) == "cteq6ll") path = "cteq6l1";
// Create the PDF set with index nset
// if (ACTIVESETS.find(nset) == ACTIVESETS.end())
ACTIVESETS[nset] = PDFSetHandler(path); ///< @todo Will be wrong if a structured path is given
CURRENTSET = nset;
}
/// Load a PDF set (non-multiset version)
void initpdfset_(const char* setpath, int setpathlength) {
int nset1 = 1;
initpdfsetm_(nset1, setpath, setpathlength);
}
/// Load a PDF set by name
void initpdfsetbynamem_(const int& nset, const char* setname, int setnamelength) {
// Truncate input to size setnamelength
string p = setname;
p.erase(setnamelength, std::string::npos);
// Strip file extension for backward compatibility
string name = LHAPDF::file_extn(p).empty() ? p : LHAPDF::file_stem(p);
// Remove trailing whitespace
name.erase( std::remove_if( name.begin(), name.end(), ::isspace ), name.end() );
/// @note We correct the misnamed CTEQ6L1/CTEQ6ll set name as a backward compatibility special case.
if (LHAPDF::to_lower(name) == "cteq6ll") name = "cteq6l1";
// Create the PDF set with index nset
// if (ACTIVESETS.find(nset) == ACTIVESETS.end())
ACTIVESETS[nset] = PDFSetHandler(name);
// Update current set focus
CURRENTSET = nset;
}
/// Load a PDF set by name (non-multiset version)
void initpdfsetbyname_(const char* setname, int setnamelength) {
int nset1 = 1;
initpdfsetbynamem_(nset1, setname, setnamelength);
}
/// Load a PDF in current set
void initpdfm_(const int& nset, const int& nmember) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
ACTIVESETS[nset].loadMember(nmember);
// Update current set focus
CURRENTSET = nset;
}
/// Load a PDF in current set (non-multiset version)
void initpdf_(const int& nmember) {
int nset1 = 1;
initpdfm_(nset1, nmember);
}
/// Get the current set number (i.e. allocation slot index)
void getnset_(int& nset) {
nset = CURRENTSET;
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
}
/// Explicitly set the current set number (i.e. allocation slot index)
void setnset_(const int& nset) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
CURRENTSET = nset;
}
/// Get the current member number in slot nset
void getnmem_(int& nset, int& nmem) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
nmem = ACTIVESETS[nset].currentmem;
// Update current set focus
CURRENTSET = nset;
}
/// Set the current member number in slot nset
void setnmem_(const int& nset, const int& nmem) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" +
LHAPDF::to_str(nset) + " but it is not initialised");
ACTIVESETS[nset].loadMember(nmem);
// Update current set focus
CURRENTSET = nset;
}
// PDF evolution functions
/// Get xf(x) values for common partons from current PDF
void evolvepdfm_(const int& nset, const double& x, const double& q, double* fxq) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
// Evaluate for the 13 LHAPDF5 standard partons (-6..6)
for (size_t i = 0; i < 13; ++i) {
try {
fxq[i] = ACTIVESETS[nset].activemember()->xfxQ(i-6, x, q);
} catch (const exception& e) {
fxq[i] = 0;
}
}
// Update current set focus
CURRENTSET = nset;
}
/// Get xf(x) values for common partons from current PDF (non-multiset version)
void evolvepdf_(const double& x, const double& q, double* fxq) {
int nset1 = 1;
evolvepdfm_(nset1, x, q, fxq);
}
/// Determine if the current PDF has a photon flavour (historically only MRST2004QED)
/// @todo Function rather than subroutine?
/// @note There is no multiset version. has_photon will respect the current set slot.
bool has_photon_() {
return ACTIVESETS[CURRENTSET].activemember()->hasFlavor(22);
}
/// Get xfx values from current PDF, including an extra photon flavour
void evolvepdfphotonm_(const int& nset, const double& x, const double& q, double* fxq, double& photonfxq) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
// First evaluate the "normal" partons
evolvepdfm_(nset, x, q, fxq);
// Then evaluate the photon flavor (historically only for MRST2004QED)
try {
photonfxq = ACTIVESETS[nset].activemember()->xfxQ(22, x, q);
} catch (const exception& e) {
photonfxq = 0;
}
// Update current set focus
CURRENTSET = nset;
}
/// Get xfx values from current PDF, including an extra photon flavour (non-multiset version)
void evolvepdfphoton_(const double& x, const double& q, double* fxq, double& photonfxq) {
int nset1 = 1;
evolvepdfphotonm_(nset1, x, q, fxq, photonfxq);
}
/// Get xf(x) values for common partons from a photon PDF
void evolvepdfpm_(const int& nset, const double& x, const double& q, const double& p2, const int& ip2, double& fxq) {
// Update current set focus
CURRENTSET = nset;
throw LHAPDF::NotImplementedError("Photon structure functions are not yet supported in LHAPDF6");
}
/// Get xf(x) values for common partons from a photon PDF (non-multiset version)
void evolvepdfp_(const double& x, const double& q, const double& p2, const int& ip2, double& fxq) {
int nset1 = 1;
evolvepdfpm_(nset1, x, q, p2, ip2, fxq);
}
// alpha_s evolution
/// Get the alpha_s order for the set
void getorderasm_(const int& nset, int& oas) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
// Set equal to the number of members for the requested set
oas = ACTIVESETS[nset].activemember()->info().get_entry_as<int>("AlphaS_OrderQCD");
// Update current set focus
CURRENTSET = nset;
}
/// Get the alpha_s order for the set (non-multiset version)
void getorderas_(int& oas) {
int nset1 = 1;
getorderasm_(nset1, oas);
}
/// Get the alpha_s(Q) value for set nset
double alphaspdfm_(const int& nset, const double& Q){
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
return ACTIVESETS[nset].activemember()->alphasQ(Q);
// Update current set focus
CURRENTSET = nset;
}
/// Get the alpha_s(Q) value for the set (non-multiset version)
double alphaspdf_(const double& Q){
int nset1 = 1;
return alphaspdfm_(nset1, Q);
}
// Metadata functions
/// Get the number of error members in the set (with special treatment for single member sets)
void numberpdfm_(const int& nset, int& numpdf) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
// Set equal to the number of members for the requested set
numpdf= ACTIVESETS[nset].activemember()->info().get_entry_as<int>("NumMembers");
// Reproduce old LHAPDF v5 behaviour, i.e. subtract 1 if more than 1 member set
if (numpdf > 1) numpdf -= 1;
// Update current set focus
CURRENTSET = nset;
}
/// Get the number of error members in the set (non-multiset version)
void numberpdf_(int& numpdf) {
int nset1 = 1;
numberpdfm_(nset1, numpdf);
}
/// Get the max number of active flavours
void getnfm_(const int& nset, int& nf) {
//nf = ACTIVESETS[nset].activemember()->info().get_entry_as<int>("AlphaS_NumFlavors");
nf = ACTIVESETS[nset].activemember()->info().get_entry_as<int>("NumFlavors");
// Update current set focus
CURRENTSET = nset;
}
/// Get the max number of active flavours (non-multiset version)
void getnf_(int& nf) {
int nset1 = 1;
getnfm_(nset1, nf);
}
/// Get nf'th quark mass
void getqmassm_(const int& nset, const int& nf, double& mass) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
if (nf*nf == 1) mass = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("MDown");
else if (nf*nf == 4) mass = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("MUp");
else if (nf*nf == 9) mass = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("MStrange");
else if (nf*nf == 16) mass = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("MCharm");
else if (nf*nf == 25) mass = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("MBottom");
else if (nf*nf == 36) mass = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("MTop");
else throw LHAPDF::UserError("Trying to get quark mass for invalid quark ID #" + LHAPDF::to_str(nf));
// Update current set focus
CURRENTSET = nset;
}
/// Get nf'th quark mass (non-multiset version)
void getqmass_(const int& nf, double& mass) {
int nset1 = 1;
getqmassm_(nset1, nf, mass);
}
/// Get the nf'th quark threshold
void getthresholdm_(const int& nset, const int& nf, double& Q) {
try {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
if (nf*nf == 1) Q = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("ThresholdDown");
else if (nf*nf == 4) Q = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("ThresholdUp");
else if (nf*nf == 9) Q = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("ThresholdStrange");
else if (nf*nf == 16) Q = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("ThresholdCharm");
else if (nf*nf == 25) Q = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("ThresholdBottom");
else if (nf*nf == 36) Q = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("ThresholdTop");
//else throw LHAPDF::UserError("Trying to get quark threshold for invalid quark ID #" + LHAPDF::to_str(nf));
} catch (...) {
getqmassm_(nset, nf, Q);
}
// Update current set focus
CURRENTSET = nset;
}
/// Get the nf'th quark threshold
void getthreshold_(const int& nf, double& Q) {
int nset1 = 1;
getthresholdm_(nset1, nf, Q);
}
/// Print PDF set's description to stdout
void getdescm_(const int& nset) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
cout << ACTIVESETS[nset].activemember()->description() << endl;
// Update current set focus
CURRENTSET = nset;
}
void getdesc_() {
int nset1 = 1;
getdescm_(nset1);
}
void getxminm_(const int& nset, const int& nmem, double& xmin) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
const int activemem = ACTIVESETS[nset].currentmem;
ACTIVESETS[nset].loadMember(nmem);
xmin = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("XMin");
ACTIVESETS[nset].loadMember(activemem);
// Update current set focus
CURRENTSET = nset;
}
void getxmin_(const int& nmem, double& xmin) {
int nset1 = 1;
getxminm_(nset1, nmem, xmin);
}
void getxmaxm_(const int& nset, const int& nmem, double& xmax) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
const int activemem = ACTIVESETS[nset].currentmem;
ACTIVESETS[nset].loadMember(nmem);
xmax = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("XMax");
ACTIVESETS[nset].loadMember(activemem);
// Update current set focus
CURRENTSET = nset;
}
void getxmax_(const int& nmem, double& xmax) {
int nset1 = 1;
getxmaxm_(nset1, nmem, xmax);
}
void getq2minm_(const int& nset, const int& nmem, double& q2min) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
const int activemem = ACTIVESETS[nset].currentmem;
ACTIVESETS[nset].loadMember(nmem);
q2min = LHAPDF::sqr(ACTIVESETS[nset].activemember()->info().get_entry_as<double>("QMin"));
ACTIVESETS[nset].loadMember(activemem);
// Update current set focus
CURRENTSET = nset;
}
void getq2min_(const int& nmem, double& q2min) {
int nset1 = 1;
getq2minm_(nset1, nmem, q2min);
}
void getq2maxm_(const int& nset, const int& nmem, double& q2max) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
const int activemem = ACTIVESETS[nset].currentmem;
ACTIVESETS[nset].loadMember(nmem);
q2max = LHAPDF::sqr(ACTIVESETS[nset].activemember()->info().get_entry_as<double>("QMax"));
ACTIVESETS[nset].loadMember(activemem);
// Update current set focus
CURRENTSET = nset;
}
void getq2max_(const int& nmem, double& q2max) {
int nset1 = 1;
getq2maxm_(nset1, nmem, q2max);
}
void getminmaxm_(const int& nset, const int& nmem, double& xmin, double& xmax, double& q2min, double& q2max) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
const int activemem = ACTIVESETS[nset].currentmem;
ACTIVESETS[nset].loadMember(nmem);
xmin = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("XMin");
xmax = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("XMax");
q2min = LHAPDF::sqr(ACTIVESETS[nset].activemember()->info().get_entry_as<double>("QMin"));
q2max = LHAPDF::sqr(ACTIVESETS[nset].activemember()->info().get_entry_as<double>("QMax"));
ACTIVESETS[nset].loadMember(activemem);
// Update current set focus
CURRENTSET = nset;
}
void getminmax_(const int& nmem, double& xmin, double& xmax, double& q2min, double& q2max) {
int nset1 = 1;
getminmaxm_(nset1, nmem, xmin, xmax, q2min, q2max);
}
void getlam4m_(const int& nset, const int& nmem, double& qcdl4) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
CURRENTSET = nset;
ACTIVESETS[nset].loadMember(nmem);
qcdl4 = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("AlphaS_Lambda4", -1.0);
}
void getlam4_(const int& nmem, double& qcdl4) {
int nset1 = 1;
getlam4m_(nset1, nmem, qcdl4);
}
void getlam5m_(const int& nset, const int& nmem, double& qcdl5) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
CURRENTSET = nset;
ACTIVESETS[nset].loadMember(nmem);
qcdl5 = ACTIVESETS[nset].activemember()->info().get_entry_as<double>("AlphaS_Lambda5", -1.0);
}
void getlam5_(const int& nmem, double& qcdl5) {
int nset1 = 1;
getlam5m_(nset1, nmem, qcdl5);
}
/// Backwards compatibility functions for LHAPDF5 calculations of
/// PDF uncertainties and PDF correlations (G. Watt, March 2014).
// subroutine GetPDFUncTypeM(nset,lMonteCarlo,lSymmetric)
void getpdfunctypem_(const int& nset, int& lmontecarlo, int& lsymmetric) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
const string errorType = ACTIVESETS[nset].activemember()->set().errorType();
if (LHAPDF::startswith(errorType, "replicas")) { // Monte Carlo PDF sets
lmontecarlo = 1;
lsymmetric = 1;
} else if (LHAPDF::startswith(errorType, "symmhessian")) { // symmetric eigenvector PDF sets
lmontecarlo = 0;
lsymmetric = 1;
} else { // default: assume asymmetric Hessian eigenvector PDF sets
lmontecarlo = 0;
lsymmetric = 0;
}
// Update current set focus
CURRENTSET = nset;
}
// subroutine GetPDFUncType(lMonteCarlo,lSymmetric)
void getpdfunctype_(int& lmontecarlo, int& lsymmetric) {
int nset1 = 1;
getpdfunctypem_(nset1, lmontecarlo, lsymmetric);
}
// subroutine GetPDFuncertaintyM(nset,values,central,errplus,errminus,errsym)
void getpdfuncertaintym_(const int& nset, const double* values, double& central, double& errplus, double& errminus, double& errsymm) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
const size_t nmem = ACTIVESETS[nset].activemember()->set().size()-1;
const vector<double> vecvalues(values, values + nmem + 1);
LHAPDF::PDFUncertainty err = ACTIVESETS[nset].activemember()->set().uncertainty(vecvalues, -1);
central = err.central;
// For a combined set, the PDF and parameter variation uncertainties will be added in quadrature.
errplus = err.errplus;
errminus = err.errminus;
errsymm = err.errsymm;
// Update current set focus
CURRENTSET = nset;
}
// subroutine GetPDFuncertainty(values,central,errplus,errminus,errsym)
void getpdfuncertainty_(const double* values, double& central, double& errplus, double& errminus, double& errsymm) {
int nset1 = 1;
getpdfuncertaintym_(nset1, values, central, errplus, errminus, errsymm);
}
// subroutine GetPDFcorrelationM(nset,valuesA,valuesB,correlation)
void getpdfcorrelationm_(const int& nset, const double* valuesA, const double* valuesB, double& correlation) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
const size_t nmem = ACTIVESETS[nset].activemember()->set().size()-1;
const vector<double> vecvaluesA(valuesA, valuesA + nmem + 1);
const vector<double> vecvaluesB(valuesB, valuesB + nmem + 1);
correlation = ACTIVESETS[nset].activemember()->set().correlation(vecvaluesA,vecvaluesB);
// Update current set focus
CURRENTSET = nset;
}
// subroutine GetPDFcorrelation(valuesA,valuesB,correlation)
void getpdfcorrelation_(const double* valuesA, const double* valuesB, double& correlation) {
int nset1 = 1;
getpdfcorrelationm_(nset1, valuesA, valuesB, correlation);
}
///////////////////////////////////////
/// REALLY OLD PDFLIB COMPATIBILITY FUNCTIONS
/// PDFLIB initialisation function
void pdfset_(const char* par, const double* value, int parlength) {
// Identify the calling program (yuck!)
string my_par(par);
if (my_par.find("NPTYPE") != string::npos) {
cout << "==== LHAPDF6 USING PYTHIA-TYPE LHAGLUE INTERFACE ====" << endl;
// Take PDF ID from value[2]
ACTIVESETS[1] = PDFSetHandler(value[2]+1000*value[1]);
} else if (my_par.find("HWLHAPDF") != string::npos) {
cout << "==== LHAPDF6 USING HERWIG-TYPE LHAGLUE INTERFACE ====" << endl;
// Take PDF ID from value[0]
ACTIVESETS[1] = PDFSetHandler(value[0]);
} else if (my_par.find("DEFAULT") != string::npos) {
cout << "==== LHAPDF6 USING DEFAULT-TYPE LHAGLUE INTERFACE ====" << endl;
// Take PDF ID from value[0]
ACTIVESETS[1] = PDFSetHandler(value[0]);
} else {
cout << "==== LHAPDF6 USING PDFLIB-TYPE LHAGLUE INTERFACE ====" << endl;
// Take PDF ID from value[2]
ACTIVESETS[1] = PDFSetHandler(value[2]+1000*value[1]);
}
CURRENTSET = 1;
// Extract parameters for common blocks (with sensible fallback values)
PDFPtr pdf = ACTIVESETS[1].activemember();
w50513_.xmin = pdf->info().get_entry_as<double>("XMin", 0.0);
w50513_.xmax = pdf->info().get_entry_as<double>("XMax", 1.0);
w50513_.q2min = LHAPDF::sqr(pdf->info().get_entry_as<double>("QMin", 1.0));
w50513_.q2max = LHAPDF::sqr(pdf->info().get_entry_as<double>("QMax", 1.0e5));
w50512_.qcdl4 = pdf->info().get_entry_as<double>("AlphaS_Lambda4", 0.0);
w50512_.qcdl5 = pdf->info().get_entry_as<double>("AlphaS_Lambda5", 0.0);
lhapdfr_.qcdlha4 = pdf->info().get_entry_as<double>("AlphaS_Lambda4", 0.0);
lhapdfr_.qcdlha5 = pdf->info().get_entry_as<double>("AlphaS_Lambda5", 0.0);
lhapdfr_.nfllha = 4;
// Activate legacy/compatibility LHAPDF5-type behaviour re. broken Lambda values
if (pdf->info().get_entry_as<bool>("Pythia6LambdaV5Compat", true)) {
w50512_.qcdl4 = 0.192;
w50512_.qcdl5 = 0.192;
lhapdfr_.qcdlha4 = 0.192;
lhapdfr_.qcdlha5 = 0.192;
}
}
/// PDFLIB nucleon structure function querying
void structm_(const double& x, const double& q,
double& upv, double& dnv, double& usea, double& dsea,
double& str, double& chm, double& bot, double& top, double& glu) {
CURRENTSET = 1;
/// Fill (partial) parton return variables
PDFPtr pdf = ACTIVESETS[1].activemember();
dsea = pdf->xfxQ(-1, x, q);
usea = pdf->xfxQ(-2, x, q);
dnv = pdf->xfxQ(1, x, q) - dsea;
upv = pdf->xfxQ(2, x, q) - usea;
str = pdf->xfxQ(3, x, q);
chm = (pdf->hasFlavor(4)) ? pdf->xfxQ(4, x, q) : 0;
bot = (pdf->hasFlavor(5)) ? pdf->xfxQ(5, x, q) : 0;
top = (pdf->hasFlavor(6)) ? pdf->xfxQ(6, x, q) : 0;
glu = pdf->xfxQ(21, x, q);
}
/// PDFLIB photon structure function querying
void structp_(const double& x, const double& q2, const double& p2, const double& ip2,
double& upv, double& dnv, double& usea, double& dsea,
double& str, double& chm, double& bot, double& top, double& glu) {
throw LHAPDF::NotImplementedError("Photon structure functions are not yet supported");
}
/// PDFLIB statistics on PDF under/overflows
void pdfsta_() {
/// @note Can't do anything...
}
}
// LHAPDF namespace C++ compatibility code
#ifdef ENABLE_LHAGLUE_CXX
void LHAPDF::setVerbosity(LHAPDF::Verbosity noiselevel) {
LHAPDF::setVerbosity((int) noiselevel);
}
void LHAPDF::setPDFPath(const string& path) {
pathsPrepend(path);
}
string LHAPDF::pdfsetsPath() {
return paths()[0];
}
int LHAPDF::numberPDF() {
int nmem;
numberpdf_(nmem);
return nmem;
}
int LHAPDF::numberPDF(int nset) {
int nmem;
numberpdfm_(nset,nmem);
return nmem;
}
void LHAPDF::initPDF( int memset) {
int nset1 = 1;
initpdfm_(nset1, memset);
}
void LHAPDF::initPDF(int nset, int memset) {
initpdfm_(nset, memset);
}
double LHAPDF::xfx(double x, double Q, int fl) {
vector<double> r(13);
evolvepdf_(x, Q, &r[0]);
return r[fl+6];
}
double LHAPDF::xfx(int nset, double x, double Q, int fl) {
vector<double> r(13);
evolvepdfm_(nset, x, Q, &r[0]);
return r[fl+6];
}
vector<double> LHAPDF::xfx(double x, double Q) {
vector<double> r(13);
evolvepdf_(x, Q, &r[0]);
return r;
}
vector<double> LHAPDF::xfx(int nset, double x, double Q) {
vector<double> r(13);
evolvepdfm_(nset, x, Q, &r[0]);
return r;
}
void LHAPDF::xfx(double x, double Q, double* results) {
evolvepdf_(x, Q, results);
}
void LHAPDF::xfx(int nset, double x, double Q, double* results) {
evolvepdfm_(nset, x, Q, results);
}
vector<double> LHAPDF::xfxphoton(double x, double Q) {
vector<double> r(13);
double mphoton;
evolvepdfphoton_(x, Q, &r[0], mphoton);
r.push_back(mphoton);
return r;
}
vector<double> LHAPDF::xfxphoton(int nset, double x, double Q) {
vector<double> r(13);
double mphoton;
evolvepdfphotonm_(nset, x, Q, &r[0], mphoton);
r.push_back(mphoton);
return r;
}
void LHAPDF::xfxphoton(double x, double Q, double* results) {
evolvepdfphoton_(x, Q, results, results[13]);
}
void LHAPDF::xfxphoton(int nset, double x, double Q, double* results) {
evolvepdfphotonm_(nset, x, Q, results, results[13]);
}
double LHAPDF::xfxphoton(double x, double Q, int fl) {
vector<double> r(13);
double mphoton;
evolvepdfphoton_(x, Q, &r[0], mphoton);
if (fl == 7) return mphoton;
return r[fl+6];
}
double LHAPDF::xfxphoton(int nset, double x, double Q, int fl) {
vector<double> r(13);
double mphoton;
evolvepdfphotonm_(nset, x, Q, &r[0], mphoton);
if ( fl == 7 ) return mphoton;
return r[fl+6];
}
void LHAPDF::initPDFSet(const string& filename, int nmem) {
initPDFSet(1,filename, nmem);
}
void LHAPDF::initPDFSet(int nset, const string& filename, int nmem) {
initPDFSetByName(nset,filename);
ACTIVESETS[nset].loadMember(nmem);
CURRENTSET = nset;
}
void LHAPDF::initPDFSet(const string& filename, SetType type, int nmem) {
// silently ignore type
initPDFSet(1,filename, nmem);
}
void LHAPDF::initPDFSet(int nset, const string& filename, SetType type, int nmem) {
// silently ignore type
initPDFSetByName(nset,filename);
ACTIVESETS[nset].loadMember(nmem);
CURRENTSET = nset;
}
void LHAPDF::initPDFSet(int nset, int setid, int nmem) {
ACTIVESETS[nset] = PDFSetHandler(setid); //
CURRENTSET = nset;
}
void LHAPDF::initPDFSet(int setid, int nmem) {
initPDFSet(1,setid,nmem);
}
#define SIZE 999
void LHAPDF::initPDFSetByName(const string& filename) {
std::cout << "initPDFSetByName: " << filename << std::endl;
char cfilename[SIZE+1];
strncpy(cfilename, filename.c_str(), SIZE);
initpdfsetbyname_(cfilename, filename.length());
}
void LHAPDF::initPDFSetByName(int nset, const string& filename) {
char cfilename[SIZE+1];
strncpy(cfilename, filename.c_str(), SIZE);
initpdfsetbynamem_(nset, cfilename, filename.length());
}
void LHAPDF::initPDFSetByName(const string& filename, SetType type) {
//silently ignore type
std::cout << "initPDFSetByName: " << filename << std::endl;
char cfilename[SIZE+1];
strncpy(cfilename, filename.c_str(), SIZE);
initpdfsetbyname_(cfilename, filename.length());
}
void LHAPDF::initPDFSetByName(int nset, const string& filename, SetType type) {
//silently ignore type
char cfilename[SIZE+1];
strncpy(cfilename, filename.c_str(), SIZE);
initpdfsetbynamem_(nset, cfilename, filename.length());
}
void LHAPDF::getDescription() {
getDescription(1);
}
void LHAPDF::getDescription(int nset) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
cout << ACTIVESETS[nset].activemember()->set().description() << endl;
}
double LHAPDF::alphasPDF(double Q) {
return LHAPDF::alphasPDF(1, Q) ;
}
double LHAPDF::alphasPDF(int nset, double Q) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
CURRENTSET = nset;
// return alphaS for the requested set
return ACTIVESETS[nset].activemember()->alphasQ(Q);
}
bool LHAPDF::hasPhoton(){
return has_photon_();
}
int LHAPDF::getOrderAlphaS() {
return LHAPDF::getOrderAlphaS(1) ;
}
int LHAPDF::getOrderAlphaS(int nset) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
CURRENTSET = nset;
// return alphaS Order for the requested set
return ACTIVESETS[nset].activemember()->info().get_entry_as<int>("AlphaS_OrderQCD", -1);
}
int LHAPDF::getOrderPDF() {
return LHAPDF::getOrderPDF(1) ;
}
int LHAPDF::getOrderPDF(int nset) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
CURRENTSET = nset;
// return PDF order for the requested set
return ACTIVESETS[nset].activemember()->info().get_entry_as<int>("OrderQCD", -1);
}
double LHAPDF::getLam4(int nmem) {
return LHAPDF::getLam4(1, nmem) ;
}
double LHAPDF::getLam4(int nset, int nmem) {
// if (ACTIVESETS.find(nset) == ACTIVESETS.end())
// throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
// CURRENTSET = nset;
// ACTIVESETS[nset].loadMember(nmem);
// return ACTIVESETS[nset].activemember()->info().get_entry_as<double>("AlphaS_Lambda4", -1.0);
double qcdl4;
getlam4m_(nset, nmem, qcdl4);
return qcdl4;
}
double LHAPDF::getLam5(int nmem) {
return LHAPDF::getLam5(1, nmem) ;
}
double LHAPDF::getLam5(int nset, int nmem) {
// if (ACTIVESETS.find(nset) == ACTIVESETS.end())
// throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
// CURRENTSET = nset;
// ACTIVESETS[nset].loadMember(nmem);
// return ACTIVESETS[nset].activemember()->info().get_entry_as<double>("AlphaS_Lambda5", -1.0);
double qcdl5;
getlam5m_(nset, nmem, qcdl5);
return qcdl5;
}
int LHAPDF::getNf() {
return LHAPDF::getNf(1) ;
}
int LHAPDF::getNf(int nset) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
CURRENTSET = nset;
// return alphaS Order for the requested set
return ACTIVESETS[nset].activemember()->info().get_entry_as<int>("NumFlavors");
}
double LHAPDF::getXmin(int nmem) {
return LHAPDF::getXmin(1, nmem) ;
}
double LHAPDF::getXmin(int nset, int nmem) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
CURRENTSET = nset;
// return alphaS Order for the requested set
ACTIVESETS[nset].loadMember(nmem);
return ACTIVESETS[nset].activemember()->info().get_entry_as<double>("XMin");
}
double LHAPDF::getXmax(int nmem) {
return LHAPDF::getXmax(1, nmem) ;
}
double LHAPDF::getXmax(int nset, int nmem) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
CURRENTSET = nset;
// return alphaS Order for the requested set
ACTIVESETS[nset].loadMember(nmem);
return ACTIVESETS[nset].activemember()->info().get_entry_as<double>("XMax");
}
double LHAPDF::getQ2min(int nmem) {
return LHAPDF::getQ2min(1, nmem) ;
}
double LHAPDF::getQ2min(int nset, int nmem) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
CURRENTSET = nset;
// return alphaS Order for the requested set
ACTIVESETS[nset].loadMember(nmem);
return pow(ACTIVESETS[nset].activemember()->info().get_entry_as<double>("QMin"),2);
}
double LHAPDF::getQ2max(int nmem) {
return LHAPDF::getQ2max(1,nmem) ;
}
double LHAPDF::getQ2max(int nset, int nmem) {
if (ACTIVESETS.find(nset) == ACTIVESETS.end())
throw LHAPDF::UserError("Trying to use LHAGLUE set #" + LHAPDF::to_str(nset) + " but it is not initialised");
CURRENTSET = nset;
// return alphaS Order for the requested set
ACTIVESETS[nset].loadMember(nmem);
return pow(ACTIVESETS[nset].activemember()->info().get_entry_as<double>("QMax"),2);
}
double LHAPDF::getQMass(int nf) {
return LHAPDF::getQMass(1, nf) ;
}
double LHAPDF::getQMass(int nset, int nf) {
double mass;
getqmassm_(nset, nf, mass);
return mass;
}
double LHAPDF::getThreshold(int nf) {
return LHAPDF::getThreshold(1, nf) ;
}
double LHAPDF::getThreshold(int nset, int nf) {
double thres;
getthresholdm_(nset, nf, thres);
return thres;
}
void LHAPDF::usePDFMember(int member) {
initpdf_(member);
}
void LHAPDF::usePDFMember(int nset, int member) {
initpdfm_(nset, member);
}
#endif // ENABLE_LHAGLUE_CXX
diff --git a/src/LogBicubicInterpolator.cc b/src/LogBicubicInterpolator.cc
--- a/src/LogBicubicInterpolator.cc
+++ b/src/LogBicubicInterpolator.cc
@@ -1,146 +1,146 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/LogBicubicInterpolator.h"
#include <iostream>
namespace LHAPDF {
namespace { // Unnamed namespace
/// One-dimensional linear interpolation for y(x)
inline double _interpolateLinear(double x, double xl, double xh, double yl, double yh) {
assert(x >= xl);
assert(xh >= x);
return yl + (x - xl) / (xh - xl) * (yh - yl);
}
/// One-dimensional cubic interpolation
inline double _interpolateCubic(double T, double VL, double VDL, double VH, double VDH) {
// Pre-calculate powers of T
const double t2 = T*T;
const double t3 = t2*T;
// Calculate left point
const double p0 = (2*t3 - 3*t2 + 1)*VL;
const double m0 = (t3 - 2*t2 + T)*VDL;
// Calculate right point
const double p1 = (-2*t3 + 3*t2)*VH;
const double m1 = (t3 - t2)*VDH;
return p0 + m0 + p1 + m1;
}
/// Calculate adjacent d(xf)/dx at all grid locations for fixed iq2
double _dxf_dlogx(const KnotArray1F& subgrid, size_t ix, size_t iq2) {
const size_t nxknots = subgrid.xs().size();
/// @todo Store pre-cached dlogxs, dlogq2s on subgrids, to replace these denominators? Any real speed gain for the extra memory?
if (ix != 0 && ix != nxknots-1) { //< If central, use the central difference
/// @note We evaluate the most likely condition first to help compiler branch prediction
const double lddx = (subgrid.xf(ix, iq2) - subgrid.xf(ix-1, iq2)) / (subgrid.logxs()[ix] - subgrid.logxs()[ix-1]);
const double rddx = (subgrid.xf(ix+1, iq2) - subgrid.xf(ix, iq2)) / (subgrid.logxs()[ix+1] - subgrid.logxs()[ix]);
return (lddx + rddx) / 2.0;
} else if (ix == 0) { //< If at leftmost edge, use forward difference
return (subgrid.xf(ix+1, iq2) - subgrid.xf(ix, iq2)) / (subgrid.logxs()[ix+1] - subgrid.logxs()[ix]);
} else if (ix == nxknots-1) { //< If at rightmost edge, use backward difference
return (subgrid.xf(ix, iq2) - subgrid.xf(ix-1, iq2)) / (subgrid.logxs()[ix] - subgrid.logxs()[ix-1]);
} else {
throw LogicError("We shouldn't be able to get here!");
}
}
}
double LogBicubicInterpolator::_interpolateXQ2(const KnotArray1F& subgrid, double x, size_t ix, double q2, size_t iq2) const {
// Raise an error if there are too few knots even for a linear fall-back
const size_t nxknots = subgrid.logxs().size();
const size_t nq2knots = subgrid.logq2s().size();
if (nxknots < 4)
throw GridError("PDF subgrids are required to have at least 4 x-knots for use with LogBicubicInterpolator");
if (nq2knots < 2)
throw GridError("PDF subgrids are required to have at least 2 Q-knots for use with LogBicubicInterpolator");
// Check x and q index ranges -- we always need i and i+1 indices to be valid
const size_t ixmax = nxknots - 1;
const size_t iq2max = nq2knots - 1;
if (ix+1 > ixmax) // also true if ix is off the end
throw GridError("Attempting to access an x-knot index past the end of the array, in linear fallback mode");
if (iq2+1 > iq2max) // also true if iq2 is off the end
throw GridError("Attempting to access an Q-knot index past the end of the array, in linear fallback mode");
const double logx = log(x);
const double logq2 = log(q2);
// Fall back to LogBilinearInterpolator if either 2 or 3 Q-knots
if (nq2knots < 4) {
// First interpolate in x
const double logx0 = subgrid.logxs()[ix];
const double logx1 = subgrid.logxs()[ix+1];
const double f_ql = _interpolateLinear(logx, logx0, logx1, subgrid.xf(ix, iq2), subgrid.xf(ix+1, iq2));
const double f_qh = _interpolateLinear(logx, logx0, logx1, subgrid.xf(ix, iq2+1), subgrid.xf(ix+1, iq2+1));
// Then interpolate in Q2, using the x-ipol results as anchor points
return _interpolateLinear(logq2, subgrid.logq2s()[iq2], subgrid.logq2s()[iq2+1], f_ql, f_qh);
}
// else proceed with cubic interpolation:
// Pre-calculate parameters
/// @todo Cache these between calls, re-using if x == x_prev and Q2 == Q2_prev
const double dlogx_1 = subgrid.logxs()[ix+1] - subgrid.logxs()[ix];
const double tlogx = (logx - subgrid.logxs()[ix]) / dlogx_1;
const double dlogq_0 = (iq2 != 0) ? subgrid.logq2s()[iq2] - subgrid.logq2s()[iq2-1] : -1; //< Don't evaluate (or use) if iq2-1 < 0
const double dlogq_1 = subgrid.logq2s()[iq2+1] - subgrid.logq2s()[iq2];
const double dlogq_2 = (iq2+1 != iq2max) ? subgrid.logq2s()[iq2+2] - subgrid.logq2s()[iq2+1] : -1; //< Don't evaluate (or use) if iq2+2 > iq2max
const double tlogq = (logq2 - subgrid.logq2s()[iq2]) / dlogq_1;
/// @todo Statically pre-compute the whole nx * nq gradiant array? I.e. _dxf_dlogx for all points in all subgrids. Memory ~doubling :-/ Could cache them as they are used...
// Points in Q2
double vl = _interpolateCubic(tlogx, subgrid.xf(ix, iq2), _dxf_dlogx(subgrid, ix, iq2) * dlogx_1,
subgrid.xf(ix+1, iq2), _dxf_dlogx(subgrid, ix+1, iq2) * dlogx_1);
double vh = _interpolateCubic(tlogx, subgrid.xf(ix, iq2+1), _dxf_dlogx(subgrid, ix, iq2+1) * dlogx_1,
subgrid.xf(ix+1, iq2+1), _dxf_dlogx(subgrid, ix+1, iq2+1) * dlogx_1);
// Derivatives in Q2
double vdl, vdh;
if (iq2 > 0 && iq2+1 < iq2max) {
// Central difference for both q
/// @note We evaluate the most likely condition first to help compiler branch prediction
double vll = _interpolateCubic(tlogx, subgrid.xf(ix, iq2-1), _dxf_dlogx(subgrid, ix, iq2-1) * dlogx_1,
subgrid.xf(ix+1, iq2-1), _dxf_dlogx(subgrid, ix+1, iq2-1) * dlogx_1);
vdl = ( (vh - vl)/dlogq_1 + (vl - vll)/dlogq_0 ) / 2.0;
double vhh = _interpolateCubic(tlogx, subgrid.xf(ix, iq2+2), _dxf_dlogx(subgrid, ix, iq2+2) * dlogx_1,
subgrid.xf(ix+1, iq2+2), _dxf_dlogx(subgrid, ix+1, iq2+2) * dlogx_1);
vdh = ( (vh - vl)/dlogq_1 + (vhh - vh)/dlogq_2 ) / 2.0;
}
else if (iq2 == 0) {
// Forward difference for lower q
vdl = (vh - vl) / dlogq_1;
// Central difference for higher q
double vhh = _interpolateCubic(tlogx, subgrid.xf(ix, iq2+2), _dxf_dlogx(subgrid, ix, iq2+2) * dlogx_1,
subgrid.xf(ix+1, iq2+2), _dxf_dlogx(subgrid, ix+1, iq2+2) * dlogx_1);
vdh = (vdl + (vhh - vh)/dlogq_2) / 2.0;
}
else if (iq2+1 == iq2max) {
// Backward difference for higher q
vdh = (vh - vl) / dlogq_1;
// Central difference for lower q
double vll = _interpolateCubic(tlogx, subgrid.xf(ix, iq2-1), _dxf_dlogx(subgrid, ix, iq2-1) * dlogx_1,
subgrid.xf(ix+1, iq2-1), _dxf_dlogx(subgrid, ix+1, iq2-1) * dlogx_1);
vdl = (vdh + (vl - vll)/dlogq_0) / 2.0;
}
else throw LogicError("We shouldn't be able to get here!");
vdl *= dlogq_1;
vdh *= dlogq_1;
return _interpolateCubic(tlogq, vl, vdl, vh, vdh);
}
}
diff --git a/src/LogBilinearInterpolator.cc b/src/LogBilinearInterpolator.cc
--- a/src/LogBilinearInterpolator.cc
+++ b/src/LogBilinearInterpolator.cc
@@ -1,39 +1,39 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/LogBilinearInterpolator.h"
namespace LHAPDF {
namespace { // Unnamed namespace
// One-dimensional linear interpolation for y(x)
inline double _interpolateLinear(double x, double xl, double xh, double yl, double yh) {
assert(x >= xl);
assert(xh >= x);
return yl + (x - xl) / (xh - xl) * (yh - yl);
}
}
double LogBilinearInterpolator::_interpolateXQ2(const KnotArray1F& subgrid, double x, size_t ix, double q2, size_t iq2) const {
if (subgrid.logxs().size() < 2)
throw GridError("PDF subgrids are required to have at least 2 x-knots for use with LogBilinearInterpolator");
if (subgrid.logq2s().size() < 2)
throw GridError("PDF subgrids are required to have at least 2 Q2-knots for use with LogBilinearInterpolator");
// First interpolate in x
const double logx = log(x);
const double logx0 = subgrid.logxs()[ix];
const double logx1 = subgrid.logxs()[ix+1];
const double f_ql = _interpolateLinear(logx, logx0, logx1, subgrid.xf(ix, iq2), subgrid.xf(ix+1, iq2));
const double f_qh = _interpolateLinear(logx, logx0, logx1, subgrid.xf(ix, iq2+1), subgrid.xf(ix+1, iq2+1));
// Then interpolate in Q2, using the x-ipol results as anchor points
return _interpolateLinear(log(q2), subgrid.logq2s()[iq2], subgrid.logq2s()[iq2+1], f_ql, f_qh);
}
}
diff --git a/src/NearestPointExtrapolator.cc b/src/NearestPointExtrapolator.cc
--- a/src/NearestPointExtrapolator.cc
+++ b/src/NearestPointExtrapolator.cc
@@ -1,39 +1,39 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/NearestPointExtrapolator.h"
#include "LHAPDF/GridPDF.h"
namespace LHAPDF {
namespace { // Unnamed namespace
// Return the value in the given list that best matches the target value
double _findClosestMatch(const vector<double>& cands, double target) {
// cout << "From NPXpol: knots = ["; for (double c : cands) cout << c << " "; cout << endl;
vector<double>::const_iterator it = lower_bound(cands.begin(), cands.end(), target);
const double upper = *it;
const double lower = (it == cands.begin()) ? upper : *(--it); //< Avoid decrementing the first entry
/// @todo Closeness in linear or log space? Hmm...
if (fabs(target - upper) < fabs(target - lower)) return upper;
return lower;
}
}
double NearestPointExtrapolator::extrapolateXQ2(int id, double x, double q2) const {
/// Find the closest valid x and Q2 points, either on- or off-grid, and use the current interpolator
// cout << "From NPXpol: x = " << x << endl;
/// @todo We should *always* interpolate x -> 1.0
const double closestX = (pdf().inRangeX(x)) ? x : _findClosestMatch(pdf().xKnots(), x);
const double closestQ2 = (pdf().inRangeQ2(q2)) ? q2 : _findClosestMatch(pdf().q2Knots(), q2);
// cout << "From NPXpol: x_closest = " << closestX << ", Q2_closest = " << closestQ2 << endl;;
return pdf().interpolator().interpolateXQ2(id, closestX, closestQ2);
}
}
diff --git a/src/PDF.cc b/src/PDF.cc
--- a/src/PDF.cc
+++ b/src/PDF.cc
@@ -1,62 +1,62 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/PDF.h"
#include "LHAPDF/PDFSet.h"
using namespace std;
namespace LHAPDF {
void PDF::print(std::ostream& os, int verbosity) const {
stringstream ss;
if (verbosity > 0) {
ss << set().name() << " PDF set, member #" << memberID()
<< ", version " << dataversion();
if (lhapdfID() > 0)
ss << "; LHAPDF ID = " << lhapdfID();
}
if (verbosity > 2 && set().description().size() > 0)
ss << "\n" << set().description();
if (verbosity > 1 && description().size() > 0)
ss << "\n" << description();
if (verbosity > 2)
ss << "\n" << "Flavor content = " << to_str(flavors());
os << ss.str() << endl;
}
int PDF::lhapdfID() const {
//return set().lhapdfID() + memberID()
/// @todo Add failure tolerance if pdfsets.index not found
try {
return lookupLHAPDFID(_setname(), memberID());
} catch (const Exception&) {
return -1; //< failure
}
}
double PDF::quarkMass(int id) const {
const unsigned int aid = std::abs(id);
if (aid == 0 || aid > 6) return -1;
const static string QNAMES[] = {"Down", "Up", "Strange", "Charm", "Bottom", "Top"}; ///< @todo Centralise?
const size_t qid = aid - 1;
const string qname = QNAMES[qid];
return info().get_entry_as<double>("M" + qname, -1);
}
double PDF::quarkThreshold(int id) const {
const unsigned int aid = std::abs(id);
if (aid == 0 || aid > 6) return -1;
const static string QNAMES[] = {"Down", "Up", "Strange", "Charm", "Bottom", "Top"}; ///< @todo Centralise?
const size_t qid = aid - 1;
const string qname = QNAMES[qid];
return info().get_entry_as<double>("Threshold" + qname, quarkMass(id));
}
}
diff --git a/src/PDFInfo.cc b/src/PDFInfo.cc
--- a/src/PDFInfo.cc
+++ b/src/PDFInfo.cc
@@ -1,65 +1,65 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/PDFInfo.h"
#include "LHAPDF/PDFSet.h"
#include "LHAPDF/Factories.h"
namespace LHAPDF {
/// Constructor from a path to a member data file.
PDFInfo::PDFInfo(const std::string& mempath) {
if (mempath.empty())
throw UserError("Empty/invalid data path given to PDFInfo constructor");
load(mempath);
// Extract the set name and member ID from the filename.
_setname = basename(dirname(mempath));
const string memname = file_stem(mempath);
assert(memname.length() > 5); // There must be more to the filename stem than just the _nnnn suffix
_member = lexical_cast<int>(memname.substr(memname.length()-4)); //< Last 4 chars should be the member number
}
/// Constructor from a set name and member ID.
PDFInfo::PDFInfo(const std::string& setname, int member) {
_setname = setname;
_member = member;
const string searchpath = findFile(pdfmempath(setname, member));
if (searchpath.empty())
throw ReadError("Couldn't find a PDF data file for " + setname + " #" + to_str(member));
load(searchpath);
}
/// Constructor from an LHAPDF ID code.
PDFInfo::PDFInfo(int lhaid) {
const pair<string,int> setname_memid = lookupPDF(lhaid);
if (setname_memid.second == -1)
throw IndexError("Can't find a PDF with LHAPDF ID = " + to_str(lhaid));
_setname = setname_memid.first; _member = setname_memid.second;
const string searchpath = pdfmempath(setname_memid.first, setname_memid.second);
if (searchpath.empty())
throw ReadError("Couldn't find a PDF data file for LHAPDF ID = " + to_str(lhaid));
load(searchpath);
}
bool PDFInfo::has_key(const string& key) const {
// cout << key << " in PDF: " << boolalpha << has_key_local(key) << endl;
// cout << key << " in Set: " << boolalpha << getPDFSet(_setname).has_key(key) << endl;
// cout << key << " in Cfg: " << boolalpha << getConfig().has_key(key) << endl;
return has_key_local(key) || getPDFSet(_setname).has_key(key);
}
const std::string& PDFInfo::get_entry(const string& key) const {
if (has_key_local(key)) return get_entry_local(key); //< value is defined locally
return getPDFSet(_setname).get_entry(key); //< fall back to the set-level info... or beyond
}
}
diff --git a/src/PDFSet.cc b/src/PDFSet.cc
--- a/src/PDFSet.cc
+++ b/src/PDFSet.cc
@@ -1,306 +1,306 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/PDFSet.h"
#include <boost/math/distributions/chi_squared.hpp>
namespace LHAPDF {
PDFSet::PDFSet(const string& setname) {
/// @todo Hmm, this relies on the standard search path system ... no specific-path API.
_setname = setname;
const string setinfopath = findpdfsetinfopath(setname);
if (!file_exists(setinfopath))
throw ReadError("Info file not found for PDF set '" + setname + "'");
// Load info file
load(setinfopath);
/// @todo Check that some mandatory metadata keys have been set: _check() function.
}
void PDFSet::print(ostream& os, int verbosity) const {
stringstream ss;
if (verbosity > 0)
ss << name() << ", version " << dataversion() << "; " << size() << " PDF members";
if (verbosity > 1)
ss << "\n" << description();
os << ss.str() << endl;
}
double PDFSet::errorConfLevel() const {
// Return -1 or similar invalid value if errorType is replicas: requires changes in uncertainty code below.
return get_entry_as<double>("ErrorConfLevel", (!startswith(errorType(), "replicas")) ? 100*erf(1/sqrt(2)) : -1);
}
PDFUncertainty PDFSet::uncertainty(const vector<double>& values, double cl, bool alternative) const {
if (values.size() != size())
throw UserError("Error in LHAPDF::PDFSet::uncertainty. Input vector must contain values for all PDF members.");
// PDF members labelled 0 to nmem, excluding possible parameter variations.
size_t nmem = size()-1;
const size_t npar = countchar(errorType(), '+');
nmem -= 2*npar;
if (nmem <= 0)
throw UserError("Error in LHAPDF::PDFSet::uncertainty. PDF set must contain more than just the central value.");
// Get set- and requested conf levels (converted from %) and check sanity (req CL = set CL if cl < 0).
// For replica sets, we internally use a nominal setCL corresponding to 1-sigma, since errorConfLevel() == -1.
const double setCL = (!startswith(errorType(), "replicas")) ? errorConfLevel() / 100.0 : erf(1/sqrt(2));
const double reqCL = (cl >= 0) ? cl / 100.0 : setCL; // convert from percentage
if (!in_range(reqCL, 0, 1) || !in_range(setCL, 0, 1))
throw UserError("Error in LHAPDF::PDFSet::uncertainty. Requested or PDF set confidence level outside [0,1] range.");
// Return value
PDFUncertainty rtn;
rtn.central = values[0];
if (alternative && startswith(errorType(), "replicas")) {
// Compute median and requested CL directly from probability distribution of replicas.
// Sort "values" into increasing order, ignoring zeroth member (average over replicas).
// Also ignore possible parameter variations included at the end of the set.
vector<double> sorted = values;
sort(sorted.begin()+1, sorted.end()-2*npar);
// Define central value to be median.
if (nmem % 2) { // odd nmem => one middle value
rtn.central = sorted[nmem/2 + 1];
} else { // even nmem => average of two middle values
rtn.central = 0.5*(sorted[nmem/2] + sorted[nmem/2 + 1]);
}
// Define uncertainties via quantiles with a CL given by reqCL.
const int upper = round(0.5*(1+reqCL)*nmem); // round to nearest integer
const int lower = 1 + round(0.5*(1-reqCL)*nmem); // round to nearest integer
rtn.errplus = sorted[upper] - rtn.central;
rtn.errminus = rtn.central - sorted[lower];
rtn.errsymm = 0.5*(rtn.errplus + rtn.errminus); // symmetrised
} else if (alternative) {
throw UserError("Error in LHAPDF::PDFSet::uncertainty. This PDF set is not in the format of replicas.");
} else if (startswith(errorType(), "replicas")) {
// Calculate the average and standard deviation using Eqs. (2.3) and (2.4) of arXiv:1106.5788v2.
double av = 0.0, sd = 0.0;
for (size_t imem = 1; imem <= nmem; imem++) {
av += values[imem];
sd += sqr(values[imem]);
}
av /= nmem; sd /= nmem;
sd = nmem/(nmem-1.0)*(sd-sqr(av));
sd = (sd > 0.0 && nmem > 1) ? sqrt(sd) : 0.0;
rtn.central = av;
rtn.errplus = rtn.errminus = rtn.errsymm = sd;
} else if (startswith(errorType(), "symmhessian")) {
double errsymm = 0;
for (size_t ieigen = 1; ieigen <= nmem; ieigen++)
errsymm += sqr(values[ieigen]-values[0]);
errsymm = sqrt(errsymm);
rtn.errplus = rtn.errminus = rtn.errsymm = errsymm;
} else if (startswith(errorType(), "hessian")) {
// Calculate the asymmetric and symmetric Hessian uncertainties
// using Eqs. (2.1), (2.2) and (2.6) of arXiv:1106.5788v2.
double errplus = 0, errminus = 0, errsymm = 0;
for (size_t ieigen = 1; ieigen <= nmem/2; ieigen++) {
errplus += sqr(max(max(values[2*ieigen-1]-values[0],values[2*ieigen]-values[0]), 0.0));
errminus += sqr(max(max(values[0]-values[2*ieigen-1],values[0]-values[2*ieigen]), 0.0));
errsymm += sqr(values[2*ieigen-1]-values[2*ieigen]);
}
rtn.errsymm = 0.5*sqrt(errsymm);
rtn.errplus = sqrt(errplus);
rtn.errminus = sqrt(errminus);
} else {
throw MetadataError("\"ErrorType: " + errorType() + "\" not supported by LHAPDF::PDFSet::uncertainty.");
}
if (setCL != reqCL) {
// Apply scaling to Hessian sets or replica sets with alternative=false.
// Calculate the qth quantile of the chi-squared distribution with one degree of freedom.
// Examples: quantile(dist, q) = {0.988946, 1, 2.70554, 3.84146, 4} for q = {0.68, 1-sigma, 0.90, 0.95, 2-sigma}.
// boost::math::chi_squared dist(1);
// double qsetCL = boost::math::quantile(dist, setCL);
// double qreqCL = boost::math::quantile(dist, reqCL);
double qsetCL = chisquared_quantile(setCL, 1);
double qreqCL = chisquared_quantile(reqCL, 1);
// Scale uncertainties from the original set CL to the requested CL.
const double scale = sqrt(qreqCL/qsetCL);
rtn.scale = scale;
if (!alternative) {
rtn.errplus *= scale;
rtn.errminus *= scale;
rtn.errsymm *= scale;
}
}
rtn.errplus_pdf = rtn.errplus;
rtn.errminus_pdf = rtn.errminus;
rtn.errsymm_pdf = rtn.errsymm;
if (npar > 0) {
// All individual parameter variation uncertainties are added in quadrature.
double err_par = 0;
for (size_t ipar = 1; ipar <= npar; ipar++) {
err_par += sqr(values[nmem+2*ipar-1]-values[nmem+2*ipar]);
}
// Calculate total uncertainty from parameter variation with same scaling as for PDF uncertainty.
rtn.err_par = rtn.scale * 0.5 * sqrt(err_par);
// Add parameter variation uncertainty in quadrature with PDF uncertainty.
rtn.errplus = sqrt( sqr(rtn.errplus_pdf) + sqr(rtn.err_par) );
rtn.errminus = sqrt( sqr(rtn.errminus_pdf) + sqr(rtn.err_par) );
rtn.errsymm = sqrt( sqr(rtn.errsymm_pdf) + sqr(rtn.err_par) );
}
return rtn;
}
double PDFSet::correlation(const vector<double>& valuesA, const vector<double>& valuesB) const {
if (valuesA.size() != size() || valuesB.size() != size())
throw UserError("Error in LHAPDF::PDFSet::correlation. Input vectors must contain values for all PDF members.");
const PDFUncertainty errA = uncertainty(valuesA, -1);
const PDFUncertainty errB = uncertainty(valuesB, -1);
// PDF members labelled 0 to nmem, excluding possible parameter variations.
size_t nmem = size()-1;
const size_t npar = countchar(errorType(), '+');
nmem -= 2*npar;
double cor = 0.0;
if (startswith(errorType(), "replicas") && nmem > 1) {
// Calculate the correlation using Eq. (2.7) of arXiv:1106.5788v2.
for (size_t imem = 1; imem <= nmem; imem++)
cor += valuesA[imem] * valuesB[imem];
cor = (cor/nmem - errA.central*errB.central) / (errA.errsymm_pdf*errB.errsymm_pdf) * nmem/(nmem-1.0);
} else if (startswith(errorType(), "symmhessian")) {
for (size_t ieigen = 1; ieigen <= nmem; ieigen++)
cor += (valuesA[ieigen]-errA.central) * (valuesB[ieigen]-errB.central);
cor /= errA.errsymm_pdf * errB.errsymm_pdf;
} else if (startswith(errorType(), "hessian")) {
// Calculate the correlation using Eq. (2.5) of arXiv:1106.5788v2.
for (size_t ieigen = 1; ieigen <= nmem/2; ieigen++)
cor += (valuesA[2*ieigen-1]-valuesA[2*ieigen]) * (valuesB[2*ieigen-1]-valuesB[2*ieigen]);
cor /= 4.0 * errA.errsymm_pdf * errB.errsymm_pdf;
}
return cor;
}
double PDFSet::randomValueFromHessian(const vector<double>& values, const vector<double>& randoms, bool symmetrise) const {
if (values.size() != size())
throw UserError("Error in LHAPDF::PDFSet::randomValueFromHessian. Input vector must contain values for all PDF members.");
double frand = 0.0;
double scale = uncertainty(values).scale;
// PDF members labelled 0 to nmem, excluding possible parameter variations.
size_t nmem = size()-1;
const size_t npar = countchar(errorType(), '+');
nmem -= 2*npar;
// Allocate number of eigenvectors based on ErrorType.
size_t neigen = 0;
if (startswith(errorType(), "hessian")) {
neigen = nmem/2;
} else if (startswith(errorType(), "symmhessian")) {
neigen = nmem;
} else {
throw UserError("Error in LHAPDF::PDFSet::randomValueFromHessian. This PDF set is not in the Hessian format.");
}
if (randoms.size() != neigen)
throw UserError("Error in LHAPDF::PDFSet::randomValueFromHessian. Input vector must contain random numbers for all eigenvectors.");
frand = values[0];
if (startswith(errorType(), "symmhessian")) {
// Loop over number of eigenvectors.
for (size_t ieigen = 1; ieigen <= neigen; ieigen++) {
double r = randoms[ieigen-1]; // Gaussian random number
frand += r*(values[ieigen]-values[0])*scale;
}
} else if (startswith(errorType(), "hessian")) {
// Use either Eq. (6.4) or corrected Eq. (6.5) of arXiv:1205.4024v2.
// Loop over number of eigenvectors.
for (size_t ieigen = 1; ieigen <= neigen; ieigen++) {
double r = randoms[ieigen-1]; // Gaussian random number
if (symmetrise) {
frand += 0.5*r*(values[2*ieigen-1]-values[2*ieigen]) * scale;
} else { // not symmetrised
if (r < 0.0) frand -= r*(values[2*ieigen]-values[0]) * scale; // negative direction
else frand += r*(values[2*ieigen-1]-values[0]) * scale; // positive direction
}
}
}
return frand;
}
void PDFSet::_checkPdfType(const std::vector<string>& pdftypes) const {
if (pdftypes.size() != size())
throw UserError("Error in LHAPDF::PDFSet::checkPdfType. Input vector must contain values for all PDF members.");
// PDF members labelled 0 to nmem, excluding possible parameter variations.
size_t nmem = size()-1;
const size_t npar = countchar(errorType(), '+');
nmem -= 2*npar;
// Check that zeroth member has "PdfType: central".
if (pdftypes[0] != "central")
throw MetadataError("Member 0, \"PdfType: " + pdftypes[0] + "\" should be \"PdfType: central\".");
// Check that PDF members have "PdfType: replica" or "PdfType: error".
if (startswith(errorType(), "replicas")) {
for (size_t imem = 1; imem <= nmem; imem++) {
if (pdftypes[imem] != "replica")
throw MetadataError("Member " + to_str(imem) + ", \"PdfType: " + pdftypes[imem] + "\" should be \"PdfType: replica\".");
}
} else if (startswith(errorType(), "symmhessian") || startswith(errorType(), "hessian")) {
for (size_t imem = 1; imem <= nmem; imem++) {
if (pdftypes[imem] != "error")
throw MetadataError("Member " + to_str(imem) + ", \"PdfType: " + pdftypes[imem] + "\" should be \"PdfType: error\".");
}
} else {
throw MetadataError("\"ErrorType: " + errorType() + "\" not supported by LHAPDF::PDFSet::checkPdfType.");
}
// Check that possible parameter variations have "PdfType: central".
for (size_t imem = nmem+1; imem <= size()-1; imem++) {
if (pdftypes[imem] != "central")
throw MetadataError("Member " + to_str(imem) + ", \"PdfType: " + pdftypes[imem] + "\" should be \"PdfType: central\".");
}
//cout << "Success: PdfType of each member matches the ErrorType of the set." << endl;
}
}
diff --git a/src/Paths.cc b/src/Paths.cc
--- a/src/Paths.cc
+++ b/src/Paths.cc
@@ -1,80 +1,80 @@
// -*- C++ -*-
//
// This file is part of LHAPDF
-// Copyright (C) 2012-2014 The LHAPDF collaboration (see AUTHORS for details)
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
//
#include "LHAPDF/Paths.h"
#include "LHAPDF/Info.h"
#include "LHAPDF/Config.h"
#include <dirent.h>
namespace LHAPDF {
std::vector<std::string> paths() {
// Use LHAPDF_DATA_PATH for all path storage
char* pathsvar = getenv("LHAPDF_DATA_PATH");
// But fall back to looking in LHAPATH if the preferred var is not defined
if (pathsvar == 0) pathsvar = getenv("LHAPATH");
const string spathsvar = (pathsvar != 0) ? pathsvar : "";
// Split the paths variable as usual
vector<string> rtn = split(spathsvar, ":");
// Look in the install prefix after other paths are exhausted, if not blocked by a trailing ::
/// @todo Move some logic to starts_with and ends_with functions in Utils.h
if (rtn.empty() || spathsvar.length() < 2 || spathsvar.substr(spathsvar.length()-2) != "::") {
const string datadir = LHAPDF_DATA_PREFIX;
rtn.push_back(datadir / "LHAPDF");
}
return rtn;
}
void setPaths(const std::string& pathstr) {
setenv("LHAPDF_DATA_PATH", pathstr.c_str(), 1);
}
string findFile(const string& target) {
if (target.empty()) return "";
for (const string& base : paths()) {
const string p = startswith(target, "/") ? target : base / target;
// if (verbosity() > 2) cout << "Trying file: " << p << endl;
if (file_exists(p)) {
// if (verbosity() > 1) cout << "Found file: " << p << endl;
return p;
}
}
return "";
}
const std::vector<std::string>& availablePDFSets() {
// Cached path list
static vector<string> rtn;
// Return cached list if valid
if (!rtn.empty()) return rtn;
// Otherwise this is the first time: populate the list
for (const string& p : paths()) {
if (!dir_exists(p)) continue;
DIR* dir;
struct dirent* ent;
if ((dir = opendir(p.c_str())) != NULL) {
while ((ent = readdir(dir)) != NULL) {
const string d = ent->d_name;
const string infopath = p / d / d + ".info";
if (file_exists(infopath)) {
if (!contains(rtn, d)) {
// cout << "@" << d << "@" << endl;
rtn.push_back(d); //< add if a set with this name isn't already known
}
}
}
closedir(dir);
}
sort(rtn.begin(), rtn.end());
}
return rtn;
}
}
diff --git a/src/Utils.cc b/src/Utils.cc
--- a/src/Utils.cc
+++ b/src/Utils.cc
@@ -1,314 +1,316 @@
-#include <iostream>
-#include <cmath>
-
-using namespace std;
+// -*- C++ -*-
+//
+// This file is part of LHAPDF
+// Copyright (C) 2012-2016 The LHAPDF collaboration (see AUTHORS for details)
+//
+#include "LHAPDF/Utils.h"
namespace LHAPDF {
namespace {
/// @todo Tidy up
static const double kMACHEP = 1.11022302462515654042363166809e-16;
static const double kMAXLOG = 709.782712893383973096206318587;
static const double kBig = 4.503599627370496e15;
static const double kBiginv = 2.22044604925031308085e-16;
double igamc(double a, double x);
double igam(double a, double x);
/// @name gamma functions from Cephes library -- http://www.netlib.org/cephes
///
/// Copyright 1985, 1987, 2000 by Stephen L. Moshier
//@{
/// @brief Incomplete gamma function (complement integral)
///
/// \f$ \gamma_c(a,x) = 1 - \gamma(a,x) \f$
/// \f$ \gamma_c(a,x) = 1/\Gamma(a) \int_x^\inf e^-t t^(a-1) dt \f$
///
/// In this implementation both arguments must be positive.
/// The integral is evaluated by either a power series or
/// continued fraction expansion, depending on the relative
/// values of a and x.
double igamc(double a, double x) {
double ans, ax, c, yc, r, t, y, z;
double pk, pkm1, pkm2, qk, qkm1, qkm2;
// LM: for negative values returns 0.0
// This is correct if a is a negative integer since Gamma(-n) = +/- inf
if (a <= 0) return 0.0;
if (x <= 0) return 1.0;
if( (x < 1.0) || (x < a) )
return( 1.0 - igam(a,x) );
ax = a * log(x) - x - lgamma(a);
if( ax < -kMAXLOG )
return( 0.0 );
ax = exp(ax);
// Continued fraction
y = 1.0 - a;
z = x + y + 1.0;
c = 0.0;
pkm2 = 1.0;
qkm2 = x;
pkm1 = x + 1.0;
qkm1 = z * x;
ans = pkm1/qkm1;
do
{
c += 1.0;
y += 1.0;
z += 2.0;
yc = y * c;
pk = pkm1 * z - pkm2 * yc;
qk = qkm1 * z - qkm2 * yc;
if(qk)
{
r = pk/qk;
t = fabs( (ans - r)/r );
ans = r;
}
else
t = 1.0;
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
if( fabs(pk) > kBig )
{
pkm2 *= kBiginv;
pkm1 *= kBiginv;
qkm2 *= kBiginv;
qkm1 *= kBiginv;
}
}
while (t > kMACHEP);
return ans*ax;
}
/// @brief Left tail of incomplete gamma function
///
/// \f$ \gamma(a,x) = x^a e^-x \sum_k=0^\inf x^k / \Gamma(a+k+1) \f$
double igam( double a, double x )
{
double ans, ax, c, r;
// LM: for negative values returns 1.0 instead of zero
// This is correct if a is a negative integer since Gamma(-n) = +/- inf
if (a <= 0) return 1.0;
if (x <= 0) return 0.0;
if( (x > 1.0) && (x > a ) )
return( 1.0 - igamc(a,x) );
/* Compute x**a * exp(-x) / gamma(a) */
ax = a * log(x) - x - lgamma(a);
if( ax < -kMAXLOG )
return( 0.0 );
ax = exp(ax);
/* power series */
r = a;
c = 1.0;
ans = 1.0;
do
{
r += 1.0;
c *= x/r;
ans += c;
}
while( c/ans > kMACHEP );
return( ans * ax/a );
}
//@}
}
/// @brief Compute quantiles for standard normal distribution N(0, 1) at probability p
///
/// ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3, 477-484.
double norm_quantile(double p) {
/// @todo Return +-inf
if (p <=0 || p >= 1) {
cerr << "norm_quantile: probability outside (0, 1)" << endl;
return 0;
}
const double a0 = 3.3871328727963666080e0;
const double a1 = 1.3314166789178437745e+2;
const double a2 = 1.9715909503065514427e+3;
const double a3 = 1.3731693765509461125e+4;
const double a4 = 4.5921953931549871457e+4;
const double a5 = 6.7265770927008700853e+4;
const double a6 = 3.3430575583588128105e+4;
const double a7 = 2.5090809287301226727e+3;
const double b1 = 4.2313330701600911252e+1;
const double b2 = 6.8718700749205790830e+2;
const double b3 = 5.3941960214247511077e+3;
const double b4 = 2.1213794301586595867e+4;
const double b5 = 3.9307895800092710610e+4;
const double b6 = 2.8729085735721942674e+4;
const double b7 = 5.2264952788528545610e+3;
const double c0 = 1.42343711074968357734e0;
const double c1 = 4.63033784615654529590e0;
const double c2 = 5.76949722146069140550e0;
const double c3 = 3.64784832476320460504e0;
const double c4 = 1.27045825245236838258e0;
const double c5 = 2.41780725177450611770e-1;
const double c6 = 2.27238449892691845833e-2;
const double c7 = 7.74545014278341407640e-4;
const double d1 = 2.05319162663775882187e0;
const double d2 = 1.67638483018380384940e0;
const double d3 = 6.89767334985100004550e-1;
const double d4 = 1.48103976427480074590e-1;
const double d5 = 1.51986665636164571966e-2;
const double d6 = 5.47593808499534494600e-4;
const double d7 = 1.05075007164441684324e-9;
const double e0 = 6.65790464350110377720e0;
const double e1 = 5.46378491116411436990e0;
const double e2 = 1.78482653991729133580e0;
const double e3 = 2.96560571828504891230e-1;
const double e4 = 2.65321895265761230930e-2;
const double e5 = 1.24266094738807843860e-3;
const double e6 = 2.71155556874348757815e-5;
const double e7 = 2.01033439929228813265e-7;
const double f1 = 5.99832206555887937690e-1;
const double f2 = 1.36929880922735805310e-1;
const double f3 = 1.48753612908506148525e-2;
const double f4 = 7.86869131145613259100e-4;
const double f5 = 1.84631831751005468180e-5;
const double f6 = 1.42151175831644588870e-7;
const double f7 = 2.04426310338993978564e-15;
const double split1 = 0.425;
const double split2=5.;
const double konst1=0.180625;
const double konst2=1.6;
double q, r, quantile;
q = p - 0.5;
if (fabs(q) < split1) {
r = konst1 - q*q;
quantile = q* (((((((a7 * r + a6) * r + a5) * r + a4) * r + a3)
* r + a2) * r + a1) * r + a0) /
(((((((b7 * r + b6) * r + b5) * r + b4) * r + b3)
* r + b2) * r + b1) * r + 1.);
} else {
r = (q < 0) ? p : 1-p;
//error case
if (r<=0)
quantile=0;
else {
r=sqrt(-log(r));
if (r<=split2) {
r=r-konst2;
quantile=(((((((c7 * r + c6) * r + c5) * r + c4) * r + c3)
* r + c2) * r + c1) * r + c0) /
(((((((d7 * r + d6) * r + d5) * r + d4) * r + d3)
* r + d2) * r + d1) * r + 1);
} else{
r=r-split2;
quantile=(((((((e7 * r + e6) * r + e5) * r + e4) * r + e3)
* r + e2) * r + e1) * r + e0) /
(((((((f7 * r + f6) * r + f5) * r + f4) * r + f3)
* r + f2) * r + f1) * r + 1);
}
if (q<0) quantile=-quantile;
}
}
return quantile;
}
/// @brief Compute quantiles of the chi-squared probability distribution function
///
/// Algorithm AS 91 Appl. Statist. (1975) Vol.24, P.35
/// implemented by Anna Kreshuk.
/// Incorporates the suggested changes in AS R85 (vol.40(1), pp.233-5, 1991)
/// Parameters:
/// @arg p - the probability value, at which the quantile is computed
/// @arg ndf - number of degrees of freedom
double chisquared_quantile(double p, double ndf) {
static const double c[] = {0, 0.01, 0.222222, 0.32, 0.4, 1.24, 2.2,
4.67, 6.66, 6.73, 13.32, 60.0, 70.0,
84.0, 105.0, 120.0, 127.0, 140.0, 175.0,
210.0, 252.0, 264.0, 294.0, 346.0, 420.0,
462.0, 606.0, 672.0, 707.0, 735.0, 889.0,
932.0, 966.0, 1141.0, 1182.0, 1278.0, 1740.0,
2520.0, 5040.0};
static const double e = 5e-7;
static const double aa = 0.6931471806;
static const int maxit = 20;
/// @todo Tidy
double ch, p1, p2, q, t, a, b, x;
double s1, s2, s3, s4, s5, s6;
if (ndf <= 0) return 0;
const double g = lgamma(0.5*ndf);
const double xx = 0.5 * ndf;
const double cp = xx - 1;
if (ndf >= log(p)*(-c[5])){
// Starting approximation for ndf less than or equal to 0.32
if (ndf > c[3]) {
x = norm_quantile(p);
// Starting approximation using Wilson and Hilferty estimate
p1 = c[2]/ndf;
ch = ndf*pow((x*sqrt(p1) + 1 - p1), 3);
if (ch > c[6]*ndf + 6)
ch = -2 * (log(1-p) - cp * log(0.5 * ch) + g);
} else {
ch = c[4];
a = log(1-p);
do {
q = ch;
p1 = 1 + ch * (c[7]+ch);
p2 = ch * (c[9] + ch * (c[8] + ch));
t = -0.5 + (c[7] + 2 * ch) / p1 - (c[9] + ch * (c[10] + 3 * ch)) / p2;
ch = ch - (1 - exp(a + g + 0.5 * ch + cp * aa) *p2 / p1) / t;
} while (fabs(q/ch - 1) > c[1]);
}
} else {
ch = pow((p * xx * exp(g + xx * aa)),(1./xx));
if (ch < e) return ch;
}
// Call to algorithm AS 239 and calculation of seven term Taylor series
for (int i = 0; i < maxit; ++i) {
q = ch;
p1 = 0.5 * ch;
p2 = p - igam(xx, p1);
t = p2 * exp(xx * aa + g + p1 - cp * log(ch));
b = t / ch;
a = 0.5 * t - b * cp;
s1 = (c[19] + a * (c[17] + a * (c[14] + a * (c[13] + a * (c[12] +c[11] * a))))) / c[24];
s2 = (c[24] + a * (c[29] + a * (c[32] + a * (c[33] + c[35] * a)))) / c[37];
s3 = (c[19] + a * (c[25] + a * (c[28] + c[31] * a))) / c[37];
s4 = (c[20] + a * (c[27] + c[34] * a) + cp * (c[22] + a * (c[30] + c[36] * a))) / c[38];
s5 = (c[13] + c[21] * a + cp * (c[18] + c[26] * a)) / c[37];
s6 = (c[15] + cp * (c[23] + c[16] * cp)) / c[38];
ch = ch + t * (1 + 0.5 * t * s1 - b * cp * (s1 - b * (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6))))));
if (fabs(q / ch - 1) > e) break;
}
return ch;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 5:42 PM (12 h, 29 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023681
Default Alt Text
(151 KB)

Event Timeline