lhapdf is hosted by Hepforge, IPPP Durham
LHAPDF 6.5.3
Loading...
Searching...
No Matches
KnotArray.h
1// -*- C++ -*-
2//
3// This file is part of LHAPDF
4// Copyright (C) 2012-2022 The LHAPDF collaboration (see AUTHORS for details)
5//
6#pragma once
7#ifndef LHAPDF_KnotArray_H
8#define LHAPDF_KnotArray_H
9
10#include "LHAPDF/Exceptions.h"
11#include "LHAPDF/Utils.h"
12
13namespace {
14
15
16 // Hide some internal functions from outside API view
17
18 // General function to find the knot below a given value
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; // can't return the last knot index
22 i -= 1; // step back to get the knot <= x behaviour
23 return i;
24 }
25
26
27 int findPidInPids(int pid, const std::vector<int>& pids) {
28 std::vector<int>::const_iterator it = std::find(pids.begin(), pids.end(), pid);
29 if (it == pids.end())
30 return -1;
31 else
32 return std::distance(pids.begin(), it);
33 }
34
35
36}
37
38namespace LHAPDF {
39
40
45 class KnotArray{
46 public:
47
49 size_t size() const { return _shape.back(); }
50
52 size_t xsize() const { return _shape[0]; }
53
55 size_t q2size() const { return _shape[1]; }
56
58 bool empty() const { return _grid.empty(); }
59
61 size_t ixbelow(double x) const { return indexbelow(x, _xs); }
62
64 size_t iq2below(double q2) const { return indexbelow(q2, _q2s); }
65
67 double xf(int ix, int iq2, int ipid) const {
68 return _grid[ix*_shape[2]*_shape[1] + iq2*_shape[2] + ipid];
69 }
70
72 const double& coeff(int ix, int iq2, int pid, int in) const {
73 return _coeffs[ix*(_shape[1])*_shape[2]*4 + iq2*_shape[2]*4 + pid*4 + in];
74 }
75
77 int lookUpPid(size_t id) const { return _lookup[id]; }
78
79 double xs(size_t id) const { return _xs[id]; }
80
81 double logxs(size_t id) const { return _logxs[id]; }
82
83 double q2s(size_t id) const { return _q2s[id]; }
84
85 double logq2s(size_t id) const { return _logq2s[id]; }
86
87 size_t shape(size_t id) const { return _shape[id]; }
88
90 bool inRangeX(double x) const {
91 if (x < _xs.front()) return false;
92 if (x > _xs.back()) return false;
93 return true;
94 }
95
97 bool inRangeQ2(double q2) const {
98 if (q2 < _q2s.front()) return false;
99 if (q2 > _q2s.back()) return false;
100 return true;
101 }
102
103 inline int get_pid(int id) const {
104 // hardcoded lookup table for particle ids
105 // -6,...,-1,21/0,1,...,6,22
106 // if id outside of this range, search in list of ids
107 if (-6 <= id && id <= 6) return _lookup[id + 6];
108 else if (id == 21) return _lookup[0 + 6];
109 else if (id == 22) return _lookup[13];
110 else return findPidInPids(id, _pids);
111 }
112
113 bool has_pid(int id) const {
114 return get_pid(id) != -1;
115 }
116
117 void initPidLookup();
118
119 void fillLogKnots();
120
122 const std::vector<double>& xs() const { return _xs; }
123
124 const std::vector<double>& logxs() const { return _logxs; }
125
126 const std::vector<double>& q2s() const { return _q2s; }
127
128 const std::vector<double>& logq2s() const { return _logq2s; }
129
131 std::vector<double>& setCoeffs() { return _coeffs; }
132
133 std::vector<double>& setGrid() { return _grid; }
134
135 std::vector<double>& setxknots() { return _xs; }
136
137 std::vector<double>& setq2knots() { return _q2s; }
138
139 std::vector<size_t>& setShape(){ return _shape; }
140
141 std::vector<int>& setPids() { return _pids; }
142
143 private:
144 // Shape of the interpolation grid
145 std::vector<size_t> _shape;
146
147 // Gridvalues
148 std::vector<double> _grid;
149
150 // Storage for the precomputed polynomial coefficients
151 std::vector<double> _coeffs;
152
153 // order the pids are filled in
154 std::vector<int> _pids;
155 std::vector<int> _lookup;
156
157 // knots
158 std::vector<double> _xs;
159 std::vector<double> _q2s;
160 std::vector<double> _logxs;
161 std::vector<double> _logq2s;
162
163 };
164
165
168 public:
169
172
175
177 AlphaSArray(const std::vector<double>& q2knots, const std::vector<double>& as)
178 : _q2s(q2knots), _as(as)
179 {
180 _syncq2s();
181 }
182
184
185
188
190 const std::vector<double>& q2s() const { return _q2s; }
191
193 const std::vector<double>& logq2s() const { return _logq2s; }
194
198 size_t iq2below(double q2) const {
199 // Test that Q2 is in the grid range
200 if (q2 < q2s().front()) throw AlphaSError("Q2 value " + to_str(q2) + " is lower than lowest-Q2 grid point at " + to_str(q2s().front()));
201 if (q2 > q2s().back()) throw AlphaSError("Q2 value " + to_str(q2) + " is higher than highest-Q2 grid point at " + to_str(q2s().back()));
203 size_t i = upper_bound(q2s().begin(), q2s().end(), q2) - q2s().begin();
204 if (i == q2s().size()) i -= 1; // can't return the last knot index
205 i -= 1; // have to step back to get the knot <= q2 behaviour
206 return i;
207 }
208
212 size_t ilogq2below(double logq2) const {
213 // Test that log(Q2) is in the grid range
214 if (logq2 < logq2s().front()) throw GridError("logQ2 value " + to_str(logq2) + " is lower than lowest-logQ2 grid point at " + to_str(logq2s().front()));
215 if (logq2 > logq2s().back()) throw GridError("logQ2 value " + to_str(logq2) + " is higher than highest-logQ2 grid point at " + to_str(logq2s().back()));
217 size_t i = upper_bound(logq2s().begin(), logq2s().end(), logq2) - logq2s().begin();
218 if (i == logq2s().size()) i -= 1; // can't return the last knot index
219 i -= 1; // have to step back to get the knot <= q2 behaviour
220 return i;
221 }
222
224
225
228
230 const std::vector<double>& alphas() const { return _as; }
231 // /// alpha_s value accessor (non-const)
232 // std::vector<double>& alphas() { return _as; }
233 // /// alpha_s value setter
234 // void setalphas(const valarray& xfs) { _as = as; }
235
237
238
241
243 double ddlogq_forward(size_t i) const {
244 return (alphas()[i+1] - alphas()[i]) / (logq2s()[i+1] - logq2s()[i]);
245 }
246
248 double ddlogq_backward(size_t i) const {
249 return (alphas()[i] - alphas()[i-1]) / (logq2s()[i] - logq2s()[i-1]);
250 }
251
253 double ddlogq_central(size_t i) const {
254 return 0.5 * (ddlogq_forward(i) + ddlogq_backward(i));
255 }
256
258
259
260 private:
261
263 void _syncq2s() {
264 _logq2s.resize(_q2s.size());
265 for (size_t i = 0; i < _q2s.size(); ++i) _logq2s[i] = log(_q2s[i]);
266 }
267
269 std::vector<double> _q2s;
271 std::vector<double> _logq2s;
273 std::vector<double> _as;
274
275 };
276}
277#endif
Internal storage class for alpha_s interpolation grids.
Definition: KnotArray.h:167
size_t iq2below(double q2) const
Definition: KnotArray.h:198
const std::vector< double > & logq2s() const
log(Q2) knot vector accessor
Definition: KnotArray.h:193
std::vector< double > _as
List of alpha_s values across the knot array.
Definition: KnotArray.h:273
AlphaSArray(const std::vector< double > &q2knots, const std::vector< double > &as)
Constructor from Q2 knot values and alpha_s values.
Definition: KnotArray.h:177
std::vector< double > _logq2s
List of log(Q2) knots.
Definition: KnotArray.h:271
AlphaSArray()
Default constructor just for std::map insertability.
Definition: KnotArray.h:174
double ddlogq_forward(size_t i) const
Forward derivative w.r.t. logQ2.
Definition: KnotArray.h:243
std::vector< double > _q2s
List of Q2 knots.
Definition: KnotArray.h:269
void _syncq2s()
Synchronise the log(Q2) array from the Q2 one.
Definition: KnotArray.h:263
const std::vector< double > & q2s() const
Q2 knot vector accessor.
Definition: KnotArray.h:190
double ddlogq_backward(size_t i) const
Backward derivative w.r.t. logQ2.
Definition: KnotArray.h:248
size_t ilogq2below(double logq2) const
Definition: KnotArray.h:212
double ddlogq_central(size_t i) const
Central (avg of forward and backward) derivative w.r.t. logQ2.
Definition: KnotArray.h:253
const std::vector< double > & alphas() const
alpha_s value accessor (const)
Definition: KnotArray.h:230
Error for general AlphaS computation problems.
Definition: Exceptions.h:94
Error for general PDF grid problems.
Definition: Exceptions.h:30
Internal storage class for PDF data point grids.
Definition: KnotArray.h:45
size_t ixbelow(double x) const
find the largest grid index below given x, such that xknots[index] < x
Definition: KnotArray.h:61
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:72
size_t q2size() const
How many q2 knots are there.
Definition: KnotArray.h:55
int lookUpPid(size_t id) const
accessor to the internal 'lookup table' for the pid's
Definition: KnotArray.h:77
bool inRangeX(double x) const
check if value within the boundaries of xknots
Definition: KnotArray.h:90
std::vector< double > & setCoeffs()
Non const accessors for programmatic filling.
Definition: KnotArray.h:131
double xf(int ix, int iq2, int ipid) const
convenient accessor to the grid values
Definition: KnotArray.h:67
const std::vector< double > & xs() const
Const accessors to the internal data container.
Definition: KnotArray.h:122
bool inRangeQ2(double q2) const
check if value within the boundaries of q2knots
Definition: KnotArray.h:97
size_t xsize() const
How many x knots are there.
Definition: KnotArray.h:52
size_t size() const
How many flavours are stored in the grid stored.
Definition: KnotArray.h:49
size_t iq2below(double q2) const
find the largest grid index below given q2, such that q2knots[index] < q2
Definition: KnotArray.h:64
bool empty() const
Is this container empty?
Definition: KnotArray.h:58
std::string to_str(const T &val)
Make a string representation of val.
Definition: Utils.h:61
Namespace for all LHAPDF functions and classes.
Definition: AlphaS.h:14