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/HcalCubicInterpolator.h"
0005 
0006 HcalCubicInterpolator::HcalCubicInterpolator() {}
0007 
0008 HcalCubicInterpolator::HcalCubicInterpolator(const std::vector<Triple>& points) {
0009   const std::size_t sz = points.size();
0010   if (sz) {
0011     std::vector<Triple> tmp(points);
0012     std::sort(tmp.begin(), tmp.end());
0013     abscissae_.reserve(sz);
0014     values_.reserve(sz);
0015     derivatives_.reserve(sz);
0016     for (std::size_t i = 0; i < sz; ++i) {
0017       const Triple& t(tmp[i]);
0018       abscissae_.push_back(std::get<0>(t));
0019       values_.push_back(std::get<1>(t));
0020       derivatives_.push_back(std::get<2>(t));
0021     }
0022     const std::size_t szm1 = sz - 1;
0023     for (std::size_t i = 0; i < szm1; ++i)
0024       if (abscissae_[i] == abscissae_[i + 1])
0025         throw cms::Exception(
0026             "In HcalCubicInterpolator constructor:"
0027             " abscissae must not coincide");
0028   }
0029 }
0030 
0031 double HcalCubicInterpolator::operator()(const double x) const {
0032   double result = 0.0;
0033   const std::size_t sz = abscissae_.size();
0034   if (sz) {
0035     if (sz > 1) {
0036       const std::size_t szm1 = sz - 1;
0037       if (x >= abscissae_[szm1])
0038         result = values_[szm1] + derivatives_[szm1] * (x - abscissae_[szm1]);
0039       else if (x <= abscissae_[0])
0040         result = values_[0] + derivatives_[0] * (x - abscissae_[0]);
0041       else {
0042         const std::size_t cell = std::upper_bound(abscissae_.begin(), abscissae_.end(), x) - abscissae_.begin() - 1;
0043         const std::size_t cellp1 = cell + 1;
0044         const double dx = abscissae_[cellp1] - abscissae_[cell];
0045         const double t = (x - abscissae_[cell]) / dx;
0046         const double onemt = 1.0 - t;
0047         const double h00 = onemt * onemt * (1.0 + 2.0 * t);
0048         const double h10 = onemt * onemt * t;
0049         const double h01 = t * t * (3.0 - 2.0 * t);
0050         const double h11 = t * t * onemt;
0051         result = h00 * values_[cell] + h10 * dx * derivatives_[cell] + h01 * values_[cellp1] -
0052                  h11 * dx * derivatives_[cellp1];
0053       }
0054     } else
0055       result = values_[0] + derivatives_[0] * (x - abscissae_[0]);
0056   }
0057   return result;
0058 }
0059 
0060 double HcalCubicInterpolator::xmin() const {
0061   double result = 0.0;
0062   if (!abscissae_.empty())
0063     result = abscissae_[0];
0064   return result;
0065 }
0066 
0067 double HcalCubicInterpolator::xmax() const {
0068   double result = 0.0;
0069   if (!abscissae_.empty())
0070     result = abscissae_.back();
0071   return result;
0072 }
0073 
0074 HcalCubicInterpolator HcalCubicInterpolator::approximateInverse() const {
0075   const bool monotonous =
0076       isStrictlyIncreasing(values_.begin(), values_.end()) || isStrictlyDecreasing(values_.begin(), values_.end());
0077   if (!monotonous)
0078     throw cms::Exception(
0079         "In HcalCubicInterpolator::inverse:"
0080         " can't invert non-monotonous functor");
0081   const std::size_t sz = abscissae_.size();
0082   std::vector<Triple> points;
0083   points.reserve(sz);
0084   for (std::size_t i = 0; i < sz; ++i) {
0085     const double dydx = derivatives_[i];
0086     if (dydx == 0.0)
0087       throw cms::Exception(
0088           "In HcalCubicInterpolator::inverse:"
0089           " can't invert functor with derivatives of 0");
0090     points.push_back(Triple(values_[i], abscissae_[i], 1.0 / dydx));
0091   }
0092   return HcalCubicInterpolator(points);
0093 }
0094 
0095 BOOST_CLASS_EXPORT_IMPLEMENT(HcalCubicInterpolator)