scripts/CLs/TConfidenceLevel.cc

// @(#)root/hist:$Id: TConfidenceLevel.cxx 27713 2009-03-07 08:12:02Z brun $
// Author: Christophe.Delaere@cern.ch   21/08/2002

///////////////////////////////////////////////////////////////////////////
//
// TConfidenceLevel
//
// Class to compute 95% CL limits
//
///////////////////////////////////////////////////////////////////////////

/*************************************************************************
 * C.Delaere                                                             *
 * adapted from the mclimit code from Tom Junk                           *
 * see http://cern.ch/thomasj/searchlimits/ecl.html                      *
 *************************************************************************/

#include "TConfidenceLevel.h"
#include "TH1F.h"
#include "TMath.h"
#include "Riostream.h"

ClassImp(TConfidenceLevel)

Double_t const TConfidenceLevel::fgMCLM2S = 0.025;
Double_t const TConfidenceLevel::fgMCLM1S = 0.16;
Double_t const TConfidenceLevel::fgMCLMED = 0.5;
Double_t const TConfidenceLevel::fgMCLP1S = 0.84;
Double_t const TConfidenceLevel::fgMCLP2S = 0.975;
// LHWG "one-sided" definition
Double_t const TConfidenceLevel::fgMCL3S1S = 2.6998E-3;
Double_t const TConfidenceLevel::fgMCL5S1S = 5.7330E-7;
// the other definition (not chosen by the LHWG)
Double_t const TConfidenceLevel::fgMCL3S2S = 1.349898E-3;
Double_t const TConfidenceLevel::fgMCL5S2S = 2.866516E-7;


//______________________________________________________________________________
TConfidenceLevel::TConfidenceLevel()
{
   // Default constructor

   fStot = 0;
   fBtot = 0;
   fDtot = 0;
   fTSD  = 0;
   fTSB  = 0;
   fTSS  = 0;
   fLRS  = 0;
   fLRB  = 0;
   fNMC  = 0;
   fNNMC = 0;
   fISS  = 0;
   fISB  = 0;
   fMCL3S = fgMCL3S1S;
   fMCL5S = fgMCL5S1S;
}


//______________________________________________________________________________
TConfidenceLevel::TConfidenceLevel(Int_t mc, bool onesided)
{
   // a constructor that fix some conventions:
   // mc is the number of Monte Carlo experiments
   // while onesided specifies if the intervals are one-sided or not.

   fStot = 0;
   fBtot = 0;
   fDtot = 0;
   fTSD  = 0;
   fTSB  = 0;
   fTSS  = 0;
   fLRS  = 0;
   fLRB  = 0;
   fNMC  = mc;
   fNNMC = mc;
   fISS  = new Int_t[mc];
   fISB  = new Int_t[mc];
   fMCL3S = onesided ? fgMCL3S1S : fgMCL3S2S;
   fMCL5S = onesided ? fgMCL5S1S : fgMCL5S2S;
}


//______________________________________________________________________________
TConfidenceLevel::~TConfidenceLevel()
{
   // The destructor

   if (fISS)
      delete[]fISS;
   if (fISB)
      delete[]fISB;
   if (fTSB)
      delete[]fTSB;
   if (fTSS)
      delete[]fTSS;
   if (fLRS)
      delete[]fLRS;
   if (fLRB)
      delete[]fLRB;
}


//______________________________________________________________________________
Double_t TConfidenceLevel::GetExpectedStatistic_b(Int_t sigma) const
{
   // Get the expected statistic value in the background only hypothesis

   switch (sigma) {
   case -2:
      return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]]) - fStot));
   case -1:
      return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]]) - fStot));
   case 0:
      return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]]) - fStot));
   case 1:
      return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]]) - fStot));
   case 2:
      return (-2 *((fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]]) - fStot));
   default:
      return 0;
   }
}


//______________________________________________________________________________
Double_t TConfidenceLevel::GetExpectedStatistic_sb(Int_t sigma) const
{
   // Get the expected statistic value in the signal plus background hypothesis

   switch (sigma) {
   case -2:
      return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]]) - fStot));
   case -1:
      return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]]) - fStot));
   case 0:
      return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]]) - fStot));
   case 1:
      return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]]) - fStot));
   case 2:
      return (-2 *((fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]]) - fStot));
   default:
      return 0;
   }
}


