NUHistoAxis.cc

Go to the documentation of this file.

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

#include "Alignment/Geners/interface/binaryIO.hh"
#include "Alignment/Geners/interface/IOException.hh"
#include "JetMETCorrections/InterpolationTables/interface/NUHistoAxis.h"

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

namespace npstat {
    NUHistoAxis::NUHistoAxis(const std::vector<double>& binEdges,
                             const char* label)
        : binEdges_(binEdges), nBins_(binEdges.size() - 1U), uniform_(false)
    {
        if (!(binEdges_.size() > 1U && binEdges_.size() < UINT_MAX/2U))
            throw npstat::NpstatInvalidArgument("In npstat::NUHistoAxis constructor: "
                                        "number of bin edges is out of range");
        std::sort(binEdges_.begin(), binEdges_.end());
        min_ = binEdges_[0];
        max_ = binEdges_[nBins_];
        if (label)
            label_ = std::string(label);
    }

    NUHistoAxis::NUHistoAxis(const unsigned nBins,
                             const double min, const double max,
                             const char* label)
        : min_(min), max_(max), nBins_(nBins), uniform_(true)
    {
        if (!(nBins_ && nBins_ < UINT_MAX/2U - 1U))
            throw npstat::NpstatInvalidArgument("In npstat::NUHistoAxis constructor: "
                                        "number of bins is out of range");
        if (min_ > max_)
            std::swap(min_, max_);
        binEdges_ = EquidistantInLinearSpace(min_, max_, nBins+1U);
        if (label)
            label_ = std::string(label);
    }

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

    bool NUHistoAxis::isClose(const NUHistoAxis& r, const double tol) const
    {
        if (!(closeWithinTolerance(min_, r.min_, tol) && 
              closeWithinTolerance(max_, r.max_, tol) &&
              label_ == r.label_ &&
              nBins_ == r.nBins_ &&
              uniform_ == r.uniform_))
            return false;
        for (unsigned i=0; i<nBins_; ++i)
            if (!closeWithinTolerance(binEdges_[i], r.binEdges_[i], tol))
                return false;
        return true;
    }

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

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

    int NUHistoAxis::binNumber(const double x) const
    {
        const int delta = std::upper_bound(binEdges_.begin(), binEdges_.end(), x) - 
                          binEdges_.begin();
        return delta - 1;
    }

    double NUHistoAxis::fltBinNumber(const double x, const bool mapLeftEdgeTo0) const
    {
        const int delta = std::upper_bound(binEdges_.begin(), binEdges_.end(), x) - 
                          binEdges_.begin();
        const int binnum = delta - 1;

        if (binnum < 0)
        {
            const double left = binEdges_[0];
            const double right = binEdges_[1];
            double bval = (x - left)/(right - left);
            if (!mapLeftEdgeTo0)
                bval -= 0.5;
            if (bval < -1.0)
                bval = -1.0;
            return bval;
        }
        else if (static_cast<unsigned>(binnum) >= nBins_)
        {
            const double left = binEdges_[nBins_ - 1U];
            const double right = binEdges_[nBins_];
            double bval = nBins_ - 1U + (x - left)/(right - left);
            if (!mapLeftEdgeTo0)
                bval -= 0.5;
            if (bval > static_cast<double>(nBins_))
                bval = nBins_;
            return bval;
        }
        else
        {
            const double left = binEdges_[binnum];
            const double right = binEdges_[delta];
            if (mapLeftEdgeTo0)
                return binnum + (x - left)/(right - left);
            else
            {
                // It is not obvious what is best to do here.
                // The following works to preserve interpolation
                // of 0th order. The commented out snippet below
                // would instead connect bin centers by straight
                // lines during interpolation.
                return binnum + (x - left)/(right - left) - 0.5;

                // Bin center is mapped to binnum.
                // Bin center of the next bin is mapped to binnum + 1.
                // Bin center of the previous bin is mapped to binnum - 1.
                //
                // const double binCenter = (left + right)/2.0;
                // if ((binnum == 0 && x <= binCenter) ||
                //     (static_cast<unsigned>(binnum) == nBins_ - 1 && x >= binCenter))
                //     return binnum + (x - left)/(right - left) - 0.5;
                // else if (x <= binCenter)
                // {
                //     const double otherBinCenter = (left + binEdges_[binnum - 1])/2.0;
                //     return binnum + (x - binCenter)/(binCenter - otherBinCenter);
                // }
                // else
                // {
                //     const double otherBinCenter = (right + binEdges_[binnum + 2])/2.0;
                //     return binnum + (x - binCenter)/(otherBinCenter - binCenter);
                // }
            }
        }
    }

    unsigned NUHistoAxis::closestValidBin(const double x) const
    {
        const int delta = std::upper_bound(binEdges_.begin(), binEdges_.end(), x) - 
                          binEdges_.begin();
        int binnum = delta - 1;
        if (binnum < 0)
            binnum = 0;
        else if (static_cast<unsigned>(binnum) >= nBins_)
            binnum = nBins_ - 1U;
        return binnum;
    }

    bool NUHistoAxis::write(std::ostream& of) const
    {
        gs::write_pod_vector(of, binEdges_);
        gs::write_pod(of, label_);
        unsigned char c = uniform_;
        gs::write_pod(of, c);
        return !of.fail();
    }

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

        std::vector<double> binEdges;
        std::string label;
        unsigned char unif;
        gs::read_pod_vector(in, &binEdges);
        gs::read_pod(in, &label);
        gs::read_pod(in, &unif);
        if (in.fail())
            throw gs::IOReadFailure("In npstat::UHistoAxis::read: "
                                    "input stream failure");
        NUHistoAxis* result = new NUHistoAxis(binEdges, label.c_str());
        result->uniform_ = unif;
        return result;
    }
}