InterpolationTables/src/NUHistoAxis.cc

0001 #include <cmath>
0002 #include <cassert>
0003 #include <climits>
0004 #include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"
0005 #include <algorithm>
0006
0007 #include "Alignment/Geners/interface/binaryIO.hh"
0008 #include "Alignment/Geners/interface/IOException.hh"
0009 #include "JetMETCorrections/InterpolationTables/interface/NUHistoAxis.h"
0010
0011 #include "JetMETCorrections/InterpolationTables/interface/EquidistantSequence.h"
0012 #include "JetMETCorrections/InterpolationTables/interface/closeWithinTolerance.h"
0013
0014 namespace npstat {
0015   NUHistoAxis::NUHistoAxis(const std::vector<double>& binEdges, const char* label)
0016       : binEdges_(binEdges), nBins_(binEdges.size() - 1U), uniform_(false) {
0017     if (!(binEdges_.size() > 1U && binEdges_.size() < UINT_MAX / 2U))
0018       throw npstat::NpstatInvalidArgument(
0019           "In npstat::NUHistoAxis constructor: "
0020           "number of bin edges is out of range");
0021     std::sort(binEdges_.begin(), binEdges_.end());
0022     min_ = binEdges_[0];
0023     max_ = binEdges_[nBins_];
0024     if (label)
0025       label_ = std::string(label);
0026   }
0027
0028   NUHistoAxis::NUHistoAxis(const unsigned nBins, const double min, const double max, const char* label)
0029       : min_(min), max_(max), nBins_(nBins), uniform_(true) {
0030     if (!(nBins_ && nBins_ < UINT_MAX / 2U - 1U))
0031       throw npstat::NpstatInvalidArgument(
0032           "In npstat::NUHistoAxis constructor: "
0033           "number of bins is out of range");
0034     if (min_ > max_)
0035       std::swap(min_, max_);
0036     binEdges_ = EquidistantInLinearSpace(min_, max_, nBins + 1U);
0037     if (label)
0038       label_ = std::string(label);
0039   }
0040
0041   NUHistoAxis NUHistoAxis::rebin(const unsigned newBins) const {
0042     return NUHistoAxis(newBins, min_, max_, label_.c_str());
0043   }
0044
0045   bool NUHistoAxis::isClose(const NUHistoAxis& r, const double tol) const {
0046     if (!(closeWithinTolerance(min_, r.min_, tol) && closeWithinTolerance(max_, r.max_, tol) && label_ == r.label_ &&
0047           nBins_ == r.nBins_ && uniform_ == r.uniform_))
0048       return false;
0049     for (unsigned i = 0; i < nBins_; ++i)
0050       if (!closeWithinTolerance(binEdges_[i], r.binEdges_[i], tol))
0051         return false;
0052     return true;
0053   }
0054
0055   bool NUHistoAxis::operator==(const NUHistoAxis& r) const {
0056     return min_ == r.min_ && max_ == r.max_ && label_ == r.label_ && nBins_ == r.nBins_ && binEdges_ == r.binEdges_ &&
0057            uniform_ == r.uniform_;
0058   }
0059
0060   bool NUHistoAxis::operator!=(const NUHistoAxis& r) const { return !(*this == r); }
0061
0062   int NUHistoAxis::binNumber(const double x) const {
0063     const int delta = std::upper_bound(binEdges_.begin(), binEdges_.end(), x) - binEdges_.begin();
0064     return delta - 1;
0065   }
0066
0067   double NUHistoAxis::fltBinNumber(const double x, const bool mapLeftEdgeTo0) const {
0068     const int delta = std::upper_bound(binEdges_.begin(), binEdges_.end(), x) - binEdges_.begin();
0069     const int binnum = delta - 1;
0070
0071     if (binnum < 0) {
0072       const double left = binEdges_[0];
0073       const double right = binEdges_[1];
0074       double bval = (x - left) / (right - left);
0075       if (!mapLeftEdgeTo0)
0076         bval -= 0.5;
0077       if (bval < -1.0)
0078         bval = -1.0;
0079       return bval;
0080     } else if (static_cast<unsigned>(binnum) >= nBins_) {
0081       const double left = binEdges_[nBins_ - 1U];
0082       const double right = binEdges_[nBins_];
0083       double bval = nBins_ - 1U + (x - left) / (right - left);
0084       if (!mapLeftEdgeTo0)
0085         bval -= 0.5;
0086       if (bval > static_cast<double>(nBins_))
0087         bval = nBins_;
0088       return bval;
0089     } else {
0090       const double left = binEdges_[binnum];
0091       const double right = binEdges_[delta];
0092       if (mapLeftEdgeTo0)
0093         return binnum + (x - left) / (right - left);
0094       else {
0095         // It is not obvious what is best to do here.
0096         // The following works to preserve interpolation
0097         // of 0th order. The commented out snippet below
0098         // would instead connect bin centers by straight
0099         // lines during interpolation.
0100         return binnum + (x - left) / (right - left) - 0.5;
0101
0102         // Bin center is mapped to binnum.
0103         // Bin center of the next bin is mapped to binnum + 1.
0104         // Bin center of the previous bin is mapped to binnum - 1.
0105         //
0106         // const double binCenter = (left + right)/2.0;
0107         // if ((binnum == 0 && x <= binCenter) ||
0108         //     (static_cast<unsigned>(binnum) == nBins_ - 1 && x >= binCenter))
0109         //     return binnum + (x - left)/(right - left) - 0.5;
0110         // else if (x <= binCenter)
0111         // {
0112         //     const double otherBinCenter = (left + binEdges_[binnum - 1])/2.0;
0113         //     return binnum + (x - binCenter)/(binCenter - otherBinCenter);
0114         // }
0115         // else
0116         // {
0117         //     const double otherBinCenter = (right + binEdges_[binnum + 2])/2.0;
0118         //     return binnum + (x - binCenter)/(otherBinCenter - binCenter);
0119         // }
0120       }
0121     }
0122   }
0123
0124   unsigned NUHistoAxis::closestValidBin(const double x) const {
0125     const int delta = std::upper_bound(binEdges_.begin(), binEdges_.end(), x) - binEdges_.begin();
0126     int binnum = delta - 1;
0127     if (binnum < 0)
0128       binnum = 0;
0129     else if (static_cast<unsigned>(binnum) >= nBins_)
0130       binnum = nBins_ - 1U;
0131     return binnum;
0132   }
0133
0134   bool NUHistoAxis::write(std::ostream& of) const {
0135     gs::write_pod_vector(of, binEdges_);
0136     gs::write_pod(of, label_);
0137     unsigned char c = uniform_;
0138     gs::write_pod(of, c);
0139     return !of.fail();
0140   }
0141
0142   NUHistoAxis* NUHistoAxis::read(const gs::ClassId& id, std::istream& in) {
0143     static const gs::ClassId current(gs::ClassId::makeId<NUHistoAxis>());
0144     current.ensureSameId(id);
0145
0146     std::vector<double> binEdges;
0147     std::string label;
0148     unsigned char unif;
0149     gs::read_pod_vector(in, &binEdges);
0150     gs::read_pod(in, &label);
0151     gs::read_pod(in, &unif);
0152     if (in.fail())
0153       throw gs::IOReadFailure(
0154           "In npstat::UHistoAxis::read: "
0155           "input stream failure");
0156     NUHistoAxis* result = new NUHistoAxis(binEdges, label.c_str());
0157     result->uniform_ = unif;
0158     return result;
0159   }
0160 }  // namespace npstat