/////////////////////////////////////////////////////////////////////////
//
// 'Hypothesis Test Inversion' RooStats tutorial macro #801
// author: Gregory Schott
// date Sep 2009
//
// This tutorial shows an example of using the HypoTestInverterOriginal class 
//
/////////////////////////////////////////////////////////////////////////

#include "RooRealVar.h"
#include "RooConstVar.h"
#include "RooProdPdf.h"
#include "RooWorkspace.h"
#include "RooDataSet.h"
#include "RooPolynomial.h"
#include "RooAddPdf.h"
#include "RooExtendPdf.h"

#include "RooStats/HypoTestInverterOriginal.h"
#include "RooStats/HypoTestInverterResult.h"
#include "RooStats/HypoTestInverterPlot.h"
#include "RooStats/HybridCalculatorOriginal.h"
#include <iostream>
#include "TGraphErrors.h"
#include "TApplication.h"
#include <sstream>
#include "TFile.h"
#include <memory>
using namespace RooFit;
using namespace RooStats;
using namespace std;

int main(int argc, char** argv)
{

  vector<double> vpoints;
  
  for (int k = 0; k<5; k++){
    //    vpoints.push_back(3 + 3.*k/10.);
    vpoints.push_back(3 + 3.*k/5.);
  }
  
  TApplication app("app", 0, 0);
  
  // prepare the model
  RooRealVar lumi("lumi","luminosity",1);
  RooRealVar r("r","cross-section ratio",3.74,0,50);
  RooFormulaVar ns("ns","1*r*lumi",RooArgList(lumi,r));
  RooRealVar nb("nb","background yield",1);
  RooRealVar x("x","dummy observable",0,1);
  RooConstVar p0(RooFit::RooConst(0));
  RooPolynomial flatPdf("flatPdf","flat PDF",x,p0);
  RooAddPdf totPdf("totPdf","S+B model",RooArgList(flatPdf,flatPdf),RooArgList(ns,nb));
  RooExtendPdf bkgPdf("bkgPdf","B-only model",flatPdf,nb);
  RooDataSet* data = totPdf.generate(x,1);

  // prepare the calculator
  HybridCalculatorOriginal myhc(*data, totPdf, bkgPdf,0,0);
  myhc.SetTestStatistic(2);
  myhc.SetNumberOfToys(200);
  myhc.UseNuisance(false);                            
  
  // run the hypothesis-test invertion
  HypoTestInverterResult* globalResult;

    //First inverter
    HypoTestInverterOriginal myInverter(myhc, r);
        myInverter.SetTestSize(0.10);
    myInverter.UseCLs(true); 
    myInverter.RunOnePoint(3);
    HypoTestInverterResult* results = static_cast<HypoTestInverterResult*>(myInverter.GetInterval());
    cout << "Results size: " << results->ArraySize() << endl;

    //Second inverter
    HypoTestInverterOriginal myInverter2(myhc, r);
    myInverter2.SetTestSize(0.10);
    myInverter2.UseCLs(true);
    myInverter2.RunOnePoint(3);
    HypoTestInverterResult* results2 = static_cast<HypoTestInverterResult*>(myInverter2.GetInterval());
    cout << "Array size II " << results2->ArraySize() << endl;

    //Third inverter
    HypoTestInverterOriginal myInverter3(myhc, r);
    myInverter3.SetTestSize(0.10);
    myInverter3.UseCLs(true);
    myInverter3.RunOnePoint(4);
    HypoTestInverterResult* results3 = static_cast<HypoTestInverterResult*>(myInverter3.GetInterval());
    cout << "Array size II " << results3->ArraySize() << endl;

    globalResult = static_cast<HypoTestInverterResult*> (results->Clone()); 
    results->Add( *(results2) );
    results->Add( *(results3) );


    
    cout << "*******************" << endl;
    cout << "Array Size: " << results->ArraySize() << endl;
    cout << "Single upper limit: " << results->UpperLimit() << endl;  
    cout << "CLs: " << results->GetLastYValue() << endl;
    
  
  return 0;
  
}