// @(#)root/roostats:$Id$

/*************************************************************************
 * 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
New contributions to this class have been written by Matthias Wolf (error estimation)

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 
using namespace std;

ClassImp(RooStats::HybridResult)

using namespace RooStats;

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

HybridResult::HybridResult( const char *name) :
   HypoTestResult(name),
   fTestStat_data(-999.),
   fComputationsNulDoneFlag(false),
   fComputationsAltDoneFlag(false),
   fSumLargerValues(false)
{
   // HybridResult default constructor (with name )
}

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

HybridResult::HybridResult( const char *name,
                            const std::vector<double>& testStat_sb_vals,
                            const std::vector<double>& testStat_b_vals,
			    bool sumLargerValues ) :
   HypoTestResult(name,0,0),
   fTestStat_data(-999.),
   fComputationsNulDoneFlag(false),
   fComputationsAltDoneFlag(false),
   fSumLargerValues(sumLargerValues)
{
   // 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()
{
   // 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;
      if (fSumLargerValues) {
	for (int iToy=0;iToy<nToys;++iToy)
	  if ( fTestStat_b[iToy] >= fTestStat_data ) ++larger_than_measured;
      } else {
	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;
      if (fSumLargerValues) {
	for (int iToy=0;iToy<nToys;++iToy)
	  if ( fTestStat_sb[iToy] >= fTestStat_data ) ++larger_than_measured;
      } else {
	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;
}

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

Double_t HybridResult::CLbError() const
{
  // Returns an estimate of the error on CLb assuming a binomial error on
  // CLb:
  // BEGIN_LATEX
  // #sigma_{CL_{b}} &=& #sqrt{CL_{b} #left( 1 - CL_{b} #right) / n_{toys}}
  // END_LATEX
  unsigned const int n = fTestStat_b.size();
  return TMath::Sqrt(CLb() * (1. - CLb()) / n);
}

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

Double_t HybridResult::CLsplusbError() const
{
  // Returns an estimate of the error on CLsplusb assuming a binomial
  // error on CLsplusb:
  // BEGIN_LATEX
  // #sigma_{CL_{s+b}} &=& #sqrt{CL_{s+b} #left( 1 - CL_{s+b} #right) / n_{toys}}
  // END_LATEX
  unsigned const int n = fTestStat_sb.size();
  return TMath::Sqrt(CLsplusb() * (1. - CLsplusb()) / n);
}

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

Double_t HybridResult::CLsError() const
{
  // Returns an estimate of the error on CLs through combination of the
  // errors on CLb and CLsplusb:
  // BEGIN_LATEX
  // #sigma_{CL_s} &=& CL_s #sqrt{#left( #frac{#sigma_{CL_{s+b}}}{CL_{s+b}} #right)^2 + #left( #frac{#sigma_{CL_{b}}}{CL_{b}} #right)^2}
  // END_LATEX
  unsigned const int n_b = fTestStat_b.size();
  unsigned const int n_sb = fTestStat_sb.size();
  
  if (CLb() == 0 || CLsplusb() == 0)
    return 0;
  
  double cl_b_err = (1. - CLb()) / (n_b * CLb());
  double cl_sb_err = (1. - CLsplusb()) / (n_sb * CLsplusb());
  
  return CLs() * TMath::Sqrt(cl_b_err + cl_sb_err);
}

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

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;
}










 HybridResult.cxx:1
 HybridResult.cxx:2
 HybridResult.cxx:3
 HybridResult.cxx:4
 HybridResult.cxx:5
 HybridResult.cxx:6
 HybridResult.cxx:7
 HybridResult.cxx:8
 HybridResult.cxx:9
 HybridResult.cxx:10
 HybridResult.cxx:11
 HybridResult.cxx:12
 HybridResult.cxx:13
 HybridResult.cxx:14
 HybridResult.cxx:15
 HybridResult.cxx:16
 HybridResult.cxx:17
 HybridResult.cxx:18
 HybridResult.cxx:19
 HybridResult.cxx:20
 HybridResult.cxx:21
 HybridResult.cxx:22
 HybridResult.cxx:23
 HybridResult.cxx:24
 HybridResult.cxx:25
 HybridResult.cxx:26
 HybridResult.cxx:27
 HybridResult.cxx:28
 HybridResult.cxx:29
 HybridResult.cxx:30
 HybridResult.cxx:31
 HybridResult.cxx:32
 HybridResult.cxx:33
 HybridResult.cxx:34
 HybridResult.cxx:35
 HybridResult.cxx:36
 HybridResult.cxx:37
 HybridResult.cxx:38
 HybridResult.cxx:39
 HybridResult.cxx:40
 HybridResult.cxx:41
 HybridResult.cxx:42
 HybridResult.cxx:43
 HybridResult.cxx:44
 HybridResult.cxx:45
 HybridResult.cxx:46
 HybridResult.cxx:47
 HybridResult.cxx:48
 HybridResult.cxx:49
 HybridResult.cxx:50
 HybridResult.cxx:51
 HybridResult.cxx:52
 HybridResult.cxx:53
 HybridResult.cxx:54
 HybridResult.cxx:55
 HybridResult.cxx:56
 HybridResult.cxx:57
 HybridResult.cxx:58
 HybridResult.cxx:59
 HybridResult.cxx:60
 HybridResult.cxx:61
 HybridResult.cxx:62
 HybridResult.cxx:63
 HybridResult.cxx:64
 HybridResult.cxx:65
 HybridResult.cxx:66
 HybridResult.cxx:67
 HybridResult.cxx:68
 HybridResult.cxx:69
 HybridResult.cxx:70
 HybridResult.cxx:71
 HybridResult.cxx:72
 HybridResult.cxx:73
 HybridResult.cxx:74
 HybridResult.cxx:75
 HybridResult.cxx:76
 HybridResult.cxx:77
 HybridResult.cxx:78
 HybridResult.cxx:79
 HybridResult.cxx:80
 HybridResult.cxx:81
 HybridResult.cxx:82
 HybridResult.cxx:83
 HybridResult.cxx:84
 HybridResult.cxx:85
 HybridResult.cxx:86
 HybridResult.cxx:87
 HybridResult.cxx:88
 HybridResult.cxx:89
 HybridResult.cxx:90
 HybridResult.cxx:91
 HybridResult.cxx:92
 HybridResult.cxx:93
 HybridResult.cxx:94
 HybridResult.cxx:95
 HybridResult.cxx:96
 HybridResult.cxx:97
 HybridResult.cxx:98
 HybridResult.cxx:99
 HybridResult.cxx:100
 HybridResult.cxx:101
 HybridResult.cxx:102
 HybridResult.cxx:103
 HybridResult.cxx:104
 HybridResult.cxx:105
 HybridResult.cxx:106
 HybridResult.cxx:107
 HybridResult.cxx:108
 HybridResult.cxx:109
 HybridResult.cxx:110
 HybridResult.cxx:111
 HybridResult.cxx:112
 HybridResult.cxx:113
 HybridResult.cxx:114
 HybridResult.cxx:115
 HybridResult.cxx:116
 HybridResult.cxx:117
 HybridResult.cxx:118
 HybridResult.cxx:119
 HybridResult.cxx:120
 HybridResult.cxx:121
 HybridResult.cxx:122
 HybridResult.cxx:123
 HybridResult.cxx:124
 HybridResult.cxx:125
 HybridResult.cxx:126
 HybridResult.cxx:127
 HybridResult.cxx:128
 HybridResult.cxx:129
 HybridResult.cxx:130
 HybridResult.cxx:131
 HybridResult.cxx:132
 HybridResult.cxx:133
 HybridResult.cxx:134
 HybridResult.cxx:135
 HybridResult.cxx:136
 HybridResult.cxx:137
 HybridResult.cxx:138
 HybridResult.cxx:139
 HybridResult.cxx:140
 HybridResult.cxx:141
 HybridResult.cxx:142
 HybridResult.cxx:143
 HybridResult.cxx:144
 HybridResult.cxx:145
 HybridResult.cxx:146
 HybridResult.cxx:147
 HybridResult.cxx:148
 HybridResult.cxx:149
 HybridResult.cxx:150
 HybridResult.cxx:151
 HybridResult.cxx:152
 HybridResult.cxx:153
 HybridResult.cxx:154
 HybridResult.cxx:155
 HybridResult.cxx:156
 HybridResult.cxx:157
 HybridResult.cxx:158
 HybridResult.cxx:159
 HybridResult.cxx:160
 HybridResult.cxx:161
 HybridResult.cxx:162
 HybridResult.cxx:163
 HybridResult.cxx:164
 HybridResult.cxx:165
 HybridResult.cxx:166
 HybridResult.cxx:167
 HybridResult.cxx:168
 HybridResult.cxx:169
 HybridResult.cxx:170
 HybridResult.cxx:171
 HybridResult.cxx:172
 HybridResult.cxx:173
 HybridResult.cxx:174
 HybridResult.cxx:175
 HybridResult.cxx:176
 HybridResult.cxx:177
 HybridResult.cxx:178
 HybridResult.cxx:179
 HybridResult.cxx:180
 HybridResult.cxx:181
 HybridResult.cxx:182
 HybridResult.cxx:183
 HybridResult.cxx:184
 HybridResult.cxx:185
 HybridResult.cxx:186
 HybridResult.cxx:187
 HybridResult.cxx:188
 HybridResult.cxx:189
 HybridResult.cxx:190
 HybridResult.cxx:191
 HybridResult.cxx:192
 HybridResult.cxx:193
 HybridResult.cxx:194
 HybridResult.cxx:195
 HybridResult.cxx:196
 HybridResult.cxx:197
 HybridResult.cxx:198
 HybridResult.cxx:199
 HybridResult.cxx:200
 HybridResult.cxx:201
 HybridResult.cxx:202
 HybridResult.cxx:203
 HybridResult.cxx:204
 HybridResult.cxx:205
 HybridResult.cxx:206
 HybridResult.cxx:207
 HybridResult.cxx:208
 HybridResult.cxx:209
 HybridResult.cxx:210
 HybridResult.cxx:211
 HybridResult.cxx:212
 HybridResult.cxx:213
 HybridResult.cxx:214
 HybridResult.cxx:215
 HybridResult.cxx:216
 HybridResult.cxx:217
 HybridResult.cxx:218
 HybridResult.cxx:219
 HybridResult.cxx:220
 HybridResult.cxx:221
 HybridResult.cxx:222
 HybridResult.cxx:223
 HybridResult.cxx:224
 HybridResult.cxx:225
 HybridResult.cxx:226
 HybridResult.cxx:227
 HybridResult.cxx:228
 HybridResult.cxx:229
 HybridResult.cxx:230
 HybridResult.cxx:231
 HybridResult.cxx:232
 HybridResult.cxx:233
 HybridResult.cxx:234
 HybridResult.cxx:235
 HybridResult.cxx:236
 HybridResult.cxx:237
 HybridResult.cxx:238
 HybridResult.cxx:239
 HybridResult.cxx:240
 HybridResult.cxx:241
 HybridResult.cxx:242
 HybridResult.cxx:243
 HybridResult.cxx:244
 HybridResult.cxx:245
 HybridResult.cxx:246
 HybridResult.cxx:247
 HybridResult.cxx:248
 HybridResult.cxx:249
 HybridResult.cxx:250
 HybridResult.cxx:251
 HybridResult.cxx:252
 HybridResult.cxx:253
 HybridResult.cxx:254
 HybridResult.cxx:255
 HybridResult.cxx:256
 HybridResult.cxx:257
 HybridResult.cxx:258
 HybridResult.cxx:259
 HybridResult.cxx:260
 HybridResult.cxx:261
 HybridResult.cxx:262
 HybridResult.cxx:263
 HybridResult.cxx:264
 HybridResult.cxx:265
 HybridResult.cxx:266
 HybridResult.cxx:267
 HybridResult.cxx:268
 HybridResult.cxx:269
 HybridResult.cxx:270
 HybridResult.cxx:271
 HybridResult.cxx:272
 HybridResult.cxx:273
 HybridResult.cxx:274
 HybridResult.cxx:275
 HybridResult.cxx:276
 HybridResult.cxx:277
 HybridResult.cxx:278
 HybridResult.cxx:279
 HybridResult.cxx:280
 HybridResult.cxx:281
 HybridResult.cxx:282
 HybridResult.cxx:283
 HybridResult.cxx:284
 HybridResult.cxx:285
 HybridResult.cxx:286
 HybridResult.cxx:287
 HybridResult.cxx:288
 HybridResult.cxx:289
 HybridResult.cxx:290
 HybridResult.cxx:291
 HybridResult.cxx:292
 HybridResult.cxx:293
 HybridResult.cxx:294
 HybridResult.cxx:295
 HybridResult.cxx:296
 HybridResult.cxx:297
 HybridResult.cxx:298
 HybridResult.cxx:299
 HybridResult.cxx:300
 HybridResult.cxx:301
 HybridResult.cxx:302
 HybridResult.cxx:303
 HybridResult.cxx:304