Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 ClassImp(TConfidenceLevel);
22 
23 Double_t const TConfidenceLevel::fgMCLM2S = 0.025;
24 Double_t const TConfidenceLevel::fgMCLM1S = 0.16;
25 Double_t const TConfidenceLevel::fgMCLMED = 0.5;
26 Double_t const TConfidenceLevel::fgMCLP1S = 0.84;
27 Double_t const TConfidenceLevel::fgMCLP2S = 0.975;
28 // LHWG "one-sided" definition
29 Double_t const TConfidenceLevel::fgMCL3S1S = 2.6998E-3;
30 Double_t const TConfidenceLevel::fgMCL5S1S = 5.7330E-7;
31 // the other definition (not chosen by the LHWG)
32 Double_t const TConfidenceLevel::fgMCL3S2S = 1.349898E-3;
33 Double_t const TConfidenceLevel::fgMCL5S2S = 2.866516E-7;
34 
35 
36 ////////////////////////////////////////////////////////////////////////////////
37 /// Default constructor
38 
39 TConfidenceLevel::TConfidenceLevel()
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;
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 
63 TConfidenceLevel::TConfidenceLevel(Int_t mc, bool onesided)
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 
85 TConfidenceLevel::~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 
105 Double_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 
127 Double_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 
149 Double_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 
168 Double_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 
188 Double_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 
201 Double_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 
250 Double_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 
299 Double_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 
346 Double_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 
361 Double_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 
376 Double_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 
392 Double_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 
414 void 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 
437 void TConfidenceLevel::SetTSB(Double_t * in)
438 {
439  fTSB = in;
440  TMath::Sort(fNNMC, fTSB, fISB, 0);
441 }
442 
443 
444 ////////////////////////////////////////////////////////////////////////////////
445 /// Set the TSS.
446 
447 void TConfidenceLevel::SetTSS(Double_t * in)
448 {
449  fTSS = in;
450  TMath::Sort(fNNMC, fTSS, fISS, 0);
451 }