#include "RooWorkspace.h"
#include "RooRealVar.h"
#include "RooPlot.h"
#include <iostream>
#include "RooDataSet.h"
#include "TFile.h"
#include "TTree.h"
#include "TList.h"
#include "TCanvas.h"
#include "RooFormulaVar.h"
#include "RooAddPdf.h"
#include "TAxis.h"
#include "RooCBShape.h"

using namespace RooFit;
using namespace std;

#define Open TFile::Open
void testScript()
{
  RooWorkspace w("w");

  Int_t myLow = 4700, myHigh = 6200, binwidth = 4;

  w.factory(Form("Lb_DTF_M_JpsiLConstr[%d,%d]",myLow,myHigh));

  RooRealVar *myVar = w.var("Lb_DTF_M_JpsiLConstr");

  myVar->setRange("simFit_window",5500,5740);

  w.var("Lb_DTF_M_JpsiLConstr")->setBins((Int_t)(myHigh-myLow)/binwidth);

  RooDataSet *ds_sim, *ds_sim_wt;
  Float_t mcNentries;
  
  w.factory("RooCBShape::Lb1_Run1(Lb_DTF_M_JpsiLConstr,mean_Run1[5619.6,5619,5621],"
	    "sigma_Run1[10.,0.,20.], alpha1_Run1[1.028, 0.8,1.3], 10.0)" );
  w.factory("RooCBShape::Lb2_Run1(Lb_DTF_M_JpsiLConstr,mean_Run1,sigma_Run1,alpha2_Run1[-1.097,-1.4,-0.7], 10.0)");
  w.factory("SUM::Lb_Run1(0.5*Lb1_Run1 , 0.5*Lb2_Run1)");

  TFile *fileIn = Open("myWS.root","READ");
  RooWorkspace *w1 = (RooWorkspace*)fileIn->Get("w");

  ds_sim_wt = (RooDataSet*)w1->data("ds_sim_wt");

  mcNentries = ds_sim_wt->sumEntries();
  cout<<"mcNentries = "<<mcNentries<<endl;

  w.import(*(ds_sim_wt));

  w.factory("Exponential::mcbkg_Run1(Lb_DTF_M_JpsiLConstr,tau_Run1[-1e-7,-1e-5,-1e-9])");
  
  w.factory(Form("mcnsig_Run1[%f,%f,%f]",mcNentries-5,mcNentries/2,mcNentries));

  w.factory(Form("mcnbkg_Run1[10.0, 1., %f]",mcNentries));

  w.factory("SUM::mcFit_Run1(mcnsig_Run1*Lb_Run1, mcnbkg_Run1*mcbkg_Run1)");
  
  (w.pdf("mcFit_Run1"))->fitTo(*(ds_sim_wt),Range("simFit_window"),Strategy(2),Extended(),SumCoefRange("simFit_window"),SumW2Error(kTRUE));

  TCanvas *sim_Run1      = new TCanvas("sim_Run1","sim_Run1");
  RooPlot *simframe_Run1 = new RooPlot(*(w.var("Lb_DTF_M_JpsiLConstr")),5500,5740,60);
  simframe_Run1->GetXaxis()->SetTitle("m_{J/#psi#Lambda}[MeV/#it{c}^{2}]");
  simframe_Run1->GetYaxis()->SetTitle(Form("Candidates/(4 MeV/#it{c}^{2})"));
  simframe_Run1->GetYaxis()->SetTitleOffset(0.75);

  (ds_sim_wt)->plotOn(simframe_Run1,Name("simdata_Run1"),DataError(RooAbsData::SumW2));
  (w.pdf("mcFit_Run1"))->plotOn(simframe_Run1);
  (w.pdf("Lb_Run1"))->plotOn(simframe_Run1,LineColor(kMagenta+2));
  (w.pdf("mcbkg_Run1"))->plotOn(simframe_Run1,LineColor(kRed));
  simframe_Run1->GetYaxis()->SetRangeUser(0.1,3000);
  simframe_Run1->Draw();
  sim_Run1->SetLogy();

  //  w.writeToFile("myWS.root",true);

}