// Author : Ajay Kumar
// Purdue University

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooChebychev.h"
#include "RooAddPdf.h"
#include "RooMCStudy.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TH2.h"
#include "RooFitResult.h"
#include "TStyle.h"
#include "TDirectory.h"

#include <RooRealVar.h>
#include <RooExponential.h>
#include <RooAbsPdf.h>
#include <RooAbsReal.h>
#include <RooExtendPdf.h>
#include <RooGenericPdf.h>
#include <RooGaussian.h>
#include <RooCBShape.h>
#include <RooVoigtian.h>
#include <RooFitResult.h>
#include <RooDataSet.h>
#include <RooArgSet.h>
#include <RooAddPdf.h>
#include <RooPlot.h>
#include <RooAbsArg.h>
#include <RooDstD0BG.h>
//#include <RooRelBreitWigner.h>
#include <RooRealIntegral.h>
#include <RooWorkspace.h>
#include <RooPrintable.h>
#include <RooFFTConvPdf.h>
#include <RooDataHist.h>

#include <RooRandom.h>
#include <RooPolynomial.h>
#include <RooConstVar.h>

#include <RooStats/HybridCalculatorOriginal.h>
#include <RooStats/HybridResult.h>
#include <RooStats/HybridPlot.h>




#include <TFile.h>
#include <TROOT.h>
#include <TCanvas.h>
#include <TLine.h>
#include <TPaveText.h>
#include <TText.h>
#include <TLegend.h>
#include <TH1.h>
#include <TMath.h>

#include <RooPowerLaw.h>
#include <RooProduct.h>
#include <RooPowerLaw.cc>
//#include <dyfit.C>

