FFTJetAlgorithms/src/LinInterpolatedTable1D.cc

0001 #include <memory>
0002
0003 #include "RecoJets/FFTJetAlgorithms/interface/LinInterpolatedTable1D.h"
0004
0005 namespace fftjetcms {
0006   LinInterpolatedTable1D::LinInterpolatedTable1D(const double c)
0007       : data_(2, c),
0008         xmin_(0.0),
0009         xmax_(1.0),
0010         binwidth_(1.0),
0011         npoints_(2U),
0012         leftExtrapolationLinear_(false),
0013         rightExtrapolationLinear_(false),
0014         monotonous_(false),
0015         monotonicityKnown_(true) {}
0016
0017   bool LinInterpolatedTable1D::operator==(const LinInterpolatedTable1D& r) const {
0018     return xmin_ == r.xmin_ && xmax_ == r.xmax_ && binwidth_ == r.binwidth_ && npoints_ == r.npoints_ &&
0019            leftExtrapolationLinear_ == r.leftExtrapolationLinear_ &&
0020            rightExtrapolationLinear_ == r.rightExtrapolationLinear_ && data_ == r.data_;
0021   }
0022
0023   bool LinInterpolatedTable1D::isMonotonous() const {
0024     if (!monotonicityKnown_) {
0025       monotonous_ = true;
0026       const double delta = data_[npoints_ - 1U] - data_[0];
0027       if (delta == 0.0)
0028         monotonous_ = false;
0029       const double sg = delta > 0.0 ? 1.0 : -1.0;
0030       for (unsigned i = 1; i < npoints_ && monotonous_; ++i)
0031         if ((data_[i] - data_[i - 1]) * sg <= 0.0)
0032           monotonous_ = false;
0033       monotonicityKnown_ = true;
0034     }
0035     return monotonous_;
0036   }
0037
0038   std::unique_ptr<LinInterpolatedTable1D> LinInterpolatedTable1D::inverse(const unsigned npoints,
0039                                                                           const bool leftExtrapolationLinear,
0040                                                                           const bool rightExtrapolationLinear) const {
0041     if (!isMonotonous())
0042       return std::unique_ptr<LinInterpolatedTable1D>(nullptr);
0043
0044     std::vector<std::pair<double, double> > points;
0045     points.reserve(npoints_);
0046
0047     if (data_[npoints_ - 1U] > data_[0]) {
0048       points.push_back(std::pair<double, double>(data_[0], xmin_));
0049       for (unsigned i = 1; i < npoints_ - 1U; ++i)
0050         points.push_back(std::pair<double, double>(data_[i], xmin_ + i * binwidth_));
0051       points.push_back(std::pair<double, double>(data_[npoints_ - 1U], xmax_));
0052     } else {
0053       points.push_back(std::pair<double, double>(data_[npoints_ - 1U], xmax_));
0054       for (unsigned i = npoints_ - 2U; i > 0; --i)
0055         points.push_back(std::pair<double, double>(data_[i], xmin_ + i * binwidth_));
0056       points.push_back(std::pair<double, double>(data_[0], xmin_));
0057     }
0058
0059     return std::make_unique<LinInterpolatedTable1D>(points, npoints, leftExtrapolationLinear, rightExtrapolationLinear);
0060   }
0061
0062   double LinInterpolatedTable1D::operator()(const double& x) const {
0063     if (x <= xmin_) {
0064       if (leftExtrapolationLinear_)
0065         return data_[0] + (data_[1] - data_[0]) * ((x - xmin_) / binwidth_);
0066       else
0067         return data_[0];
0068     } else if (x >= xmax_) {
0069       if (rightExtrapolationLinear_)
0070         return data_[npoints_ - 1U] - (data_[npoints_ - 2U] - data_[npoints_ - 1U]) * ((x - xmax_) / binwidth_);
0071       else
0072         return data_[npoints_ - 1U];
0073     } else {
0074       const unsigned ux = static_cast<unsigned>((x - xmin_) / binwidth_);
0075       if (ux >= npoints_ - 1U)
0076         return data_[npoints_ - 1U];
0077       const double delta = x - (ux * binwidth_ + xmin_);
0078       return data_[ux] + (data_[ux + 1U] - data_[ux]) * delta / binwidth_;
0079     }
0080   }
0081 }  // namespace fftjetcms