//______________________________________________________________________________
Double_t TConfidenceLevel::CLb(bool use_sMC) const
{
   // Get the Confidence Level for the background only

   Double_t result = 0;
   switch (use_sMC) {
   case kFALSE:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSB[fISB[i]] < fTSD)
               result = (Double_t(i + 1)) / fNMC;
         return result;
      }
   case kTRUE:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSS[fISS[i]] < fTSD)
               result += (1 / (fLRS[fISS[i]] * fNMC));
         return result;
      }
   }
   return result;
}


//______________________________________________________________________________
Double_t TConfidenceLevel::CLsb(bool use_sMC) const
{
   // Get the Confidence Level for the signal plus background hypothesis

   Double_t result = 0;
   switch (use_sMC) {
   case kFALSE:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSB[fISB[i]] <= fTSD)
               result += (fLRB[fISB[i]]) / fNMC;
         return result;
      }
   case kTRUE:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSS[fISS[i]] <= fTSD)
               result = i / fNMC;
         return result;
      }
   }
   return result;
}


//______________________________________________________________________________
Double_t TConfidenceLevel::CLs(bool use_sMC) const
{
   // Get the Confidence Level defined by CLs = CLsb/CLb.
   // This quantity is stable w.r.t. background fluctuations.

   Double_t clb = CLb(kFALSE);
   Double_t clsb = CLsb(use_sMC);
   if(clb==0) { cout << "Warning: clb = 0 !" << endl; return 0;}
   else return clsb/clb;
}


//______________________________________________________________________________
Double_t TConfidenceLevel::GetExpectedCLsb_b(Int_t sigma) const
{
   // Get the expected Confidence Level for the signal plus background hypothesis
   // if there is only background.

   Double_t result = 0;
   switch (sigma) {
   case -2:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
               result += fLRB[fISB[i]] / fNMC;
         return result;
      }
   case -1:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
               result += fLRB[fISB[i]] / fNMC;
         return result;
      }
   case 0:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
               result += fLRB[fISB[i]] / fNMC;
         return result;
      }
   case 1:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
               result += fLRB[fISB[i]] / fNMC;
         return result;
      }
   case 2:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
               result += fLRB[fISB[i]] / fNMC;
         return result;
      }
   default:
      return 0;
   }
}


//______________________________________________________________________________
Double_t TConfidenceLevel::GetExpectedCLb_sb(Int_t sigma) const
{
   // Get the expected Confidence Level for the background only
   // if there is signal and background.

   Double_t result = 0;
   switch (sigma) {
   case 2:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
               result += fLRS[fISS[i]] / fNMC;
         return result;
      }
   case 1:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
               result += fLRS[fISS[i]] / fNMC;
         return result;
      }
   case 0:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
               result += fLRS[fISS[i]] / fNMC;
         return result;
      }
   case -1:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
               result += fLRS[fISS[i]] / fNMC;
         return result;
      }
   case -2:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSS[fISS[i]] <= fTSS[fISS[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
               result += fLRS[fISS[i]] / fNMC;
         return result;
      }
   default:
      return 0;
   }
}


//______________________________________________________________________________
Double_t TConfidenceLevel::GetExpectedCLb_b(Int_t sigma) const
{
   // Get the expected Confidence Level for the background only
   // if there is only background.

   Double_t result = 0;
   switch (sigma) {
   case 2:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM2S)))]])
               result = (i + 1) / double (fNMC);
         return result;
      }
   case 1:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLM1S)))]])
               result = (i + 1) / double (fNMC);
         return result;
      }
   case 0:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLMED)))]])
               result = (i + 1) / double (fNMC);
         return result;
      }
   case -1:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP1S)))]])
               result = (i + 1) / double (fNMC);
         return result;
      }
   case -2:
      {
         for (Int_t i = 0; i < fNMC; i++)
            if (fTSB[fISB[i]] <= fTSB[fISB[TMath::Min((Int_t) fNMC,(Int_t) TMath::Max((Int_t) 1,(Int_t) (fNMC * fgMCLP2S)))]])
               result = (i + 1) / double (fNMC);
         return result;
      }
   }
   return result;
}


