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