lhapdf is hosted by Hepforge, IPPP Durham
LHAPDF  6.5.4
KnotArray.h
1 // -*- C++ -*-
2 //
3 // This file is part of LHAPDF
4 // Copyright (C) 2012-2023 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 
13 namespace {
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 
39 namespace LHAPDF {
40 
41 
46  class KnotArray{
47  public:
48 
50  size_t size() const { return _shape.back(); }
51 
53  size_t xsize() const { return _shape[0]; }
54 
56  size_t q2size() const { return _shape[1]; }
57 
59  bool empty() const { return _grid.empty(); }
60 
62  size_t ixbelow(double x) const { return indexbelow(x, _xs); }
63 
65  size_t iq2below(double q2) const { return indexbelow(q2, _q2s); }
66 
68  double xf(int ix, int iq2, int ipid) const {
69  return _grid[ix*_shape[2]*_shape[1] + iq2*_shape[2] + ipid];
70  }
71 
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 
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 
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 
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 
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 
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 
168  class AlphaSArray {
169  public:
170 
173 
176 
178  AlphaSArray(const std::vector<double>& q2knots, const std::vector<double>& as)
179  : _q2s(q2knots), _as(as)
180  {
181  _syncq2s();
182  }
183 
185 
186 
189 
191  const std::vector<double>& q2s() const { return _q2s; }
192 
194  const std::vector<double>& logq2s() const { return _logq2s; }
195 
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()));
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 
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()));
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 
225 
226 
229 
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 
238 
239 
242 
244  double ddlogq_forward(size_t i) const {
245  return (alphas()[i+1] - alphas()[i]) / (logq2s()[i+1] - logq2s()[i]);
246  }
247 
249  double ddlogq_backward(size_t i) const {
250  return (alphas()[i] - alphas()[i-1]) / (logq2s()[i] - logq2s()[i-1]);
251  }
252 
254  double ddlogq_central(size_t i) const {
255  return 0.5 * (ddlogq_forward(i) + ddlogq_backward(i));
256  }
257 
259 
260 
261  private:
262 
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 
270  std::vector<double> _q2s;
272  std::vector<double> _logq2s;
274  std::vector<double> _as;
275 
276  };
277 }
278 #endif
bool inRangeQ2(double q2) const
check if value within the boundaries of q2knots
Definition: KnotArray.h:98
int lookUpPid(size_t id) const
accessor to the internal &#39;lookup table&#39; for the pid&#39;s
Definition: KnotArray.h:78
const std::vector< double > & q2s() const
Q2 knot vector accessor.
Definition: KnotArray.h:191
Internal storage class for alpha_s interpolation grids.
Definition: KnotArray.h:168
std::vector< double > _q2s
List of Q2 knots.
Definition: KnotArray.h:270
bool inRangeX(double x) const
check if value within the boundaries of xknots
Definition: KnotArray.h:91
std::vector< double > _logq2s
List of log(Q2) knots.
Definition: KnotArray.h:272
size_t size() const
How many flavours are stored in the grid stored.
Definition: KnotArray.h:50
Internal storage class for PDF data point grids.
Definition: KnotArray.h:46
size_t xsize() const
How many x knots are there.
Definition: KnotArray.h:53
size_t ilogq2below(double logq2) const
Definition: KnotArray.h:213
Error for general AlphaS computation problems.
Definition: Exceptions.h:94
size_t iq2below(double q2) const
Definition: KnotArray.h:199
const std::vector< double > & xs() const
Const accessors to the internal data container.
Definition: KnotArray.h:123
bool empty() const
Is this container empty?
Definition: KnotArray.h:59
size_t iq2below(double q2) const
find the largest grid index below given q2, such that q2knots[index] &lt; q2
Definition: KnotArray.h:65
Error for general PDF grid problems.
Definition: Exceptions.h:30
void _syncq2s()
Synchronise the log(Q2) array from the Q2 one.
Definition: KnotArray.h:264
std::vector< double > & setCoeffs()
Non const accessors for programmatic filling.
Definition: KnotArray.h:132
std::string to_str(const T &val)
Make a string representation of val.
Definition: Utils.h:61
size_t q2size() const
How many q2 knots are there.
Definition: KnotArray.h:56
size_t ixbelow(double x) const
find the largest grid index below given x, such that xknots[index] &lt; x
Definition: KnotArray.h:62
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 > & logq2s() const
log(Q2) knot vector accessor
Definition: KnotArray.h:194
double ddlogq_backward(size_t i) const
Backward derivative w.r.t. logQ2.
Definition: KnotArray.h:249
double xf(int ix, int iq2, int ipid) const
convenient accessor to the grid values
Definition: KnotArray.h:68
std::vector< double > _as
List of alpha_s values across the knot array.
Definition: KnotArray.h:274
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
AlphaSArray()
Default constructor just for std::map insertability.
Definition: KnotArray.h:175
const std::vector< double > & alphas() const
alpha_s value accessor (const)
Definition: KnotArray.h:231
double ddlogq_forward(size_t i) const
Forward derivative w.r.t. logQ2.
Definition: KnotArray.h:244
AlphaSArray(const std::vector< double > &q2knots, const std::vector< double > &as)
Constructor from Q2 knot values and alpha_s values.
Definition: KnotArray.h:178