#include "ostream"
#include "TFile.h"
#include "TTreeReader.h"
#include "TTreeReaderArray.h"
#include "TLorentzVector.h"
#include "TROOT.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TPostScript.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TPaveStats.h"
#include "TPaveLabel.h"
#include "TPaveText.h"
#include "TF1.h"
#include "TMath.h"
#include "TAxis.h"
#include "TGraphErrors.h"
#include "TMultiGraph.h"
#include "TLegend.h"
#include "TLatex.h"
#include <math.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include "TTree.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TLorentzVector.h"
#include "TVector3.h"
#include "TChain.h"
#include "Fit/FitResult.h"
#include "TFile.h"
#include "TTree.h"
#include "TNtuple.h"
#include "TCanvas.h"
#include "TMath.h"
#include "RooRealVar.h"
#include "RooGaussian.h"
#include "RooAddPdf.h"
#include "RooDataSet.h"
#include "RooChebychev.h"
#include "RooPlot.h"
#include "RooDataHist.h"

using namespace RooFit;
void HistoSum_Zroofit(){
 gROOT->Reset();

//  gStyle->SetOptTitle(kFALSE);
//  gStyle->SetOptStat(000000000);

 TFile *f1 = TFile::Open("hist_data.root");

 gDirectory->ls();
 TH1F* d_z = (TH1F*)f1->Get("k11");
 
 RooRealVar di_mu("di_mu", "di_mu",  60,120);

// Assigning Histogram dataset to di_mu object
//---------------------------------------------

RooDataHist data("data", "dataset",RooArgSet(di_mu) ,d_z);

// Signal Model and Parameters, here we define two Gaussians PDFs for signal region
// ------------------------------------------------------------------------------
   RooRealVar mean("mean","mean of gaussians",91) ;
   RooRealVar sigma1("sigma1","width of gaussians",1.5) ;
   RooRealVar sigma2("sigma2","width of gaussians",2.0) ;

   RooGaussian sig1("sig1","Signal component 1",di_mu ,mean,sigma1) ;  
   RooGaussian sig2("sig2","Signal component 2",di_mu, mean,sigma2) ;  

// Background Model and parameters (Chebychev polynomial PDF)
//-----------------------------------------------------------

   RooRealVar a0("a0","a0",0.5,0.,1.) ;
   RooRealVar a1("a1","a1",0.6,0.,1.) ;
   RooChebychev bkg("bkg","Background",di_mu,RooArgSet(a0,a1)) ;
   //RooChebychev bkg("bkg","Background",di_mu,a0) ;


// Sum the signal components into a composite signal PDF
//---------------------------------------------------------

   RooRealVar sig1frac("sig1frac","fraction of component 1 in signal",0.8,0.,1.) ;
   RooAddPdf sig("sig","Signal",RooArgList(sig1,sig2),sig1frac) ;

// Define signal range in which events counts are to be defined
//-------------------------------------------------------------

di_mu.setRange("signalRange",80,100) ;

//  Associated nsig/nbkg as expected number of events with sig & bkg in the range  "signalRange" 
//---------------------------------------------------------------------------------------------
   RooRealVar nsig("nsig","number of signal events in signalRange",5000,0.,1000000) ;
   RooRealVar nbkg("nbkg","number of background events in signalRange",3000,0,1000000) ;
   RooRealVar bkgfrac("bkgfrac","fraction of background",0.5,0.,1.) ;
// Use AddPdf to extend the model:
//--------------------------------

RooAddPdf  model("model","(g1+g2)+a", RooArgList(bkg,sig), RooArgList(nbkg,nsig)) ;

// Need to clone these models because the interpretation of normalisation coefficients changes when different ranges are used: 
//---------------------------------------------------------------------------------------------------------------------------
 
   RooAddPdf model2(model);
   RooAddPdf model3(model);

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// <|>  S a m p l e   d a t a ,   F i t   m o d e l  <|>
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

 

 TCanvas *c5 = new TCanvas("c5", "c5");


//  Fit full range
// ---------------------
 c5->cd();
RooFitResult* r = model.fitTo(data,Save()) ;
   r->Print() ;
// Plotting the full range
// -----------------------

RooPlot * frame = di_mu.frame(Title("Full range fitted"));
   data.plotOn(frame);
   model.plotOn(frame, VisualizeError(*r));
   model.plotOn(frame);
   model.paramOn(frame);
   data.statOn(frame);
   frame->Draw();
c5->Draw();

c5->SaveAs("Z1roofit_mass.png");

TCanvas *c6 = new TCanvas("c6", "c6");

   // Fit in the left and right sideband regions
  // -------------------------------------------

c6->cd();
   di_mu.setRange("left", 60., 80.);
   di_mu.setRange("right", 100., 120.);
   
   RooFitResult* r2 = model2.fitTo(data,Range("left,right"),Save()) ;
   r2->Print();
   
// Plotting the left and right sidebands
// -------------------------------------
   
   RooPlot * frame2 = di_mu.frame(Title("Fit in left/right sideband"));
   data.plotOn(frame2);
   model2.plotOn(frame2, VisualizeError(*r2));
   model2.plotOn(frame2);
   model2.paramOn(frame2);
   data.statOn(frame2);
   frame2->Draw();
   

c6->Draw(); 
 
c6->SaveAs("Z2roofit_mass.png");
}