Utilities/interface/HistoPdf.h

0001 #ifndef PhysicsTools_Utilities_HistoPdf_h
0002 #define PhysicsTools_Utilities_HistoPdf_h
0003
0004 #include "FWCore/Utilities/interface/EDMException.h"
0005
0006 #include <iostream>
0007 #include "TH1.h"
0008
0009 namespace funct {
0010   class HistoPdf {
0011   public:
0012     template <typename Iterator>
0013     HistoPdf(double xMin, double xMax, const Iterator& begin, const Iterator& end)
0014         : xMin_(xMin), xMax_(xMax), delta_(xMax - xMin), binSize_(delta_ / (end - begin)), y_(end - begin) {
0015       double s = 0;
0016       unsigned int i = 0;
0017       for (Iterator it = begin; it != end; ++it)
0018         s += (y_[i++] = *it);
0019       for (std::vector<double>::iterator i = y_.begin(); i != y_.end(); ++i)
0020         *i /= s;
0021     }
0022     HistoPdf() {}
0023     template <typename Iterator>
0024     void init(double xMin, double xMax, const Iterator& begin, const Iterator& end) {
0025       xMin_ = xMin;
0026       xMax_ = xMax;
0027       delta_ = xMax - xMin;
0028       unsigned int n = end - begin;
0029       binSize_ = delta_ / n;
0030       y_.resize(n);
0031       double s = 0;
0032       unsigned int i = 0;
0033       for (Iterator it = begin; it != end; ++it)
0034         s += (y_[i++] = *it);
0035       for (std::vector<double>::iterator i = y_.begin(); i != y_.end(); ++i)
0036         *i /= s;
0037     }
0038     double operator()(double x) const {
0039       if (x < xMin_ || x > xMax_)
0040         return 0;
0041       double pdf = y_[static_cast<unsigned int>(((x - xMin_) / delta_) * y_.size())] / binSize_;
0042       return pdf;
0043     }
0044     void rebin(unsigned int r) {
0045       if (y_.size() % r != 0)
0046         throw edm::Exception(edm::errors::Configuration)
0047             << "HistoPdf: can't rebin histogram of " << y_.size() << " entries by " << r << "\n";
0048       unsigned int n = y_.size() / r;
0049       std::vector<double> y(n, 0);
0050       for (unsigned int i = 0, j = 0; i < n; ++i)
0051         for (unsigned int k = 0; k < r; ++k)
0052           y[i] += y_[j++];
0053       y_ = y;
0054       binSize_ *= r;
0055     }
0056     void dump() {
0057       std::cout << ">>> range: [" << xMin_ << ", " << xMax_ << "], bin size: " << delta_ << "/" << y_.size() << " = "
0058                 << binSize_ << std::endl;
0059       double s = 0;
0060       for (unsigned int i = 0; i != y_.size(); ++i) {
0061         double x = xMin_ + (0.5 + i) * binSize_;
0062         double y = operator()(x);
0063         std::cout << ">>> pdf(" << x << ") = " << y << std::endl;
0064         s += y * binSize_;
0065       }
0066       std::cout << ">>>: PDF normalization is " << s << std::endl;
0067     }
0068
0069   private:
0070     double xMin_, xMax_, delta_, binSize_;
0071     std::vector<double> y_;
0072   };
0073
0074   class RootHistoPdf : public HistoPdf {
0075   public:
0076     explicit RootHistoPdf(const TH1& histo, double fMin, double fMax) {
0077       unsigned int nBins = histo.GetNbinsX();
0078       std::vector<double> y;
0079       y.reserve(nBins);
0080       double xMin = histo.GetXaxis()->GetXmin();
0081       double xMax = histo.GetXaxis()->GetXmax();
0082       double deltaX = (xMax - xMin) / nBins;
0083       for (unsigned int i = 0; i != nBins; ++i) {
0084         double x = xMin + (i + .5) * deltaX;
0085         if (x > fMin && x < fMax) {
0086           y.push_back(histo.GetBinContent(i + 1));
0087         }
0088       }
0089       init(fMin, fMax, y.begin(), y.end());
0090     }
0091   };
0092
0093 }  // namespace funct
0094
0095 #endif