//______________________________________________________________________________
Double_t TConfidenceLevel::GetAverageCLsb() const
{
   // Get average CLsb.

   Double_t result = 0;
   Double_t psumsb = 0;
   for (Int_t i = 0; i < fNMC; i++) {
      psumsb += fLRB[fISB[i]] / fNMC;
      result += psumsb / fNMC;
   }
   return result;
}


//______________________________________________________________________________
Double_t TConfidenceLevel::GetAverageCLs() const
{
   // Get average CLs.

   Double_t result = 0;
   Double_t psumsb = 0;
   for (Int_t i = 0; i < fNMC; i++) {
      psumsb += fLRB[fISB[i]] / fNMC;
      result += ((psumsb / fNMC) / ((i + 1) / fNMC));
   }
   return result;
}


//______________________________________________________________________________
Double_t TConfidenceLevel::Get3sProbability() const
{
   // Get 3s probability.

   Double_t result = 0;
   Double_t psumbs = 0;
   for (Int_t i = 0; i < fNMC; i++) {
      psumbs += 1 / (Double_t) (fLRS[(fISS[(Int_t) (fNMC - i)])] * fNMC);
      if (psumbs <= fMCL3S)
         result = i / fNMC;
   }
   return result;
}


//______________________________________________________________________________
Double_t TConfidenceLevel::Get5sProbability() const
{
   // Get 5s probability.

   Double_t result = 0;
   Double_t psumbs = 0;
   for (Int_t i = 0; i < fNMC; i++) {
      psumbs += 1 / (Double_t) (fLRS[(fISS[(Int_t) (fNMC - i)])] * fNMC);
      if (psumbs <= fMCL5S)
         result = i / fNMC;
   }
   return result;
}

//______________________________________________________________________________
void  TConfidenceLevel::Draw(const Option_t*)
{
   // Display sort of a "canonical" -2lnQ plot.
   // This results in a plot with 2 elements:
   // - The histogram of -2lnQ for background hypothesis (full)
   // - The histogram of -2lnQ for signal and background hypothesis (dashed)
   // The 2 histograms are respectively named b_hist and sb_hist.

   TH1F h("TConfidenceLevel_Draw","",50,0,0);
   Int_t i;
   for (i=0; i<fNMC; i++) {
      h.Fill(-2*(fTSB[i]-fStot));
      h.Fill(-2*(fTSS[i]-fStot));
   }
   TH1F* b_hist  = new TH1F("b_hist", "-2lnQ",50,h.GetXaxis()->GetXmin(),h.GetXaxis()->GetXmax());
   TH1F* sb_hist = new TH1F("sb_hist","-2lnQ",50,h.GetXaxis()->GetXmin(),h.GetXaxis()->GetXmax());
   for (i=0; i<fNMC; i++) {
      b_hist->Fill(-2*(fTSB[i]-fStot));
      sb_hist->Fill(-2*(fTSS[i]-fStot));
   }
   b_hist->Draw();
   sb_hist->Draw("Same");
   sb_hist->SetLineStyle(3);
}


//______________________________________________________________________________
void  TConfidenceLevel::SetTSB(Double_t * in)
{
   // Set the TSB.
   fTSB = in; 
   TMath::Sort(fNNMC, fTSB, fISB, 0); 
}


//______________________________________________________________________________
void  TConfidenceLevel::SetTSS(Double_t * in)
{
   // Set the TSS.
   fTSS = in; 
   TMath::Sort(fNNMC, fTSS, fISS, 0); 
}
Revision:	1.1
Committed:	Thu Apr 1 15:14:10 2010 UTC (15 years, 1 month ago) by auterman
Content type:	text/plain
Branch point for:	CLs, MAIN
Log Message:	Initial revision
#	User	Rev	Content
1	auterman	1.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			}