File indexing completed on 2024-04-06 12:02:15
0001 #include <algorithm>
0002 #include "FWCore/Utilities/interface/Exception.h"
0003
0004 #include "CondFormats/HcalObjects/interface/HcalCubicInterpolator.h"
0005
0006 HcalCubicInterpolator::HcalCubicInterpolator() {}
0007
0008 HcalCubicInterpolator::HcalCubicInterpolator(const std::vector<Triple>& points) {
0009 const std::size_t sz = points.size();
0010 if (sz) {
0011 std::vector<Triple> tmp(points);
0012 std::sort(tmp.begin(), tmp.end());
0013 abscissae_.reserve(sz);
0014 values_.reserve(sz);
0015 derivatives_.reserve(sz);
0016 for (std::size_t i = 0; i < sz; ++i) {
0017 const Triple& t(tmp[i]);
0018 abscissae_.push_back(std::get<0>(t));
0019 values_.push_back(std::get<1>(t));
0020 derivatives_.push_back(std::get<2>(t));
0021 }
0022 const std::size_t szm1 = sz - 1;
0023 for (std::size_t i = 0; i < szm1; ++i)
0024 if (abscissae_[i] == abscissae_[i + 1])
0025 throw cms::Exception(
0026 "In HcalCubicInterpolator constructor:"
0027 " abscissae must not coincide");
0028 }
0029 }
0030
0031 double HcalCubicInterpolator::operator()(const double x) const {
0032 double result = 0.0;
0033 const std::size_t sz = abscissae_.size();
0034 if (sz) {
0035 if (sz > 1) {
0036 const std::size_t szm1 = sz - 1;
0037 if (x >= abscissae_[szm1])
0038 result = values_[szm1] + derivatives_[szm1] * (x - abscissae_[szm1]);
0039 else if (x <= abscissae_[0])
0040 result = values_[0] + derivatives_[0] * (x - abscissae_[0]);
0041 else {
0042 const std::size_t cell = std::upper_bound(abscissae_.begin(), abscissae_.end(), x) - abscissae_.begin() - 1;
0043 const std::size_t cellp1 = cell + 1;
0044 const double dx = abscissae_[cellp1] - abscissae_[cell];
0045 const double t = (x - abscissae_[cell]) / dx;
0046 const double onemt = 1.0 - t;
0047 const double h00 = onemt * onemt * (1.0 + 2.0 * t);
0048 const double h10 = onemt * onemt * t;
0049 const double h01 = t * t * (3.0 - 2.0 * t);
0050 const double h11 = t * t * onemt;
0051 result = h00 * values_[cell] + h10 * dx * derivatives_[cell] + h01 * values_[cellp1] -
0052 h11 * dx * derivatives_[cellp1];
0053 }
0054 } else
0055 result = values_[0] + derivatives_[0] * (x - abscissae_[0]);
0056 }
0057 return result;
0058 }
0059
0060 double HcalCubicInterpolator::xmin() const {
0061 double result = 0.0;
0062 if (!abscissae_.empty())
0063 result = abscissae_[0];
0064 return result;
0065 }
0066
0067 double HcalCubicInterpolator::xmax() const {
0068 double result = 0.0;
0069 if (!abscissae_.empty())
0070 result = abscissae_.back();
0071 return result;
0072 }
0073
0074 HcalCubicInterpolator HcalCubicInterpolator::approximateInverse() const {
0075 const bool monotonous =
0076 isStrictlyIncreasing(values_.begin(), values_.end()) || isStrictlyDecreasing(values_.begin(), values_.end());
0077 if (!monotonous)
0078 throw cms::Exception(
0079 "In HcalCubicInterpolator::inverse:"
0080 " can't invert non-monotonous functor");
0081 const std::size_t sz = abscissae_.size();
0082 std::vector<Triple> points;
0083 points.reserve(sz);
0084 for (std::size_t i = 0; i < sz; ++i) {
0085 const double dydx = derivatives_[i];
0086 if (dydx == 0.0)
0087 throw cms::Exception(
0088 "In HcalCubicInterpolator::inverse:"
0089 " can't invert functor with derivatives of 0");
0090 points.push_back(Triple(values_[i], abscissae_[i], 1.0 / dydx));
0091 }
0092 return HcalCubicInterpolator(points);
0093 }
0094
0095 BOOST_CLASS_EXPORT_IMPLEMENT(HcalCubicInterpolator)
