// @(#)root/hist:$Name:  $:$Id: TLimit.cxx,v 1.5 2003/03/21 14:53:49 brun Exp $
// Author: Christophe.Delaere@cern.ch   21/08/2002

///////////////////////////////////////////////////////////////////////////
//
// TLimit
//
// 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 "TLimit.h"
#include "TArrayD.h"
#include "TOrdCollection.h"
#include "TConfidenceLevel.h"
#include "TLimitDataSource.h"
#include "TRandom3.h"
#include "TH1.h"
#include "TObjArray.h"
#include "TMath.h"
#include "TIterator.h"
#include "TObjString.h"
#include "TClassTable.h"
#include "Riostream.h"

ClassImp(TLimit)

TArrayD *TLimit::fgTable = new TArrayD(0);
TOrdCollection *TLimit::fgSystNames = new TOrdCollection();

 TConfidenceLevel *TLimit::ComputeLimit(TLimitDataSource * data,
                                       Int_t nmc, TRandom * generator,
                                       Double_t(*statistic) (Double_t,
                                                             Double_t,
                                                             Double_t))
{
   // class TLimit
   // ------------
   //
   // Algorithm to compute 95% C.L. limits using the Likelihood ratio
   // semi-bayesian method.
   // It takes signal, background and data histograms wrapped in a
   // TLimitDataSource as input and runs a set of Monte Carlo experiments in
   // order to compute the limits. If needed, inputs are fluctuated according
   // to systematics. The output is a TConfidenceLevel.
   //
   // class TLimitDataSource
   // ----------------------
   //
   // Takes the signal, background and data histograms as well as different
   // systematics sources to form the TLimit input.
   //
   //  class TConfidenceLevel
   //  ----------------------
   //
   // Final result of the TLimit algorithm. It is created just after the
   // time-consuming part and can be stored in a TFile for further processing.
   // It contains light methods to return CLs, CLb and other interesting
   // quantities.
   //
   // The actual algorithm...
   // From an input (TLimitDataSource) it produces an output TConfidenceLevel.
   // For this, nmc Monte Carlo experiments are performed.
   // As usual, the larger this number, the longer the compute time,
   // but the better the result.
   //
   /*
   
   
Supposing that there is a plotfile.root file containing 3 histograms
           (signal, background and data), you can imagine doing things like:
   

   

    TFile* infile=new TFile("plotfile.root","READ");
    infile->cd();
    TH1D* sh=(TH1D*)infile->Get("signal");
    TH1D* bh=(TH1D*)infile->Get("background");
    TH1D* dh=(TH1D*)infile->Get("data");
    TLimitDataSource* mydatasource = new TLimitDataSource(sh,bh,dh);
    TConfidenceLevel *myconfidence = TLimit::ComputeLimit(mydatasource,50000);
    cout << "  CLs    : " << myconfidence->CLs()  << endl;
    cout << "  CLsb   : " << myconfidence->CLsb() << endl;
    cout << "  CLb    : " << myconfidence->CLb()  << endl;
    cout << "< CLs >  : " << myconfidence->GetExpectedCLs_b()  << endl;
    cout << "< CLsb > : " << myconfidence->GetExpectedCLsb_b() << endl;
    cout << "< CLb >  : " << myconfidence->GetExpectedCLb_b()  << endl;
    delete myconfidence;
    delete mydatasource;
    infile->Close();