HcalInterpolatedTableFunctor.cc

Go to the documentation of this file.

#include <algorithm>
#include "FWCore/Utilities/interface/Exception.h"

#include "CondFormats/HcalObjects/interface/HcalInterpolatedTableFunctor.h"

inline static double interpolateStep(
    const double x0, const double step, const double y0, const double y1, const double x) {
  return y0 + (y1 - y0) * ((x - x0) / step);
}

HcalInterpolatedTableFunctor::HcalInterpolatedTableFunctor()
    : xmin_(0.), xmax_(0.), leftExtrapolationLinear_(false), rightExtrapolationLinear_(false) {}

HcalInterpolatedTableFunctor::HcalInterpolatedTableFunctor(const std::vector<double>& values,
                                                           const double ixmin,
                                                           const double ixmax,
                                                           const bool leftExtrapolationLinear,
                                                           const bool rightExtrapolationLinear)
    : values_(values),
      xmin_(ixmin),
      xmax_(ixmax),
      leftExtrapolationLinear_(leftExtrapolationLinear),
      rightExtrapolationLinear_(rightExtrapolationLinear) {
  if (values_.size() < 2)
    throw cms::Exception(
        "In HcalInterpolatedTableFunctor constructor:"
        " insufficient number of points");
  if (xmin_ >= xmax_)
    throw cms::Exception(
        "In HcalInterpolatedTableFunctor constructor:"
        " invalid min and/or max coordinates");
}

double HcalInterpolatedTableFunctor::operator()(const double x) const {
  double result = 0.0;
  const std::size_t sz = values_.size();
  const std::size_t szm1 = sz - 1;
  const double step = (xmax_ - xmin_) / szm1;

  if (x >= xmax_) {
    if (rightExtrapolationLinear_)
      result = interpolateStep(xmax_ - step, step, values_[sz - 2], values_[szm1], x);
    else
      result = values_[szm1];
  } else if (x <= xmin_) {
    if (leftExtrapolationLinear_)
      result = interpolateStep(xmin_, step, values_[0], values_[1], x);
    else
      result = values_[0];
  } else {
    const std::size_t ux = static_cast<std::size_t>((x - xmin_) / step);
    if (ux >= szm1)
      return values_[szm1];
    result = interpolateStep(ux * step + xmin_, step, values_[ux], values_[ux + 1], x);
  }

  return result;
}

bool HcalInterpolatedTableFunctor::isStrictlyMonotonous() const {
  return isStrictlyIncreasing(values_.begin(), values_.end()) || isStrictlyDecreasing(values_.begin(), values_.end());
}

HcalPiecewiseLinearFunctor HcalInterpolatedTableFunctor::inverse() const {
  if (!isStrictlyMonotonous())
    throw cms::Exception(
        "In HcalInterpolatedTableFunctor::inverse:"
        " can't invert non-monotonous functor");
  const std::size_t sz = values_.size();
  const std::size_t szm1 = sz - 1;
  const double step = (xmax_ - xmin_) / szm1;
  std::vector<std::pair<double, double> > points;
  points.reserve(sz);
  for (std::size_t i = 0; i < sz; ++i) {
    const double x = (i == szm1 ? xmax_ : xmin_ + step * i);
    points.push_back(std::make_pair(values_[i], x));
  }
  bool l = leftExtrapolationLinear_;
  bool r = rightExtrapolationLinear_;
  if (values_[0] > values_[sz - 1])
    std::swap(l, r);
  return HcalPiecewiseLinearFunctor(points, l, r);
}

BOOST_CLASS_EXPORT_IMPLEMENT(HcalInterpolatedTableFunctor)