lhapdf is hosted by Hepforge, IPPP Durham
LHAPDF 6.5.5
Loading...
Searching...
No Matches
KnotArray.h
1// -*- C++ -*-
2//
3// This file is part of LHAPDF
4// Copyright (C) 2012-2024 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 // save to cast, given small number of pid's
33 return static_cast<int>(std::distance(pids.begin(), it));
34 }
35
36
37}
38
39namespace LHAPDF {
40
41
42 /// @brief Internal storage class for PDF data point grids
43 ///
44 /// We use "array" to refer to the "raw" knot grid, while "grid" means a grid-based PDF.
45 /// The "1F" means that this is a single-flavour array
46 class KnotArray{
47 public:
48
49 /// How many flavours are stored in the grid stored
50 size_t size() const { return _shape.back(); }
51
52 /// How many x knots are there
53 size_t xsize() const { return _shape[0]; }
54
55 /// How many q2 knots are there
56 size_t q2size() const { return _shape[1]; }
57
58 /// Is this container empty?
59 bool empty() const { return _grid.empty(); }
60
61 /// find the largest grid index below given x, such that xknots[index] < x
62 size_t ixbelow(double x) const { return indexbelow(x, _xs); }
63
64 /// find the largest grid index below given q2, such that q2knots[index] < q2
65 size_t iq2below(double q2) const { return indexbelow(q2, _q2s); }
66
67 /// convenient accessor to the grid values
68 double xf(int ix, int iq2, int ipid) const {
69 return _grid[ix*_shape[2]*_shape[1] + iq2*_shape[2] + ipid];
70 }
71
72 /// convenient accessor to the polynomial coefficients, returns reference rather than value, to be able to read multiple adjacent at once
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];
75 }
76
77 /// accessor to the internal 'lookup table' for the pid's
78 int lookUpPid(size_t id) const { return _lookup[id]; }
79
80 double xs(size_t id) const { return _xs[id]; }
81
82 double logxs(size_t id) const { return _logxs[id]; }
83
84 double q2s(size_t id) const { return _q2s[id]; }
85
86 double logq2s(size_t id) const { return _logq2s[id]; }
87
88 size_t shape(size_t id) const { return _shape[id]; }
89
90 /// check if value within the boundaries of xknots
91 bool inRangeX(double x) const {
92 if (x < _xs.front()) return false;
93 if (x > _xs.back()) return false;
94 return true;
95 }
96
97 /// check if value within the boundaries of q2knots
98 bool inRangeQ2(double q2) const {
99 if (q2 < _q2s.front()) return false;
100 if (q2 > _q2s.back()) return false;
101 return true;
102 }
103
104 inline int get_pid(int id) const {
105 // hardcoded lookup table for particle ids
106 // -6,...,-1,21/0,1,...,6,22
107 // if id outside of this range, search in list of ids
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);
112 }
113
114 bool has_pid(int id) const {
115 return get_pid(id) != -1;
116 }
117
118 void initPidLookup();
119
120 void fillLogKnots();
121
122 /// Const accessors to the internal data container
123 const std::vector<double>& xs() const { return _xs; }
124
125 const std::vector<double>& logxs() const { return _logxs; }
126
127 const std::vector<double>& q2s() const { return _q2s; }
128
129 const std::vector<double>& logq2s() const { return _logq2s; }
130
131 /// Non const accessors for programmatic filling
132 std::vector<double>& setCoeffs() { return _coeffs; }
133
134 std::vector<double>& setGrid() { return _grid; }
135
136 std::vector<double>& setxknots() { return _xs; }
137
138 std::vector<double>& setq2knots() { return _q2s; }
139
140 std::vector<size_t>& setShape(){ return _shape; }
141
142 std::vector<int>& setPids() { return _pids; }
143
144 private:
145 // Shape of the interpolation grid
146 std::vector<size_t> _shape;
147
148 // Gridvalues
149 std::vector<double> _grid;
150
151 // Storage for the precomputed polynomial coefficients
152 std::vector<double> _coeffs;
153
154 // order the pids are filled in
155 std::vector<int> _pids;
156 std::vector<int> _lookup;
157
158 // knots
159 std::vector<double> _xs;
160 std::vector<double> _q2s;
161 std::vector<double> _logxs;
162 std::vector<double> _logq2s;
163
164 };
165
166
167 /// Internal storage class for alpha_s interpolation grids
169 public:
170
171 /// @name Construction etc.
172 ///@{
173
174 /// Default constructor just for std::map insertability
176
177 /// Constructor from Q2 knot values and alpha_s values
178 AlphaSArray(const std::vector<double>& q2knots, const std::vector<double>& as)
179 : _q2s(q2knots), _as(as)
180 {
182 }
183
184 ///@}
185
186
187 /// @name Q2 stuff
188 ///@{
189
190 /// Q2 knot vector accessor
191 const std::vector<double>& q2s() const { return _q2s; }
192
193 /// log(Q2) knot vector accessor
194 const std::vector<double>& logq2s() const { return _logq2s; }
195
196 /// Get the index of the closest Q2 knot row <= q2
197 ///
198 /// If the value is >= q2_max, return i_max-1 (for polynomial spine construction)
199 size_t iq2below(double q2) const {
200 // Test that Q2 is in the grid range
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()));
203 /// Find the closest knot below the requested value
204 size_t i = upper_bound(q2s().begin(), q2s().end(), q2) - q2s().begin();
205 if (i == q2s().size()) i -= 1; // can't return the last knot index
206 i -= 1; // have to step back to get the knot <= q2 behaviour
207 return i;
208 }
209
210 /// Get the index of the closest logQ2 knot row <= logq2
211 ///
212 /// If the value is >= q2_max, return i_max-1 (for polynomial spine construction)
213 size_t ilogq2below(double logq2) const {
214 // Test that log(Q2) is in the grid range
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()));
217 /// Find the closest knot below the requested value
218 size_t i = upper_bound(logq2s().begin(), logq2s().end(), logq2) - logq2s().begin();
219 if (i == logq2s().size()) i -= 1; // can't return the last knot index
220 i -= 1; // have to step back to get the knot <= q2 behaviour
221 return i;
222 }
223
224 ///@}
225
226
227 /// @name alpha_s values at Q2 points
228 ///@{
229
230 /// alpha_s value accessor (const)
231 const std::vector<double>& alphas() const { return _as; }
232 // /// alpha_s value accessor (non-const)
233 // std::vector<double>& alphas() { return _as; }
234 // /// alpha_s value setter
235 // void setalphas(const valarray& xfs) { _as = as; }
236
237 ///@}
238
239
240 /// @name alpha_s derivatives vs (log)Q2, useful for interpolation
241 ///@{
242
243 /// Forward derivative w.r.t. logQ2
244 double ddlogq_forward(size_t i) const {
245 return (alphas()[i+1] - alphas()[i]) / (logq2s()[i+1] - logq2s()[i]);
246 }
247
248 /// Backward derivative w.r.t. logQ2
249 double ddlogq_backward(size_t i) const {
250 return (alphas()[i] - alphas()[i-1]) / (logq2s()[i] - logq2s()[i-1]);
251 }
252
253 /// Central (avg of forward and backward) derivative w.r.t. logQ2
254 double ddlogq_central(size_t i) const {
255 return 0.5 * (ddlogq_forward(i) + ddlogq_backward(i));
256 }
257
258 ///@}
259
260
261 private:
262
263 /// Synchronise the log(Q2) array from the Q2 one
264 void _syncq2s() {
265 _logq2s.resize(_q2s.size());
266 for (size_t i = 0; i < _q2s.size(); ++i) _logq2s[i] = log(_q2s[i]);
267 }
268
269 /// List of Q2 knots
270 std::vector<double> _q2s;
271 /// List of log(Q2) knots
273 /// List of alpha_s values across the knot array
274 std::vector<double> _as;
275
276 };
277}
278#endif
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