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
21ClassImp(TConfidenceLevel);
22
23Double_t const TConfidenceLevel::fgMCLM2S = 0.025;
24Double_t const TConfidenceLevel::fgMCLM1S = 0.16;
25Double_t const TConfidenceLevel::fgMCLMED = 0.5;
26Double_t const TConfidenceLevel::fgMCLP1S = 0.84;
27Double_t const TConfidenceLevel::fgMCLP2S = 0.975;
28// LHWG "one-sided" definition
29Double_t const TConfidenceLevel::fgMCL3S1S = 2.6998E-3;
30Double_t const TConfidenceLevel::fgMCL5S1S = 5.7330E-7;
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
39TConfidenceLevel::TConfidenceLevel()
40{
41 fStot = 0;
42 fBtot = 0;
43 fDtot = 0;
44 fTSD = 0;
45 fTSB = nullptr;
46 fTSS = nullptr;
47 fLRS = nullptr;
48 fLRB = nullptr;
49 fNMC = 0;
50 fNNMC = 0;
51 fISS = nullptr;
52 fISB = nullptr;
53 fMCL3S = fgMCL3S1S;
54 fMCL5S = fgMCL5S1S;
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
63TConfidenceLevel::TConfidenceLevel(Int_t mc, bool onesided)
64{
65 fStot = 0;
66 fBtot = 0;
67 fDtot = 0;
68 fTSD = 0;
69 fTSB = nullptr;
70 fTSS = nullptr;
71 fLRS = nullptr;
72 fLRB = nullptr;
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
85TConfidenceLevel::~TConfidenceLevel()
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
105Double_t TConfidenceLevel::GetExpectedStatistic_b(Int_t sigma) const
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
127Double_t TConfidenceLevel::GetExpectedStatistic_sb(Int_t sigma) const
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
149Double_t TConfidenceLevel::CLb(bool use_sMC) const
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
168Double_t TConfidenceLevel::CLsb(bool use_sMC) const
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
188Double_t TConfidenceLevel::CLs(bool use_sMC) const
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
201Double_t TConfidenceLevel::GetExpectedCLsb_b(Int_t sigma) const
202{
203 Double_t result = 0;
204 switch (sigma) {
205 case -2:
206 {
207 for (Int_t i = 0; i < fNMC; i++)
208 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
209 result += fLRB[fISB[i]] / fNMC;
210 return result;
211 }
212 case -1:
213 {
214 for (Int_t i = 0; i < fNMC; i++)
215 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
216 result += fLRB[fISB[i]] / fNMC;
217 return result;
218 }
219 case 0:
220 {
221 for (Int_t i = 0; i < fNMC; i++)
222 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
223 result += fLRB[fISB[i]] / fNMC;
224 return result;
225 }
226 case 1:
227 {
228 for (Int_t i = 0; i < fNMC; i++)
229 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
230 result += fLRB[fISB[i]] / fNMC;
231 return result;
232 }
233 case 2:
234 {
235 for (Int_t i = 0; i < fNMC; i++)
236 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
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
250Double_t TConfidenceLevel::GetExpectedCLb_sb(Int_t sigma) const
251{
252 Double_t result = 0;
253 switch (sigma) {
254 case 2:
255 {
256 for (Int_t i = 0; i < fNMC; i++)
257 if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
258 result += fLRS[fISS[i]] / fNMC;
259 return result;
260 }
261 case 1:
262 {
263 for (Int_t i = 0; i < fNMC; i++)
264 if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
265 result += fLRS[fISS[i]] / fNMC;
266 return result;
267 }
268 case 0:
269 {
270 for (Int_t i = 0; i < fNMC; i++)
271 if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
272 result += fLRS[fISS[i]] / fNMC;
273 return result;
274 }
275 case -1:
276 {
277 for (Int_t i = 0; i < fNMC; i++)
278 if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
279 result += fLRS[fISS[i]] / fNMC;
280 return result;
281 }
282 case -2:
283 {
284 for (Int_t i = 0; i < fNMC; i++)
285 if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
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
299Double_t TConfidenceLevel::GetExpectedCLb_b(Int_t sigma) const
300{
301 Double_t result = 0;
302 switch (sigma) {
303 case 2:
304 {
305 for (Int_t i = 0; i < fNMC; i++)
306 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
307 result = (i + 1) / double (fNMC);
308 return result;
309 }
310 case 1:
311 {
312 for (Int_t i = 0; i < fNMC; i++)
313 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
314 result = (i + 1) / double (fNMC);
315 return result;
316 }
317 case 0:
318 {
319 for (Int_t i = 0; i < fNMC; i++)
320 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
321 result = (i + 1) / double (fNMC);
322 return result;
323 }
324 case -1:
325 {
326 for (Int_t i = 0; i < fNMC; i++)
327 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
328 result = (i + 1) / double (fNMC);
329 return result;
330 }
331 case -2:
332 {
333 for (Int_t i = 0; i < fNMC; i++)
334 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
335 result = (i + 1) / double (fNMC);
336 return result;
337 }
338 }
339 return result;
340}
341
342
343////////////////////////////////////////////////////////////////////////////////
344/// Get average CLsb.
345
346Double_t TConfidenceLevel::GetAverageCLsb() const
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
361Double_t TConfidenceLevel::GetAverageCLs() const
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
376Double_t TConfidenceLevel::Get3sProbability() const
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
392Double_t TConfidenceLevel::Get5sProbability() const
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
414void TConfidenceLevel::Draw(const Option_t*)
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
437void TConfidenceLevel::SetTSB(Double_t * in)
438{
439 fTSB = in;
440 TMath::Sort(fNNMC, fTSB, fISB, false);
441}
442
443
444////////////////////////////////////////////////////////////////////////////////
445/// Set the TSS.
446
447void TConfidenceLevel::SetTSS(Double_t * in)
448{
449 fTSS = in;
450 TMath::Sort(fNNMC, fTSS, fISS, false);
451}
h
#define h(i)
Definition RSha256.hxx:106
kFALSE
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
Double_t
double Double_t
Definition RtypesCore.h:59
Option_t
const char Option_t
Definition RtypesCore.h:66
ClassImp
#define ClassImp(name)
Definition Rtypes.h:377
result
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
Definition TGWin32VirtualXProxy.cxx:174
TAttLine::SetLineStyle
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition TAttLine.h:42
TConfidenceLevel
Class to compute 95% CL limits.
Definition TConfidenceLevel.h:20
TConfidenceLevel::fgMCL5S2S
static const Double_t fgMCL5S2S
Definition TConfidenceLevel.h:76
TConfidenceLevel::SetTSS
void SetTSS(Double_t *in)
Set the TSS.
Definition TConfidenceLevel.cxx:447
TConfidenceLevel::fgMCLM2S
static const Double_t fgMCLM2S
Definition TConfidenceLevel.h:68
TConfidenceLevel::GetExpectedCLb_sb
Double_t GetExpectedCLb_sb(Int_t sigma=0) const
Get the expected Confidence Level for the background only if there is signal and background.
Definition TConfidenceLevel.cxx:250
TConfidenceLevel::SetTSB
void SetTSB(Double_t *in)
Set the TSB.
Definition TConfidenceLevel.cxx:437
TConfidenceLevel::Get3sProbability
Double_t Get3sProbability() const
Get 3s probability.
Definition TConfidenceLevel.cxx:376
TConfidenceLevel::fgMCLP2S
static const Double_t fgMCLP2S
Definition TConfidenceLevel.h:72
TConfidenceLevel::Draw
void Draw(const Option_t *option="") override
Display sort of a "canonical" -2lnQ plot.
Definition TConfidenceLevel.cxx:414
TConfidenceLevel::fgMCL3S2S
static const Double_t fgMCL3S2S
Definition TConfidenceLevel.h:75
TConfidenceLevel::GetExpectedStatistic_b
Double_t GetExpectedStatistic_b(Int_t sigma=0) const
Get the expected statistic value in the background only hypothesis.
Definition TConfidenceLevel.cxx:105
TConfidenceLevel::fgMCLP1S
static const Double_t fgMCLP1S
Definition TConfidenceLevel.h:71
TConfidenceLevel::fgMCLMED
static const Double_t fgMCLMED
Definition TConfidenceLevel.h:70
TConfidenceLevel::TConfidenceLevel
TConfidenceLevel()
Default constructor.
Definition TConfidenceLevel.cxx:39
TConfidenceLevel::GetExpectedStatistic_sb
Double_t GetExpectedStatistic_sb(Int_t sigma=0) const
Get the expected statistic value in the signal plus background hypothesis.
Definition TConfidenceLevel.cxx:127
TConfidenceLevel::GetAverageCLs
Double_t GetAverageCLs() const
Get average CLs.
Definition TConfidenceLevel.cxx:361
TConfidenceLevel::GetExpectedCLb_b
Double_t GetExpectedCLb_b(Int_t sigma=0) const
Get the expected Confidence Level for the background only if there is only background.
Definition TConfidenceLevel.cxx:299
TConfidenceLevel::fgMCLM1S
static const Double_t fgMCLM1S
Definition TConfidenceLevel.h:69
TConfidenceLevel::~TConfidenceLevel
~TConfidenceLevel() override
The destructor.
Definition TConfidenceLevel.cxx:85
TConfidenceLevel::CLsb
Double_t CLsb(bool use_sMC=kFALSE) const
Get the Confidence Level for the signal plus background hypothesis.
Definition TConfidenceLevel.cxx:168
TConfidenceLevel::fgMCL3S1S
static const Double_t fgMCL3S1S
Definition TConfidenceLevel.h:73
TConfidenceLevel::Get5sProbability
Double_t Get5sProbability() const
Get 5s probability.
Definition TConfidenceLevel.cxx:392
TConfidenceLevel::CLb
Double_t CLb(bool use_sMC=kFALSE) const
Get the Confidence Level for the background only.
Definition TConfidenceLevel.cxx:149
TConfidenceLevel::fgMCL5S1S
static const Double_t fgMCL5S1S
Definition TConfidenceLevel.h:74
TConfidenceLevel::GetAverageCLsb
Double_t GetAverageCLsb() const
Get average CLsb.
Definition TConfidenceLevel.cxx:346
TConfidenceLevel::GetExpectedCLsb_b
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...
Definition TConfidenceLevel.cxx:201
TConfidenceLevel::CLs
Double_t CLs(bool use_sMC=kFALSE) const
Get the Confidence Level defined by CLs = CLsb/CLb.
Definition TConfidenceLevel.cxx:188
TH1F
1-D histogram with a float per channel (see TH1 documentation)}
Definition TH1.h:577
TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3345
TH1::Draw
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3067
sigma
const Double_t sigma
Definition h1analysisProxy.h:11
TMath::Max
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:250
TMath::Min
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:198
TMath::Sort
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Sort the n elements of the array a of generic templated type Element.
Definition TMathBase.h:431