ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/auterman/scripts/CLs/TConfidenceLevel.cc
Revision: 1.1.1.1 (vendor branch)
Committed: Thu Apr 1 15:14:10 2010 UTC (15 years, 1 month ago) by auterman
Content type: text/plain
Branch: CLs, MAIN
CVS Tags: start, HEAD
Changes since 1.1: +0 -0 lines
Error occurred while calculating annotation data.
Log Message:
some cls demo plots

File Contents

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