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
21
22Double_t const TConfidenceLevel::fgMCLM2S = 0.025;
23Double_t const TConfidenceLevel::fgMCLM1S = 0.16;
24Double_t const TConfidenceLevel::fgMCLMED = 0.5;
25Double_t const TConfidenceLevel::fgMCLP1S = 0.84;
26Double_t const TConfidenceLevel::fgMCLP2S = 0.975;
27// LHWG "one-sided" definition
28Double_t const TConfidenceLevel::fgMCL3S1S = 2.6998E-3;
29Double_t const TConfidenceLevel::fgMCL5S1S = 5.7330E-7;
30// the other definition (not chosen by the LHWG)
31Double_t const TConfidenceLevel::fgMCL3S2S = 1.349898E-3;
32Double_t const TConfidenceLevel::fgMCL5S2S = 2.866516E-7;
33
34
35////////////////////////////////////////////////////////////////////////////////
36/// Default constructor
37
38TConfidenceLevel::TConfidenceLevel()
39{
40 fStot = 0;
41 fBtot = 0;
42 fDtot = 0;
43 fTSD = 0;
44 fTSB = nullptr;
45 fTSS = nullptr;
46 fLRS = nullptr;
47 fLRB = nullptr;
48 fNMC = 0;
49 fNNMC = 0;
50 fISS = nullptr;
51 fISB = nullptr;
52 fMCL3S = fgMCL3S1S;
53 fMCL5S = fgMCL5S1S;
54}
55
56
57////////////////////////////////////////////////////////////////////////////////
58/// Constructor that fix some conventions
59/// \param mc is the number of Monte Carlo experiments
60/// \param onesided specifies if the intervals are one-sided or not.
61
62TConfidenceLevel::TConfidenceLevel(Int_t mc, bool onesided)
63{
64 fStot = 0;
65 fBtot = 0;
66 fDtot = 0;
67 fTSD = 0;
68 fTSB = nullptr;
69 fTSS = nullptr;
70 fLRS = nullptr;
71 fLRB = nullptr;
72 fNMC = mc;
73 fNNMC = mc;
74 fISS = new Int_t[mc];
75 fISB = new Int_t[mc];
76 fMCL3S = onesided ? fgMCL3S1S : fgMCL3S2S;
77 fMCL5S = onesided ? fgMCL5S1S : fgMCL5S2S;
78}
79
80
81////////////////////////////////////////////////////////////////////////////////
82/// The destructor
83
84TConfidenceLevel::~TConfidenceLevel()
85{
86 if (fISS)
87 delete[]fISS;
88 if (fISB)
89 delete[]fISB;
90 if (fTSB)
91 delete[]fTSB;
92 if (fTSS)
93 delete[]fTSS;
94 if (fLRS)
95 delete[]fLRS;
96 if (fLRB)
97 delete[]fLRB;
98}
99
100
101////////////////////////////////////////////////////////////////////////////////
102/// Get the expected statistic value in the background only hypothesis
103
104Double_t TConfidenceLevel::GetExpectedStatistic_b(Int_t sigma) const
105{
106 switch (sigma) {
107 case -2:
108 return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]]) - fStot));
109 case -1:
110 return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]]) - fStot));
111 case 0:
112 return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]]) - fStot));
113 case 1:
114 return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]]) - fStot));
115 case 2:
116 return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]]) - fStot));
117 default:
118 return 0;
119 }
120}
121
122
123////////////////////////////////////////////////////////////////////////////////
124/// Get the expected statistic value in the signal plus background hypothesis
125
126Double_t TConfidenceLevel::GetExpectedStatistic_sb(Int_t sigma) const
127{
128 switch (sigma) {
129 case -2:
130 return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]]) - fStot));
131 case -1:
132 return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]]) - fStot));
133 case 0:
134 return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]]) - fStot));
135 case 1:
136 return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]]) - fStot));
137 case 2:
138 return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]]) - fStot));
139 default:
140 return 0;
141 }
142}
143
144
145////////////////////////////////////////////////////////////////////////////////
146/// Get the Confidence Level for the background only
147
148Double_t TConfidenceLevel::CLb(bool use_sMC) const
149{
150 Double_t result = 0;
151 if (use_sMC) {
152 for (Int_t i = 0; i < fNMC; i++)
153 if (fTSS[fISS[i]] < fTSD)
154 result += (1 / (fLRS[fISS[i]] * fNMC));
155 } else {
156 for (Int_t i = 0; i < fNMC; i++)
157 if (fTSB[fISB[i]] < fTSD)
158 result = (Double_t(i + 1)) / fNMC;
159 }
160 return result;
161}
162
163
164////////////////////////////////////////////////////////////////////////////////
165/// Get the Confidence Level for the signal plus background hypothesis
166
167Double_t TConfidenceLevel::CLsb(bool use_sMC) const
168{
169 Double_t result = 0;
170 if (use_sMC) {
171 for (Int_t i = 0; i < fNMC; i++)
172 if (fTSS[fISS[i]] <= fTSD)
173 result = i / fNMC;
174 } else {
175 for (Int_t i = 0; i < fNMC; i++)
176 if (fTSB[fISB[i]] <= fTSD)
177 result += (fLRB[fISB[i]]) / fNMC;
178 }
179 return result;
180}
181
182
183////////////////////////////////////////////////////////////////////////////////
184/// Get the Confidence Level defined by CLs = CLsb/CLb.
185/// This quantity is stable w.r.t. background fluctuations.
186
187Double_t TConfidenceLevel::CLs(bool use_sMC) const
188{
189 Double_t clb = CLb(kFALSE);
190 Double_t clsb = CLsb(use_sMC);
191 if(clb==0) { std::cout << "Warning: clb = 0 !" << std::endl; return 0;}
192 else return clsb/clb;
193}
194
195
196////////////////////////////////////////////////////////////////////////////////
197/// Get the expected Confidence Level for the signal plus background hypothesis
198/// if there is only background.
199
200Double_t TConfidenceLevel::GetExpectedCLsb_b(Int_t sigma) const
201{
202 Double_t result = 0;
203 switch (sigma) {
204 case -2:
205 {
206 for (Int_t i = 0; i < fNMC; i++)
207 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
208 result += fLRB[fISB[i]] / fNMC;
209 return result;
210 }
211 case -1:
212 {
213 for (Int_t i = 0; i < fNMC; i++)
214 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
215 result += fLRB[fISB[i]] / fNMC;
216 return result;
217 }
218 case 0:
219 {
220 for (Int_t i = 0; i < fNMC; i++)
221 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
222 result += fLRB[fISB[i]] / fNMC;
223 return result;
224 }
225 case 1:
226 {
227 for (Int_t i = 0; i < fNMC; i++)
228 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
229 result += fLRB[fISB[i]] / fNMC;
230 return result;
231 }
232 case 2:
233 {
234 for (Int_t i = 0; i < fNMC; i++)
235 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
236 result += fLRB[fISB[i]] / fNMC;
237 return result;
238 }
239 default:
240 return 0;
241 }
242}
243
244
245////////////////////////////////////////////////////////////////////////////////
246/// Get the expected Confidence Level for the background only
247/// if there is signal and background.
248
249Double_t TConfidenceLevel::GetExpectedCLb_sb(Int_t sigma) const
250{
251 Double_t result = 0;
252 switch (sigma) {
253 case 2:
254 {
255 for (Int_t i = 0; i < fNMC; i++)
256 if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
257 result += fLRS[fISS[i]] / fNMC;
258 return result;
259 }
260 case 1:
261 {
262 for (Int_t i = 0; i < fNMC; i++)
263 if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
264 result += fLRS[fISS[i]] / fNMC;
265 return result;
266 }
267 case 0:
268 {
269 for (Int_t i = 0; i < fNMC; i++)
270 if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
271 result += fLRS[fISS[i]] / fNMC;
272 return result;
273 }
274 case -1:
275 {
276 for (Int_t i = 0; i < fNMC; i++)
277 if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
278 result += fLRS[fISS[i]] / fNMC;
279 return result;
280 }
281 case -2:
282 {
283 for (Int_t i = 0; i < fNMC; i++)
284 if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
285 result += fLRS[fISS[i]] / fNMC;
286 return result;
287 }
288 default:
289 return 0;
290 }
291}
292
293
294////////////////////////////////////////////////////////////////////////////////
295/// Get the expected Confidence Level for the background only
296/// if there is only background.
297
298Double_t TConfidenceLevel::GetExpectedCLb_b(Int_t sigma) const
299{
300 Double_t result = 0;
301 switch (sigma) {
302 case 2:
303 {
304 for (Int_t i = 0; i < fNMC; i++)
305 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
306 result = (i + 1) / double (fNMC);
307 return result;
308 }
309 case 1:
310 {
311 for (Int_t i = 0; i < fNMC; i++)
312 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
313 result = (i + 1) / double (fNMC);
314 return result;
315 }
316 case 0:
317 {
318 for (Int_t i = 0; i < fNMC; i++)
319 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
320 result = (i + 1) / double (fNMC);
321 return result;
322 }
323 case -1:
324 {
325 for (Int_t i = 0; i < fNMC; i++)
326 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
327 result = (i + 1) / double (fNMC);
328 return result;
329 }
330 case -2:
331 {
332 for (Int_t i = 0; i < fNMC; i++)
333 if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
334 result = (i + 1) / double (fNMC);
335 return result;
336 }
337 }
338 return result;
339}
340
341
342////////////////////////////////////////////////////////////////////////////////
343/// Get average CLsb.
344
345Double_t TConfidenceLevel::GetAverageCLsb() const
346{
347 Double_t result = 0;
348 Double_t psumsb = 0;
349 for (Int_t i = 0; i < fNMC; i++) {
350 psumsb += fLRB[fISB[i]] / fNMC;
351 result += psumsb / fNMC;
352 }
353 return result;
354}
355
356
357////////////////////////////////////////////////////////////////////////////////
358/// Get average CLs.
359
360Double_t TConfidenceLevel::GetAverageCLs() const
361{
362 Double_t result = 0;
363 Double_t psumsb = 0;
364 for (Int_t i = 0; i < fNMC; i++) {
365 psumsb += fLRB[fISB[i]] / fNMC;
366 result += ((psumsb / fNMC) / ((i + 1) / fNMC));
367 }
368 return result;
369}
370
371
372////////////////////////////////////////////////////////////////////////////////
373/// Get 3s probability.
374
375Double_t TConfidenceLevel::Get3sProbability() const
376{
377 Double_t result = 0;
378 Double_t psumbs = 0;
379 for (Int_t i = 0; i < fNMC; i++) {
380 psumbs += 1 / (Double_t) (fLRS[(fISS[(Int_t) (fNMC - i)])] * fNMC);
381 if (psumbs <= fMCL3S)
382 result = i / fNMC;
383 }
384 return result;
385}
386
387
388////////////////////////////////////////////////////////////////////////////////
389/// Get 5s probability.
390
391Double_t TConfidenceLevel::Get5sProbability() const
392{
393 Double_t result = 0;
394 Double_t psumbs = 0;
395 for (Int_t i = 0; i < fNMC; i++) {
396 psumbs += 1 / (Double_t) (fLRS[(fISS[(Int_t) (fNMC - i)])] * fNMC);
397 if (psumbs <= fMCL5S)
398 result = i / fNMC;
399 }
400 return result;
401}
402
403
404////////////////////////////////////////////////////////////////////////////////
405/// Display sort of a "canonical" -2lnQ plot.
406/// This results in a plot with 2 elements:
407///
408/// - The histogram of -2lnQ for background hypothesis (full)
409/// - The histogram of -2lnQ for signal and background hypothesis (dashed)
410///
411/// The 2 histograms are respectively named b_hist and sb_hist.
412
413void TConfidenceLevel::Draw(const Option_t*)
414{
415 TH1F h("TConfidenceLevel_Draw","",50,0,0);
416 Int_t i;
417 for (i=0; i<fNMC; i++) {
418 h.Fill(-2*(fTSB[i]-fStot));
419 h.Fill(-2*(fTSS[i]-fStot));
420 }
421 TH1F* b_hist = new TH1F("b_hist", "-2lnQ",50,h.GetXaxis()->GetXmin(),h.GetXaxis()->GetXmax());
422 TH1F* sb_hist = new TH1F("sb_hist","-2lnQ",50,h.GetXaxis()->GetXmin(),h.GetXaxis()->GetXmax());
423 for (i=0; i<fNMC; i++) {
424 b_hist->Fill(-2*(fTSB[i]-fStot));
425 sb_hist->Fill(-2*(fTSS[i]-fStot));
426 }
427 b_hist->Draw();
428 sb_hist->Draw("Same");
429 sb_hist->SetLineStyle(3);
430}
431
432
433////////////////////////////////////////////////////////////////////////////////
434/// Set the TSB.
435
436void TConfidenceLevel::SetTSB(Double_t * in)
437{
438 fTSB = in;
439 TMath::Sort(fNNMC, fTSB, fISB, false);
440}
441
442
443////////////////////////////////////////////////////////////////////////////////
444/// Set the TSS.
445
446void TConfidenceLevel::SetTSS(Double_t * in)
447{
448 fTSS = in;
449 TMath::Sort(fNNMC, fTSS, fISS, false);
450}
h
#define h(i)
Definition RSha256.hxx:106
kFALSE
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
Double_t
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
Option_t
const char Option_t
Option string (const char)
Definition RtypesCore.h:80
TRangeDynCast
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Definition TCollection.h:358
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
TConfidenceLevel::fgMCL5S2S
static const Double_t fgMCL5S2S
Definition TConfidenceLevel.h:76
TConfidenceLevel::SetTSS
void SetTSS(Double_t *in)
Set the TSS.
Definition TConfidenceLevel.cxx:446
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:249
TConfidenceLevel::SetTSB
void SetTSB(Double_t *in)
Set the TSB.
Definition TConfidenceLevel.cxx:436
TConfidenceLevel::Get3sProbability
Double_t Get3sProbability() const
Get 3s probability.
Definition TConfidenceLevel.cxx:375
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:413
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:104
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:38
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:126
TConfidenceLevel::GetAverageCLs
Double_t GetAverageCLs() const
Get average CLs.
Definition TConfidenceLevel.cxx:360
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:298
TConfidenceLevel::fgMCLM1S
static const Double_t fgMCLM1S
Definition TConfidenceLevel.h:69
TConfidenceLevel::~TConfidenceLevel
~TConfidenceLevel() override
The destructor.
Definition TConfidenceLevel.cxx:84
TConfidenceLevel::CLsb
Double_t CLsb(bool use_sMC=kFALSE) const
Get the Confidence Level for the signal plus background hypothesis.
Definition TConfidenceLevel.cxx:167
TConfidenceLevel::fgMCL3S1S
static const Double_t fgMCL3S1S
Definition TConfidenceLevel.h:73
TConfidenceLevel::Get5sProbability
Double_t Get5sProbability() const
Get 5s probability.
Definition TConfidenceLevel.cxx:391
TConfidenceLevel::CLb
Double_t CLb(bool use_sMC=kFALSE) const
Get the Confidence Level for the background only.
Definition TConfidenceLevel.cxx:148
TConfidenceLevel::fgMCL5S1S
static const Double_t fgMCL5S1S
Definition TConfidenceLevel.h:74
TConfidenceLevel::GetAverageCLsb
Double_t GetAverageCLsb() const
Get average CLsb.
Definition TConfidenceLevel.cxx:345
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:200
TConfidenceLevel::CLs
Double_t CLs(bool use_sMC=kFALSE) const
Get the Confidence Level defined by CLs = CLsb/CLb.
Definition TConfidenceLevel.cxx:187
TH1F
1-D histogram with a float per channel (see TH1 documentation)
Definition TH1.h:879
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:251
TMath::Min
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:199
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:432