UniformAxis.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/UniformAxis.h"
#include "JetMETCorrections/InterpolationTables/interface/closeWithinTolerance.h"

namespace npstat {
    UniformAxis::UniformAxis(const unsigned nCoords,
                             const double min, const double max,
                             const char* label)
        : min_(min), max_(max), label_(label ? label : ""), npt_(nCoords)
    {
        if (!(npt_ > 1U && npt_ < UINT_MAX/2U - 1U))
            throw npstat::NpstatInvalidArgument("In npstat::UniformAxis constructor: "
                                        "number of points is out of range");
        if (min_ > max_)
            std::swap(min_, max_);
        bw_ = (max_ - min_)/(npt_ - 1U);
        if (max_ == min_)
            throw npstat::NpstatInvalidArgument(
                "In npstat::UniformAxis constructor: "
                "minimum and maximum must be distinct");
    }

    std::pair<unsigned,double> UniformAxis::getInterval(const double x) const
    {
        if (x <= min_)
            return std::pair<unsigned,double>(0U, 1.0);
        else if (x >= max_)
            return std::pair<unsigned,double>(npt_ - 2U, 0.0);
        else
        {
            unsigned binnum = static_cast<unsigned>(floor((x - min_)/bw_));
            if (binnum > npt_ - 2U)
                binnum = npt_ - 2U;
            double w = binnum + 1.0 - (x - min_)/bw_;
            if (w < 0.0)
                w = 0.0;
            else if (w > 1.0)
                w = 1.0;
            return std::pair<unsigned,double>(binnum, w);
        }
    }

    std::pair<unsigned,double> UniformAxis::linearInterval(const double x) const
    {
        if (x <= min_)
            return std::pair<unsigned,double>(0U, 1.0 - (x - min_)/bw_);
        else if (x >= max_)
            return std::pair<unsigned,double>(npt_ - 2U, (max_ - x)/bw_);
        else
        {
            unsigned binnum = static_cast<unsigned>(floor((x - min_)/bw_));
            if (binnum > npt_ - 2U)
                binnum = npt_ - 2U;
            double w = binnum + 1.0 - (x - min_)/bw_;
            if (w < 0.0)
                w = 0.0;
            else if (w > 1.0)
                w = 1.0;
            return std::pair<unsigned,double>(binnum, w);
        }
    }

    std::vector<double> UniformAxis::coords() const
    {
        std::vector<double> vec;
        vec.reserve(npt_);
        const unsigned nptm1 = npt_ - 1U;
        for (unsigned i=0; i<nptm1; ++i)
            vec.push_back(min_ + bw_*i);
        vec.push_back(max_);
        return vec;
    }

    double UniformAxis::coordinate(const unsigned i) const
    {
        if (i >= npt_)
            throw npstat::NpstatOutOfRange(
                "In npstat::UniformAxis::coordinate: index out of range");
        if (i == npt_ - 1U)
            return max_;
        else
            return min_ + bw_*i;
    }

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

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

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

    UniformAxis* UniformAxis::read(const gs::ClassId& id, std::istream& in)
    {
        static const gs::ClassId current(gs::ClassId::makeId<UniformAxis>());
        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 UniformAxis(nBins, min, max, label.c_str());
        else
            throw gs::IOReadFailure("In npstat::UniformAxis::read: "
                                    "input stream failure");
    }
}