File indexing completed on 2024-04-06 12:24:18
0001 #ifndef PhysicsTools_Utilities_Likelihood_h
0002 #define PhysicsTools_Utilities_Likelihood_h
0003 #include "PhysicsTools/Utilities/interface/RootMinuitResultPrinter.h"
0004 #include "PhysicsTools/Utilities/interface/RootMinuitFuncEvaluator.h"
0005 #include <cmath>
0006 #include "TMath.h"
0007
0008 namespace fit {
0009 template <typename PDF, typename Tuple>
0010 struct LikelihoodEvaluator {
0011 static double evaluate(const PDF& pdf, const Tuple& tuple) { return pdf(tuple); }
0012 };
0013
0014
0015 template <typename PDF>
0016 struct LikelihoodEvaluator<PDF, double> {
0017 static double evaluate(const PDF& pdf, const double val) { return pdf(val); }
0018 };
0019
0020 struct NoExtendedLikelihood {};
0021
0022 template <typename Sample, typename PDF, typename Yield = NoExtendedLikelihood>
0023 class Likelihood {
0024 public:
0025 Likelihood() {}
0026 Likelihood(const Sample& sample, PDF& pdf, Yield& yield) : pdf_(&pdf), yield_(&yield), sample_(sample) {}
0027 double operator()() const { return log(); }
0028 double log() const {
0029 double l = -(*yield_)();
0030 for (typename Sample::const_iterator i = sample_.begin(); i != sample_.end(); ++i) {
0031 double p = Evaluator::evaluate(*pdf_, *i);
0032 l += std::log(p);
0033 }
0034 sampleSize_ = sample_.size();
0035 return l;
0036 }
0037 double logNFactorial() const { return std::log(TMath::Factorial(sampleSize_)); }
0038 double absoluteLog() const { return log() - logNFactorial(); }
0039 PDF& pdf() { return *pdf_; }
0040 const PDF& pdf() const { return *pdf_; }
0041 Yield& yield() { return *yield_; }
0042 const Yield& yield() const { return *yield_; }
0043 unsigned int sampleSize() const { return sampleSize_; }
0044
0045 private:
0046 typedef LikelihoodEvaluator<PDF, typename Sample::value_type> Evaluator;
0047 PDF* pdf_;
0048 Yield* yield_;
0049 Sample sample_;
0050 mutable unsigned int sampleSize_ = 0u;
0051 };
0052
0053 template <typename Sample, typename PDF>
0054 class Likelihood<Sample, PDF, NoExtendedLikelihood> {
0055 public:
0056 Likelihood() {}
0057 Likelihood(const Sample& sample, PDF& pdf) : pdf_(&pdf), sample_(sample) {}
0058 double operator()() const { return log(); }
0059 double log() const {
0060 double l = 0;
0061 for (typename Sample::const_iterator i = sample_.begin(); i != sample_.end(); ++i) {
0062 l += std::log(Evaluator::evaluate(*pdf_, *i));
0063 }
0064 return l;
0065 }
0066 PDF& pdf() { return *pdf_; }
0067 const PDF& pdf() const { return *pdf_; }
0068
0069 private:
0070 typedef LikelihoodEvaluator<PDF, typename Sample::value_type> Evaluator;
0071 PDF* pdf_;
0072 Sample sample_;
0073 };
0074
0075 template <typename Sample, typename PDF, typename Yield>
0076 struct RootMinuitResultPrinter<Likelihood<Sample, PDF, Yield> > {
0077 static void print(double amin, unsigned int numberOfFreeParameters, const Likelihood<Sample, PDF, Yield>& f) {
0078 std::cout << "-2 log(maximum-likelihood) = " << amin << ", free parameters = " << numberOfFreeParameters
0079 << std::endl;
0080 }
0081 };
0082
0083 template <typename Sample, typename PDF, typename Yield>
0084 struct RootMinuitFuncEvaluator<Likelihood<Sample, PDF, Yield> > {
0085 static double evaluate(const Likelihood<Sample, PDF, Yield>& f) { return -2 * f(); }
0086 };
0087
0088 }
0089
0090 #endif