/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitCore                                                       *
 *    File: $Id: RooHistError.h,v 1.14 2007/05/11 09:11:30 verkerke Exp $
 * Authors:                                                                  *
 *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
 *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
 *                                                                           *
 * Copyright (c) 2000-2005, Regents of the University of California          *
 *                          and Stanford University. All rights reserved.    *
 *                                                                           *
 * Redistribution and use in source and binary forms,                        *
 * with or without modification, are permitted according to the terms        *
 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
 *****************************************************************************/
#ifndef ROO_HIST_ERROR
#define ROO_HIST_ERROR

#include "Rtypes.h"
#include "RooNumber.h"
#include "RooAbsFunc.h"
#include <math.h>
#include <iostream>

class RooHistError {
public:
  static const RooHistError &instance();
  virtual ~RooHistError() {} ;

  Bool_t getPoissonInterval(Int_t n, Double_t &mu1, Double_t &mu2, Double_t nSigma= 1) const;
  Bool_t getBinomialIntervalAsym(Int_t n, Int_t m, Double_t &a1, Double_t &a2, Double_t nSigma= 1) const;
  Bool_t getBinomialIntervalEff(Int_t n, Int_t m, Double_t &a1, Double_t &a2, Double_t nSigma= 1) const;
  Bool_t getInterval(const RooAbsFunc *Qu, const RooAbsFunc *Ql, Double_t pointEstimate, Double_t stepSize,
		     Double_t &lo, Double_t &hi, Double_t nSigma) const;
  
  static RooAbsFunc *createPoissonSum(Int_t n) ;
  static RooAbsFunc *createBinomialSum(Int_t n, Int_t m, Bool_t eff) ;

private:


  Bool_t getPoissonIntervalCalc(Int_t n, Double_t &mu1, Double_t &mu2, Double_t nSigma= 1) const;
  Double_t _poissonLoLUT[1000] ;
  Double_t _poissonHiLUT[1000] ;

  RooHistError();
  Double_t seek(const RooAbsFunc &f, Double_t startAt, Double_t step, Double_t value) const;

  // -----------------------------------------------------------
  // Define a 1-dim RooAbsFunc of mu that evaluates the sum:
  //
  //  Q(n|mu) = Sum_{k=0}^{n} P(k|mu)
  //
  // where P(n|mu) = exp(-mu) mu**n / n! is the Poisson PDF.
  // -----------------------------------------------------------
  class PoissonSum : public RooAbsFunc {
  public:
    inline PoissonSum(Int_t n) : RooAbsFunc(1), _n(n) { }
    inline Double_t operator()(const Double_t xvec[]) const {
      Double_t mu(xvec[0]),result(1),factorial(1);
      for(Int_t k= 1; k <= _n; k++) {
	factorial*= k;
	result+= pow(mu,k)/factorial;
      }
      return exp(-mu)*result;
    };
    inline Double_t getMinLimit(UInt_t /*index*/) const { return 0; }
    inline Double_t getMaxLimit(UInt_t /*index*/) const { return RooNumber::infinity() ; }
  private:
    Int_t _n;
  };

  // -----------------------------------------------------------
  // Define a 1-dim RooAbsFunc of a that evaluates the sum:
  //
  //  Q(n|n+m,a) = Sum_{k=0}^{n} B(k|n+m,a)
  //
  // where B(n|n+m,a) = (n+m)!/(n!m!) ((1+a)/2)**n ((1-a)/2)**m
  // is the Binomial PDF.
  // -----------------------------------------------------------
  class BinomialSumAsym : public RooAbsFunc {
  public:
    BinomialSumAsym(Int_t n, Int_t m) : RooAbsFunc(1), _n1(n), _N1(n+m) { 
    }
    inline Double_t operator()(const Double_t xvec[]) const 
      {
	Double_t p1(0.5*(1+xvec[0])),p2(1-p1),result(0),fact1(1),fact2(1);
	for(Int_t k= 0; k <= _n1; k++) {
	  if(k > 0) { fact2*= k; fact1*= _N1-k+1; }
	  result+= fact1/fact2*pow(p1,k)*pow(p2,_N1-k);
	}
	return result;
      };

    inline Double_t getMinLimit(UInt_t /*index*/) const { return -1; }
    inline Double_t getMaxLimit(UInt_t /*index*/) const { return +1; }

  private:
    Int_t _n1 ; // WVE Solaris CC5 doesn't want _n or _N here (likely compiler bug)
    Int_t _N1 ;
  } ;


