HistoAxis.cc

Go to the documentation of this file.

#include <cmath>
#include <climits>
#include <algorithm>
#include "JetMETCorrections/InterpolationTables/interface/NpstatException.h"

#include "Alignment/Geners/interface/binaryIO.hh"
#include "Alignment/Geners/interface/IOException.hh"

#include "JetMETCorrections/InterpolationTables/interface/closeWithinTolerance.h"
#include "JetMETCorrections/InterpolationTables/interface/HistoAxis.h"

namespace npstat {
  HistoAxis::HistoAxis(const unsigned nbins, const double min, const double max, const char* label)
      : min_(min), max_(max), label_(label ? label : ""), nBins_(nbins) {
    if (!(nBins_ && nBins_ < UINT_MAX / 2U - 1U))
      throw npstat::NpstatInvalidArgument(
          "In npstat::HistoAxis constructor: "
          "number of bins is out of range");
    if (min_ > max_)
      std::swap(min_, max_);
    bw_ = (max_ - min_) / nBins_;
  }

  HistoAxis HistoAxis::rebin(const unsigned nbins) const { return HistoAxis(nbins, min_, max_, label_.c_str()); }

  bool HistoAxis::isClose(const HistoAxis& r, const double tol) const {
    return closeWithinTolerance(min_, r.min_, tol) && closeWithinTolerance(max_, r.max_, tol) && label_ == r.label_ &&
           nBins_ == r.nBins_;
  }

  bool HistoAxis::operator==(const HistoAxis& r) const {
    return min_ == r.min_ && max_ == r.max_ && label_ == r.label_ && nBins_ == r.nBins_;
  }

  bool HistoAxis::operator!=(const HistoAxis& r) const { return !(*this == r); }

  int HistoAxis::binNumber(const double x) const {
    if (bw_) {
      int binnum = static_cast<int>(floor((x - min_) / bw_));
      if (x >= max_) {
        if (binnum < static_cast<int>(nBins_))
          binnum = nBins_;
      } else {
        if (binnum >= static_cast<int>(nBins_))
          binnum = nBins_ - 1U;
      }
      return binnum;
    } else {
      if (x < min_)
        return -1;
      else
        return nBins_;
    }
  }

  unsigned HistoAxis::closestValidBin(const double x) const {
    if (x <= min_)
      return 0U;
    else if (bw_ && x < max_) {
      const unsigned binnum = static_cast<unsigned>(floor((x - min_) / bw_));
      if (binnum < nBins_)
        return binnum;
    }
    return nBins_ - 1U;
  }

  LinearMapper1d HistoAxis::binNumberMapper(const bool mapLeftEdgeTo0) const {
    if (!bw_)
      throw npstat::NpstatDomainError(
          "In npstat::HistoAxis::binNumberMapper: "
          "bin width is zero. Mapper can not be constructed.");
    const double base = mapLeftEdgeTo0 ? min_ / bw_ : min_ / bw_ + 0.5;
    return LinearMapper1d(1.0 / bw_, -base);
  }

  CircularMapper1d HistoAxis::kernelScanMapper(const bool doubleRange) const {
    if (!bw_)
      throw npstat::NpstatDomainError(
          "In npstat::HistoAxis::kernelScanMapper: "
          "bin width is zero. Mapper can not be constructed.");
    double range = max_ - min_;
    if (doubleRange)
      range *= 2.0;
    return CircularMapper1d(bw_, 0.0, range);
  }

  unsigned HistoAxis::overflowIndexWeighted(const double x, unsigned* binNumber, double* weight) const {
    if (x < min_)
      return 0U;
    else if (x >= max_)
      return 2U;
    else {
      if (nBins_ <= 1U)
        throw npstat::NpstatInvalidArgument(
            "In npstat::HistoAxis::overflowIndexWeighted: "
            "must have more than one bin");
      const double dbin = (x - min_) / bw_;
      if (dbin <= 0.5) {
        *binNumber = 0;
        *weight = 1.0;
      } else if (dbin >= nBins_ - 0.5) {
        *binNumber = nBins_ - 2;
        *weight = 0.0;
      } else {
        const unsigned bin = static_cast<unsigned>(dbin - 0.5);
        *binNumber = bin >= nBins_ - 1U ? nBins_ - 2U : bin;
        *weight = 1.0 - (dbin - 0.5 - *binNumber);
      }
      return 1U;
    }
  }

  bool HistoAxis::write(std::ostream& of) const {
    gs::write_pod(of, min_);
    gs::write_pod(of, max_);
    gs::write_pod(of, label_);
    gs::write_pod(of, nBins_);
    return !of.fail();
  }

  HistoAxis* HistoAxis::read(const gs::ClassId& id, std::istream& in) {
    static const gs::ClassId current(gs::ClassId::makeId<HistoAxis>());
    current.ensureSameId(id);

    double min = 0.0, max = 0.0;
    std::string label;
    unsigned nBins = 0;

    gs::read_pod(in, &min);
    gs::read_pod(in, &max);
    gs::read_pod(in, &label);
    gs::read_pod(in, &nBins);

    if (!in.fail())
      return new HistoAxis(nBins, min, max, label.c_str());
    else
      throw gs::IOReadFailure(
          "In npstat::HistoAxis::read: "
          "input stream failure");
  }
}  // namespace npstat