/// \file
/// \ingroup tutorial_roostats
/// \notebook
/// Neutrino Oscillation Example from Feldman & Cousins
///
/// This tutorial shows a more complex example using the FeldmanCousins utility
/// to create a confidence interval for a toy neutrino oscillation experiment.
/// The example attempts to faithfully reproduce the toy example described in Feldman & Cousins'
/// original paper, Phys.Rev.D57:3873-3889,1998.
///
/// \macro_image
/// \macro_output
/// \macro_code
///
/// \author Kyle Cranmer

#include "RooGlobalFunc.h"
#include "RooStats/ConfInterval.h"
#include "RooStats/FeldmanCousins.h"
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/MCMCCalculator.h"
#include "RooStats/UniformProposal.h"
#include "RooStats/LikelihoodIntervalPlot.h"
#include "RooStats/MCMCIntervalPlot.h"
#include "RooStats/MCMCInterval.h"

#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooRealVar.h"
#include "RooConstVar.h"
#include "RooAddition.h"
#include "RooProduct.h"
#include "RooProdPdf.h"
#include "RooAddPdf.h"

#include "TROOT.h"
#include "RooPolynomial.h"
#include "RooRandom.h"

#include "RooNLLVar.h"
#include "RooProfileLL.h"

#include "RooPlot.h"

#include "TCanvas.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TTree.h"
#include "TMarker.h"
#include "TStopwatch.h"
#include "Riostream.h"

#include "RooAbsReal.h"
#include "RooAbsCategory.h"
#include <math.h>
#include "TMath.h"

#include <iostream>

using namespace RooFit;
using namespace RooStats;

void FC_2D_test_8_ws(bool doFeldmanCousins = true)

