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/HcalInterpolatedTableFunctor.h"
0005
0006 inline static double interpolateStep(
0007 const double x0, const double step, const double y0, const double y1, const double x) {
0008 return y0 + (y1 - y0) * ((x - x0) / step);
0009 }
0010
0011 HcalInterpolatedTableFunctor::HcalInterpolatedTableFunctor()
0012 : xmin_(0.), xmax_(0.), leftExtrapolationLinear_(false), rightExtrapolationLinear_(false) {}
0013
0014 HcalInterpolatedTableFunctor::HcalInterpolatedTableFunctor(const std::vector<double>& values,
0015 const double ixmin,
0016 const double ixmax,
0017 const bool leftExtrapolationLinear,
0018 const bool rightExtrapolationLinear)
0019 : values_(values),
0020 xmin_(ixmin),
0021 xmax_(ixmax),
0022 leftExtrapolationLinear_(leftExtrapolationLinear),
0023 rightExtrapolationLinear_(rightExtrapolationLinear) {
0024 if (values_.size() < 2)
0025 throw cms::Exception(
0026 "In HcalInterpolatedTableFunctor constructor:"
0027 " insufficient number of points");
0028 if (xmin_ >= xmax_)
0029 throw cms::Exception(
0030 "In HcalInterpolatedTableFunctor constructor:"
0031 " invalid min and/or max coordinates");
0032 }
0033
0034 double HcalInterpolatedTableFunctor::operator()(const double x) const {
0035 double result = 0.0;
0036 const std::size_t sz = values_.size();
0037 const std::size_t szm1 = sz - 1;
0038 const double step = (xmax_ - xmin_) / szm1;
0039
0040 if (x >= xmax_) {
0041 if (rightExtrapolationLinear_)
0042 result = interpolateStep(xmax_ - step, step, values_[sz - 2], values_[szm1], x);
0043 else
0044 result = values_[szm1];
0045 } else if (x <= xmin_) {
0046 if (leftExtrapolationLinear_)
0047 result = interpolateStep(xmin_, step, values_[0], values_[1], x);
0048 else
0049 result = values_[0];
0050 } else {
0051 const std::size_t ux = static_cast<std::size_t>((x - xmin_) / step);
0052 if (ux >= szm1)
0053 return values_[szm1];
0054 result = interpolateStep(ux * step + xmin_, step, values_[ux], values_[ux + 1], x);
0055 }
0056
0057 return result;
0058 }
0059
0060 bool HcalInterpolatedTableFunctor::isStrictlyMonotonous() const {
0061 return isStrictlyIncreasing(values_.begin(), values_.end()) || isStrictlyDecreasing(values_.begin(), values_.end());
0062 }
0063
0064 HcalPiecewiseLinearFunctor HcalInterpolatedTableFunctor::inverse() const {
0065 if (!isStrictlyMonotonous())
0066 throw cms::Exception(
0067 "In HcalInterpolatedTableFunctor::inverse:"
0068 " can't invert non-monotonous functor");
0069 const std::size_t sz = values_.size();
0070 const std::size_t szm1 = sz - 1;
0071 const double step = (xmax_ - xmin_) / szm1;
0072 std::vector<std::pair<double, double> > points;
0073 points.reserve(sz);
0074 for (std::size_t i = 0; i < sz; ++i) {
0075 const double x = (i == szm1 ? xmax_ : xmin_ + step * i);
0076 points.push_back(std::make_pair(values_[i], x));
0077 }
0078 bool l = leftExtrapolationLinear_;
0079 bool r = rightExtrapolationLinear_;
0080 if (values_[0] > values_[sz - 1])
0081 std::swap(l, r);
0082 return HcalPiecewiseLinearFunctor(points, l, r);
0083 }
0084
0085 BOOST_CLASS_EXPORT_IMPLEMENT(HcalInterpolatedTableFunctor)
