7#ifndef LHAPDF_KnotArray_H
8#define LHAPDF_KnotArray_H
10#include "LHAPDF/Exceptions.h"
11#include "LHAPDF/Utils.h"
19 size_t indexbelow(
double value,
const std::vector<
double>& knots) {
20 size_t i = upper_bound(knots.begin(), knots.end(), value) - knots.begin();
21 if (i == knots.size()) i -= 1;
27 int findPidInPids(
int pid,
const std::vector<
int>& pids) {
28 std::vector<
int>::const_iterator it = std::find(pids.begin(), pids.end(), pid);
33 return static_cast<
int>(std::distance(pids.begin(), it));
59 bool empty()
const {
return _grid.empty(); }
68 double xf(
int ix,
int iq2,
int ipid)
const {
69 return _grid[ix*_shape[2]*_shape[1] + iq2*_shape[2] + ipid];
73 const double&
coeff(
int ix,
int iq2,
int pid,
int in)
const {
74 return _coeffs[ix*(_shape[1])*_shape[2]*4 + iq2*_shape[2]*4 + pid*4 + in];
78 int lookUpPid(size_t id)
const {
return _lookup[id]; }
80 double xs(size_t id)
const {
return _xs[id]; }
82 double logxs(size_t id)
const {
return _logxs[id]; }
84 double q2s(size_t id)
const {
return _q2s[id]; }
86 double logq2s(size_t id)
const {
return _logq2s[id]; }
88 size_t shape(size_t id)
const {
return _shape[id]; }
92 if (x < _xs.front())
return false;
93 if (x > _xs.back())
return false;
99 if (q2 < _q2s.front())
return false;
100 if (q2 > _q2s.back())
return false;
104 inline int get_pid(
int id)
const {
108 if (-6 <= id && id <= 6)
return _lookup[id + 6];
109 else if (id == 21)
return _lookup[0 + 6];
110 else if (id == 22)
return _lookup[13];
111 else return findPidInPids(id, _pids);
114 bool has_pid(
int id)
const {
115 return get_pid(id) != -1;
118 void initPidLookup();
125 const std::vector<
double>& logxs()
const {
return _logxs; }
127 const std::vector<
double>& q2s()
const {
return _q2s; }
129 const std::vector<
double>& logq2s()
const {
return _logq2s; }
134 std::vector<
double>& setGrid() {
return _grid; }
136 std::vector<
double>& setxknots() {
return _xs; }
138 std::vector<
double>& setq2knots() {
return _q2s; }
140 std::vector<size_t>& setShape(){
return _shape; }
142 std::vector<
int>& setPids() {
return _pids; }
146 std::vector<size_t> _shape;
149 std::vector<
double> _grid;
152 std::vector<
double> _coeffs;
155 std::vector<
int> _pids;
156 std::vector<
int> _lookup;
159 std::vector<
double> _xs;
160 std::vector<
double> _q2s;
161 std::vector<
double> _logxs;
162 std::vector<
double> _logq2s;
178 AlphaSArray(
const std::vector<
double>& q2knots,
const std::vector<
double>& as)
201 if (q2 < q2s().front())
throw AlphaSError(
"Q2 value " + to_str(q2) +
" is lower than lowest-Q2 grid point at " + to_str(q2s().front()));
202 if (q2 > q2s().back())
throw AlphaSError(
"Q2 value " + to_str(q2) +
" is higher than highest-Q2 grid point at " + to_str(q2s().back()));
204 size_t i = upper_bound(q2s().begin(), q2s().end(), q2) - q2s().begin();
205 if (i == q2s().size()) i -= 1;
215 if (logq2 < logq2s().front())
throw GridError(
"logQ2 value " + to_str(logq2) +
" is lower than lowest-logQ2 grid point at " + to_str(logq2s().front()));
216 if (logq2 > logq2s().back())
throw GridError(
"logQ2 value " + to_str(logq2) +
" is higher than highest-logQ2 grid point at " + to_str(logq2s().back()));
218 size_t i = upper_bound(logq2s().begin(), logq2s().end(), logq2) - logq2s().begin();
219 if (i == logq2s().size()) i -= 1;
245 return (alphas()[i+1] - alphas()[i]) / (logq2s()[i+1] - logq2s()[i]);
250 return (alphas()[i] - alphas()[i-1]) / (logq2s()[i] - logq2s()[i-1]);
255 return 0.5 * (ddlogq_forward(i) + ddlogq_backward(i));
265 _logq2s.resize(_q2s.size());
266 for (size_t i = 0; i < _q2s.size(); ++i) _logq2s[i] = log(_q2s[i]);
Internal storage class for alpha_s interpolation grids.
Definition KnotArray.h:168
size_t iq2below(double q2) const
Definition KnotArray.h:199
const std::vector< double > & logq2s() const
log(Q2) knot vector accessor
Definition KnotArray.h:194
std::vector< double > _as
List of alpha_s values across the knot array.
Definition KnotArray.h:274
AlphaSArray(const std::vector< double > &q2knots, const std::vector< double > &as)
Constructor from Q2 knot values and alpha_s values.
Definition KnotArray.h:178
std::vector< double > _logq2s
List of log(Q2) knots.
Definition KnotArray.h:272
AlphaSArray()
Default constructor just for std::map insertability.
Definition KnotArray.h:175
double ddlogq_forward(size_t i) const
Forward derivative w.r.t. logQ2.
Definition KnotArray.h:244
std::vector< double > _q2s
List of Q2 knots.
Definition KnotArray.h:270
void _syncq2s()
Synchronise the log(Q2) array from the Q2 one.
Definition KnotArray.h:264
const std::vector< double > & q2s() const
Q2 knot vector accessor.
Definition KnotArray.h:191
double ddlogq_backward(size_t i) const
Backward derivative w.r.t. logQ2.
Definition KnotArray.h:249
size_t ilogq2below(double logq2) const
Definition KnotArray.h:213
double ddlogq_central(size_t i) const
Central (avg of forward and backward) derivative w.r.t. logQ2.
Definition KnotArray.h:254
const std::vector< double > & alphas() const
alpha_s value accessor (const)
Definition KnotArray.h:231
Internal storage class for PDF data point grids.
Definition KnotArray.h:46
size_t ixbelow(double x) const
find the largest grid index below given x, such that xknots[index] < x
Definition KnotArray.h:62
const double & coeff(int ix, int iq2, int pid, int in) const
convenient accessor to the polynomial coefficients, returns reference rather than value,...
Definition KnotArray.h:73
size_t q2size() const
How many q2 knots are there.
Definition KnotArray.h:56
int lookUpPid(size_t id) const
accessor to the internal 'lookup table' for the pid's
Definition KnotArray.h:78
bool inRangeX(double x) const
check if value within the boundaries of xknots
Definition KnotArray.h:91
std::vector< double > & setCoeffs()
Non const accessors for programmatic filling.
Definition KnotArray.h:132
double xf(int ix, int iq2, int ipid) const
convenient accessor to the grid values
Definition KnotArray.h:68
const std::vector< double > & xs() const
Const accessors to the internal data container.
Definition KnotArray.h:123
bool inRangeQ2(double q2) const
check if value within the boundaries of q2knots
Definition KnotArray.h:98
size_t xsize() const
How many x knots are there.
Definition KnotArray.h:53
size_t size() const
How many flavours are stored in the grid stored.
Definition KnotArray.h:50
size_t iq2below(double q2) const
find the largest grid index below given q2, such that q2knots[index] < q2
Definition KnotArray.h:65
bool empty() const
Is this container empty?
Definition KnotArray.h:59
Namespace for all LHAPDF functions and classes.
Definition AlphaS.h:14