using namespace RooFit;

#include "RooRealVar.h"
#include "RooDataSet.h"
#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"



void test_conv_bug()
{
  // define x1 shape
  
  RooRealVar x1("x1", "mass of x1", 400, 1000);
  RooRealVar x1_mean("x1_mean", "mean mass", 780, 760, 850);
  RooRealVar x1_width("x1_width", "width of x1 resonance", 150);

  RooBreitWigner x1_bw("x1_bw", "nonrel bw", x1, x1_mean, x1_width);

  // define gaussian describing resolution of x1

  RooRealVar mg("mg","mg",0) ; 
  RooRealVar sg("sg","sg",35, 10, 60) ; 
  RooGaussian resolution("resolution","gauss",x1,mg,sg) ; 
  


  //-----------------------------------------------------------
  //-----------------------------------------------------------
  // Construct convolution

  // 1) use RooVoigtian
  RooVoigtian x1_bw_conv("x1_bw_conv", "Voigtian", x1,  x1_mean, x1_width, sg);

  // 2) use RooFFTConvPdf

  //RooFFTConvPdf x1_bw_conv("x1_bw_conv", " BW (X) gauss", x1, x1_bw, resolution);
  //-----------------------------------------------------------
  //-----------------------------------------------------------



  // define parameters x2
  RooRealVar x2("x2", "mass of x2", 5000, 5600);
  RooRealVar x2_mean("x2_mean", "mean of Bmass from mc",5279);
  RooRealVar x2_sigma("x2_sigma", "sigma of gaussian 1",13, 5, 20);

  RooGaussian gauss("gauss","gauss",x2,x2_mean,x2_sigma);

  // create 2 dimensional pdf for x1 and x2
  RooProdPdf *sig = new RooProdPdf("sig", "signal from Bplus and x1plus", RooArgList(x1_bw_conv, gauss));
  

  // generate dataset from pdfs

  RooDataSet *data = sig->generate(RooArgSet(x2, x1), 8000);
  
  // fit data
  RooFitResult * result = sig->fitTo(*data);
  
  /////////////////////////////////////////////////////
  // plotting
  // plot result for x1
  TCanvas* canvas1 = new TCanvas("canvas1","canvas1");

  RooPlot * frame1 = x1.frame(Title("X1"));//, Range("full")
  data->plotOn(frame1, Binning(60));
  sig->plotOn(frame1, LineColor(kBlack));
  sig->paramOn(frame1, data);
   
  frame1->Draw();

  // plot result for x2
  TCanvas* canvas2 = new TCanvas("canvas2","canvas2");


  RooPlot * frame22 =  x2.frame(Title("X2"));//, Range("full")
  data->plotOn(frame22, Binning(60));
  sig->plotOn(frame22, LineColor(kBlack));
  sig->paramOn(frame22, data);
   
  frame22->Draw();
   

}