LHAPDF is hosted by Hepforge, IPPP Durham
LHAPDF  6.1.6
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 #include "boost/multi_array.hpp"
12 
13 namespace LHAPDF {
14 
15 
20  class KnotArray1F {
21  public:
22 
24  typedef boost::multi_array<double, 2> valarray;
25 
26 
28 
29 
32 
34  KnotArray1F(const std::vector<double>& xknots, const std::vector<double>& q2knots, const valarray& xfs)
35  : _xs(xknots), _q2s(q2knots), _xfs(xfs)
36  {
37  assert(_xfs.shape()[0] == xknots.size() && _xfs.shape()[1] == q2knots.size());
38  _sync();
39  }
40 
42  KnotArray1F(const std::vector<double>& xknots, const std::vector<double>& q2knots)
43  : _xs(xknots), _q2s(q2knots),
44  _xfs(boost::extents[xknots.size()][q2knots.size()])
45  {
46  assert(_xfs.shape()[0] == xknots.size() && _xfs.shape()[1] == q2knots.size());
47  _sync();
48  }
49 
53  KnotArray1F(const KnotArray1F& other)
54  : _xs(other._xs), _q2s(other._q2s),
55  _logxs(other._logxs), _logq2s(other._logq2s)
56  {
57  _xfs.resize(boost::extents[other._xfs.shape()[0]][other._xfs.shape()[1]]);
58  _xfs = other._xfs;
59  assert(_xfs.shape()[0] == _xs.size() && _xfs.shape()[1] == _q2s.size());
60  }
61 
64  _xs = other._xs;
65  _q2s = other._q2s;
66  _logxs = other._logxs;
67  _logq2s = other._logq2s;
68  _xfs.resize(boost::extents[other._xfs.shape()[0]][other._xfs.shape()[1]]);
69  _xfs = other._xfs;
70  assert(_xfs.shape()[0] == _xs.size() && _xfs.shape()[1] == _q2s.size());
71  return *this;
72  }
73 
75 
76 
78 
79 
81  void setxs(const std::vector<double>& xs) { _xs = xs; _xfs.resize(boost::extents[_xs.size()][_q2s.size()]); }
82 
84  const std::vector<double>& xs() const { return _xs; }
85 
87  const std::vector<double>& logxs() const { return _logxs; }
88 
92  size_t ixbelow(double x) const {
93  // Test that x is in the grid range
94  if (x < xs().front()) throw GridError("x value " + to_str(x) + " is lower than lowest-x grid point at " + to_str(xs().front()));
95  if (x > xs().back()) throw GridError("x value " + to_str(x) + " is higher than highest-x grid point at " + to_str(xs().back()));
96  // Find the closest knot below the requested value
97  size_t i = upper_bound(xs().begin(), xs().end(), x) - xs().begin();
98  if (i == xs().size()) i -= 1; // can't return the last knot index
99  i -= 1; // have to step back to get the knot <= x behaviour
100  return i;
101  }
102 
104 
105 
107 
108 
110  void setq2s(const std::vector<double>& q2s) { _q2s = q2s; _xfs.resize(boost::extents[_xs.size()][_q2s.size()]); }
111 
113  const std::vector<double>& q2s() const { return _q2s; }
114 
116  const std::vector<double>& logq2s() const { return _logq2s; }
117 
121  size_t iq2below(double q2) const {
122  // Test that Q2 is in the grid range
123  if (q2 < q2s().front()) throw GridError("Q2 value " + to_str(q2) + " is lower than lowest-Q2 grid point at " + to_str(q2s().front()));
124  if (q2 > q2s().back()) throw GridError("Q2 value " + to_str(q2) + " is higher than highest-Q2 grid point at " + to_str(q2s().back()));
126  size_t i = upper_bound(q2s().begin(), q2s().end(), q2) - q2s().begin();
127  if (i == q2s().size()) i -= 1; // can't return the last knot index
128  i -= 1; // have to step back to get the knot <= q2 behaviour
129  return i;
130  }
131 
133 
134 
136 
137 
139  const valarray& xfs() const { return _xfs; }
141  valarray& xfs() { return _xfs; }
143  void setxfs(const valarray& xfs) { _xfs = xfs; }
144 
147  const double& xf(size_t ix, size_t iq2) const { return _xfs[ix][iq2]; }
148 
150 
151 
152  private:
153 
155  void _sync() {
156  _logxs.resize(_xs.size());
157  _logq2s.resize(_q2s.size());
158  for (size_t i = 0; i < _xs.size(); ++i) _logxs[i] = log(_xs[i]);
159  for (size_t i = 0; i < _q2s.size(); ++i) _logq2s[i] = log(_q2s[i]);
160  }
161 
163  std::vector<double> _xs;
165  std::vector<double> _q2s;
167  std::vector<double> _logxs;
169  std::vector<double> _logq2s;
171  valarray _xfs;
172 
173  };
174 
175 
179  class KnotArrayNF {
180  public:
181 
183  size_t size() const { return _map.size(); }
184 
186  bool empty() const { return _map.empty(); }
187 
189  bool has_pid(int id) const {
190  return _map.find(id) != _map.end();
191  }
192 
194  const KnotArray1F& get_pid(int id) const {
195  if (!has_pid(id)) throw FlavorError("Undefined particle ID requested: " + to_str(id));
196  return _map.find(id)->second;
197  }
198 
200  const KnotArray1F& get_first() const {
201  if (empty()) throw GridError("Tried to access grid indices when no flavour grids were loaded");
202  return _map.begin()->second;
203  }
204 
206  void set_pid(int id, const KnotArray1F& ka) {
207  _map[id] = ka;
208  }
209 
211  KnotArray1F& operator[](int id) { return _map[id]; }
212 
214  const std::vector<double>& xs() const { return get_first().xs(); }
216  const std::vector<double>& logxs() const { return get_first().logxs(); }
218  size_t ixbelow(double x) const { return get_first().ixbelow(x); }
219 
221  const std::vector<double>& q2s() const { return get_first().q2s(); }
223  const std::vector<double>& logq2s() const { return get_first().logq2s(); }
225  size_t iq2below(double q2) const { return get_first().iq2below(q2); }
226 
227  private:
228 
230  std::map<int, KnotArray1F> _map;
231 
232  };
233 
234 
235 
237  class AlphaSArray {
238  public:
239 
241 
242 
245 
247  AlphaSArray(const std::vector<double>& q2knots, const std::vector<double>& as)
248  : _q2s(q2knots), _as(as)
249  {
250  _sync();
251  }
252 
254 
255 
257 
258 
259  // /// Q2 knot setter
260  // void setq2s(const std::vector<double>& q2s) { _q2s = q2s; _xfs.resize(boost::extents[_xs.size()][_q2s.size()]); }
261 
263  const std::vector<double>& q2s() const { return _q2s; }
264 
265  // /// Get the Q2 value at a particular indexed Q2 knot
266  // const double& q2(size_t iq2) const { return _q2s[iq2]; }
267 
269  const std::vector<double>& logq2s() const { return _logq2s; }
270 
271  // /// Get the log(Q2) value at a particular indexed Q2 knot
272  // const double& logq2(size_t iq2) const { return _logq2s[iq2]; }
273 
277  size_t iq2below(double q2) const {
278  // Test that Q2 is in the grid range
279  if (q2 < q2s().front()) throw AlphaSError("Q2 value " + to_str(q2) + " is lower than lowest-Q2 grid point at " + to_str(q2s().front()));
280  if (q2 > q2s().back()) throw AlphaSError("Q2 value " + to_str(q2) + " is higher than highest-Q2 grid point at " + to_str(q2s().back()));
282  size_t i = upper_bound(q2s().begin(), q2s().end(), q2) - q2s().begin();
283  if (i == q2s().size()) i -= 1; // can't return the last knot index
284  i -= 1; // have to step back to get the knot <= q2 behaviour
285  return i;
286  }
287 
291  size_t ilogq2below(double logq2) const {
292  // Test that log(Q2) is in the grid range
293  if (logq2 < logq2s().front()) throw GridError("logQ2 value " + to_str(logq2) + " is lower than lowest-logQ2 grid point at " + to_str(logq2s().front()));
294  if (logq2 > logq2s().back()) throw GridError("logQ2 value " + to_str(logq2) + " is higher than highest-logQ2 grid point at " + to_str(logq2s().back()));
296  size_t i = upper_bound(logq2s().begin(), logq2s().end(), logq2) - logq2s().begin();
297  if (i == logq2s().size()) i -= 1; // can't return the last knot index
298  i -= 1; // have to step back to get the knot <= q2 behaviour
299  return i;
300  }
301 
303 
304 
306 
307 
309  const std::vector<double>& alphas() const { return _as; }
310  // /// alpha_s value accessor (non-const)
311  // std::vector<double>& alphas() { return _as; }
312  // /// alpha_s value setter
313  // void setalphas(const valarray& xfs) { _as = as; }
314 
315  // /// Get the alpha_s value at a particular indexed Q2 knot
316  // const double& alpha(size_t iq2) const { return _as[iq2]; }
317 
319 
320 
322 
323 
325  double ddlogq_forward(size_t i) const {
326  return (alphas()[i+1] - alphas()[i]) / (logq2s()[i+1] - logq2s()[i]);
327  }
328 
330  double ddlogq_backward(size_t i) const {
331  return (alphas()[i] - alphas()[i-1]) / (logq2s()[i] - logq2s()[i-1]);
332  }
333 
335  double ddlogq_central(size_t i) const {
336  return 0.5 * (ddlogq_forward(i) + ddlogq_backward(i));
337  }
338 
340 
341 
342  private:
343 
345  void _sync() {
346  _logq2s.resize(_q2s.size());
347  for (size_t i = 0; i < _q2s.size(); ++i) _logq2s[i] = log(_q2s[i]);
348  }
349 
351  std::vector<double> _q2s;
353  std::vector<double> _logq2s;
355  std::vector<double> _as;
356 
357  };
358 
359 
360 }
361 #endif
size_t ixbelow(double x) const
Get the index of the closest x knot row <= x.
Definition: KnotArray.h:92
const std::vector< double > & logxs() const
Access the log(x)s array.
Definition: KnotArray.h:216
std::string to_str(const T &val)
Make a string representation of val.
Definition: Utils.h:61
KnotArray1F(const KnotArray1F &other)
Definition: KnotArray.h:53
size_t ixbelow(double x) const
Get the index of the closest x knot column <= x (see KnotArray1F)
Definition: KnotArray.h:218
std::vector< double > _logq2s
List of log(Q2) knots.
Definition: KnotArray.h:169
std::vector< double > _logxs
List of log(x) knots.
Definition: KnotArray.h:167
KnotArray1F()
Default constructor just for std::map insertability.
Definition: KnotArray.h:31
const std::vector< double > & q2s() const
Q2 knot vector accessor.
Definition: KnotArray.h:263
Internal storage class for alpha_s interpolation grids.
Definition: KnotArray.h:237
const std::vector< double > & logxs() const
log(x) knot accessor
Definition: KnotArray.h:87
Error for requests for unsupported/invalid flavour PIDs.
Definition: Exceptions.h:70
KnotArray1F & operator[](int id)
Indexing operator (non-const)
Definition: KnotArray.h:211
void set_pid(int id, const KnotArray1F &ka)
Get the KnotArray1F for PID code id.
Definition: KnotArray.h:206
std::vector< double > _q2s
List of Q2 knots.
Definition: KnotArray.h:351
const std::vector< double > & xs() const
x knot accessor
Definition: KnotArray.h:84
std::vector< double > _logq2s
List of log(Q2) knots.
Definition: KnotArray.h:353
A collection of {KnotArray1F}s accessed by PID code.
Definition: KnotArray.h:179
void setq2s(const std::vector< double > &q2s)
Q2 knot setter.
Definition: KnotArray.h:110
valarray _xfs
List of xf values across the knot array.
Definition: KnotArray.h:171
size_t iq2below(double q2) const
Get the index of the closest Q2 knot row <= q2 (see KnotArray1F)
Definition: KnotArray.h:225
size_t iq2below(double q2) const
Definition: KnotArray.h:121
void setxfs(const valarray &xfs)
xf value setter
Definition: KnotArray.h:143
boost::multi_array< double, 2 > valarray
Use the Boost multi_array for efficiency and ease of indexing.
Definition: KnotArray.h:24
std::vector< double > _xs
List of x knots.
Definition: KnotArray.h:163
const std::vector< double > & q2s() const
Q2 knot accessor.
Definition: KnotArray.h:113
size_t ilogq2below(double logq2) const
Definition: KnotArray.h:291
Error for general AlphaS computation problems.
Definition: Exceptions.h:94
std::vector< double > _q2s
List of Q2 knots.
Definition: KnotArray.h:165
const std::vector< double > & q2s() const
Access the Q2s array.
Definition: KnotArray.h:221
const std::vector< double > & logq2s() const
log(Q2) knot accessor
Definition: KnotArray.h:116
size_t iq2below(double q2) const
Definition: KnotArray.h:277
const KnotArray1F & get_first() const
Convenience accessor for any valid subgrid, to get access to the x/Q2/etc. arrays.
Definition: KnotArray.h:200
const std::vector< double > & xs() const
Access the xs array.
Definition: KnotArray.h:214
const valarray & xfs() const
xf value accessor (const)
Definition: KnotArray.h:139
const KnotArray1F & get_pid(int id) const
Get the KnotArray1F for PID code id.
Definition: KnotArray.h:194
Error for general PDF grid problems.
Definition: Exceptions.h:30
valarray & xfs()
xf value accessor (non-const)
Definition: KnotArray.h:141
void _sync()
Synchronise the log(Q2) array from the Q2 one.
Definition: KnotArray.h:345
bool empty() const
Is this container empty?
Definition: KnotArray.h:186
void _sync()
Synchronise log(x) and log(Q2) arrays from the x and Q2 ones.
Definition: KnotArray.h:155
Namespace for all LHAPDF functions and classes.
Definition: AlphaS.h:14
std::map< int, KnotArray1F > _map
Storage.
Definition: KnotArray.h:230
KnotArray1F & operator=(const KnotArray1F &other)
An explicit operator= is needed due to the Boost multi_array copy semantics.
Definition: KnotArray.h:63
const std::vector< double > & logq2s() const
Access the log(Q2)s array.
Definition: KnotArray.h:223
void setxs(const std::vector< double > &xs)
x knot setter
Definition: KnotArray.h:81
double ddlogq_central(size_t i) const
Central (avg of forward and backward) derivative w.r.t. logQ2.
Definition: KnotArray.h:335
const std::vector< double > & logq2s() const
log(Q2) knot vector accessor
Definition: KnotArray.h:269
double ddlogq_backward(size_t i) const
Backward derivative w.r.t. logQ2.
Definition: KnotArray.h:330
bool has_pid(int id) const
Does this contain a KnotArray1F for PID code id?
Definition: KnotArray.h:189
Internal storage class for PDF data point grids.
Definition: KnotArray.h:20
const double & xf(size_t ix, size_t iq2) const
Definition: KnotArray.h:147
KnotArray1F(const std::vector< double > &xknots, const std::vector< double > &q2knots, const valarray &xfs)
Constructor from x and Q2 knot values, and an xf value grid.
Definition: KnotArray.h:34
KnotArray1F(const std::vector< double > &xknots, const std::vector< double > &q2knots)
Constructor from x and Q2 knot values.
Definition: KnotArray.h:42
std::vector< double > _as
List of alpha_s values across the knot array.
Definition: KnotArray.h:355
AlphaSArray()
Default constructor just for std::map insertability.
Definition: KnotArray.h:244
const std::vector< double > & alphas() const
alpha_s value accessor (const)
Definition: KnotArray.h:309
double ddlogq_forward(size_t i) const
Forward derivative w.r.t. logQ2.
Definition: KnotArray.h:325
AlphaSArray(const std::vector< double > &q2knots, const std::vector< double > &as)
Constructor from Q2 knot values and alpha_s values.
Definition: KnotArray.h:247
size_t size() const
How many {KnotArray1F}s are stored in this container?
Definition: KnotArray.h:183