Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:51:11

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)