InterpolationTables/src/HistoAxis.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/closeWithinTolerance.h"
0010 #include "JetMETCorrections/InterpolationTables/interface/HistoAxis.h"
0011
0012 namespace npstat {
0013   HistoAxis::HistoAxis(const unsigned nbins, const double min, const double max, const char* label)
0014       : min_(min), max_(max), label_(label ? label : ""), nBins_(nbins) {
0015     if (!(nBins_ && nBins_ < UINT_MAX / 2U - 1U))
0016       throw npstat::NpstatInvalidArgument(
0017           "In npstat::HistoAxis constructor: "
0018           "number of bins is out of range");
0019     if (min_ > max_)
0020       std::swap(min_, max_);
0021     bw_ = (max_ - min_) / nBins_;
0022   }
0023
0024   HistoAxis HistoAxis::rebin(const unsigned nbins) const { return HistoAxis(nbins, min_, max_, label_.c_str()); }
0025
0026   bool HistoAxis::isClose(const HistoAxis& r, const double tol) const {
0027     return closeWithinTolerance(min_, r.min_, tol) && closeWithinTolerance(max_, r.max_, tol) && label_ == r.label_ &&
0028            nBins_ == r.nBins_;
0029   }
0030
0031   bool HistoAxis::operator==(const HistoAxis& r) const {
0032     return min_ == r.min_ && max_ == r.max_ && label_ == r.label_ && nBins_ == r.nBins_;
0033   }
0034
0035   bool HistoAxis::operator!=(const HistoAxis& r) const { return !(*this == r); }
0036
0037   int HistoAxis::binNumber(const double x) const {
0038     if (bw_) {
0039       int binnum = static_cast<int>(floor((x - min_) / bw_));
0040       if (x >= max_) {
0041         if (binnum < static_cast<int>(nBins_))
0042           binnum = nBins_;
0043       } else {
0044         if (binnum >= static_cast<int>(nBins_))
0045           binnum = nBins_ - 1U;
0046       }
0047       return binnum;
0048     } else {
0049       if (x < min_)
0050         return -1;
0051       else
0052         return nBins_;
0053     }
0054   }
0055
0056   unsigned HistoAxis::closestValidBin(const double x) const {
0057     if (x <= min_)
0058       return 0U;
0059     else if (bw_ && x < max_) {
0060       const unsigned binnum = static_cast<unsigned>(floor((x - min_) / bw_));
0061       if (binnum < nBins_)
0062         return binnum;
0063     }
0064     return nBins_ - 1U;
0065   }
0066
0067   LinearMapper1d HistoAxis::binNumberMapper(const bool mapLeftEdgeTo0) const {
0068     if (!bw_)
0069       throw npstat::NpstatDomainError(
0070           "In npstat::HistoAxis::binNumberMapper: "
0071           "bin width is zero. Mapper can not be constructed.");
0072     const double base = mapLeftEdgeTo0 ? min_ / bw_ : min_ / bw_ + 0.5;
0073     return LinearMapper1d(1.0 / bw_, -base);
0074   }
0075
0076   CircularMapper1d HistoAxis::kernelScanMapper(const bool doubleRange) const {
0077     if (!bw_)
0078       throw npstat::NpstatDomainError(
0079           "In npstat::HistoAxis::kernelScanMapper: "
0080           "bin width is zero. Mapper can not be constructed.");
0081     double range = max_ - min_;
0082     if (doubleRange)
0083       range *= 2.0;
0084     return CircularMapper1d(bw_, 0.0, range);
0085   }
0086
0087   unsigned HistoAxis::overflowIndexWeighted(const double x, unsigned* binNumber, double* weight) const {
0088     if (x < min_)
0089       return 0U;
0090     else if (x >= max_)
0091       return 2U;
0092     else {
0093       if (nBins_ <= 1U)
0094         throw npstat::NpstatInvalidArgument(
0095             "In npstat::HistoAxis::overflowIndexWeighted: "
0096             "must have more than one bin");
0097       const double dbin = (x - min_) / bw_;
0098       if (dbin <= 0.5) {
0099         *binNumber = 0;
0100         *weight = 1.0;
0101       } else if (dbin >= nBins_ - 0.5) {
0102         *binNumber = nBins_ - 2;
0103         *weight = 0.0;
0104       } else {
0105         const unsigned bin = static_cast<unsigned>(dbin - 0.5);
0106         *binNumber = bin >= nBins_ - 1U ? nBins_ - 2U : bin;
0107         *weight = 1.0 - (dbin - 0.5 - *binNumber);
0108       }
0109       return 1U;
0110     }
0111   }
0112
0113   bool HistoAxis::write(std::ostream& of) const {
0114     gs::write_pod(of, min_);
0115     gs::write_pod(of, max_);
0116     gs::write_pod(of, label_);
0117     gs::write_pod(of, nBins_);
0118     return !of.fail();
0119   }
0120
0121   HistoAxis* HistoAxis::read(const gs::ClassId& id, std::istream& in) {
0122     static const gs::ClassId current(gs::ClassId::makeId<HistoAxis>());
0123     current.ensureSameId(id);
0124
0125     double min = 0.0, max = 0.0;
0126     std::string label;
0127     unsigned nBins = 0;
0128
0129     gs::read_pod(in, &min);
0130     gs::read_pod(in, &max);
0131     gs::read_pod(in, &label);
0132     gs::read_pod(in, &nBins);
0133
0134     if (!in.fail())
0135       return new HistoAxis(nBins, min, max, label.c_str());
0136     else
0137       throw gs::IOReadFailure(
0138           "In npstat::HistoAxis::read: "
0139           "input stream failure");
0140   }
0141 }  // namespace npstat