Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TConfidenceLevel.cxx
Go to the documentation of this file.
1// @(#)root/hist:$Id$
2// Author: Christophe.Delaere@cern.ch 21/08/2002
3
4////////////////////////////////////////////////////////////////////////////////
5/** \class TConfidenceLevel
6 \ingroup Hist
7 \brief Class to compute 95% CL limits
8*///////////////////////////////////////////////////////////////////////////////
9
10/*************************************************************************
11 * C.Delaere *
12 * adapted from the mclimit code from Tom Junk *
13 * see http://cern.ch/thomasj/searchlimits/ecl.html *
14 *************************************************************************/
15
16#include "TConfidenceLevel.h"
17#include "TH1F.h"
18#include "TMath.h"
19#include "Riostream.h"
20
22
28// LHWG "one-sided" definition
31// the other definition (not chosen by the LHWG)
32Double_t const TConfidenceLevel::fgMCL3S2S = 1.349898E-3;
33Double_t const TConfidenceLevel::fgMCL5S2S = 2.866516E-7;
34
35
36////////////////////////////////////////////////////////////////////////////////
37/// Default constructor
38
40{
41 fStot = 0;
42 fBtot = 0;
43 fDtot = 0;
44 fTSD = 0;
45 fTSB = 0;
46 fTSS = 0;
47 fLRS = 0;
48 fLRB = 0;
49 fNMC = 0;
50 fNNMC = 0;
51 fISS = 0;
52 fISB = 0;
55}
56
57
58////////////////////////////////////////////////////////////////////////////////
59/// Constructor that fix some conventions
60/// \param mc is the number of Monte Carlo experiments
61/// \param onesided specifies if the intervals are one-sided or not.
62
64{
65 fStot = 0;
66 fBtot = 0;
67 fDtot = 0;
68 fTSD = 0;
69 fTSB = 0;
70 fTSS = 0;
71 fLRS = 0;
72 fLRB = 0;
73 fNMC = mc;
74 fNNMC = mc;
75 fISS = new Int_t[mc];
76 fISB = new Int_t[mc];
77 fMCL3S = onesided ? fgMCL3S1S : fgMCL3S2S;
78 fMCL5S = onesided ? fgMCL5S1S : fgMCL5S2S;
79}
80
81
82////////////////////////////////////////////////////////////////////////////////
83/// The destructor
84
86{
87 if (fISS)
88 delete[]fISS;
89 if (fISB)
90 delete[]fISB;
91 if (fTSB)
92 delete[]fTSB;
93 if (fTSS)
94 delete[]fTSS;
95 if (fLRS)
96 delete[]fLRS;
97 if (fLRB)
98 delete[]fLRB;
99}
100
101
102////////////////////////////////////////////////////////////////////////////////
103/// Get the expected statistic value in the background only hypothesis
104
106{
107 switch (sigma) {
108 case -2:
109 return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]]) - fStot));
110 case -1:
111 return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]]) - fStot));
112 case 0:
113 return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]]) - fStot));
114 case 1:
115 return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]]) - fStot));
116 case 2:
117 return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]]) - fStot));
118 default:
119 return 0;
120 }
121}
122
123
124////////////////////////////////////////////////////////////////////////////////
125/// Get the expected statistic value in the signal plus background hypothesis
126
128{
129 switch (sigma) {
130 case -2:
131 return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]]) - fStot));
132 case -1:
133 return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]]) - fStot));
134 case 0:
135 return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]]) - fStot));
136 case 1:
137 return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]]) - fStot));
138 case 2:
139 return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]]) - fStot));
140 default:
141 return 0;
142 }
143}
144
145
146////////////////////////////////////////////////////////////////////////////////
147/// Get the Confidence Level for the background only
148
150{
151 Double_t result = 0;
152 if (use_sMC) {
153 for (Int_t i = 0; i < fNMC; i++)
154 if (fTSS[fISS[i]] < fTSD)
155 result += (1 / (fLRS[fISS[i]] * fNMC));
156 } else {
157 for (Int_t i = 0; i < fNMC; i++)
158 if (fTSB[fISB[i]] < fTSD)
159 result = (Double_t(i + 1)) / fNMC;
160 }
161 return result;
162}
163
164
165////////////////////////////////////////////////////////////////////////////////
166/// Get the Confidence Level for the signal plus background hypothesis
167
169{
170 Double_t result = 0;
171 if (use_sMC) {
172 for (Int_t i = 0; i < fNMC; i++)
173 if (fTSS[fISS[i]] <= fTSD)
174 result = i / fNMC;
175 } else {
176 for (Int_t i = 0; i < fNMC; i++)
177 if (fTSB[fISB[i]] <= fTSD)
178 result += (fLRB[fISB[i]]) / fNMC;
179 }
180 return result;
181}
182
183
184////////////////////////////////////////////////////////////////////////////////
185/// Get the Confidence Level defined by CLs = CLsb/CLb.
186/// This quantity is stable w.r.t. background fluctuations.
187
189{
190 Double_t clb = CLb(kFALSE);
191 Double_t clsb = CLsb(use_sMC);
192 if(clb==0) { std::cout << "Warning: clb = 0 !" << std::endl; return 0;}
193 else return clsb/clb;
194}
195
196
197////////////////////////////////////////////////////////////////////////////////
198/// Get the expected Confidence Level for the signal plus background hypothesis
199/// if there is only background.
200
202{
203 Double_t result = 0;
204 switch (sigma) {
205 case -2:
206 {
207 for (Int_t i = 0; i < fNMC; i++)
209 result += fLRB[fISB[i]] / fNMC;
210 return result;
211 }
212 case -1:
213 {
214 for (Int_t i = 0; i < fNMC; i++)
216 result += fLRB[fISB[i]] / fNMC;
217 return result;
218 }
219 case 0:
220 {
221 for (Int_t i = 0; i < fNMC; i++)
223 result += fLRB[fISB[i]] / fNMC;
224 return result;
225 }
226 case 1:
227 {
228 for (Int_t i = 0; i < fNMC; i++)
230 result += fLRB[fISB[i]] / fNMC;
231 return result;
232 }
233 case 2:
234 {
235 for (Int_t i = 0; i < fNMC; i++)
237 result += fLRB[fISB[i]] / fNMC;
238 return result;
239 }
240 default:
241 return 0;
242 }
243}
244
245
246////////////////////////////////////////////////////////////////////////////////
247/// Get the expected Confidence Level for the background only
248/// if there is signal and background.
249
251{
252 Double_t result = 0;
253 switch (sigma) {
254 case 2:
255 {
256 for (Int_t i = 0; i < fNMC; i++)
258 result += fLRS[fISS[i]] / fNMC;
259 return result;
260 }
261 case 1:
262 {
263 for (Int_t i = 0; i < fNMC; i++)
265 result += fLRS[fISS[i]] / fNMC;
266 return result;
267 }
268 case 0:
269 {
270 for (Int_t i = 0; i < fNMC; i++)
272 result += fLRS[fISS[i]] / fNMC;
273 return result;
274 }
275 case -1:
276 {
277 for (Int_t i = 0; i < fNMC; i++)
279 result += fLRS[fISS[i]] / fNMC;
280 return result;
281 }
282 case -2:
283 {
284 for (Int_t i = 0; i < fNMC; i++)
286 result += fLRS[fISS[i]] / fNMC;
287 return result;
288 }
289 default:
290 return 0;
291 }
292}
293
294
295////////////////////////////////////////////////////////////////////////////////
296/// Get the expected Confidence Level for the background only
297/// if there is only background.
298
300{
301 Double_t result = 0;
302 switch (sigma) {
303 case 2:
304 {
305 for (Int_t i = 0; i < fNMC; i++)
307 result = (i + 1) / double (fNMC);
308 return result;
309 }
310 case 1:
311 {
312 for (Int_t i = 0; i < fNMC; i++)
314 result = (i + 1) / double (fNMC);
315 return result;
316 }
317 case 0:
318 {
319 for (Int_t i = 0; i < fNMC; i++)
321 result = (i + 1) / double (fNMC);
322 return result;
323 }
324 case -1:
325 {
326 for (Int_t i = 0; i < fNMC; i++)
328 result = (i + 1) / double (fNMC);
329 return result;
330 }
331 case -2:
332 {
333 for (Int_t i = 0; i < fNMC; i++)
335 result = (i + 1) / double (fNMC);
336 return result;
337 }
338 }
339 return result;
340}
341
342
343////////////////////////////////////////////////////////////////////////////////
344/// Get average CLsb.
345
347{
348 Double_t result = 0;
349 Double_t psumsb = 0;
350 for (Int_t i = 0; i < fNMC; i++) {
351 psumsb += fLRB[fISB[i]] / fNMC;
352 result += psumsb / fNMC;
353 }
354 return result;
355}
356
357
358////////////////////////////////////////////////////////////////////////////////
359/// Get average CLs.
360
362{
363 Double_t result = 0;
364 Double_t psumsb = 0;
365 for (Int_t i = 0; i < fNMC; i++) {
366 psumsb += fLRB[fISB[i]] / fNMC;
367 result += ((psumsb / fNMC) / ((i + 1) / fNMC));
368 }
369 return result;
370}
371
372
373////////////////////////////////////////////////////////////////////////////////
374/// Get 3s probability.
375
377{
378 Double_t result = 0;
379 Double_t psumbs = 0;
380 for (Int_t i = 0; i < fNMC; i++) {
381 psumbs += 1 / (Double_t) (fLRS[(fISS[(Int_t) (fNMC - i)])] * fNMC);
382 if (psumbs <= fMCL3S)
383 result = i / fNMC;
384 }
385 return result;
386}
387
388
389////////////////////////////////////////////////////////////////////////////////
390/// Get 5s probability.
391
393{
394 Double_t result = 0;
395 Double_t psumbs = 0;
396 for (Int_t i = 0; i < fNMC; i++) {
397 psumbs += 1 / (Double_t) (fLRS[(fISS[(Int_t) (fNMC - i)])] * fNMC);
398 if (psumbs <= fMCL5S)
399 result = i / fNMC;
400 }
401 return result;
402}
403
404
405////////////////////////////////////////////////////////////////////////////////
406/// Display sort of a "canonical" -2lnQ plot.
407/// This results in a plot with 2 elements:
408///
409/// - The histogram of -2lnQ for background hypothesis (full)
410/// - The histogram of -2lnQ for signal and background hypothesis (dashed)
411///
412/// The 2 histograms are respectively named b_hist and sb_hist.
413
415{
416 TH1F h("TConfidenceLevel_Draw","",50,0,0);
417 Int_t i;
418 for (i=0; i<fNMC; i++) {
419 h.Fill(-2*(fTSB[i]-fStot));
420 h.Fill(-2*(fTSS[i]-fStot));
421 }
422 TH1F* b_hist = new TH1F("b_hist", "-2lnQ",50,h.GetXaxis()->GetXmin(),h.GetXaxis()->GetXmax());
423 TH1F* sb_hist = new TH1F("sb_hist","-2lnQ",50,h.GetXaxis()->GetXmin(),h.GetXaxis()->GetXmax());
424 for (i=0; i<fNMC; i++) {
425 b_hist->Fill(-2*(fTSB[i]-fStot));
426 sb_hist->Fill(-2*(fTSS[i]-fStot));
427 }
428 b_hist->Draw();
429 sb_hist->Draw("Same");
430 sb_hist->SetLineStyle(3);
431}
432
433
434////////////////////////////////////////////////////////////////////////////////
435/// Set the TSB.
436
438{
439 fTSB = in;
441}
442
443
444////////////////////////////////////////////////////////////////////////////////
445/// Set the TSS.
446
448{
449 fTSS = in;
451}
double
#define h(i)
Definition RSha256.hxx:106
const Bool_t kFALSE
Definition RtypesCore.h:92
double Double_t
Definition RtypesCore.h:59
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:364
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition TAttLine.h:42
Class to compute 95% CL limits.
static const Double_t fgMCL5S2S
void SetTSS(Double_t *in)
Set the TSS.
static const Double_t fgMCLM2S
Double_t GetExpectedCLb_sb(Int_t sigma=0) const
Get the expected Confidence Level for the background only if there is signal and background.
void SetTSB(Double_t *in)
Set the TSB.
Double_t Get3sProbability() const
Get 3s probability.
static const Double_t fgMCLP2S
static const Double_t fgMCL3S2S
Double_t GetExpectedStatistic_b(Int_t sigma=0) const
Get the expected statistic value in the background only hypothesis.
static const Double_t fgMCLP1S
static const Double_t fgMCLMED
TConfidenceLevel()
Default constructor.
Double_t GetExpectedStatistic_sb(Int_t sigma=0) const
Get the expected statistic value in the signal plus background hypothesis.
Double_t GetAverageCLs() const
Get average CLs.
Double_t GetExpectedCLb_b(Int_t sigma=0) const
Get the expected Confidence Level for the background only if there is only background.
static const Double_t fgMCLM1S
Double_t CLsb(bool use_sMC=kFALSE) const
Get the Confidence Level for the signal plus background hypothesis.
static const Double_t fgMCL3S1S
Double_t Get5sProbability() const
Get 5s probability.
virtual ~TConfidenceLevel()
The destructor.
Double_t CLb(bool use_sMC=kFALSE) const
Get the Confidence Level for the background only.
static const Double_t fgMCL5S1S
Double_t GetAverageCLsb() const
Get average CLsb.
Double_t GetExpectedCLsb_b(Int_t sigma=0) const
Get the expected Confidence Level for the signal plus background hypothesis if there is only backgrou...
Double_t CLs(bool use_sMC=kFALSE) const
Get the Confidence Level defined by CLs = CLsb/CLb.
void Draw(const Option_t *option="")
Display sort of a "canonical" -2lnQ plot.
1-D histogram with a float per channel (see TH1 documentation)}
Definition TH1.h:575
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3350
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition TH1.cxx:3073
const Double_t sigma
Short_t Max(Short_t a, Short_t b)
Definition TMathBase.h:212
Short_t Min(Short_t a, Short_t b)
Definition TMathBase.h:180
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition TMathBase.h:362