#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooPlot.h"
#include "RooAbsPdf.h"
#include "TROOT.h"

using namespace RooFit;

int fit(){

  TFile *f = new TFile("data.root","READ");
  TTree *t = (TTree*)f->Get("data");
  

  RooRealVar var1("var1","var1",5.25,5.287);
  RooRealVar var2("var2","var2",-0.2,0.2);
  RooRealVar var3("var3","var3",-3,1.5);


  RooDataSet *data = new RooDataSet("data","data",t, 
				    RooArgSet(var1,var2,var3));

  
  //var1 PDF signal

  RooRealVar var1m1_sig("var1m1_sig","var1m1_sig",5.279602);
  RooRealVar var1m2_sig("var1m2_sig","var1m2_sig",5.2774);
  RooRealVar var1s1_sig("var1s1_sig","var1s1_sig",0.002610);
  RooRealVar var1s2_sig("var1s2_sig","var1s2_sig",0.0096);
  RooRealVar var1f_sig("var1f_sig","var1f_sig",0.931);

  RooGaussian var1g1_sig("var1g1_sig","var1g1_sig",var1,var1m1_sig,var1s1_sig);
  RooGaussian var1g2_sig("var1g2_sig","var1g2_sig",var1,var1m2_sig,var1s2_sig);

  RooAddPdf var1pdf_sig("var1pdf_sig","var1pdf_sig",
		       RooArgList(var1g1_sig,var1g2_sig),
		       RooArgList(var1f_sig)); 


  //var2 PDF signal

  RooRealVar var2m1_sig("var2m1_sig","var2m1_sig",-0.005933);
  RooRealVar var2m2_sig("var2m2_sig","var2m2_sig",-0.02994);
  RooRealVar var2s1_sig("var2s1_sig","var2s1_sig",0.02673);
  RooRealVar var2s2_sig("var2s2_sig","var2s2_sig",0.0928);
  RooRealVar var2f_sig("var2f_sig","var2f_sig",0.843);

  RooGaussian var2g1_sig("var2g1_sig","var2g1_sig",var2,var2m1_sig,var2s1_sig);
  RooGaussian var2g2_sig("var2g2_sig","var2g2_sig",var2,var2m2_sig,var2s2_sig);

  RooAddPdf var2pdf_sig("var2pdf_sig","var2pdf_sig",
		       RooArgList(var2g1_sig,var2g2_sig),
		       RooArgList(var2f_sig)); 

  //var3 PDF signal

  RooRealVar var3m1_sig("var3m1_sig","var3m1_sig",-1.013);
  RooRealVar var3m3_sig("var3m2_sig","var3m2_sig",-0.4930);
  RooRealVar var3s1_sig("var3s1_sig","var3s1_sig",0.271);
  RooRealVar var3s2_sig("var3s2_sig","var3s2_sig",0.525);
  RooRealVar var3s3_sig("var3s3_sig","var3s3_sig",0.313);
  RooRealVar var3f_sig("var3f_sig","var3f_sig",0.14);

  RooBifurGauss var3g1_sig("var3g1_sig","var3g1_sig",var3,var3m1_sig,
			  var3s1_sig,var3s2_sig);
  RooGaussian var3g2_sig("var3g2_sig","var3g2_sig",var3,var3m3_sig,var3s3_sig);

  RooAddPdf var3pdf_sig("var3pdf_sig","var3pdf_sig",
		       RooArgList(var3g1_sig,var3g2_sig),
		       RooArgList(var3f_sig)); 

  //var1 PDF background

  RooRealVar var1xi_bkg("var1xi_bkg","var1xi_bkg",-24.8170);
  RooRealVar var1ep_bkg("var1ep_bkg","var1ep_bkg",5.2893);

  RooArgusBG var1pdf_bkg("var1pdf_bkg","var1pdf_bkg",var1,
			 var1ep_bkg,var1xi_bkg);

  //var2 PDF background

  RooRealVar var2c1_bkg("var2c1_bkg","var2c1_bkg",-0.3119);
  
  RooChebychev var2pdf_bkg("var2pdf_bkg","var2pdf_bkg",var2,
			   RooArgList(var2c1_bkg));

  //var3 PDF background

  RooRealVar var3m1_bkg("var3m1_bkg","var3m1_bkg",-1.1117);
  RooRealVar var3s1_bkg("var3s1_bkg","var3s1_bkg",0.3225);

  RooGaussian var3pdf_bkg("var3pdf_bkg","var3pdf_bkg",
			  var3,var3m1_bkg,var3s1_bkg);


  //product pdfs for signal and background

  RooProdPdf pdf_sig("pdf_sig","pdf_sig",RooArgList(var1pdf_sig,
						    var2pdf_sig,
						    var3pdf_sig));

  RooProdPdf pdf_bkg("pdf_bkg","pdf_bkg",RooArgList(var1pdf_bkg,
						    var2pdf_bkg,
						    var3pdf_bkg));

  //yields declariation and fit pdf

  RooRealVar nsig("nsig","nsig",500,0,5000);
  RooRealVar nbkg("nbkg","nbkg",500,0,5000);

  RooAddPdf *fitPDF = new RooAddPdf("fitPDF","fitPDF",
				    RooArgList(pdf_sig,pdf_bkg),
				    RooArgList(nsig,nbkg));

  //fit and print results
  fitPDF->fitTo(*data,Extended(kTRUE),Minos(kTRUE));
  fitPDF->Print("t");

  RooArgSet *params = fitPDF->getVariables();
  params->Print("v");

  cout <<endl;
  cout <<"number of sig is " <<nsig.getVal() <<"+/-" <<nsig.getError()
       <<"( "<<nsig.getAsymErrorLo() <<" +" <<nsig.getAsymErrorHi() 
       <<")" <<endl;
  cout <<"number of bkg is " <<nbkg.getVal() <<"+/-" <<nbkg.getError()
       <<"( "<<nbkg.getAsymErrorLo() <<" +" <<nbkg.getAsymErrorHi() 
       <<")" <<endl;


  //draw projections
  TCanvas *c1 = new TCanvas("c1","c1",1);
  c1->Divide(3);
  c1->Draw();

  c1->cd(1);
  RooPlot *var1frame = var1.frame();
  data->plotOn(var1frame,Binning(25));
  fitPDF->plotOn(var1frame);
  fitPDF->plotOn(var1frame,Components(pdf_bkg),
		 LineStyle(kDashed));
  var1frame->Draw();

  c1->cd(2);
  RooPlot *var2frame = var2.frame();
  data->plotOn(var2frame,Binning(25));
  fitPDF->plotOn(var2frame);
  fitPDF->plotOn(var2frame,Components(pdf_bkg),
		 LineStyle(kDashed));
  var2frame->Draw();

  c1->cd(3);
  TPad *pad = gPad;
  pad->SetLogy();
  RooPlot *var3frame = var3.frame();
  data->plotOn(var3frame,Binning(25));
  fitPDF->plotOn(var3frame);
  fitPDF->plotOn(var3frame,Components(pdf_bkg),
		 LineStyle(kDashed));
  var3frame->Draw();

}