/*
 * =====================================================================================
 *    Description:  
 *
 *        Version:  1.0
 *        Created:  12/27/2011 10:47:16 AM
 * =====================================================================================
 */

#include "RooCFunction1Binding.h"
#include "RooRealVar.h"
#include "RooDataHist.h"
#include "RooDecay.h"
#include "RooGaussModel.h"
#include "RooFFTConvPdf.h"
#include "RooPlot.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TStopwatch.h"
using namespace RooFit;
double exp_fun(double, double, double);
double gaus(double, double,double);

// define new class
class MyExpF : public RooAbsPdf {
public:
  MyExpF(const char *name, const char *title,
	      RooAbsReal& _x,
	      RooAbsReal& _E,
	      RooAbsReal& _tau,
	      RooAbsReal& _t0,
	      RooAbsReal& _s1);
  MyExpF(const MyExpF& other, const char* name=0) ; 
  virtual TObject* clone(const char* newname) const { return new MyExpF(*this,newname); } 
  inline virtual ~MyExpF() { }

protected:

  RooRealProxy x ;
  RooRealProxy E ;
  RooRealProxy tau ;
  RooRealProxy t0 ;
  RooRealProxy s1 ;
  Double_t evaluate() const ;

private:

 ClassDef(MyExpF,1) // Your description goes here...

};

 MyExpF::MyExpF(const char *name, const char *title, 
                        RooAbsReal& _x,
                        RooAbsReal& _E,
                        RooAbsReal& _tau,
                        RooAbsReal& _t0,
			RooAbsReal& _s1) :
   RooAbsPdf(name,title), 
   x("x","x",this,_x),
   E("E","E",this,_E),
   tau("tau","",this,_tau),
   t0("t0","",this,_t0),
   s1("s1","",this,_s1)
 { 
 } 

  MyExpF::MyExpF(const MyExpF& other, const char* name) :  
   RooAbsPdf(other,name), 
   x("x",this,other.x),
   E("E",this,other.E),
   tau("tau",this,other.tau),
   t0("t0",this,other.t0),
   s1("s1",this,other.s1)
 { 
}  //constructor called from clone function

 Double_t MyExpF::evaluate() const 
 { 
   // ENTER EXPRESSION IN TERMS OF VARIABLE ARGUMENTS HERE 
   return exp_fun(x,tau,t0)*gaus(E,0,s1) ; 
 } 

// end define new class


void interp_prot_2d1()
{
	RooRealVar t("t","time",-2,10);
	RooRealVar E("E","Energy",-5,5);
	RooRealVar tau("tau","decay const.",1.2);
	RooRealVar t0("t0","",0.);
	RooRealVar c_t("c_t","",0.25);
	RooRealVar c_E("c_E","",0.25);
	RooGaussModel g_t("g_t","",t,t0,c_t) ;
	
	RooPlot* fr1 = t.frame() ;
	RooPlot* fr2 = E.frame() ;


//make binding

	MyExpF mydecay("mydecay","",t,E,tau,t0,c_E);
//==== I want to convolve this only along t axis
//==== and here begins the strange ====
	RooFFTConvPdf exg("exg","exp (X) gauss",t,mydecay,g_t) ;
//==== Here ends the strange ====
	TStopwatch w;
	w.Start();	
	RooDataHist * data = exg.generateBinned(RooArgSet(t,E),10000) ;
	w.Stop();
	cout<<"generate binned time= "<<w.RealTime()<<endl;
//	fr1->Draw() ;
//  TH1* hh_data = data->createHistogram("t,E",50,50) ;
        TCanvas* c1 = new TCanvas("c1","preview",800,400) ;
//  c1->SetLogz();
//  hh_data->Draw("lego"); 
//	cout<<"generated ok!"<<endl;
	c1->Divide(2,1) ;
  	w.Start();	
        data->plotOn(fr1) ;
//	exg.plotOn(fr1); // takes too long
	mydecay.plotOn(fr1); // almost the same as 'exg'
	data->plotOn(fr2) ;
//	exg.plotOn(fr2); // takes too long
	mydecay.plotOn(fr2); // almost the same as 'exg'
	w.Stop();
	cout<<"plot time= "<<w.RealTime()<<endl;
	c1->cd(1) ; /* c2->cd(1)->SetLogy();*/ fr1->Draw() ;
	c1->cd(2) ; /* c2->cd(2)->SetLogy();*/ fr2->Draw() ;
	return;
}

double gaus(double E, double E0,double sig)
{
	return 1./sqrt(2.*M_PI)/sig*exp(-1.*(E-E0)*(E-E0)/(2.*sig*sig));
}


double exp_fun(double x, double tau, double t0)
{
	if(x>t0){return exp(-(x-t0)/tau);}
	else{return 0.;}
}