// @(#)root/roostats:$Id$
// Author: Kyle Cranmer   28/07/2008

/*************************************************************************
 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/


#ifndef RooStats_NumberCountingUtils
#include "RooStats/NumberCountingUtils.h"
#endif

#ifndef RooStats_RooStatsUtils
#include "RooStats/RooStatsUtils.h"
#endif

// // Without this macro the THtml doc  can not be generated
// #if !defined(R__ALPHA) && !defined(R__SOLARIS) && !defined(R__ACC) && !defined(R__FBSD)
// NamespaceImp(RooStats)
// //NamespaceImp(NumberCountingUtils)
// #endif

//using namespace RooStats;

Double_t RooStats::NumberCountingUtils::BinomialExpP(Double_t signalExp, Double_t backgroundExp, Double_t relativeBkgUncert){
  // Expected P-value for s=0 in a ratio of Poisson means.  
  // Here the background and its uncertainty are provided directly and 
  // assumed to be from the double Poisson counting setup described in the 
  // BinomialWithTau functions.  
  // Normally one would know tau directly, but here it is determiend from
  // the background uncertainty.  This is not strictly correct, but a useful 
  // approximation.


  //SIDE BAND EXAMPLE
  //See Eqn. (19) of Cranmer and pp. 36-37 of Linnemann.
  //150 total events in signalExp region, 100 in sideband of equal size
  Double_t mainInf = signalExp+backgroundExp;  //Given
  Double_t tau = 1./backgroundExp/(relativeBkgUncert*relativeBkgUncert);
  Double_t auxiliaryInf = backgroundExp*tau;  //Given
  
  Double_t P_Bi = TMath::BetaIncomplete(1./(1.+tau),mainInf,auxiliaryInf+1);
  return P_Bi;
  
/*
Now, if instead the mean background level b in the signal region is
specified, along with Gaussian rms sigb, then one can fake a Poisson
sideband (see Linnemann, p. 35, converted to Cranmer's notation) by
letting tau = b/(sigb*sigb) and y = b*tau.  Thus, for example, if one
has x=150 and b = 100 +/- 10, one then derives tau and y.  Then one
has the same two lines of ROOT calling BetaIncomplete and ErfInverse.
Since I chose these numbers to revert to the previous example, we get
the same answer:
*/
/*
//GAUSSIAN FAKED AS POISSON EXAMPLE
x = 150.    //Given
b = 100.    //Given
sigb = 10.  //Given
tau = b/(sigb*sigb)
y = tau*b   
Z_Bi = TMath::BetaIncomplete(1./(1.+tau),x,y+1)    
S = sqrt(2)*TMath::ErfInverse(1 - 2*Z_Bi)     

*/

}


Double_t RooStats::NumberCountingUtils::BinomialWithTauExpP(Double_t signalExp, Double_t backgroundExp, Double_t tau){
  // Expected P-value for s=0 in a ratio of Poisson means.  
  // Based on two expectations, a main measurement that might have signal
  // and an auxiliarly measurement for the background that is signal free.
  // The expected background in the auxiliary measurement is a factor
  // tau larger than in the main measurement.

  Double_t mainInf = signalExp+backgroundExp;  //Given
  Double_t auxiliaryInf = backgroundExp*tau;  //Given
  
  Double_t P_Bi = TMath::BetaIncomplete(1./(1.+tau),mainInf,auxiliaryInf+1);
  
  return P_Bi;
  
}

Double_t RooStats::NumberCountingUtils::BinomialObsP(Double_t mainObs, Double_t backgroundObs, Double_t relativeBkgUncert){
  // P-value for s=0 in a ratio of Poisson means.  
  // Here the background and its uncertainty are provided directly and 
  // assumed to be from the double Poisson counting setup.  
  // Normally one would know tau directly, but here it is determiend from
  // the background uncertainty.  This is not strictly correct, but a useful 
  // approximation.
  
  Double_t tau = 1./backgroundObs/(relativeBkgUncert*relativeBkgUncert);
  Double_t auxiliaryInf = backgroundObs*tau;  //Given
  
    
  //SIDE BAND EXAMPLE
  //See Eqn. (19) of Cranmer and pp. 36-37 of Linnemann.
  Double_t P_Bi = TMath::BetaIncomplete(1./(1.+tau),mainObs,auxiliaryInf+1);
  
  return P_Bi;

}


Double_t RooStats::NumberCountingUtils::BinomialWithTauObsP(Double_t mainObs, Double_t auxiliaryObs, Double_t tau){
  // P-value for s=0 in a ratio of Poisson means.  
  // Based on two observations, a main measurement that might have signal
  // and an auxiliarly measurement for the background that is signal free.
  // The expected background in the auxiliary measurement is a factor
  // tau larger than in the main measurement.

  //SIDE BAND EXAMPLE
  //See Eqn. (19) of Cranmer and pp. 36-37 of Linnemann.
  Double_t P_Bi = TMath::BetaIncomplete(1./(1.+tau),mainObs,auxiliaryObs+1);
  
  return P_Bi;
  
}

Double_t RooStats::NumberCountingUtils::BinomialExpZ(Double_t signalExp, Double_t backgroundExp, Double_t relativeBkgUncert) {    
  // See BinomialExpP
  return RooStats::PValueToSignificance( BinomialExpP(signalExp,backgroundExp,relativeBkgUncert) ) ;
  }

Double_t RooStats::NumberCountingUtils::BinomialWithTauExpZ(Double_t signalExp, Double_t backgroundExp, Double_t tau){
  // See BinomialWithTauExpP
  return RooStats::PValueToSignificance( BinomialWithTauExpP(signalExp,backgroundExp,tau) ) ;
}


Double_t RooStats::NumberCountingUtils::BinomialObsZ(Double_t mainObs, Double_t backgroundObs, Double_t relativeBkgUncert){
  // See BinomialObsP
  return RooStats::PValueToSignificance( BinomialObsP(mainObs,backgroundObs,relativeBkgUncert) ) ;
}

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