Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:33:27

0001 #ifndef PhysicsTools_Utilities_HistoPoissonLikelihoodRatio_h
0002 #define PhysicsTools_Utilities_HistoPoissonLikelihoodRatio_h
0003 #include "PhysicsTools/Utilities/interface/RootMinuitResultPrinter.h"
0004 #include <vector>
0005 #include <cmath>
0006 #include "TH1.h"
0007 #include "TMath.h"
0008 
0009 namespace fit {
0010   template <typename T>
0011   class HistoPoissonLikelihoodRatio {
0012   public:
0013     HistoPoissonLikelihoodRatio() {}
0014     HistoPoissonLikelihoodRatio(T &t, TH1 *histo, double rangeMin, double rangeMax)
0015         : t_(&t), rangeMin_(rangeMin), rangeMax_(rangeMax) {
0016       nBins_ = histo->GetNbinsX();
0017       xMin_ = histo->GetXaxis()->GetXmin();
0018       xMax_ = histo->GetXaxis()->GetXmax();
0019       deltaX_ = (xMax_ - xMin_) / nBins_;
0020       for (size_t i = 0; i < nBins_; ++i) {
0021         cont_.push_back(histo->GetBinContent(i + 1));
0022       }
0023     }
0024     double operator()() const {
0025       double chi2lambda = 0;
0026       for (size_t i = 0; i < nBins_; ++i) {
0027         double x = xMin_ + (i + .5) * deltaX_;
0028         if ((x > rangeMin_) && (x < rangeMax_)) {
0029           double nu = (*t_)(x);
0030           if (nu > 0 && cont_[i] > 0)
0031             chi2lambda += nu - cont_[i] + cont_[i] * log(cont_[i] / nu);
0032         }
0033       }
0034       chi2lambda *= 2;
0035       return chi2lambda;
0036     }
0037     void setHistos(TH1 *histo) {
0038       nBins_ = histo->GetNbinsX();
0039       xMin_ = histo->GetXaxis()->GetXmin();
0040       xMax_ = histo->GetXaxis()->GetXmax();
0041       deltaX_ = (xMax_ - xMin_) / nBins_;
0042     }
0043     size_t numberOfBins() const {
0044       size_t fullBins = 0;
0045       for (size_t i = 0; i < nBins_; ++i) {
0046         double x = xMin_ + (i + .5) * deltaX_;
0047         if ((x > rangeMin_) && (x < rangeMax_))
0048           fullBins++;
0049       }
0050       return fullBins;
0051     }
0052     T &function() { return *t_; }
0053     const T &function() const { return *t_; }
0054 
0055   private:
0056     T *t_;
0057     double rangeMin_, rangeMax_;
0058     size_t nBins_;
0059     double xMin_, xMax_, deltaX_;
0060     std::vector<double> cont_;
0061   };
0062 
0063   template <typename T>
0064   struct RootMinuitResultPrinter<HistoPoissonLikelihoodRatio<T> > {
0065     static void print(double amin, unsigned int numberOfFreeParameters, const HistoPoissonLikelihoodRatio<T> &f) {
0066       unsigned int ndof = f.numberOfBins() - numberOfFreeParameters;
0067       std::cout << "chi-squared/n.d.o.f. = " << amin << "/" << ndof << " = " << amin / ndof
0068                 << "; prob: " << TMath::Prob(amin, ndof) << std::endl;
0069     }
0070   };
0071 }  // namespace fit
0072 
0073 #endif