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/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 }