LHAPDF  6.3.0
KnotArray.h
1 // -*- C++ -*-
2 //
3 // This file is part of LHAPDF
4 // Copyright (C) 2012-2020 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 LHAPDF {
14 
15 
20  class KnotArray1F {
21  public:
22 
25 
27  KnotArray1F(const std::vector<double>& xknots, const std::vector<double>& q2knots, const std::vector<double>& xfs)
28  : _xs(xknots), _q2s(q2knots), _xfs(xfs)
29  {
30  assert(_xfs.size() == size());
31  _syncx();
32  _syncq2();
33  }
34 
36  KnotArray1F(const std::vector<double>& xknots, const std::vector<double>& q2knots)
37  : _xs(xknots), _q2s(q2knots),
38  _xfs(size(), 0.0)
39  {
40  assert(_xfs.size() == size());
41  _syncx();
42  _syncq2();
43  }
44 
45 
48 
52  void setxs(const std::vector<double>& xs) {
53  _xs = xs;
54  _syncx();
55  _xfs = std::vector<double>(size(), 0.0);
56  }
57 
59  size_t xsize() const { return _xs.size(); }
60 
62  const std::vector<double>& xs() const { return _xs; }
63 
65  const std::vector<double>& logxs() const { return _logxs; }
66 
68  bool samexs(const KnotArray1F& other) const { return _xgridhash == other._xgridhash; }
69  size_t xhash() const { return _xgridhash; }
70 
74  size_t ixbelow(double x) const {
75  // Test that x is in the grid range
76  if (x < xs().front()) throw GridError("x value " + to_str(x) + " is lower than lowest-x grid point at " + to_str(xs().front()));
77  if (x > xs().back()) throw GridError("x value " + to_str(x) + " is higher than highest-x grid point at " + to_str(xs().back()));
78  // Find the closest knot below the requested value
79  size_t i = upper_bound(xs().begin(), xs().end(), x) - xs().begin();
80  if (i == xs().size()) i -= 1; // can't return the last knot index
81  i -= 1; // have to step back to get the knot <= x behaviour
82  return i;
83  }
84 
86 
87 
90 
94  void setq2s(const std::vector<double>& q2s) {
95  _q2s = q2s;
96  _syncq2();
97  _xfs = std::vector<double>(size(), 0.0);
98  }
99 
101  size_t q2size() const { return _q2s.size(); }
102 
104  const std::vector<double>& q2s() const { return _q2s; }
105 
107  const std::vector<double>& logq2s() const { return _logq2s; }
108 
110  bool sameq2s(const KnotArray1F& other) const { return _q2gridhash == other._q2gridhash; }
111  size_t q2hash() const { return _q2gridhash; }
112 
116  size_t iq2below(double q2) const {
117  // Test that Q2 is in the grid range
118  if (q2 < q2s().front()) throw GridError("Q2 value " + to_str(q2) + " is lower than lowest-Q2 grid point at " + to_str(q2s().front()));
119  if (q2 > q2s().back()) throw GridError("Q2 value " + to_str(q2) + " is higher than highest-Q2 grid point at " + to_str(q2s().back()));
121  size_t i = upper_bound(q2s().begin(), q2s().end(), q2) - q2s().begin();
122  if (i == q2s().size()) i -= 1; // can't return the last knot index
123  i -= 1; // have to step back to get the knot <= q2 behaviour
124  return i;
125  }
126 
128 
129 
132 
134  size_t size() const { return xsize()*q2size(); }
135 
137  const std::vector<double>& xfs() const { return _xfs; }
139  std::vector<double>& xfs() { return _xfs; }
141  void setxfs(const std::vector<double>& xfs) { _xfs = xfs; }
142 
144  const double& xf(size_t ix, size_t iq2) const { return _xfs[ix*q2size() + iq2]; }
145 
147 
148 
149  private:
150 
152  void _syncx() {
153  _logxs.resize(_xs.size());
154  for (size_t i = 0; i < _xs.size(); ++i) _logxs[i] = log(_xs[i]);
156  }
157 
158 
160  void _syncq2() {
161  _logq2s.resize(_q2s.size());
162  for (size_t i = 0; i < _q2s.size(); ++i) _logq2s[i] = log(_q2s[i]);
164  }
165 
166 
168  size_t _mkhash(const std::vector<double>& xx) const;
169 
171  std::vector<double> _xs;
173  std::vector<double> _logxs;
175  size_t _xgridhash = 0;
176 
178  std::vector<double> _q2s;
180  std::vector<double> _logq2s;
182  size_t _q2gridhash = 0;
183 
185  std::vector<double> _xfs;
186 
187 
188  };
189 
190 
194  class KnotArrayNF {
195  public:
196 
198  size_t size() const { return _map.size(); }
199 
201  bool empty() const { return _map.empty(); }
202 
204  bool has_pid(int id) const {
205  return _map.find(id) != _map.end();
206  }
207 
209  const KnotArray1F& get_pid(int id) const {
210  if (!has_pid(id)) throw FlavorError("Undefined particle ID requested: " + to_str(id));
211  return _map.find(id)->second;
212  }
213 
215  const KnotArray1F& get_first() const {
216  if (empty()) throw GridError("Tried to access grid indices when no flavour grids were loaded");
217  return _map.begin()->second;
218  }
219 
221  void set_pid(int id, const KnotArray1F& ka) {
222  _map[id] = ka;
223  }
224 
226  KnotArray1F& operator[](int id) { return _map[id]; }
227 
229  const std::vector<double>& xs() const { return get_first().xs(); }
231  const std::vector<double>& logxs() const { return get_first().logxs(); }
233  size_t ixbelow(double x) const { return get_first().ixbelow(x); }
234 
236  const std::vector<double>& q2s() const { return get_first().q2s(); }
238  const std::vector<double>& logq2s() const { return get_first().logq2s(); }
240  size_t iq2below(double q2) const { return get_first().iq2below(q2); }
241 
242  private:
243 
245  std::map<int, KnotArray1F> _map;
246 
247  };
248 
249 
250 
252  class AlphaSArray {
253  public:
254 
257 
260 
262  AlphaSArray(const std::vector<double>& q2knots, const std::vector<double>& as)
263  : _q2s(q2knots), _as(as)
264  {
265  _syncq2s();
266  }
267 
269 
270 
273 
275  const std::vector<double>& q2s() const { return _q2s; }
276 
278  const std::vector<double>& logq2s() const { return _logq2s; }
279 
283  size_t iq2below(double q2) const {
284  // Test that Q2 is in the grid range
285  if (q2 < q2s().front()) throw AlphaSError("Q2 value " + to_str(q2) + " is lower than lowest-Q2 grid point at " + to_str(q2s().front()));
286  if (q2 > q2s().back()) throw AlphaSError("Q2 value " + to_str(q2) + " is higher than highest-Q2 grid point at " + to_str(q2s().back()));
288  size_t i = upper_bound(q2s().begin(), q2s().end(), q2) - q2s().begin();
289  if (i == q2s().size()) i -= 1; // can't return the last knot index
290  i -= 1; // have to step back to get the knot <= q2 behaviour
291  return i;
292  }
293 
297  size_t ilogq2below(double logq2) const {
298  // Test that log(Q2) is in the grid range
299  if (logq2 < logq2s().front()) throw GridError("logQ2 value " + to_str(logq2) + " is lower than lowest-logQ2 grid point at " + to_str(logq2s().front()));
300  if (logq2 > logq2s().back()) throw GridError("logQ2 value " + to_str(logq2) + " is higher than highest-logQ2 grid point at " + to_str(logq2s().back()));
302  size_t i = upper_bound(logq2s().begin(), logq2s().end(), logq2) - logq2s().begin();
303  if (i == logq2s().size()) i -= 1; // can't return the last knot index
304  i -= 1; // have to step back to get the knot <= q2 behaviour
305  return i;
306  }
307 
309 
310 
313 
315  const std::vector<double>& alphas() const { return _as; }
316  // /// alpha_s value accessor (non-const)
317  // std::vector<double>& alphas() { return _as; }
318  // /// alpha_s value setter
319  // void setalphas(const valarray& xfs) { _as = as; }
320 
322 
323 
326 
328  double ddlogq_forward(size_t i) const {
329  return (alphas()[i+1] - alphas()[i]) / (logq2s()[i+1] - logq2s()[i]);
330  }
331 
333  double ddlogq_backward(size_t i) const {
334  return (alphas()[i] - alphas()[i-1]) / (logq2s()[i] - logq2s()[i-1]);
335  }
336 
338  double ddlogq_central(size_t i) const {
339  return 0.5 * (ddlogq_forward(i) + ddlogq_backward(i));
340  }
341 
343 
344 
345  private:
346 
348  void _syncq2s() {
349  _logq2s.resize(_q2s.size());
350  for (size_t i = 0; i < _q2s.size(); ++i) _logq2s[i] = log(_q2s[i]);
351  }
352 
354  std::vector<double> _q2s;
356  std::vector<double> _logq2s;
358  std::vector<double> _as;
359 
360  };
361 
362 
363 }
364 #endif
Internal storage class for alpha_s interpolation grids.
Definition: KnotArray.h:252
size_t iq2below(double q2) const
Definition: KnotArray.h:283
const std::vector< double > & alphas() const
alpha_s value accessor (const)
Definition: KnotArray.h:315
std::vector< double > _as
List of alpha_s values across the knot array.
Definition: KnotArray.h:358
AlphaSArray(const std::vector< double > &q2knots, const std::vector< double > &as)
Constructor from Q2 knot values and alpha_s values.
Definition: KnotArray.h:262
std::vector< double > _logq2s
List of log(Q2) knots.
Definition: KnotArray.h:356
AlphaSArray()
Default constructor just for std::map insertability.
Definition: KnotArray.h:259
double ddlogq_forward(size_t i) const
Forward derivative w.r.t. logQ2.
Definition: KnotArray.h:328
std::vector< double > _q2s
List of Q2 knots.
Definition: KnotArray.h:354
const std::vector< double > & logq2s() const
log(Q2) knot vector accessor
Definition: KnotArray.h:278
void _syncq2s()
Synchronise the log(Q2) array from the Q2 one.
Definition: KnotArray.h:348
double ddlogq_backward(size_t i) const
Backward derivative w.r.t. logQ2.
Definition: KnotArray.h:333
const std::vector< double > & q2s() const
Q2 knot vector accessor.
Definition: KnotArray.h:275
size_t ilogq2below(double logq2) const
Definition: KnotArray.h:297
double ddlogq_central(size_t i) const
Central (avg of forward and backward) derivative w.r.t. logQ2.
Definition: KnotArray.h:338
Error for general AlphaS computation problems.
Definition: Exceptions.h:94
Error for requests for unsupported/invalid flavour PIDs.
Definition: Exceptions.h:70
Error for general PDF grid problems.
Definition: Exceptions.h:30
Internal storage class for PDF data point grids.
Definition: KnotArray.h:20
size_t _q2gridhash
Hash for this set of Q2 knots.
Definition: KnotArray.h:182
size_t size() const
Number of x knots.
Definition: KnotArray.h:134
std::vector< double > _q2s
List of Q2 knots.
Definition: KnotArray.h:178
void _syncq2()
Synchronise log(x) and log(Q2) arrays from the x and Q2 ones.
Definition: KnotArray.h:160
size_t q2size() const
Number of Q2 knots.
Definition: KnotArray.h:101
std::vector< double > _logq2s
List of log(Q2) knots, precomputed for efficiency.
Definition: KnotArray.h:180
KnotArray1F(const std::vector< double > &xknots, const std::vector< double > &q2knots)
Constructor of a zero-valued array from x and Q2 knot values.
Definition: KnotArray.h:36
KnotArray1F()
Default constructor just for std::map insertability.
Definition: KnotArray.h:24
std::vector< double > _xfs
List of xf values across the 2D knot array, stored as a strided [ix][iQ2] 1D array.
Definition: KnotArray.h:185
void setq2s(const std::vector< double > &q2s)
Definition: KnotArray.h:94
size_t iq2below(double q2) const
Definition: KnotArray.h:116
const std::vector< double > & xs() const
x knot accessor
Definition: KnotArray.h:62
std::vector< double > _xs
List of x knots.
Definition: KnotArray.h:171
size_t _xgridhash
Hash for this set of x knots.
Definition: KnotArray.h:175
const std::vector< double > & xfs() const
xf value accessor (const)
Definition: KnotArray.h:137
std::vector< double > & xfs()
xf value accessor (non-const)
Definition: KnotArray.h:139
size_t ixbelow(double x) const
Get the index of the closest x knot row <= x.
Definition: KnotArray.h:74
KnotArray1F(const std::vector< double > &xknots, const std::vector< double > &q2knots, const std::vector< double > &xfs)
Constructor from x and Q2 knot values, and an xf value grid as strided list.
Definition: KnotArray.h:27
const std::vector< double > & logxs() const
log(x) knot accessor
Definition: KnotArray.h:65
size_t _mkhash(const std::vector< double > &xx) const
Utility function for making a hash code from a vector<double>
const double & xf(size_t ix, size_t iq2) const
Get the xf value at a particular indexed x,Q2 knot.
Definition: KnotArray.h:144
void setxfs(const std::vector< double > &xfs)
xf value setter
Definition: KnotArray.h:141
const std::vector< double > & q2s() const
Q2 knot accessor.
Definition: KnotArray.h:104
bool samexs(const KnotArray1F &other) const
Hash comparator.
Definition: KnotArray.h:68
bool sameq2s(const KnotArray1F &other) const
Hash comparator for Q2 knots.
Definition: KnotArray.h:110
void _syncx()
Synchronise log(x) array and hash from the x array.
Definition: KnotArray.h:152
std::vector< double > _logxs
List of log(x) knots, precomputed for efficiency.
Definition: KnotArray.h:173
void setxs(const std::vector< double > &xs)
Definition: KnotArray.h:52
const std::vector< double > & logq2s() const
log(Q2) knot accessor
Definition: KnotArray.h:107
size_t xsize() const
Number of x knots.
Definition: KnotArray.h:59
A collection of {KnotArray1F}s accessed by PID code.
Definition: KnotArray.h:194
std::map< int, KnotArray1F > _map
Storage.
Definition: KnotArray.h:245
bool empty() const
Is this container empty?
Definition: KnotArray.h:201
const std::vector< double > & logxs() const
Access the log(x)s array.
Definition: KnotArray.h:231
size_t size() const
How many {KnotArray1F}s are stored in this container?
Definition: KnotArray.h:198
void set_pid(int id, const KnotArray1F &ka)
Get the KnotArray1F for PID code id.
Definition: KnotArray.h:221
KnotArray1F & operator[](int id)
Indexing operator (non-const)
Definition: KnotArray.h:226
size_t ixbelow(double x) const
Get the index of the closest x knot column <= x (see KnotArray1F)
Definition: KnotArray.h:233
const std::vector< double > & q2s() const
Access the Q2s array.
Definition: KnotArray.h:236
const std::vector< double > & xs() const
Access the xs array.
Definition: KnotArray.h:229
const KnotArray1F & get_first() const
Convenience accessor for any valid subgrid, to get access to the x/Q2/etc. arrays.
Definition: KnotArray.h:215
const std::vector< double > & logq2s() const
Access the log(Q2)s array.
Definition: KnotArray.h:238
const KnotArray1F & get_pid(int id) const
Get the KnotArray1F for PID code id.
Definition: KnotArray.h:209
bool has_pid(int id) const
Does this contain a KnotArray1F for PID code id?
Definition: KnotArray.h:204
size_t iq2below(double q2) const
Get the index of the closest Q2 knot row <= q2 (see KnotArray1F)
Definition: KnotArray.h:240
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