Logo ROOT   6.10/09
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 
16 Class to calculate the CL upper limit using
17 the Feldman-Cousins method as described in PRD V57 #7, p3873-3889
18 
19 The default confidence interval calculated using this method is 90%
20 This is set either by having a default the constructor, or using the
21 appropriate fraction when instantiating an object of this class (e.g. 0.9)
22 
23 The simple extension to a gaussian resolution function bounded at zero
24 has not been addressed as yet -> `time is of the essence' as they write
25 on 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,
29 fMuMax, fMuStep and fNMax are those used in the PRD:
30 ~~~ {.cpp}
31  fMuMin = 0.0
32  fMuMax = 50.0
33  fMuStep= 0.005
34 ~~~
35 but there is total flexibility in changing this should you desire.
36 
37 
38 see example of use in $ROOTSYS/tutorials/math/FeldmanCousins.C
39 
40 see note about: "Should I use TRolke, TFeldmanCousins, TLimit?"
41  in the TRolke class description.
42 
43 \author: Adrian Bevan, Liverpool University
44 
45 Copyright 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 }
static double B[]
void SetMuStep(Double_t newMuStep=0.005)
Set the step in signal to use when generating tables.
static long int sum(long int i)
Definition: Factory.cxx:2162
Double_t fNbackground
#define N
virtual ~TFeldmanCousins()
Basic string class.
Definition: TString.h:129
int Int_t
Definition: RtypesCore.h:41
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition: TMath.h:1151
TRandom2 r(17)
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...
Double_t Poisson(Double_t x, Double_t par)
Compute the Poisson distribution function for (x,par) The Poisson PDF is implemented by means of Eule...
Definition: TMath.cxx:571
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...
void SetMuMax(Double_t newMax=50.0)
Set maximum value of signal to use in calculating the tables.
#define ClassImp(name)
Definition: Rtypes.h:336
double Double_t
Definition: RtypesCore.h:55
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...
Class to calculate the CL upper limit using the Feldman-Cousins method as described in PRD V57 #7...
Double_t fLowerLimit
Double_t fUpperLimit
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...