// @(#)root/hist:$Name:  $:$Id: TConfidenceLevel.cxx,v 1.8 2007/02/02 19:08:39 brun Exp $
// 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); 
}