void ToyMyModel(int ntoys = 1000)
{
  using namespace RooFit;
  using namespace RooStats;
   // set RooFit random seed
   RooRandom::randomGenerator()->SetSeed(3007);
   // build the models for background and signal+background
   RooRealVar x("x","",-3,3);
   RooRealVar mass("mass","",200,2000);
   RooArgList observables(x); // variables to be generated
  
      // Create signal (=voigt)
   RooRealVar mean("mean","mean of voigt",1000) ;
   RooRealVar phi_width("phi_width","width of phi",27.0) ;
   //RooRealVar sigma_reso("sigma_reso","reso",2.0,0.1,10.) ;
  //Barrel Barrel reso used
   RooRealVar p0("p0","p0",0.01291);
   RooRealVar p1("p1","p1",1.835e-09);
   RooRealVar p2("p2","p2",-2.733e-09);
   RooPolyVar fm("fm","fm",mass,RooArgSet(p0,p1,p2)) ;


   //RooFormulaVar sigma_reso("sigma_reso","reso",p0 + p1*x + p2 * x * x, RooArgList(p0,p1,p2,x));

   //RooVoigtian sig_pdf("sig_pdf" , "", x, mean, phi_width, sigma_reso, kFALSE)   ;
   RooVoigtian sig_pdf("sig_pdf" , "", x, mean, phi_width, fm, kFALSE) ;
   RooRealVar sig_yield("sig_yield","",20,0,300);

  //bkg pdf

   RooRealVar N("N","N", 600,500,8.96732e+07);
   RooRealVar a1("a1","a1",0);
   RooRealVar a("a","a",-8.59843e-03);
   RooRealVar b("b","b",2.86870e-06);
   RooRealVar k("k","k",-2.05193e+00);

   RooPolyVar fp("fp","fp",mass,RooArgSet(a1,a,b)) ;
   RooExponential expo("expo", "", x, fp);
   RooPowerLaw powerfac("powerfac", "", mass, k);

  //RooAddPdf bkg_pdf("bkg_pdf","",RooArgList(expo,powerfac),1/N,1-1/N) ;

 // RooAddPdf bkg_pdf("bkg_pdf","",RooArgList(expo,powerfac),N);
    RooProdPdf bkg_pdf("bkg_pdf","",RooArgSet(expo,powerfac)) ;


   //RooProdPdf bkg_pdf("" ,   ,RooArgSet(px,py)) ;  

   //RooProdPdf bkg_pdfl("bkg_pdf","",RooArgSet(expo, powerfac));

   //RooFormulaVar fitval2("fival2",pow_func, TMath::Power(x,k),RooArgList(x,k));
   //RooFormulaVar fitval1("fitval1",func, N * exp(x*(-a)  +  x * x *b), RooArgList(N,x,a,b));
  //RooGenericPdf bkg_pdf("bkg_pdf","", N * exp(  +  mass * mass *b) ,RooArgList(N,mass,a,b);

  // RooPolynomial bkg_pdf("bkg_pdf","", x, RooConst(0));
 // RooPolynomial bkg_pdf("bkg_pdf","", x, fp);
  RooRealVar bkg_yield("bkg_yield","",500,0,2000);

  RooExtendPdf bkg_ext_pdf("bkg_ext_pdf","",bkg_pdf,bkg_yield);

  //bkg_yield.setConstant(kTRUE);
  sig_yield.setConstant(kTRUE);

  // total sig+bkg (extended PDF)
  RooAddPdf tot_pdf("tot_pdf","",RooArgList(sig_pdf,bkg_pdf),RooArgList(sig_yield,bkg_yield));

  /// build the prior PDF on the parameters to be integrated
  //  // gaussian contraint on the background yield ( N_B = 40 +/- 10  ie. 25% )

  RooGaussian bkg_yield_prior("bkg_yield_prior","",bkg_yield,RooConst(bkg_yield.getVal()),RooConst(10.));
  RooArgSet nuisance_parameters(bkg_yield); // variables to be integrated
  // generate a data sample
  RooDataSet* data = tot_pdf.generate(observables,RooFit::Extended());
  // run HybridCalculator on those inputs
  // use interface from HypoTest calculator by default
  HybridCalculatorOriginal myHybridCalc(*data, tot_pdf , bkg_ext_pdf ,
                                   &nuisance_parameters, &bkg_yield_prior);
  // here I use the default test statistics: 2*lnQ (optional)
  myHybridCalc.SetTestStatistic(1);


  myHybridCalc.SetNumberOfToys(ntoys);
  myHybridCalc.UseNuisance(true);
  // for speed up generation (do binned data)
  myHybridCalc.SetGenerateBinned(false);
  // calculate by running ntoys for the S+B and B hypothesis and retrieve the result
  HybridResult* myHybridResult = myHybridCalc.GetHypoTest();
  if (! myHybridResult) {
     std::cerr << "\nError returned from Hypothesis test" << std::endl;
     return;
  }

  /// nice plot of the results  
  HybridPlot* myHybridPlot = myHybridResult->GetPlot("myHybridPlot","Plot of results with HybridCalculatorOriginal",100);
  myHybridPlot->Draw(); 
  /// recover and display the results
  double clsb_data = myHybridResult->CLsplusb();
  double clb_data = myHybridResult->CLb();
  double cls_data = myHybridResult->CLs();
  double data_significance = myHybridResult->Significance();
  double min2lnQ_data = myHybridResult->GetTestStat_data();



  /// compute the mean expected significance from toys
  //

  double mean_sb_toys_test_stat = myHybridPlot->GetSBmean();
  myHybridResult->SetDataTestStatistics(mean_sb_toys_test_stat);
  double toys_significance = myHybridResult->Significance();

  std::cout << "Completed HybridCalculatorOriginal example:\n";
  std::cout << " - -2lnQ = " << min2lnQ_data << endl;
  std::cout << " - CL_sb = " << clsb_data << std::endl;
  std::cout << " - CL_b  = " << clb_data << std::endl;
  std::cout << " - CL_s  = " << cls_data << std::endl;
  std::cout << " - significance of data  = " << data_significance << std::endl;
  std::cout << " - mean significance of toys  = " << toys_significance << std::endl;

}