#include "TH1F.h"
#include "TFile.h"

using namespace RooFit;
using namespace RooStats;
using namespace std;

void ComputeCLS(){

  
  TString filename;
  filename="file.root"; //with add plots
  
  float yieldSig_input = 0.;
  float yieldBkg_input = 0.;

  TFile * filechannel = new TFile(filename);

 ...............
      
// Here I build an arbitrary distorded histo corresponding to a +/- one sigma systematic effect affecting the shape & normaliastion
// test for a tiny scaling effect only  

      for (unsigned int k = 1; k < histo_Sig->GetNbinsX()+1; k++) {
        double shiftedbin_minus = 0;
	shiftedbin_minus = 0.9999*histo_Sig->GetBinContent(k);
        double shiftedbin_plus = 0;
	shiftedbin_plus = 1.0001*histo_Sig->GetBinContent(k);
        histo_Sig_MinusOneSigma->SetBinContent(k,shiftedbin_minus );
        histo_Sig_PlusOneSigma->SetBinContent(k,shiftedbin_plus );
      }


  /// set RooFit random seed
  RooRandom::randomGenerator()->SetSeed(1234); 

  RooRealVar the_Obs("the_Observable", "Observable", 0., 200.);  // 

  //create RooDataHisto for the data 
  RooDataHist*  rooData_data        = new RooDataHist("Obs_data",       "Obs_data",        the_Obs,  histo_Data      );
  RooHistPdf    rooPDF_data    ("rooPDF_data", "rooPDF_data",      the_Obs,  *rooData_data      , 0);
  
  //create RooDataHisto for signal
  RooDataHist*  rooData_signal      = new RooDataHist("Obs_signal",     "Obs_signal",      the_Obs,  histo_Sig    );
  RooHistPdf    rooPDF_signal    ("rooPDF_signal", "rooPDF_signal",    the_Obs,  *rooData_signal    , 0);
  
  //create RooDataHisto for background
  RooDataHist*   rooData_background = new RooDataHist("Obs_background", "Obs_background",  the_Obs, histo_Bkg );
  RooHistPdf  rooPDF_background ("rooPDF_background", "rooPDF_background",the_Obs,  *rooData_background, 0);
  RooAbsPdf*  rooPDF_bgd = dynamic_cast<RooAbsPdf*>( rooPDF_background );
  
  //create RooDataHisto for signal with systematic inducing distorsion
  RooDataHist*  rooData_signal_minus      = new RooDataHist("Obs_signal_minus",     "Obs_signal_minus",      the_Obs,  histo_Sig_MinusOneSigma    );
  RooHistPdf    rooPDF_signal_minus    ("rooPDF_signal_minus", "rooPDF_signal_minus",    the_Obs,  *rooData_signal_minus   );
  RooDataHist*  rooData_signal_plus      = new RooDataHist("Obs_signal_plus",     "Obs_signal_plus",      the_Obs,  histo_Sig_PlusOneSigma    );
  RooHistPdf    rooPDF_signal_plus    ("rooPDF_signal_plus", "rooPDF_signal_plus",    the_Obs,  *rooData_signal_plus   );
  
// Construct interpolating pdf in (the_Obs,alpha), represent histo_Sig_MinusOneSigma at min and histo_Sig_PlusOneSigma at max
  RooArgList pdfList;
  pdfList.add(rooPDF_signal);
  pdfList.add(rooPDF_signal_plus);
  pdfList.add(rooPDF_signal_minus);

  TVectorD vec(3);
  vec(0) = 0.0; vec(1) = 1.0; vec(2) = -1.0;
  RooRealVar alpha("alpha","alpha",0.,-5,5);
  RooMomentMorph rooPDF_signal_interpol("rooPDF_signal_interpol","rooPDF_signal_interpol",alpha,RooArgList(the_Obs), pdfList,vec);
  rooPDF_signal_interpol.Print();

// // Test to see if interpol is ok
//   rooPDF_signal.Print();
//   RooPlot* frame = the_Obs.frame();
//   rooPDF_signal.plotOn(frame,LineColor(2));
// 
//   frame ->Draw();
//   
//   alpha.setVal(1.0);
// //   rooPDF_signal_interpol.plotOn(frame,LineColor(8)); 
// //   rooPDF_signal_plus.plotOn(frame,LineStyle(2),LineColor(1)); 
//   alpha.setVal(-1.0);
//   rooPDF_signal_interpol.plotOn(frame,LineColor(46)); 
//   rooPDF_signal_minus.plotOn(frame,LineStyle(2),LineColor(1)); 
//   
//   frame ->Draw();
//   

  // Set yields

  yieldSig_input = histo_Sig->Integral();
  yieldBkg_input = histo_Bkg->Integral();
  RooRealVar yieldSig("yieldSig", "", yieldSig_input, 0,10000000000);
  RooRealVar yieldBkg("yieldBkg", "", yieldBkg_input, 0,10000000000);
  
  // Define the model (ie the set of PDF) 
  RooAddPdf tot_pdf("tot_pdf","",RooArgList(rooPDF_signal, rooPDF_background), RooArgList(yieldSig,yieldBkg));

////***********************************************************************//

////   HybridCalculator w/wo systematics

////***********************************************************************//

      RooWorkspace* ws = new RooWorkspace("ws","ws");
      
//// Importation
      ws->import(rooPDF_data);
      ws->import(rooPDF_background);
      ws->import(rooPDF_signal_interpol);
////  

/// Parameters definition
      RooRealVar yieldBkg1("yieldBkg1", "yieldBkg1",yieldBkg_input , 0.,1000000000, "events");
      ws->import(yieldBkg1);
      RooRealVar yieldSig1("yieldSig1", "yieldSig1",yieldSig_input , 0.,1000000000, "events");
      ws->import(yieldSig1);
      
//// Pdf definitions 
//// s+b model   
// Here I use alternatively pdfs with morphing (and a nuisance parameter) or a reference pdf without any nuisance                       
//       RooAbsPdf* sbPDF = (RooAbsPdf*) ws->factory("SUM::sbPDF( yieldSig1*rooPDF_signal ,yieldBkg1*rooPDF_background )");  //  w/o morphing
       RooAbsPdf* sbPDF = (RooAbsPdf*) ws->factory("SUM::sbPDF( yieldSig1*rooPDF_signal_interpol ,yieldBkg1*rooPDF_background )");  // w morphing 
////      

// Definition of ModefConfig (sb and b) 
// sbModel
       ModelConfig* sbModel = new ModelConfig("sbModel");
       sbModel->SetWorkspace(*ws);
       sbModel->SetPdf("sbPDF");
       sbModel->SetObservables(the_Obs);
      sbModel->SetParametersOfInterest(yieldSig1); 
       sbModel->SetSnapshot(yieldSig1);

// bModel (= sbModel with s=0)
       ModelConfig* bModel = new ModelConfig("bModel");
       bModel->SetWorkspace(*ws);
       bModel->SetPdf("sbPDF");  // Here we take sb 
       bModel->SetObservables(the_Obs);
       yieldSig1.setVal(0.);
       bModel->SetParametersOfInterest(yieldSig1); 
       bModel->SetSnapshot(yieldSig1);

// here I use the simplest test statistics (LEP)
   SimpleLikelihoodRatioTestStat slrts(*bModel->GetPdf(),*sbModel->GetPdf());
   slrts.SetNullParameters(*bModel->GetSnapshot());
   slrts.SetAltParameters(*sbModel->GetSnapshot());
   
   
    HybridCalculator* hc = new HybridCalculator(*rooData_background, *sbModel, *bModel);
    
// Here 2 nuisance paramters: bkg scaling and morphing of signal 
    RooRealVar bkg0("bkg0", "bkg0", yieldBkg_input); 
    bkg0.setConstant(kTRUE);         
    ws->import(bkg0);
    RooRealVar sigma_bkg0("sigma_bkg0", "sigma_bkg0", 0.00000000000001);  // no error
    ws->import(sigma_bkg0);
    RooAbsPdf* nuisPdf_bkg  = (RooAbsPdf*) ws->factory("Gaussian::nuisPdf_bkg(yieldBkg1,bkg0,sigma_bkg0)");  

    RooAbsPdf* nuisPdf_alpha  = (RooAbsPdf*) ws->factory("Gaussian::nuisPdf_alpha(alpha,0.,1.)");
    RooAbsPdf* nuisPdfs =  (RooAbsPdf*) ws->factory("PROD::nuisPdfs(nuisPdf_bkg,nuisPdf_alpha)");
    hc->ForcePriorNuisanceAlt(*nuisPdfs);
    hc->ForcePriorNuisanceNull(*nuisPdfs);
   
//// Define Toy Sampler
    ToyMCSampler *toymcs = (ToyMCSampler*)hc->GetTestStatSampler();
    toymcs->SetTestStatistic(&slrts);
    toymcs->SetGenerateBinned(true); // to speed up 

    hc->SetToys(5000,5000);
        
//// calculate by running ntoys for the s+b and b hypothesis and retrieve the result
    HypoTestResult* myHypoTestResult = hc->GetHypoTest();
    myHypoTestResult->Print();
    
// // Inversion
///////////////

    std::cout<<" !! HypoTestInversion !!"<<std::endl;

// Create first an HypoTestCalculator with Null hypothesis = s+b model
    HybridCalculator* hcinv = new HybridCalculator(*rooData_background, *bModel, *sbModel);

    hcinv->ForcePriorNuisanceAlt(*nuisPdfs);
    hcinv->ForcePriorNuisanceNull(*nuisPdfs);
    
    
    ToyMCSampler *toymcsinv = (ToyMCSampler*)hcinv->GetTestStatSampler();
    SimpleLikelihoodRatioTestStat slrts1(*sbModel->GetPdf(),*bModel->GetPdf());
    slrts1.SetNullParameters(*sbModel->GetSnapshot());
    slrts1.SetAltParameters(*bModel->GetSnapshot());
   
    toymcsinv->SetTestStatistic(&slrts1);
    toymcsinv->SetGenerateBinned(true); // to speed up 
    hcinv->SetToys(1000,1000);
     
    HypoTestInverter InvHypo(*hcinv);

    InvHypo.SetConfidenceLevel(0.95);
  
  
    InvHypo.UseCLs(true);
    InvHypo.SetVerbose(true);

// configure and run the scan
    InvHypo.SetFixedScan(5,10.,100.);  
    HypoTestInverterResult* results = InvHypo.GetInterval();
   
   
// get result and plot it
   double upperlimit = results->UpperLimit();
   std::cout<<"UpperLimit = "<<upperlimit<<std::endl;
   HypoTestInverterPlot hcinvplot("hcinvplot","",results);
   hcinvplot.Draw();
   
   


}