File indexing completed on 2024-04-06 12:19:22
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 }
