#include "TStopwatch.h"
#include "TCanvas.h"
#include "TRandom3.h"
#include "TROOT.h"
#include "RooPlot.h"
#include "RooAbsPdf.h"
#include "RooWorkspace.h"
#include "RooDataSet.h"
#include "RooGlobalFunc.h"
#include "RooFitResult.h"
#include "RooRandom.h"
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/LikelihoodInterval.h"
#include "RooStats/LikelihoodIntervalPlot.h"
#include "RooStats/BayesianCalculator.h"
#include "RooStats/MCMCCalculator.h"
#include "RooStats/MCMCInterval.h"
#include "RooStats/MCMCIntervalPlot.h"
#include "RooStats/ProposalHelper.h"
#include "RooStats/SimpleInterval.h"
#include "RooStats/FeldmanCousins.h"
#include "RooStats/PointSetInterval.h"
#include "RooStats/HypoTestResult.h"
#include "RooStats/HypoTestResult.h"
#include "RooStats/SamplingDistribution.h"
#include "TF1.h"

using namespace RooFit;
using namespace RooStats;


/*
 * Problem: Have a known background of b=noffObs events and observe nonObs=s+b events. 
 * Test with the profilelikelihoodratio for the null hypothesis s=0 and for the alternative hypothesis s>0 
 * -> One Parameter (s)
 * -> Likelihood function for Null hyptohesis L0=TMath::Poisson(nonObs,b)
 * -> Likelihood function for alternative hypothesis L1=TMath::Poisson(nonObs,b+s) which is maximal at s=noffObs-nonObs.
 *
 * For the concrete example of nonObs=150. and noffObs=100. this means 
 * -> Null hypothesis -log(L0)=-log(TMath::Poisson(150.,100.)) ~ 14.24
 * -> Alternative hypothesis -log(L1)=-log(TMath::Poisson(150.,150)) ~ 3.22
 *
 * The problem is that in the GetHypoTest() function below (implemented in ProfileLikelihoodCalculator.cxx) the following values are calculated:
 * 
 * Variables in ProfileLikelihoodCalculator::GetHypoTest():
 * Double_t NLLatMLE = fFitResult->minNll(); (this variable should hold -log(L1) as it is the maximum likelihood for the best fit of the likelihood function to the data leaving all variables floating)
 * What I get is however ~18.56 instead of the expected 3.22
 *
 * For the variable 
 * Double_t NLLatCondMLE = fit2->minNll(); 
 * holding the -log(L0) for the null hypothesis 
 * in GetHypoTest() I get the value 40.4 instead of the expected value 14.24 for the null hypthesis likelihood
 *
 * What is going on here?
*/
void test(double nonObs=150., double noffObs=100.){
 
  RooWorkspace* wspace = new RooWorkspace("wspace");

  RooRealVar* non = (RooRealVar*) wspace->factory("non[0., 400.]");

  RooRealVar* b=(RooRealVar*) wspace->factory("b[0., 400.]");
  RooRealVar* s=(RooRealVar*) wspace->factory("s[0., 100.]"); 

  // make model
  wspace->factory("Poisson::on(non, sum::splusb(s,b))");

  wspace->factory("PROD::model(on)");

  wspace->defineSet("obs","non");
  wspace->defineSet("poi","s");
  wspace->defineSet("nuis","b");

  non->setVal(nonObs);
  b->setVal(noffObs);

  RooArgSet* ArgSet = new RooArgSet("args");   
  ArgSet->add(*non);
  ArgSet->add(*b);

  RooDataSet *data = new RooDataSet("modelData", "modelData", *ArgSet); 
  data->add(*non);
  data->add(*b);

  wspace->import(*data);

  wspace->var("non")->setVal(nonObs);
  wspace->var("non")->setConstant(true);

  wspace->var("b")->setVal(noffObs);
  wspace->var("b")->setConstant(true); 

  ModelConfig* modelConfig = new ModelConfig("FourBins");
  modelConfig->SetWorkspace(*wspace);
  modelConfig->SetPdf(*wspace->pdf("model"));
  modelConfig->SetParametersOfInterest(*wspace->set("poi"));
  modelConfig->SetNuisanceParameters(*wspace->set("nuis"));
  wspace->import(*modelConfig);

  // use ProfileLikelihood
  ProfileLikelihoodCalculator plc(*data, *modelConfig); 

  RooRealVar *xs=wspace->var("s");

  RooArgSet *nullpar = new RooArgSet("nullParams");
  nullpar->addClone(*xs);
  
  nullpar->setRealValue("s",0.);
  plc.SetNullParameters(*nullpar);

  HypoTestResult *htr = plc.GetHypoTest();
  
  double pval0=htr->NullPValue(); 

  std::cout<<"nullpval="<<pval0<<" sig="<<PValueToSignificance(pval0)<<std::endl;



}