#include "TROOT.h"
#include <RooFit.h>
#include <RooRealVar.h>
#include <RooArgSet.h>
#include <RooGaussian.h>
#include <RooExponential.h>
#include <RooProdPdf.h>
#include <RooAbsPdf.h>
#include <RooDataSet.h>
#include <RooPlot.h>
#include <RooCmdArg.h>
#include <RooNLLVar.h>
#include <RooMinuit.h>

#include "TCanvas.h"

int roofitTest()
{
  RooRealVar * m_m = new RooRealVar("mass","mass",5200,5600,"MeV/c^{2}");
  RooRealVar * m_t = new RooRealVar("t","t",-0.5,10,"ps");

  RooRealVar *m_meanB = new RooRealVar("mean_B","mean_B",0.,-0.02,0.02);
  RooRealVar *m_sigB = new RooRealVar("sig_B","sig_B",0.04,0.02,0.06);
  RooGaussian  *m_B_decay = new RooGaussian("m_B_decay","m_B_decay",
                                        *m_t,*m_meanB,*m_sigB);

  RooRealVar *m_slope = new RooRealVar("slope_B","slope_B",-0.5,-1,-0.6);
  RooExponential  *m_B_mass = new RooExponential("m_B_mass","m_B_mass",
                                            *m_m,*m_slope);
  
  RooAbsPdf  *m_B_pdf = NULL;
  
  m_B_pdf =  new RooProdPdf("m_B_pdf","m_B_pdf",
                          RooArgSet(*m_B_decay,*m_B_mass));
  //m_B_pdf =m_B_mass;
  
  RooArgSet  *m_observables = new RooArgSet(*m_t,*m_m);
  //RooArgSet  *m_observables = new RooArgSet(*m_m);
  
  RooDataSet *m_data = m_B_pdf->generate(*m_observables,10000);
  

  RooNLLVar  *m_bkg_fit = new RooNLLVar("Bkg fit","Bkg fit",*m_B_pdf,
                                        *m_data);
                                        //, RooFit::NumCPU(2));
  RooMinuit minuit_Bkg(*m_bkg_fit);
  minuit_Bkg.setStrategy(1);
  minuit_Bkg.optimizeConst(kTRUE);
  minuit_Bkg.setProfile(kTRUE);
  minuit_Bkg.setVerbose(kFALSE);
  minuit_Bkg.setEps(1e-20);
  minuit_Bkg.migrad();
  minuit_Bkg.hesse();

  TCanvas* c2 = new TCanvas("c2","sidebands");
  c2->cd();
  c2->SetLogy();
  RooPlot* m_t_frame_side = m_t->frame(-0.5,10.,100);
  m_data->plotOn(m_t_frame_side);
  m_B_pdf->plotOn(m_t_frame_side);
  m_t_frame_side->SetName("sidefit");
  m_t_frame_side->Draw();

  TCanvas* c3 = new TCanvas("c3","sidebands");
  c3->cd();
  RooPlot* m_m_frame_side = m_m->frame();
  m_data->plotOn(m_m_frame_side);
  m_B_pdf->plotOn(m_m_frame_side);
  m_m_frame_side->SetName("msidefit");
  m_m_frame_side->Draw();

  return 0;
}