// @(#)root/roostats:$Id: HybridResult.cxx 28963 2009-06-12 15:47:45Z wouter $

/*************************************************************************
 * Project: RooStats                                                     *
 * Package: RooFit/RooStats                                              *
 * Authors:                                                              *
 *   Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke       *
 *************************************************************************
 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

//____________________________________________________________________
/*
HybridResult class: this class is a fresh rewrite in RooStats of
	RooStatsCms/LimitResults developped by D. Piparo and G. Schott

The objects of this class store and access with lightweight methods the
information calculated by LimitResults through a Lent calculation using
MC toy experiments.
In some ways can be considered an extended and extensible implementation of the
TConfidenceLevel class (http://root.cern.ch/root/html/TConfidenceLevel.html).
*/

#include "RooDataHist.h"
#include "RooDataSet.h"
#include "RooGlobalFunc.h" // for RooFit::Extended()
#include "RooNLLVar.h"
#include "RooRealVar.h"
#include "RooAbsData.h"

#include "RooStats/HybridResult.h"
#include "RooStats/HybridPlot.h"


/// ClassImp for building the THtml documentation of the class 
ClassImp(RooStats::HybridResult)

using namespace RooStats;

///////////////////////////////////////////////////////////////////////////

HybridResult::HybridResult( const char *name, const char *title,
                            std::vector<double>& testStat_sb_vals,
                            std::vector<double>& testStat_b_vals) :
  //TNamed(name,title),
   HypoTestResult(name,title,0,0),
   fTestStat_data(-999.),
   fComputationsNulDoneFlag(false),
   fComputationsAltDoneFlag(false)
{
   // HybridResult constructor (with name, title and vectors of S+B and B values)

   int vector_size_sb = testStat_sb_vals.size();
   assert(vector_size_sb>0);

   int vector_size_b = testStat_b_vals.size();
   assert(vector_size_b>0);

   fTestStat_sb.reserve(vector_size_sb);
   for (int i=0;i<vector_size_sb;++i)
      fTestStat_sb.push_back(testStat_sb_vals[i]);

   fTestStat_b.reserve(vector_size_b);
   for (int i=0;i<vector_size_b;++i)
      fTestStat_b.push_back(testStat_b_vals[i]);
}

///////////////////////////////////////////////////////////////////////////

HybridResult::HybridResult( const char *name, const char *title) :
   //TNamed(name,title),
   HypoTestResult(name,title,0,0),
   fTestStat_data(-999.),
   fComputationsNulDoneFlag(false),
   fComputationsAltDoneFlag(false)
{
   // HybridResult constructor (with name and title)
}

///////////////////////////////////////////////////////////////////////////

HybridResult::HybridResult( ) :
   HypoTestResult("HybridResult_DefaultName","HybridResult",0,0),
//   TNamed("HybridResult_DefaultName","HybridResult"),
   fTestStat_data(-999.),
   fComputationsNulDoneFlag(false),
   fComputationsAltDoneFlag(false)
{
   // HybridResult default constructor
}

///////////////////////////////////////////////////////////////////////////

HybridResult::~HybridResult()
{
   // HybridResult destructor

   fTestStat_sb.clear();
   fTestStat_b.clear();
}

///////////////////////////////////////////////////////////////////////////

void HybridResult::SetDataTestStatistics(double testStat_data_val)
{
   // set the value of the test statistics on data

   fComputationsAltDoneFlag = false;
   fComputationsNulDoneFlag = false;
   fTestStat_data = testStat_data_val;
   return;
}

///////////////////////////////////////////////////////////////////////////

double HybridResult::NullPValue() const
{
   // return 1-CL_b : the B p-value

   if (fComputationsNulDoneFlag==false) {
      int nToys = fTestStat_b.size();
      if (nToys==0) {
         std::cout << "Error: no toy data present. Returning -1.\n";
         return -1;
      }

      double larger_than_measured=0;
      for (int iToy=0;iToy<nToys;++iToy)
         if ( fTestStat_b[iToy] > fTestStat_data ) ++larger_than_measured;

      if (larger_than_measured==0) std::cout << "Warning: CLb = 0 ... maybe more toys are needed!\n";

      fComputationsNulDoneFlag = true;
      fNullPValue = 1-larger_than_measured/nToys;
   }

   return fNullPValue;
}

///////////////////////////////////////////////////////////////////////////

double HybridResult::AlternatePValue() const
{
   // return CL_s+b : the S+B p-value

   if (fComputationsAltDoneFlag==false) {
      int nToys = fTestStat_b.size();
      if (nToys==0) {
         std::cout << "Error: no toy data present. Returning -1.\n";
         return -1;
      }

      double larger_than_measured=0;
      for (int iToy=0;iToy<nToys;++iToy)
         if ( fTestStat_sb[iToy] > fTestStat_data ) ++larger_than_measured;

      if (larger_than_measured==0) std::cout << "Warning: CLsb = 0 ... maybe more toys are needed!\n";

      fComputationsAltDoneFlag = true;
      fAlternatePValue = larger_than_measured/nToys;
   }

   return fAlternatePValue;
}

///////////////////////////////////////////////////////////////////////////

void HybridResult::Add(HybridResult* other)
{
   // add additional toy-MC experiments to the current results
   // use the data test statistics of the added object if none is already present (otherwise, ignore the new one)

   int other_size_sb = other->GetTestStat_sb().size();
   for (int i=0;i<other_size_sb;++i)
      fTestStat_sb.push_back(other->GetTestStat_sb()[i]);

   int other_size_b = other->GetTestStat_b().size();
   for (int i=0;i<other_size_b;++i)
      fTestStat_b.push_back(other->GetTestStat_b()[i]);

   // if no data is present use the other's HybridResult's data
   if (fTestStat_data==-999.)
      fTestStat_data = other->GetTestStat_data();

   fComputationsAltDoneFlag = false;
   fComputationsNulDoneFlag = false;

   return;
}

///////////////////////////////////////////////////////////////////////////

HybridPlot* HybridResult::GetPlot(const char* name,const char* title, int n_bins)
{
   // prepare a plot showing a result and return a pointer to a HybridPlot object
   // the needed arguments are: an object name, a title and the number of bins in the plot

   // default plot name
   TString plot_name;
   if ( TString(name)=="" ) {
      plot_name += GetName();
      plot_name += "_plot";
   } else plot_name = name;

   // default plot title
   TString plot_title;
   if ( TString(title)=="" ) {
      plot_title += GetTitle();
      plot_title += "_plot (";
      plot_title += fTestStat_b.size();
      plot_title += " toys)";
   } else plot_title = title;

   HybridPlot* plot = new HybridPlot( plot_name.Data(),
                                      plot_title.Data(),
                                      fTestStat_sb,
                                      fTestStat_b,
                                      fTestStat_data,
                                      n_bins,
                                      true );
   return plot;
}

///////////////////////////////////////////////////////////////////////////

void HybridResult::PrintMore(const char* /*options */)
{
   // Print out some information about the results

   std::cout << "\nResults " << GetName() << ":\n"
             << " - Number of S+B toys: " << fTestStat_b.size() << std::endl
             << " - Number of B toys: " << fTestStat_sb.size() << std::endl
             << " - test statistics evaluated on data: " << fTestStat_data << std::endl
             << " - CL_b " << CLb() << std::endl
             << " - CL_s+b " << CLsplusb() << std::endl
             << " - CL_s " << CLs() << std::endl;

   return;
}