// @(#)root/roostats:$Id: NumberCountingUtils.cxx 26324 2008-11-20 17:17:32Z moneta$
// Author: Kyle Cranmer   28/07/2008

/*************************************************************************
*                                                                       *
* For the licensing terms see $ROOTSYS/LICENSE. * * For the list of contributors see$ROOTSYS/README/CREDITS.             *
*************************************************************************/

/////////////////////////////////////////
// NumberCountingUtils
//
// Encapsulates common number counting utilities
/////////////////////////////////////////
///////////////////////////////////
// Standalone Functions.
// Naming conventions:
//  Exp = Expected
//  Obs = Observed
//  P   = p-value
//  Z   = Z-value or significance in Sigma (one-sided convention)
//////////////////////////////////

#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.
//
//
// This is based on code and comments from Bob Cousins
//  based on the following papers:
//
// Statistical Challenges for Searches for New Physics at the LHC
// Authors: Kyle Cranmer
// http://arxiv.org/abs/physics/0511028
//
//  Measures of Significance in HEP and Astrophysics
//  Authors: J. T. Linnemann
//  http://arxiv.org/abs/physics/0312059
//
// In short, this is the exact frequentist solution to the problem of
// a main measurement x distributed as a Poisson around s+b and a sideband or
// auxiliary measurement y distributed as a Poisson around \taub.  Eg.
// L(x,y|s,b,\tau) = Pois(x|s+b)Pois(y|\tau b)

//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
*/
/*
//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 BinomialObsZ
return RooStats::PValueToSignificance( BinomialObsP(mainObs,backgroundObs,relativeBkgUncert) ) ;
}

Double_t RooStats::NumberCountingUtils::BinomialWithTauObsZ(Double_t mainObs, Double_t auxiliaryObs, Double_t tau){
// See BinomialWithTauObsZ
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
NumberCountingUtils.cxx:146
NumberCountingUtils.cxx:147
NumberCountingUtils.cxx:148
NumberCountingUtils.cxx:149
NumberCountingUtils.cxx:150
NumberCountingUtils.cxx:151
NumberCountingUtils.cxx:152
NumberCountingUtils.cxx:153
NumberCountingUtils.cxx:154
NumberCountingUtils.cxx:155
NumberCountingUtils.cxx:156
NumberCountingUtils.cxx:157
NumberCountingUtils.cxx:158
NumberCountingUtils.cxx:159
NumberCountingUtils.cxx:160
NumberCountingUtils.cxx:161
NumberCountingUtils.cxx:162
NumberCountingUtils.cxx:163
NumberCountingUtils.cxx:164
NumberCountingUtils.cxx:165
NumberCountingUtils.cxx:166
NumberCountingUtils.cxx:167
NumberCountingUtils.cxx:168
NumberCountingUtils.cxx:169
NumberCountingUtils.cxx:170
NumberCountingUtils.cxx:171
NumberCountingUtils.cxx:172
NumberCountingUtils.cxx:173
NumberCountingUtils.cxx:174
NumberCountingUtils.cxx:175