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
0041
0042   bool getPoissonIntervalCalc(Int_t n, double &mu1, double &mu2, double nSigma= 1) const;
0043   double _poissonLoLUT[1000] ;
0044   double _poissonHiLUT[1000] ;
0045
0046   RooHistError();
0047   double seek(const RooAbsFunc &f, double startAt, double step, double value) const;
0048
0049   // -----------------------------------------------------------
0050   // Define a 1-dim RooAbsFunc of mu that evaluates the sum:
0051   //
0052   //  Q(n|mu) = Sum_{k=nullptr}^{n} P(k|mu)
0053   //
0054   // where P(n|mu) = exp(-mu) mu**n / n! is the Poisson PDF.
0055   // -----------------------------------------------------------
0056   class PoissonSum : public RooAbsFunc {
0057   public:
0058     inline PoissonSum(Int_t n) : RooAbsFunc(1), _n(n) { }
0059     inline double operator()(const double xvec[]) const override {
0060        double mu(xvec[0]);
0061        double result(1);
0062        double factorial(1);
0063        for (Int_t k = 1; k <= _n; k++) {
0064           factorial *= k;
0065           result += pow(mu, k) / factorial;
0066       }
0067       return exp(-mu)*result;
0068     };
0069     inline double getMinLimit(UInt_t /*index*/) const override { return 0; }
0070     inline double getMaxLimit(UInt_t /*index*/) const override { return RooNumber::infinity() ; }
0071   private:
0072     Int_t _n;
0073   };
0074
0075   // -----------------------------------------------------------
0076   // Define a 1-dim RooAbsFunc of a that evaluates the sum:
0077   //
0078   //  Q(n|n+m,a) = Sum_{k=nullptr}^{n} B(k|n+m,a)
0079   //
0080   // where B(n|n+m,a) = (n+m)!/(n!m!) ((1+a)/2)**n ((1-a)/2)**m
0081   // is the Binomial PDF.
0082   // -----------------------------------------------------------
0083   class BinomialSumAsym : public RooAbsFunc {
0084   public:
0085     BinomialSumAsym(Int_t n, Int_t m) : RooAbsFunc(1), _n1(n), _N1(n+m) {
0086     }
0087     inline double operator()(const double xvec[]) const override
0088       {
0089       double p1(0.5 * (1 + xvec[0]));
0090       double p2(1 - p1);
0091       double result(0);
0092       double fact1(1);
0093       double fact2(1);
0094       for (Int_t k = 0; k <= _n1; k++) {
0095          if (k > 0) {
0096             fact2 *= k;
0097             fact1 *= _N1 - k + 1;
0098          }
0099          result += fact1 / fact2 * pow(p1, k) * pow(p2, _N1 - k);
0100    }
0101    return result;
0102       };
0103
0104     inline double getMinLimit(UInt_t /*index*/) const override { return -1; }
0105     inline double getMaxLimit(UInt_t /*index*/) const override { return +1; }
0106
0107   private:
0108     Int_t _n1 ; ///< WVE Solaris CC5 doesn't want _n or _N here (likely compiler bug)
0109     Int_t _N1 ;
0110   } ;
0111
0112
0113   // -----------------------------------------------------------
0114   // Define a 1-dim RooAbsFunc of a that evaluates the sum:
0115   //
0116   //  Q(n|n+m,a) = Sum_{k=nullptr}^{n} B(k|n+m,a)
0117   //
0118   // where B(n|n+m,a) = (n+m)!/(n!m!) ((1+a)/2)**n ((1-a)/2)**m
0119   // is the Binomial PDF.
0120   // -----------------------------------------------------------
0121   class BinomialSumEff : public RooAbsFunc {
0122   public:
0123     BinomialSumEff(Int_t n, Int_t m) : RooAbsFunc(1), _n1(n), _N1(n+m) {
0124     }
0125     inline double operator()(const double xvec[]) const override
0126       {
0127    double p1(xvec[0]);
0128    double p2(1 - p1);
0129    double result(0);
0130    double fact1(1);
0131    double fact2(1);
0132    for(Int_t k= 0; k <= _n1; k++) {
0133      if(k > 0) { fact2*= k; fact1*= _N1-k+1; }
0134      result+= fact1/fact2*pow(p1,k)*pow(p2,_N1-k);
0135    }
0136    return result;
0137       };
0138
0139     inline double getMinLimit(UInt_t /*index*/) const override { return  0; }
0140     inline double getMaxLimit(UInt_t /*index*/) const override { return +1; }
0141
0142   private:
0143     Int_t _n1 ; ///< WVE Solaris CC5 doesn't want _n or _N here (likely compiler bug)
0144     Int_t _N1 ;
0145   } ;
0146
0147 };
0148
0149 #endif