{

   //gROOT->SetBatch(1);

   // to time the macro
   TStopwatch t;
   t.Start();

   // R e a d   w o r k s p a c e   f r o m   f i l e
   // -----------------------------------------------

   // Open input file with workspace 
   TFile *f = new TFile("wspace_summaryLowQ2.root");

   // Retrieve workspace from file
   RooWorkspace *w = (RooWorkspace *)f->Get("wspace.summaryLowQ2");

   // R e t r i e v e   p d f ,   d a t a   f r o m   w o r k s p a c e
   // -----------------------------------------------------------------

   // Retrieve variables, model and data from workspace

   RooRealVar *Bmass = w->var("Bmass");

   RooRealVar *CosThetaK = w->var("CosThetaK");

   RooRealVar *CosThetaL = w->var("CosThetaL");

   RooArgSet variables(*Bmass, *CosThetaK, *CosThetaL);

 ///////////////////////////////POIs////////////////////////////////////////////////

   RooRealVar *unboundAfb = w->var("unboundAfb");

   RooRealVar *unboundFl = w->var("unboundFl");

   RooAbsReal *afb = w->function("afb");

   RooAbsReal *fl = w->function("fl");

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

   RooRealVar *nSig = w->var("nSig");
   nSig->setVal(122.725481956) ;
   nSig->setConstant(kTRUE) ;

   RooRealVar *nBkgComb = w->var("nBkgComb");
   nBkgComb->setVal(266.314687526) ;
   nBkgComb->setConstant(kTRUE) ;

   RooRealVar *PeakFrac = w->var("PeakFrac");
   PeakFrac->setVal(0.177102920935) ;
   PeakFrac->setConstant(kTRUE) ;

/////////////////////////////MC FITTED PARAMETERS///////////////////////////////////

   w->var("cbs1_sigma")->setVal(0.0218718001495) ;
   w->var("cbs1_sigma")->setConstant(kTRUE) ;

   w->var("cbs2_sigma")->setVal(0.041364965623) ;
   w->var("cbs2_sigma")->setConstant(kTRUE) ;

   w->var("cbs1_alpha")->setVal(2.30354878805) ;
   w->var("cbs1_alpha")->setConstant(kTRUE) ;

   w->var("cbs2_alpha")->setVal(-1.69224048367) ;
   w->var("cbs2_alpha")->setConstant(kTRUE) ;

   w->var("cbs1_n")->setVal(0.0911597192092) ;
   w->var("cbs1_n")->setConstant(kTRUE) ;

   w->var("cbs2_n")->setVal(10.0) ;
   w->var("cbs2_n")->setConstant(kTRUE) ;

   w->var("sigMDCB_frac")->setVal(0.407336326511) ;
   w->var("sigMDCB_frac")->setConstant(kTRUE) ;


   w->var("cbs_mean_pb")->setVal(5.38266552896) ;
   w->var("cbs_mean_pb")->setConstant(kTRUE) ;

   w->var("cbs1_sigma_pb")->setVal(0.0703155900621) ;
   w->var("cbs1_sigma_pb")->setConstant(kTRUE) ;

   w->var("cbs2_sigma_pb")->setVal(0.0341134981745) ;
   w->var("cbs2_sigma_pb")->setConstant(kTRUE) ;

   w->var("cbs1_alpha_pb")->setVal(1.68570454302) ;
   w->var("cbs1_alpha_pb")->setConstant(kTRUE) ;

   w->var("cbs2_alpha_pb")->setVal(-0.320785330534) ;
   w->var("cbs2_alpha_pb")->setConstant(kTRUE) ;

   w->var("cbs1_n_pb")->setVal(1.49648286393e-07) ;
   w->var("cbs1_n_pb")->setConstant(kTRUE) ;

   w->var("cbs2_n_pb")->setVal(119.412743734) ;
   w->var("cbs2_n_pb")->setConstant(kTRUE) ;

   w->var("MDCB_frac_pb")->setVal(0.477642358082) ;
   w->var("MDCB_frac_pb")->setConstant(kTRUE) ;


   w->var("u1")->setVal(0.0241661006756) ;
   w->var("u1")->setConstant(kTRUE) ;

   w->var("u2")->setVal(0.837791050259) ;
   w->var("u2")->setConstant(kTRUE) ;

   w->var("u3")->setVal(-0.0479622573662) ;
   w->var("u3")->setConstant(kTRUE) ;

   w->var("u4")->setVal(-0.569056480793) ;
   w->var("u4")->setConstant(kTRUE) ;

   w->var("v1")->setVal(0.000143574127913) ;
   w->var("v1")->setConstant(kTRUE) ;

   w->var("v2")->setVal(-0.992878066968) ;
   w->var("v2")->setConstant(kTRUE) ;

   w->var("v3")->setVal(-0.000153869514138) ;
   w->var("v3")->setConstant(kTRUE) ;

   w->var("v4")->setVal(-0.00710422751506) ;
   w->var("v4")->setConstant(kTRUE) ;


   w->var("effi_norm")->setVal(0.0244680409599) ;
   w->var("effi_norm")->setConstant(kTRUE) ;

   w->var("z0")->setVal(-0.957584369326) ;
   w->var("z0")->setConstant(kTRUE) ;

   w->var("z1")->setVal(5.81156425241e-05) ;
   w->var("z1")->setConstant(kTRUE) ;

   w->var("z2")->setVal(-0.00790898695679) ;
   w->var("z2")->setConstant(kTRUE) ;

   w->var("z3")->setVal(5.55490118899e-05) ;
   w->var("z3")->setConstant(kTRUE) ;

   w->var("z4")->setVal(0.00534212060178) ;
   w->var("z4")->setConstant(kTRUE) ;

   w->var("z5")->setVal(0.000361181052416) ;
   w->var("z5")->setConstant(kTRUE) ;

   w->var("z6")->setVal(0.00115440071328) ;
   w->var("z6")->setConstant(kTRUE) ;

   w->var("z7")->setVal(-0.000119737359302) ;
   w->var("z7")->setConstant(kTRUE) ;

   w->var("z8")->setVal(-0.000179452769173) ;
   w->var("z8")->setConstant(kTRUE) ;

   w->var("z9")->setVal(0.00248231849833) ;
   w->var("z9")->setConstant(kTRUE) ;

   w->var("z10")->setVal(0.0291778527321) ;
   w->var("z10")->setConstant(kTRUE) ;

   w->var("z11")->setVal(-6.3185579279e-05) ;
   w->var("z11")->setConstant(kTRUE) ;

   w->var("z12")->setVal(-0.00368004872762) ;
   w->var("z12")->setConstant(kTRUE) ;

   w->var("z13")->setVal(-0.0016219301288) ;
   w->var("z13")->setConstant(kTRUE) ;

   w->var("z14")->setVal(-0.00419420157757) ;
   w->var("z14")->setConstant(kTRUE) ;

   w->var("z15")->setVal(-2.51346462257e-05) ;
   w->var("z15")->setConstant(kTRUE) ;
   
   w->var("z16")->setVal(-0.00255093697933) ;
   w->var("z16")->setConstant(kTRUE) ;

   w->var("z17")->setVal(-0.000562549977301) ;
   w->var("z17")->setConstant(kTRUE) ;

   w->var("z18")->setVal(0.00241213396072) ;
   w->var("z18")->setConstant(kTRUE) ;

   w->var("z19")->setVal(-0.0058014926427) ;
   w->var("z19")->setConstant(kTRUE) ;

   w->var("z20")->setVal(-0.0639044613716) ;
   w->var("z20")->setConstant(kTRUE) ;

   w->var("z21")->setVal(-0.000563336289174) ;
   w->var("z21")->setConstant(kTRUE) ;

   w->var("z22")->setVal(0.00875756548499) ;
   w->var("z22")->setConstant(kTRUE) ;

   w->var("z23")->setVal(0.00300495062685) ;
   w->var("z23")->setConstant(kTRUE) ;

   w->var("z24")->setVal(0.00030525198677) ;
   w->var("z24")->setConstant(kTRUE) ;

   w->var("pb_K_c1")->setVal(-0.117775974278) ;
   w->var("pb_K_c1")->setConstant(kTRUE) ;

   w->var("pb_K_c2")->setVal(1.40945886463) ;
   w->var("pb_K_c2")->setConstant(kTRUE) ;

   w->var("pb_K_c3")->setVal(0.242942860423) ;
   w->var("pb_K_c3")->setConstant(kTRUE) ;
  
   w->var("pb_K_c4")->setVal(0.63259832486) ;
   w->var("pb_K_c4")->setConstant(kTRUE) ;

   w->var("pb_L_c1")->setVal(-0.00252408213876) ;
   w->var("pb_L_c1")->setConstant(kTRUE) ;

   w->var("pb_L_c2")->setVal(0.384953389845) ;
   w->var("pb_L_c2")->setConstant(kTRUE) ;

   w->var("pb_L_c3")->setVal(-1.22728157281) ;
   w->var("pb_L_c3")->setConstant(kTRUE) ;

   w->var("pb_L_c4")->setVal(0.0100004361491) ;
   w->var("pb_L_c4")->setConstant(kTRUE) ;

   w->var("pb_L_c5")->setVal(1.62255369868) ;
   w->var("pb_L_c5")->setConstant(kTRUE) ;


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

   RooAbsPdf *sigG = w->pdf("f_sig3DAltM");

   RooAbsPdf *BKG = w->pdf("f_bkgComb");

   RooAbsPdf *BKG_pb = w->pdf("f_pb_Kstar");

   RooAbsPdf *model = w->pdf("f_final_WithKstar");

   //RooAbsData *data = w->data("data");


   // Print structure of composite p.d.f.
   model->Print("t");

   // for debugging, check model tree
  model->printCompactTree();
  model->graphVizTree("model.dot");

   // turn off some messages
   RooMsgService::instance().setStreamStatus(0, kFALSE);
   RooMsgService::instance().setStreamStatus(1, kFALSE);
   RooMsgService::instance().setStreamStatus(2, kFALSE);

   // needed to remove some sporious messages from RooAddPdf when generating events
   //model->fixCoefNormalization(*Bmass);   

   // ----------------------------------------------------------------- //

   // n events in data to data, simply sum of sig+bkg
   Int_t nEventsData = nBkgComb->getVal() + nSig->getVal() + nSig->getVal()*PeakFrac->getVal() ;
   cout << "generate toy data with nEvents = " << nEventsData << endl;

   // adjust random seed to get a toy dataset similar to one in paper.
   // Found by trial and error (3 trials, so not very "fine tuned")

   RooRandom::randomGenerator()->SetSeed(3);

   // create a toy dataset
   //RooDataSet *data = model->generate(RooArgSet(*Bmass), nEventsData);
   RooDataSet *data = model->generate(variables, nEventsData);

   // ---------------------------------------------------------------- //

   //------------------------- make some plots------------------------//

   TCanvas *c = new TCanvas("c","c",1000,800);
   c->Divide(2, 2);


   // ------------------plot the data with the best fit------------ //

   c->cd(1);
   RooPlot *Bframe = Bmass->frame(Title("B_{s} invariant mass : toy data with best fit model"),Bins(30)); //50
   data->plotOn(Bframe);
   model->fitTo(*data, Extended());
   model->plotOn(Bframe);
   model->plotOn(Bframe, RooFit::Name("f_sig3DAltM"), Components("f_sig3DAltM"), /*FillStyle(3013), FillColor(kGreen), DrawOption("F")*/LineColor(kGreen) ); //pdf for only signal
   model->plotOn(Bframe, RooFit::Name("f_pb_Kstar"), Components("f_pb_Kstar"), /*FillStyle(3013), FillColor(kOrange), DrawOption("F"),*/ LineColor(kOrange)); //pdf for only signal
   model->plotOn(Bframe, RooFit::Name("f_bkgComb"), Components("f_bkgComb"),LineStyle(9), LineColor(kRed),LineWidth(5) );
   model->plotOn(Bframe);
   Bframe->Draw();

   c->cd(2);
   RooPlot *Kframe = CosThetaK->frame(Title("Cos#theta_{K} : toy data with best fit model"),Bins(20)); //50
   data->plotOn(Kframe);
   model->fitTo(*data, Extended());
   model->plotOn(Kframe);
   model->plotOn(Kframe, RooFit::Name("f_sig3DAltM"), Components("f_sig3DAltM"), /*FillStyle(3013), FillColor(kGreen), DrawOption("F")*/LineColor(kGreen) ); //pdf for only signal
   model->plotOn(Kframe, RooFit::Name("f_pb_Kstar"), Components("f_pb_Kstar"), /*FillStyle(3013), FillColor(kOrange), DrawOption("F"),*/ LineColor(kOrange)); //pdf for only signal
   model->plotOn(Kframe, RooFit::Name("f_bkgComb"), Components("f_bkgComb"),LineStyle(9), LineColor(kRed),LineWidth(5) );
   model->plotOn(Kframe);
   Kframe->Draw();

   c->cd(3);
   RooPlot *Lframe = CosThetaL->frame(Title("Cos#theta_{L} : toy data with best fit model"),Bins(20)); //50
   data->plotOn(Lframe);
   model->fitTo(*data, Extended());
   model->plotOn(Lframe);
   model->plotOn(Lframe, RooFit::Name("f_sig3DAltM"), Components("f_sig3DAltM"), /*FillStyle(3013), FillColor(kGreen), DrawOption("F")*/LineColor(kGreen) ); //pdf for only signal
   model->plotOn(Lframe, RooFit::Name("f_pb_Kstar"), Components("f_pb_Kstar"), /*FillStyle(3013), FillColor(kOrange), DrawOption("F"),*/ LineColor(kOrange)); //pdf for only signal
   model->plotOn(Lframe, RooFit::Name("f_bkgComb"), Components("f_bkgComb"),LineStyle(9), LineColor(kRed),LineWidth(5) );
   model->plotOn(Lframe);
   Lframe->Draw();

   //RooFormulaVar fl("fl", "0.5+TMath::ATan(*unboundFl)/TMath::Pi()", RooArgList(*unboundFl));

   //RooFormulaVar afb("afb", "2.0*(1.-fl)*TMath::ATan(*unboundAfb)/TMath::Pi()", RooArgList(*unboundAfb,fl));


   // ------------------plot the likelihood function--------------- //

/*   c->cd(4);
   RooNLLVar nll("nll", "nll", *model, *data, Extended());
   RooProfileLL pll("pll", "", nll, RooArgSet(*unboundAfb, *unboundFl));
   TH1 *hhh = pll.createHistogram("hhh", *unboundAfb, Binning(20), YVar(*unboundFl, Binning(20)), Scaling(kFALSE));
   hhh->SetLineColor(kViolet);
   hhh->SetTitle("Likelihood Function");
   hhh->Draw("surf"); */

   c->cd(4);
   RooNLLVar nll("nll", "nll", *model, *data, Extended());
   RooProfileLL pll("pll", "", nll, RooArgSet(*afb, *fl));
   TH1 *hhh = pll.createHistogram("hhh", *afb, Binning(20), YVar(*fl, Binning(20)), Scaling(kFALSE));
   hhh->SetLineColor(kViolet);
   hhh->SetTitle("Likelihood Function");
   hhh->Draw("surf");

}