  // -----------------------------------------------------------
  // Define a 1-dim RooAbsFunc of a that evaluates the sum:
  //
  //  Q(n|n+m,a) = Sum_{k=0}^{n} B(k|n+m,a)
  //
  // where B(n|n+m,a) = (n+m)!/(n!m!) ((1+a)/2)**n ((1-a)/2)**m
  // is the Binomial PDF.
  // -----------------------------------------------------------
  class BinomialSumEff : public RooAbsFunc {
  public:
    BinomialSumEff(Int_t n, Int_t m) : RooAbsFunc(1), _n1(n), _N1(n+m) { 
    }
    inline Double_t operator()(const Double_t xvec[]) const 
      {
	Double_t p1(xvec[0]),p2(1-p1),result(0),fact1(1),fact2(1);
	for(Int_t k= 0; k <= _n1; k++) {
	  if(k > 0) { fact2*= k; fact1*= _N1-k+1; }
	  result+= fact1/fact2*pow(p1,k)*pow(p2,_N1-k);
	}
	return result;
      };

    inline Double_t getMinLimit(UInt_t /*index*/) const { return  0; }
    inline Double_t getMaxLimit(UInt_t /*index*/) const { return +1; }

  private:
    Int_t _n1 ; // WVE Solaris CC5 doesn't want _n or _N here (likely compiler bug)
    Int_t _N1 ;
  } ;

  ClassDef(RooHistError,1) // Utility class for calculating histogram errors
};

#endif
 RooHistError.h:1
 RooHistError.h:2
 RooHistError.h:3
 RooHistError.h:4
 RooHistError.h:5
 RooHistError.h:6
 RooHistError.h:7
 RooHistError.h:8
 RooHistError.h:9
 RooHistError.h:10
 RooHistError.h:11
 RooHistError.h:12
 RooHistError.h:13
 RooHistError.h:14
 RooHistError.h:15
 RooHistError.h:16
 RooHistError.h:17
 RooHistError.h:18
 RooHistError.h:19
 RooHistError.h:20
 RooHistError.h:21
 RooHistError.h:22
 RooHistError.h:23
 RooHistError.h:24
 RooHistError.h:25
 RooHistError.h:26
 RooHistError.h:27
 RooHistError.h:28
 RooHistError.h:29
 RooHistError.h:30
 RooHistError.h:31
 RooHistError.h:32
 RooHistError.h:33
 RooHistError.h:34
 RooHistError.h:35
 RooHistError.h:36
 RooHistError.h:37
 RooHistError.h:38
 RooHistError.h:39
 RooHistError.h:40
 RooHistError.h:41
 RooHistError.h:42
 RooHistError.h:43
 RooHistError.h:44
 RooHistError.h:45
 RooHistError.h:46
 RooHistError.h:47
 RooHistError.h:48
 RooHistError.h:49
 RooHistError.h:50
 RooHistError.h:51
 RooHistError.h:52
 RooHistError.h:53
 RooHistError.h:54
 RooHistError.h:55
 RooHistError.h:56
 RooHistError.h:57
 RooHistError.h:58
 RooHistError.h:59
 RooHistError.h:60
 RooHistError.h:61
 RooHistError.h:62
 RooHistError.h:63
 RooHistError.h:64
 RooHistError.h:65
 RooHistError.h:66
 RooHistError.h:67
 RooHistError.h:68
 RooHistError.h:69
 RooHistError.h:70
 RooHistError.h:71
 RooHistError.h:72
 RooHistError.h:73
 RooHistError.h:74
 RooHistError.h:75
 RooHistError.h:76
 RooHistError.h:77
 RooHistError.h:78
 RooHistError.h:79
 RooHistError.h:80
 RooHistError.h:81
 RooHistError.h:82
 RooHistError.h:83
 RooHistError.h:84
 RooHistError.h:85
 RooHistError.h:86
 RooHistError.h:87
 RooHistError.h:88
 RooHistError.h:89
 RooHistError.h:90
 RooHistError.h:91
 RooHistError.h:92
 RooHistError.h:93
 RooHistError.h:94
 RooHistError.h:95
 RooHistError.h:96
 RooHistError.h:97
 RooHistError.h:98
 RooHistError.h:99
 RooHistError.h:100
 RooHistError.h:101
 RooHistError.h:102
 RooHistError.h:103
 RooHistError.h:104
 RooHistError.h:105
 RooHistError.h:106
 RooHistError.h:107
 RooHistError.h:108
 RooHistError.h:109
 RooHistError.h:110
 RooHistError.h:111
 RooHistError.h:112
 RooHistError.h:113
 RooHistError.h:114
 RooHistError.h:115
 RooHistError.h:116
 RooHistError.h:117
 RooHistError.h:118
 RooHistError.h:119
 RooHistError.h:120
 RooHistError.h:121
 RooHistError.h:122
 RooHistError.h:123
 RooHistError.h:124
 RooHistError.h:125
 RooHistError.h:126
 RooHistError.h:127
 RooHistError.h:128
 RooHistError.h:129
 RooHistError.h:130
 RooHistError.h:131
 RooHistError.h:132
 RooHistError.h:133
 RooHistError.h:134
 RooHistError.h:135
 RooHistError.h:136
 RooHistError.h:137