#include "TStopwatch.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TMath.h"
#include "RooPlot.h"
#include "RooAbsPdf.h"
#include "RooWorkspace.h"
#include "RooDataSet.h"
#include "RooGlobalFunc.h"
#include "RooFitResult.h"
#include "RooRandom.h"
#include "RooPoisson.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/ToyMCSampler.h"
#include "RooStats/ProfileLikelihoodTestStat.h"
#include <vector>
#include <string>

using namespace RooFit;
using namespace RooStats;


std::vector<double> runPL(RooWorkspace wspace, RooDataSet* data, const char *pdfName,std::vector<RooRealVar*> theVars);

//Set the confidence level
double CL = 0.6827;

/*
  ======
  TO DO: (24 Nov. 2013)
  ======

  - add total cross-section measurement (now only fiducial)  
  - add nuisance parameters for all systematics (the loop is ready)
  - add estimation of lower systematic uncertainty (now only upper)
  - make code run for many channels, now only one (maybe create histFactory with results from each channel and combine with another script?)


 */

void myPL(){

  gROOT->ForceStyle();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(000);

  TH1F *h;

  //************************************************************************
  //Set true for fiducial or total x section measurement
  //************************************************************************
  bool doFiducialXsection = true; //if false, extract total cross-section
  //// bool totalXsection = false;


  //ZZ -> 4mu

  //Set the number of observed events
  double NObs = 85;
  
  std::vector<string> nuisPar_name;
  std::vector<string> nuisPar_title;
  std::vector<double> nuisPar_val;
  std::vector<double> nuisPar_staterr;
  std::vector<double> nuisPar_systerr;

  //
  nuisPar_name.push_back("b");
  nuisPar_title.push_back("background");
  nuisPar_val.push_back(1.1);
  nuisPar_staterr.push_back(1.4);
  nuisPar_systerr.push_back(0.5);

  //
  nuisPar_name.push_back("Czz");
  nuisPar_title.push_back("Czz");
  nuisPar_val.push_back(0.83);
  nuisPar_staterr.push_back(0);
  nuisPar_systerr.push_back(0.03);

  //
  nuisPar_name.push_back("Lumi");
  nuisPar_title.push_back("luminosity");
  nuisPar_val.push_back(20.3);
  nuisPar_staterr.push_back(0);
  nuisPar_systerr.push_back(20.3*0.028);
  

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

  //create the workspace
  RooWorkspace* wspace = new RooWorkspace("wspace");
  /// variables of this workspace: n, b, Czz, Azz, Lumi, br

  //create parameters, set bounds according to errors or to dummy values
  RooRealVar* n = (RooRealVar*) wspace->factory("n[0., 1000.]"); //number of observed events
  n->setBins(1);
  n->setVal(NObs);
  // First create empty data set, set data value and put into workspace
  RooDataSet* data = new RooDataSet("data", "data", RooArgSet(*n));  // empty 
  data->add(*n);
  wspace->import(*data);

  std::vector<RooRealVar* > vars;
  for (unsigned int i=0; i<nuisPar_name.size(); i++) {
    double x=nuisPar_val[i];
    double xerr=sqrt(pow(nuisPar_staterr[i],2)+pow(nuisPar_systerr[i],2));
    double xmin=x-xerr;
    double xmax=x+xerr;
    if (xmin<0) {
      cout << "WARNING! lower bound for parameter " <<  nuisPar_name[i] << " is negative (" << xmin << ")" << endl;
    }
    vars.push_back( 
		   (RooRealVar*) wspace->factory(Form("%s[%f,%f,%f]",nuisPar_name[i].c_str(),x,xmin,xmax))
		    );
    
  }
  cout << "============================ " << endl;
  cout << "Created nuisance parameters: " << endl;
  cout << "============================ " << endl;
 
  for (unsigned int i=0; i<nuisPar_name.size(); i++) {
    cout << vars[i]->GetName() << " " << vars[i]->getVal() << " " << vars[i]->getMin() << " " << vars[i]->getMax() << endl;;
  }

}