Logo ROOT  
Reference Guide
TFeldmanCousins.cxx
Go to the documentation of this file.
1// @(#)root/physics:$Id$
2// Author: Adrian Bevan 2001
3
4/*************************************************************************
5 * Copyright (C) 1995-2004, Rene Brun and Fons Rademakers. *
6 * Copyright (C) 2001, Liverpool University. *
7 * All rights reserved. *
8 * *
9 * For the licensing terms see $ROOTSYS/LICENSE. *
10 * For the list of contributors see $ROOTSYS/README/CREDITS. *
11 *************************************************************************/
12
13/** \class TFeldmanCousins
14 \ingroup Physics
15
16Class to calculate the CL upper limit using
17the Feldman-Cousins method as described in PRD V57 #7, p3873-3889
18
19The default confidence interval calculated using this method is 90%
20This is set either by having a default the constructor, or using the
21appropriate fraction when instantiating an object of this class (e.g. 0.9)
22
23The simple extension to a gaussian resolution function bounded at zero
24has not been addressed as yet -> `time is of the essence' as they write
25on the wall of the maze in that classic game ...
26
27#### VARIABLES THAT CAN BE ALTERED
28=> depending on your desired precision: The initial values of fMuMin,
29fMuMax, fMuStep and fNMax are those used in the PRD:
30~~~ {.cpp}
31 fMuMin = 0.0
32 fMuMax = 50.0
33 fMuStep= 0.005
34~~~
35but there is total flexibility in changing this should you desire.
36
37
38see example of use in $ROOTSYS/tutorials/math/FeldmanCousins.C
39
40see note about: "Should I use TRolke, TFeldmanCousins, TLimit?"
41 in the TRolke class description.
42
43\author: Adrian Bevan, Liverpool University
44
45Copyright Liverpool University 2001 bevan@slac.stanford.edu
46*/
47
48#include <iostream>
49#include "TMath.h"
50#include "TFeldmanCousins.h"
51
53
54////////////////////////////////////////////////////////////////////////////////
55/// Constructor.
56
58{
59 fCL = newFC;
60 fUpperLimit = 0.0;
61 fLowerLimit = 0.0;
62 fNobserved = 0.0;
63 fNbackground = 0.0;
64 options.ToLower();
65 if (options.Contains("q")) fQUICK = 1;
66 else fQUICK = 0;
67
68 fNMax = 50;
69 fMuStep = 0.005;
70 SetMuMin();
71 SetMuMax();
72 SetMuStep();
73}
74
75
76////////////////////////////////////////////////////////////////////////////////
77// Destructor.
78
80{
81}
82
83////////////////////////////////////////////////////////////////////////////////
84/// given Nobserved and Nbackground, try different values of mu that give lower limits that
85/// are consistent with Nobserved. The closed interval (plus any stragglers) corresponds
86/// to the F&C interval
87
89{
90 CalculateUpperLimit(Nobserved, Nbackground);
91 return fLowerLimit;
92}
93
94
95////////////////////////////////////////////////////////////////////////////////
96/// given Nobserved and Nbackground, try different values of mu that give upper limits that
97/// are consistent with Nobserved. The closed interval (plus any stragglers) corresponds
98/// to the F&C interval
99
101{
102 fNobserved = Nobserved;
103 fNbackground = Nbackground;
104
105 Double_t mu = 0.0;
106
107 // for each mu construct the ranked table of probabilities and test the
108 // observed number of events with the upper limit
109 Double_t min = -999.0;
110 Double_t max = 0;
111 Int_t iLower = 0;
112
113 Int_t i;
114 for(i = 0; i <= fNMuStep; i++) {
115 mu = fMuMin + (Double_t)i*fMuStep;
116 Int_t goodChoice = FindLimitsFromTable( mu );
117 if( goodChoice ) {
118 min = mu;
119 iLower = i;
120 break;
121 }
122 }
123
124 //==================================================================
125 // For quicker evaluation, assume that you get the same results when
126 // you expect the uppper limit to be > Nobserved-Nbackground.
127 // This is certainly true for all of the published tables in the PRD
128 // and is a reasonable assumption in any case.
129 //==================================================================
130
131 Double_t quickJump = 0.0;
132 if (fQUICK) quickJump = Nobserved-Nbackground-fMuMin;
133 if (quickJump < 0.0) quickJump = 0.0;
134
135 for(i = iLower+1; i <= fNMuStep; i++) {
136 mu = fMuMin + (Double_t)i*fMuStep + quickJump;
137 Int_t goodChoice = FindLimitsFromTable( mu );
138 if( !goodChoice ) {
139 max = mu;
140 break;
141 }
142 }
143
144 fUpperLimit = max;
145 fLowerLimit = min;
146
147 return max;
148}
149
150////////////////////////////////////////////////////////////////////////////////
151/// calculate the probability table for a given mu for n = 0, NMAX
152/// and return 1 if the number of observed events is consistent
153/// with the CL bad
154
156{
157 Double_t *p = new Double_t[fNMax]; //the array of probabilities in the interval MUMIN-MUMAX
158 Double_t *r = new Double_t[fNMax]; //the ratio of likliehoods = P(Mu|Nobserved)/P(MuBest|Nobserved)
159 Int_t *rank = new Int_t[fNMax]; //the ranked array corresponding to R (largest first)
160 Double_t *muBest = new Double_t[fNMax];
161 Double_t *probMuBest = new Double_t[fNMax];
162
163 //calculate P(i | mu) and P(i | mu)/P(i | mubest)
164 Int_t i;
165 for(i = 0; i < fNMax; i++) {
166 muBest[i] = (Double_t)(i - fNbackground);
167 if(muBest[i]<0.0) muBest[i] = 0.0;
168 probMuBest[i] = Prob(i, muBest[i], fNbackground);
169 p[i] = Prob(i, mu, fNbackground);
170 if(probMuBest[i] == 0.0) r[i] = 0.0;
171 else r[i] = p[i]/probMuBest[i];
172 }
173
174 //rank the likelihood ratio
175 TMath::Sort(fNMax, r, rank);
176
177 //search through the probability table and get the i for the CL
178 Double_t sum = 0.0;
179 Int_t iMax = rank[0];
180 Int_t iMin = rank[0];
181 for(i = 0; i < fNMax; i++) {
182 sum += p[rank[i]];
183 if(iMax < rank[i]) iMax = rank[i];
184 if(iMin > rank[i]) iMin = rank[i];
185 if(sum >= fCL) break;
186 }
187
188 delete [] p;
189 delete [] r;
190 delete [] rank;
191 delete [] muBest;
192 delete [] probMuBest;
193
194 if((fNobserved <= iMax) && (fNobserved >= iMin)) return 1;
195 else return 0;
196}
197
198////////////////////////////////////////////////////////////////////////////////
199/// Calculate the poissonian probability for a mean of mu+B events with a variance of N.
200
202{
203 return TMath::Poisson( N, mu+B);
204}
205
206////////////////////////////////////////////////////////////////////////////////
207/// Set maximum value of signal to use in calculating the tables.
208
210{
211 fMuMax = newMax;
212 fNMax = (Int_t)newMax;
214}
215
216////////////////////////////////////////////////////////////////////////////////
217/// Set the step in signal to use when generating tables.
218
220{
221 if(newMuStep == 0.0) {
222 std::cout << "TFeldmanCousins::SetMuStep ERROR New step size is zero - unable to change value"<< std::endl;
223 return;
224 } else {
225 fMuStep = newMuStep;
227 }
228}
ROOT::R::TRInterface & r
Definition: Object.C:4
int Int_t
Definition: RtypesCore.h:43
double Double_t
Definition: RtypesCore.h:57
#define ClassImp(name)
Definition: Rtypes.h:361
#define N
Class to calculate the CL upper limit using the Feldman-Cousins method as described in PRD V57 #7,...
void SetMuMin(Double_t newMin=0.0)
Double_t fUpperLimit
Double_t fNbackground
Double_t CalculateLowerLimit(Double_t Nobserved, Double_t Nbackground)
given Nobserved and Nbackground, try different values of mu that give lower limits that are consisten...
Double_t Prob(Int_t N, Double_t mu, Double_t B)
Calculate the poissonian probability for a mean of mu+B events with a variance of N.
Int_t FindLimitsFromTable(Double_t mu)
calculate the probability table for a given mu for n = 0, NMAX and return 1 if the number of observed...
void SetMuStep(Double_t newMuStep=0.005)
Set the step in signal to use when generating tables.
Double_t fLowerLimit
void SetMuMax(Double_t newMax=50.0)
Set maximum value of signal to use in calculating the tables.
Double_t CalculateUpperLimit(Double_t Nobserved, Double_t Nbackground)
given Nobserved and Nbackground, try different values of mu that give upper limits that are consisten...
TFeldmanCousins(Double_t newCL=0.9, TString options="")
Constructor.
virtual ~TFeldmanCousins()
Basic string class.
Definition: TString.h:131
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
static double B[]
Double_t Poisson(Double_t x, Double_t par)
Compute the Poisson distribution function for (x,par).
Definition: TMath.cxx:564
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition: TMathBase.h:362
static long int sum(long int i)
Definition: Factory.cxx:2275