HcalPiecewiseLinearFunctor.cc

Go to the documentation of this file.

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

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

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

HcalPiecewiseLinearFunctor::HcalPiecewiseLinearFunctor()
    : leftExtrapolationLinear_(false), rightExtrapolationLinear_(false) {}

HcalPiecewiseLinearFunctor::HcalPiecewiseLinearFunctor(const std::vector<std::pair<double, double> >& points,
                                                       const bool leftExtrapolationLinear,
                                                       const bool rightExtrapolationLinear)
    : leftExtrapolationLinear_(leftExtrapolationLinear), rightExtrapolationLinear_(rightExtrapolationLinear) {
  const std::size_t sz = points.size();
  if (sz) {
    std::vector<std::pair<double, double> > tmp(points);
    std::sort(tmp.begin(), tmp.end());
    abscissae_.reserve(sz);
    values_.reserve(sz);
    for (std::size_t i = 0; i < sz; ++i) {
      abscissae_.push_back(tmp[i].first);
      values_.push_back(tmp[i].second);
    }
    if (!isStrictlyIncreasing(abscissae_.begin(), abscissae_.end()))
      throw cms::Exception(
          "In HcalPiecewiseLinearFunctor constructor:"
          " abscissae must not coincide");
  }
}

double HcalPiecewiseLinearFunctor::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]) {
        if (rightExtrapolationLinear_)
          result = interpolateSimple(abscissae_[sz - 2], abscissae_[szm1], values_[sz - 2], values_[szm1], x);
        else
          result = values_[szm1];
      } else if (x <= abscissae_[0]) {
        if (leftExtrapolationLinear_)
          result = interpolateSimple(abscissae_[0], abscissae_[1], values_[0], values_[1], x);
        else
          result = values_[0];
      } else {
        const std::size_t cell = std::upper_bound(abscissae_.begin(), abscissae_.end(), x) - abscissae_.begin() - 1;
        result = interpolateSimple(abscissae_[cell], abscissae_[cell + 1], values_[cell], values_[cell + 1], x);
      }
    } else
      result = values_[0];
  }
  return result;
}

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

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

double HcalPiecewiseLinearFunctor::xmin() const {
  double result = 0.0;
  if (!abscissae_.empty())
    result = abscissae_[0];
  return result;
}

double HcalPiecewiseLinearFunctor::xmax() const {
  double result = 0.0;
  if (!abscissae_.empty())
    result = abscissae_.back();
  return result;
}

BOOST_CLASS_EXPORT_IMPLEMENT(HcalPiecewiseLinearFunctor)