InterpolationTables/src/UniformAxis.cc

0001 #include <cmath>
0002 #include <climits>
0003 #include <algorithm>
0004 #include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"
0005
0006 #include "Alignment/Geners/interface/binaryIO.hh"
0007 #include "Alignment/Geners/interface/IOException.hh"
0008
0009 #include "JetMETCorrections/InterpolationTables/interface/UniformAxis.h"
0010 #include "JetMETCorrections/InterpolationTables/interface/closeWithinTolerance.h"
0011
0012 namespace npstat {
0013   UniformAxis::UniformAxis(const unsigned nCoords, const double min, const double max, const char* label)
0014       : min_(min), max_(max), label_(label ? label : ""), npt_(nCoords) {
0015     if (!(npt_ > 1U && npt_ < UINT_MAX / 2U - 1U))
0016       throw npstat::NpstatInvalidArgument(
0017           "In npstat::UniformAxis constructor: "
0018           "number of points is out of range");
0019     if (min_ > max_)
0020       std::swap(min_, max_);
0021     bw_ = (max_ - min_) / (npt_ - 1U);
0022     if (max_ == min_)
0023       throw npstat::NpstatInvalidArgument(
0024           "In npstat::UniformAxis constructor: "
0025           "minimum and maximum must be distinct");
0026   }
0027
0028   std::pair<unsigned, double> UniformAxis::getInterval(const double x) const {
0029     if (x <= min_)
0030       return std::pair<unsigned, double>(0U, 1.0);
0031     else if (x >= max_)
0032       return std::pair<unsigned, double>(npt_ - 2U, 0.0);
0033     else {
0034       unsigned binnum = static_cast<unsigned>(floor((x - min_) / bw_));
0035       if (binnum > npt_ - 2U)
0036         binnum = npt_ - 2U;
0037       double w = binnum + 1.0 - (x - min_) / bw_;
0038       if (w < 0.0)
0039         w = 0.0;
0040       else if (w > 1.0)
0041         w = 1.0;
0042       return std::pair<unsigned, double>(binnum, w);
0043     }
0044   }
0045
0046   std::pair<unsigned, double> UniformAxis::linearInterval(const double x) const {
0047     if (x <= min_)
0048       return std::pair<unsigned, double>(0U, 1.0 - (x - min_) / bw_);
0049     else if (x >= max_)
0050       return std::pair<unsigned, double>(npt_ - 2U, (max_ - x) / bw_);
0051     else {
0052       unsigned binnum = static_cast<unsigned>(floor((x - min_) / bw_));
0053       if (binnum > npt_ - 2U)
0054         binnum = npt_ - 2U;
0055       double w = binnum + 1.0 - (x - min_) / bw_;
0056       if (w < 0.0)
0057         w = 0.0;
0058       else if (w > 1.0)
0059         w = 1.0;
0060       return std::pair<unsigned, double>(binnum, w);
0061     }
0062   }
0063
0064   std::vector<double> UniformAxis::coords() const {
0065     std::vector<double> vec;
0066     vec.reserve(npt_);
0067     const unsigned nptm1 = npt_ - 1U;
0068     for (unsigned i = 0; i < nptm1; ++i)
0069       vec.push_back(min_ + bw_ * i);
0070     vec.push_back(max_);
0071     return vec;
0072   }
0073
0074   double UniformAxis::coordinate(const unsigned i) const {
0075     if (i >= npt_)
0076       throw npstat::NpstatOutOfRange("In npstat::UniformAxis::coordinate: index out of range");
0077     if (i == npt_ - 1U)
0078       return max_;
0079     else
0080       return min_ + bw_ * i;
0081   }
0082
0083   bool UniformAxis::isClose(const UniformAxis& r, const double tol) const {
0084     return closeWithinTolerance(min_, r.min_, tol) && closeWithinTolerance(max_, r.max_, tol) && label_ == r.label_ &&
0085            npt_ == r.npt_;
0086   }
0087
0088   bool UniformAxis::operator==(const UniformAxis& r) const {
0089     return min_ == r.min_ && max_ == r.max_ && label_ == r.label_ && npt_ == r.npt_;
0090   }
0091
0092   bool UniformAxis::write(std::ostream& of) const {
0093     gs::write_pod(of, min_);
0094     gs::write_pod(of, max_);
0095     gs::write_pod(of, label_);
0096     gs::write_pod(of, npt_);
0097     return !of.fail();
0098   }
0099
0100   UniformAxis* UniformAxis::read(const gs::ClassId& id, std::istream& in) {
0101     static const gs::ClassId current(gs::ClassId::makeId<UniformAxis>());
0102     current.ensureSameId(id);
0103
0104     double min = 0.0, max = 0.0;
0105     std::string label;
0106     unsigned nBins = 0;
0107
0108     gs::read_pod(in, &min);
0109     gs::read_pod(in, &max);
0110     gs::read_pod(in, &label);
0111     gs::read_pod(in, &nBins);
0112
0113     if (!in.fail())
0114       return new UniformAxis(nBins, min, max, label.c_str());
0115     else
0116       throw gs::IOReadFailure(
0117           "In npstat::UniformAxis::read: "
0118           "input stream failure");
0119   }
0120 }  // namespace npstat