Utilities/interface/HistoChiSquare.h

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