#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussModel.h"
#include "RooAddModel.h"
#include "RooTruthModel.h"
#include "TPaveText.h"
#include "RooDecay.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TH1.h"
using namespace RooFit ;

void testdecay()
{
	RooRealVar t("t","t",100,200) ;
	RooRealVar tau("tau","tau",10,0,100) ;
	RooRealVar tau2("tau2","tau2",25) ;
	RooRealVar t_shift("t_shift","",110.0,100,200) ;
	RooRealVar t_shift2("t_shift2","",150.0) ;
	RooRealVar c_1("c_1","c1",15/2.35,1,100) ;
	RooRealVar c_2("c_2","c2",5/2.35) ;
	RooGaussModel gauss1("gauss1","res. gauss 1",t,t_shift,c_1) ;
	RooDecay dec_1("dec_1","decay",t,tau,gauss1,RooDecay::SingleSided) ;
	RooDecay dec_2("dec_2","decay2",t,tau2,gauss1,RooDecay::SingleSided) ;
	RooGaussModel gauss2("gauss2","something to the right of the peak",t,t_shift2,c_2) ;
	RooRealVar bg("bg","background",0.01,0,1);
	RooRealVar I2("I2","I2",0.01);
	RooRealVar Idec2("Idec2","Idec2",0.2);
	RooPolynomial bkg("bkg","bkg",t) ;
	RooAddPdf dec_bkg("dec_bkg","decay.sum w.bg",RooArgList(bkg,dec_1),bg); 
	RooAddPdf dec_bkg2("dec_bkg2","decay.sum w.bg",RooArgList(bkg,gauss2,dec_2,dec_1),RooArgList(bg,I2,Idec2)); 
	RooDataSet* modelData = dec_bkg2.generate(t,100000) ;
	RooDataHist* dh = modelData->binnedClone() ;
        t.setRange("sig1",100,128) ;  
        t.setRange("sig2",130,200) ;  

	RooChi2Var chi2("chi2","chi2",dec_bkg,*dh,RooFit::Range("sig1,sig2"));
  // Use RooMinuit interface to minimize chi^2
	RooMinuit m(chi2) ;
  	m.migrad() ;
  	m.improve();
  	m.hesse() ;

//  chi2.chi2FitTo(*dh);//,Range("sig1,sig2")); //doesnt work at all

	RooPlot* frame = t.frame() ;
	dh->plotOn(frame);
//	dec_bkg.plotOn(frame,RooFit::Range("sig1,sig2"));
	dec_bkg.plotOn(frame,RooFit::Range("sig1,sig2"),RooFit::NormRange("sig1,sig2"),RooFit::LineColor(kGreen));
	dec_bkg.plotOn(frame,RooFit::LineColor(kRed));
//	frame->Draw() ;
	RooPlot* frame2 = t.frame(RooFit::Title("chi2 residuals")) ;
	RooHist* hresid = frame->pullHist();
        frame2->addPlotable(hresid,"P");
        TCanvas* c = new TCanvas("rf","rf",850, 500);
        c->Divide(2);
        c->cd(1);frame->Draw();
	c->cd(2);frame2->Draw();

	return;
}