|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
|
#include <algorithm>
#include "FWCore/Utilities/interface/Exception.h"
#include "CondFormats/HcalObjects/interface/HcalCubicInterpolator.h"
HcalCubicInterpolator::HcalCubicInterpolator() {}
HcalCubicInterpolator::HcalCubicInterpolator(const std::vector<Triple>& points) {
const std::size_t sz = points.size();
if (sz) {
std::vector<Triple> tmp(points);
std::sort(tmp.begin(), tmp.end());
abscissae_.reserve(sz);
values_.reserve(sz);
derivatives_.reserve(sz);
for (std::size_t i = 0; i < sz; ++i) {
const Triple& t(tmp[i]);
abscissae_.push_back(std::get<0>(t));
values_.push_back(std::get<1>(t));
derivatives_.push_back(std::get<2>(t));
}
const std::size_t szm1 = sz - 1;
for (std::size_t i = 0; i < szm1; ++i)
if (abscissae_[i] == abscissae_[i + 1])
throw cms::Exception(
"In HcalCubicInterpolator constructor:"
" abscissae must not coincide");
}
}
double HcalCubicInterpolator::operator()(const double x) const {
double result = 0.0;
const std::size_t sz = abscissae_.size();
if (sz) {
if (sz > 1) {
const std::size_t szm1 = sz - 1;
if (x >= abscissae_[szm1])
result = values_[szm1] + derivatives_[szm1] * (x - abscissae_[szm1]);
else if (x <= abscissae_[0])
result = values_[0] + derivatives_[0] * (x - abscissae_[0]);
else {
const std::size_t cell = std::upper_bound(abscissae_.begin(), abscissae_.end(), x) - abscissae_.begin() - 1;
const std::size_t cellp1 = cell + 1;
const double dx = abscissae_[cellp1] - abscissae_[cell];
const double t = (x - abscissae_[cell]) / dx;
const double onemt = 1.0 - t;
const double h00 = onemt * onemt * (1.0 + 2.0 * t);
const double h10 = onemt * onemt * t;
const double h01 = t * t * (3.0 - 2.0 * t);
const double h11 = t * t * onemt;
result = h00 * values_[cell] + h10 * dx * derivatives_[cell] + h01 * values_[cellp1] -
h11 * dx * derivatives_[cellp1];
}
} else
result = values_[0] + derivatives_[0] * (x - abscissae_[0]);
}
return result;
}
double HcalCubicInterpolator::xmin() const {
double result = 0.0;
if (!abscissae_.empty())
result = abscissae_[0];
return result;
}
double HcalCubicInterpolator::xmax() const {
double result = 0.0;
if (!abscissae_.empty())
result = abscissae_.back();
return result;
}
HcalCubicInterpolator HcalCubicInterpolator::approximateInverse() const {
const bool monotonous =
isStrictlyIncreasing(values_.begin(), values_.end()) || isStrictlyDecreasing(values_.begin(), values_.end());
if (!monotonous)
throw cms::Exception(
"In HcalCubicInterpolator::inverse:"
" can't invert non-monotonous functor");
const std::size_t sz = abscissae_.size();
std::vector<Triple> points;
points.reserve(sz);
for (std::size_t i = 0; i < sz; ++i) {
const double dydx = derivatives_[i];
if (dydx == 0.0)
throw cms::Exception(
"In HcalCubicInterpolator::inverse:"
" can't invert functor with derivatives of 0");
points.push_back(Triple(values_[i], abscissae_[i], 1.0 / dydx));
}
return HcalCubicInterpolator(points);
}
BOOST_CLASS_EXPORT_IMPLEMENT(HcalCubicInterpolator)
|