include/root/RooHistError.h

0001 /*****************************************************************************
0002  * Project: RooFit                                                           *
0003  * Package: RooFitCore                                                       *
0004  *    File: $Id: RooHistError.h,v 1.14 2007/05/11 09:11:30 verkerke Exp $
0005  * Authors:                                                                  *
0006  *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
0007  *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
0008  *                                                                           *
0009  * Copyright (c) 2000-2005, Regents of the University of California          *
0010  *                          and Stanford University. All rights reserved.    *
0011  *                                                                           *
0012  * Redistribution and use in source and binary forms,                        *
0013  * with or without modification, are permitted according to the terms        *
0014  * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
0015  *****************************************************************************/
0016 #ifndef ROO_HIST_ERROR
0017 #define ROO_HIST_ERROR
0018
0019 #include "Rtypes.h"
0020 #include "RooNumber.h"
0021 #include "RooAbsFunc.h"
0022 #include <cmath>
0023 #include <iostream>
0024
0025 class RooHistError {
0026 public:
0027   static const RooHistError &instance();
0028   virtual ~RooHistError() {} ;
0029
0030   bool getPoissonInterval(Int_t n, double &mu1, double &mu2, double nSigma= 1) const;
0031   bool getBinomialIntervalAsym(Int_t n, Int_t m, double &a1, double &a2, double nSigma= 1) const;
0032   bool getBinomialIntervalEff(Int_t n, Int_t m, double &a1, double &a2, double nSigma= 1) const;
0033   bool getInterval(const RooAbsFunc *Qu, const RooAbsFunc *Ql, double pointEstimate, double stepSize,
0034            double &lo, double &hi, double nSigma) const;
0035
0036   static RooAbsFunc *createPoissonSum(Int_t n) ;
0037   static RooAbsFunc *createBinomialSum(Int_t n, Int_t m, bool eff) ;
0038
0039 private:
0040   double seek(const RooAbsFunc &f, double startAt, double step, double value) const;
0041
0042   // -----------------------------------------------------------
0043   // Define a 1-dim RooAbsFunc of mu that evaluates the sum:
0044   //
0045   //  Q(n|mu) = Sum_{k=nullptr}^{n} P(k|mu)
0046   //
0047   // where P(n|mu) = exp(-mu) mu**n / n! is the Poisson PDF.
0048   // -----------------------------------------------------------
0049   class PoissonSum : public RooAbsFunc {
0050   public:
0051     inline PoissonSum(Int_t n) : RooAbsFunc(1), _n(n) { }
0052     inline double operator()(const double xvec[]) const override {
0053        double mu(xvec[0]);
0054        double result(1);
0055        double factorial(1);
0056        for (Int_t k = 1; k <= _n; k++) {
0057           factorial *= k;
0058           result += pow(mu, k) / factorial;
0059       }
0060       return exp(-mu)*result;
0061     };
0062     inline double getMinLimit(UInt_t /*index*/) const override { return 0; }
0063     inline double getMaxLimit(UInt_t /*index*/) const override { return RooNumber::infinity() ; }
0064   private:
0065     Int_t _n;
0066   };
0067
0068   // -----------------------------------------------------------
0069   // Define a 1-dim RooAbsFunc of a that evaluates the sum:
0070   //
0071   //  Q(n|n+m,a) = Sum_{k=nullptr}^{n} B(k|n+m,a)
0072   //
0073   // where B(n|n+m,a) = (n+m)!/(n!m!) ((1+a)/2)**n ((1-a)/2)**m
0074   // is the Binomial PDF.
0075   // -----------------------------------------------------------
0076   class BinomialSumAsym : public RooAbsFunc {
0077   public:
0078     BinomialSumAsym(Int_t n, Int_t m) : RooAbsFunc(1), _n1(n), _N1(n+m) {
0079     }
0080     inline double operator()(const double xvec[]) const override
0081       {
0082       double p1(0.5 * (1 + xvec[0]));
0083       double p2(1 - p1);
0084       double result(0);
0085       double fact1(1);
0086       double fact2(1);
0087       for (Int_t k = 0; k <= _n1; k++) {
0088          if (k > 0) {
0089             fact2 *= k;
0090             fact1 *= _N1 - k + 1;
0091          }
0092          result += fact1 / fact2 * pow(p1, k) * pow(p2, _N1 - k);
0093    }
0094    return result;
0095       };
0096
0097     inline double getMinLimit(UInt_t /*index*/) const override { return -1; }
0098     inline double getMaxLimit(UInt_t /*index*/) const override { return +1; }
0099
0100   private:
0101     Int_t _n1 ; ///< WVE Solaris CC5 doesn't want _n or _N here (likely compiler bug)
0102     Int_t _N1 ;
0103   } ;
0104
0105
0106   // -----------------------------------------------------------
0107   // Define a 1-dim RooAbsFunc of a that evaluates the sum:
0108   //
0109   //  Q(n|n+m,a) = Sum_{k=nullptr}^{n} B(k|n+m,a)
0110   //
0111   // where B(n|n+m,a) = (n+m)!/(n!m!) ((1+a)/2)**n ((1-a)/2)**m
0112   // is the Binomial PDF.
0113   // -----------------------------------------------------------
0114   class BinomialSumEff : public RooAbsFunc {
0115   public:
0116     BinomialSumEff(Int_t n, Int_t m) : RooAbsFunc(1), _n1(n), _N1(n+m) {
0117     }
0118     inline double operator()(const double xvec[]) const override
0119       {
0120    double p1(xvec[0]);
0121    double p2(1 - p1);
0122    double result(0);
0123    double fact1(1);
0124    double fact2(1);
0125    for(Int_t k= 0; k <= _n1; k++) {
0126      if(k > 0) { fact2*= k; fact1*= _N1-k+1; }
0127      result+= fact1/fact2*pow(p1,k)*pow(p2,_N1-k);
0128    }
0129    return result;
0130       };
0131
0132     inline double getMinLimit(UInt_t /*index*/) const override { return  0; }
0133     inline double getMaxLimit(UInt_t /*index*/) const override { return +1; }
0134
0135   private:
0136     Int_t _n1 ; ///< WVE Solaris CC5 doesn't want _n or _N here (likely compiler bug)
0137     Int_t _N1 ;
0138   } ;
0139
0140 };
0141
0142 #endif