LHAPDF is hosted by Hepforge, IPPP Durham
LHAPDF  6.2.1
KnotArray.h
1 // -*- C++ -*-
2 //
3 // This file is part of LHAPDF
4 // Copyright (C) 2012-2016 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 
12 namespace LHAPDF {
13 
14 
19  class KnotArray1F {
20  public:
21 
23 
24 
27 
29  KnotArray1F(const std::vector<double>& xknots, const std::vector<double>& q2knots, const std::vector<double>& xfs)
30  : _xs(xknots), _q2s(q2knots), _xfs(xfs)
31  {
32  assert(_xfs.size() == size());
33  _synclogs();
34  }
35 
37  KnotArray1F(const std::vector<double>& xknots, const std::vector<double>& q2knots)
38  : _xs(xknots), _q2s(q2knots),
39  _xfs(size(), 0.0)
40  {
41  assert(_xfs.size() == size());
42  _synclogs();
43  }
44 
46 
47 
49 
50 
53  void setxs(const std::vector<double>& xs) {
54  _xs = xs;
55  _synclogs();
56  _xfs = std::vector<double>(size(), 0.0);
57  }
58 
60  const size_t xsize() const { return _xs.size(); }
61 
63  const std::vector<double>& xs() const { return _xs; }
64 
66  const std::vector<double>& logxs() const { return _logxs; }
67 
71  size_t ixbelow(double x) const {
72  // Test that x is in the grid range
73  if (x < xs().front()) throw GridError("x value " + to_str(x) + " is lower than lowest-x grid point at " + to_str(xs().front()));
74  if (x > xs().back()) throw GridError("x value " + to_str(x) + " is higher than highest-x grid point at " + to_str(xs().back()));
75  // Find the closest knot below the requested value
76  size_t i = upper_bound(xs().begin(), xs().end(), x) - xs().begin();
77  if (i == xs().size()) i -= 1; // can't return the last knot index
78  i -= 1; // have to step back to get the knot <= x behaviour
79  return i;
80  }
81 
83 
84 
86 
87 
90  void setq2s(const std::vector<double>& q2s) {
91  _q2s = q2s;
92  _synclogs();
93  _xfs = std::vector<double>(size(), 0.0);
94  }
95 
97  const size_t q2size() const { return _q2s.size(); }
98 
100  const std::vector<double>& q2s() const { return _q2s; }
101 
103  const std::vector<double>& logq2s() const { return _logq2s; }
104 
108  size_t iq2below(double q2) const {
109  // Test that Q2 is in the grid range
110  if (q2 < q2s().front()) throw GridError("Q2 value " + to_str(q2) + " is lower than lowest-Q2 grid point at " + to_str(q2s().front()));
111  if (q2 > q2s().back()) throw GridError("Q2 value " + to_str(q2) + " is higher than highest-Q2 grid point at " + to_str(q2s().back()));
113  size_t i = upper_bound(q2s().begin(), q2s().end(), q2) - q2s().begin();
114  if (i == q2s().size()) i -= 1; // can't return the last knot index
115  i -= 1; // have to step back to get the knot <= q2 behaviour
116  return i;
117  }
118 
120 
121 
123 
124 
126  const size_t size() const { return xsize()*q2size(); }
127 
129  const std::vector<double>& xfs() const { return _xfs; }
131  std::vector<double>& xfs() { return _xfs; }
133  void setxfs(const std::vector<double>& xfs) { _xfs = xfs; }
134 
136  const double& xf(size_t ix, size_t iq2) const { return _xfs[ix*q2size() + iq2]; }
137 
139 
140 
141  private:
142 
144  void _synclogs() {
145  _logxs.resize(_xs.size());
146  _logq2s.resize(_q2s.size());
147  for (size_t i = 0; i < _xs.size(); ++i) _logxs[i] = log(_xs[i]);
148  for (size_t i = 0; i < _q2s.size(); ++i) _logq2s[i] = log(_q2s[i]);
149  }
150 
152  std::vector<double> _xs;
154  std::vector<double> _q2s;
156  std::vector<double> _logxs;
158  std::vector<double> _logq2s;
160  std::vector<double> _xfs;
161 
162  };
163 
164 
168  class KnotArrayNF {
169  public:
170 
172  size_t size() const { return _map.size(); }
173 
175  bool empty() const { return _map.empty(); }
176 
178  bool has_pid(int id) const {
179  return _map.find(id) != _map.end();
180  }
181 
183  const KnotArray1F& get_pid(int id) const {
184  if (!has_pid(id)) throw FlavorError("Undefined particle ID requested: " + to_str(id));
185  return _map.find(id)->second;
186  }
187 
189  const KnotArray1F& get_first() const {
190  if (empty()) throw GridError("Tried to access grid indices when no flavour grids were loaded");
191  return _map.begin()->second;
192  }
193 
195  void set_pid(int id, const KnotArray1F& ka) {
196  _map[id] = ka;
197  }
198 
200  KnotArray1F& operator[](int id) { return _map[id]; }
201 
203  const std::vector<double>& xs() const { return get_first().xs(); }
205  const std::vector<double>& logxs() const { return get_first().logxs(); }
207  size_t ixbelow(double x) const { return get_first().ixbelow(x); }
208 
210  const std::vector<double>& q2s() const { return get_first().q2s(); }
212  const std::vector<double>& logq2s() const { return get_first().logq2s(); }
214  size_t iq2below(double q2) const { return get_first().iq2below(q2); }
215 
216  private:
217 
219  std::map<int, KnotArray1F> _map;
220 
221  };
222 
223 
224 
226  class AlphaSArray {
227  public:
228 
230 
231 
234 
236  AlphaSArray(const std::vector<double>& q2knots, const std::vector<double>& as)
237  : _q2s(q2knots), _as(as)
238  {
239  _synclogs();
240  }
241 
243 
244 
246 
247 
249  const std::vector<double>& q2s() const { return _q2s; }
250 
252  const std::vector<double>& logq2s() const { return _logq2s; }
253 
257  size_t iq2below(double q2) const {
258  // Test that Q2 is in the grid range
259  if (q2 < q2s().front()) throw AlphaSError("Q2 value " + to_str(q2) + " is lower than lowest-Q2 grid point at " + to_str(q2s().front()));
260  if (q2 > q2s().back()) throw AlphaSError("Q2 value " + to_str(q2) + " is higher than highest-Q2 grid point at " + to_str(q2s().back()));
262  size_t i = upper_bound(q2s().begin(), q2s().end(), q2) - q2s().begin();
263  if (i == q2s().size()) i -= 1; // can't return the last knot index
264  i -= 1; // have to step back to get the knot <= q2 behaviour
265  return i;
266  }
267 
271  size_t ilogq2below(double logq2) const {
272  // Test that log(Q2) is in the grid range
273  if (logq2 < logq2s().front()) throw GridError("logQ2 value " + to_str(logq2) + " is lower than lowest-logQ2 grid point at " + to_str(logq2s().front()));
274  if (logq2 > logq2s().back()) throw GridError("logQ2 value " + to_str(logq2) + " is higher than highest-logQ2 grid point at " + to_str(logq2s().back()));
276  size_t i = upper_bound(logq2s().begin(), logq2s().end(), logq2) - logq2s().begin();
277  if (i == logq2s().size()) i -= 1; // can't return the last knot index
278  i -= 1; // have to step back to get the knot <= q2 behaviour
279  return i;
280  }
281 
283 
284 
286 
287 
289  const std::vector<double>& alphas() const { return _as; }
290  // /// alpha_s value accessor (non-const)
291  // std::vector<double>& alphas() { return _as; }
292  // /// alpha_s value setter
293  // void setalphas(const valarray& xfs) { _as = as; }
294 
296 
297 
299 
300 
302  double ddlogq_forward(size_t i) const {
303  return (alphas()[i+1] - alphas()[i]) / (logq2s()[i+1] - logq2s()[i]);
304  }
305 
307  double ddlogq_backward(size_t i) const {
308  return (alphas()[i] - alphas()[i-1]) / (logq2s()[i] - logq2s()[i-1]);
309  }
310 
312  double ddlogq_central(size_t i) const {
313  return 0.5 * (ddlogq_forward(i) + ddlogq_backward(i));
314  }
315 
317 
318 
319  private:
320 
322  void _synclogs() {
323  _logq2s.resize(_q2s.size());
324  for (size_t i = 0; i < _q2s.size(); ++i) _logq2s[i] = log(_q2s[i]);
325  }
326 
328  std::vector<double> _q2s;
330  std::vector<double> _logq2s;
332  std::vector<double> _as;
333 
334  };
335 
336 
337 }
338 #endif
const std::vector< double > & xs() const
Access the xs array.
Definition: KnotArray.h:203
const size_t xsize() const
Number of x knots.
Definition: KnotArray.h:60
std::string to_str(const T &val)
Make a string representation of val.
Definition: Utils.h:60
const KnotArray1F & get_pid(int id) const
Get the KnotArray1F for PID code id.
Definition: KnotArray.h:183
size_t ixbelow(double x) const
Get the index of the closest x knot column <= x (see KnotArray1F)
Definition: KnotArray.h:207
size_t iq2below(double q2) const
Get the index of the closest Q2 knot row <= q2 (see KnotArray1F)
Definition: KnotArray.h:214
std::vector< double > _logq2s
List of log(Q2) knots.
Definition: KnotArray.h:158
std::vector< double > _logxs
List of log(x) knots.
Definition: KnotArray.h:156
bool empty() const
Is this container empty?
Definition: KnotArray.h:175
KnotArray1F()
Default constructor just for std::map insertability.
Definition: KnotArray.h:26
Internal storage class for alpha_s interpolation grids.
Definition: KnotArray.h:226
const std::vector< double > & logxs() const
Access the log(x)s array.
Definition: KnotArray.h:205
Error for requests for unsupported/invalid flavour PIDs.
Definition: Exceptions.h:70
KnotArray1F & operator[](int id)
Indexing operator (non-const)
Definition: KnotArray.h:200
const std::vector< double > & q2s() const
Q2 knot accessor.
Definition: KnotArray.h:100
size_t size() const
How many {KnotArray1F}s are stored in this container?
Definition: KnotArray.h:172
void _synclogs()
Synchronise the log(Q2) array from the Q2 one.
Definition: KnotArray.h:322
const std::vector< double > & logq2s() const
log(Q2) knot accessor
Definition: KnotArray.h:103
void _synclogs()
Synchronise log(x) and log(Q2) arrays from the x and Q2 ones.
Definition: KnotArray.h:144
void set_pid(int id, const KnotArray1F &ka)
Get the KnotArray1F for PID code id.
Definition: KnotArray.h:195
std::vector< double > _q2s
List of Q2 knots.
Definition: KnotArray.h:328
std::vector< double > _logq2s
List of log(Q2) knots.
Definition: KnotArray.h:330
A collection of {KnotArray1F}s accessed by PID code.
Definition: KnotArray.h:168
void setq2s(const std::vector< double > &q2s)
Definition: KnotArray.h:90
const std::vector< double > & xfs() const
xf value accessor (const)
Definition: KnotArray.h:129
std::vector< double > _xs
List of x knots.
Definition: KnotArray.h:152
std::vector< double > _xfs
List of xf values across the 2D knot array, stored as a strided [ix][iQ2] 1D array.
Definition: KnotArray.h:160
const std::vector< double > & logq2s() const
log(Q2) knot vector accessor
Definition: KnotArray.h:252
const std::vector< double > & alphas() const
alpha_s value accessor (const)
Definition: KnotArray.h:289
const std::vector< double > & q2s() const
Access the Q2s array.
Definition: KnotArray.h:210
Error for general AlphaS computation problems.
Definition: Exceptions.h:94
std::vector< double > _q2s
List of Q2 knots.
Definition: KnotArray.h:154
void setxfs(const std::vector< double > &xfs)
xf value setter
Definition: KnotArray.h:133
std::vector< double > & xfs()
xf value accessor (non-const)
Definition: KnotArray.h:131
double ddlogq_forward(size_t i) const
Forward derivative w.r.t. logQ2.
Definition: KnotArray.h:302
Error for general PDF grid problems.
Definition: Exceptions.h:30
size_t ilogq2below(double logq2) const
Definition: KnotArray.h:271
const std::vector< double > & xs() const
x knot accessor
Definition: KnotArray.h:63
size_t iq2below(double q2) const
Definition: KnotArray.h:108
double ddlogq_backward(size_t i) const
Backward derivative w.r.t. logQ2.
Definition: KnotArray.h:307
Namespace for all LHAPDF functions and classes.
Definition: AlphaS.h:14
std::map< int, KnotArray1F > _map
Storage.
Definition: KnotArray.h:219
const std::vector< double > & logq2s() const
Access the log(Q2)s array.
Definition: KnotArray.h:212
void setxs(const std::vector< double > &xs)
Definition: KnotArray.h:53
size_t iq2below(double q2) const
Definition: KnotArray.h:257
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:29
const std::vector< double > & logxs() const
log(x) knot accessor
Definition: KnotArray.h:66
const KnotArray1F & get_first() const
Convenience accessor for any valid subgrid, to get access to the x/Q2/etc. arrays.
Definition: KnotArray.h:189
size_t ixbelow(double x) const
Get the index of the closest x knot row <= x.
Definition: KnotArray.h:71
double ddlogq_central(size_t i) const
Central (avg of forward and backward) derivative w.r.t. logQ2.
Definition: KnotArray.h:312
const double & xf(size_t ix, size_t iq2) const
Get the xf value at a particular indexed x,Q2 knot.
Definition: KnotArray.h:136
const size_t q2size() const
Number of Q2 knots.
Definition: KnotArray.h:97
Internal storage class for PDF data point grids.
Definition: KnotArray.h:19
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:37
std::vector< double > _as
List of alpha_s values across the knot array.
Definition: KnotArray.h:332
AlphaSArray()
Default constructor just for std::map insertability.
Definition: KnotArray.h:233
bool has_pid(int id) const
Does this contain a KnotArray1F for PID code id?
Definition: KnotArray.h:178
const std::vector< double > & q2s() const
Q2 knot vector accessor.
Definition: KnotArray.h:249
const size_t size() const
Number of x knots.
Definition: KnotArray.h:126
AlphaSArray(const std::vector< double > &q2knots, const std::vector< double > &as)
Constructor from Q2 knot values and alpha_s values.
Definition: KnotArray.h:236