#include "TROOT.h"
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "TGraph.h"
#include "TStyle.h"
#include "TMath.h"
#include "RooAbsPdf.h"
#include "RooGaussian.h"
#include "RooFFTConvPdf.h"
#include "RooNumConvPdf.h"
#include "RooFormulaVar.h"
#include "RooPlot.h"

using namespace RooFit ;

void prodConv()
{
  Double_t E_res = 0.075;

  RooRealVar Enu("Enu","Enu",1.9,1.0,10.0);

  //Spectrum PDF
  RooGenericPdf spect("spect","spect","exp((-@0/1.54)-pow(@0/6.05,2)-pow(@0/7.73,10))",RooArgList(Enu));

  //Box function
  RooRealVar box_mid("box_mid","box_mid",5.0) ; //location of box middle
  RooRealVar box_low("box_low","box_low",4.0) ; //location of box low edge
  RooGenericPdf box("box","box","abs(@0-@1)<(@1-@2)",RooArgList(Enu,box_mid,box_low)) ;

  //Product of GenericPdf with UnoscPdf
  RooProdPdf prod("prod","prod",RooArgList(box,spect));
  
  //Gaussian function with Enu-dependent width
  RooRealVar res_mean("res_mean","res_mean",0.0);
  RooFormulaVar res_sigma("res_sigma","res_sigma",Form("%f*pow(abs(@0)+0.0001*(@0==0),-0.5)",E_res),RooArgList(Enu));
  RooGaussian res_gaus("res_gaus","Gaussian resolution DRF",Enu,res_mean,res_sigma);
  
  //Convolve box*spectrum with a gaussian
  Enu.setBins(10000,"cache"); 
  RooNumConvPdf conv("conv","conv",Enu,res_gaus,prod);//works
  //RooFFTConvPdf conv("conv","conv",Enu,res_gaus,prod);//doesn't work

  //Plot convolution
  RooPlot *Enuframe = Enu.frame();
  conv.plotOn(Enuframe);
  Enuframe->Draw();

}