#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooLandau.h"
#include "RooFFTConvPdf.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TH1.h"
using namespace RooFit;

TH1 *makeTH1();

void gexp_conv()
{
   // ---------------------------------------
   // ---------------------------------------------------
   // I m p o r t i n g   R O O T   h i s t o g r a m s
   // ===================================================

   // I m p o r t   T H 1   i n t o   a   R o o D a t a H i s t
   // ---------------------------------------------------------

   // Create a ROOT TH1 histogram
   TH1 *hh = makeTH1();

   // Declare observable x
   //RooRealVar x("x", "x", -10, 10);
   RooRealVar x("x", "x", 0, 800);

   // Create a binned dataset that imports contents of TH1 and associates its contents to observable 'x'
   RooDataHist dh("dh", "dh", x, Import(*hh));

      // P l o t   a n d   f i t   a   R o o D a t a H i s t
   // ---------------------------------------------------

   RooAbsData::PlotOpt opts;
   opts.correctForBinWidth = true;
   opts.etype = RooAbsData::None;

   // Make plot of binned dataset showing Poisson error bars (RooFit default)
   RooPlot *frame = x.frame(Title("Imported TH1 for fitting"));
   //dh.plotOn(frame, DrawOption("hist"), DataError(RooAbsData::None), XErrorSize(1));
   dh.plotOn(frame, DataError(RooAbsData::None), DrawOption("HIST"));
   //dh.plotOn(frame,opts);

   // Fit a Gaussian pdf to the data
   RooRealVar mean("mean", "mean", 400., 300., 700.);
   RooRealVar sigma("sigma", "sigma", 50., 15., 100.);
   RooGaussian gauss("gauss", "gauss", x, mean, sigma);

   RooRealVar sl("sl", "slope", -0.01, -1, 0.001);
   RooRealVar sl1("sl1", "slope1", -0.03, -0.0001, 10.);
   RooExponential exp("exp", "exp", x, sl);
   RooExponential exp1("exp1", "exp1", x, sl1);

   RooGenericPdf step_func("step_func","step_func","@0 > 0.0 ? @1 : 0.0",RooArgSet(x, exp));

   RooFFTConvPdf gxe("gxe", "gauss (X) exponential", x, gauss, step_func);

   //gxe.fitTo(dh,Range("signal"));
   //gxe.plotOn(frame);

   x.setRange("signal",315,640);
   std::unique_ptr<RooFitResult> fitResult{gxe.fitTo(dh,Range("signal"),NormRange("signal"), PrintLevel(-1), Save(true))};
   fitResult->Print();

   gxe.plotOn(frame, Range("signal"), NormRange("signal"));

   //x.setRange("signal",360,465);
   //gauss.fitTo(dh,Range("signal"));
   gauss.plotOn(frame, LineColor(kRed));

   //x.setRange("signal",465,660);
   //exp.fitTo(dh,Range("signal"),NormRange("signal"));
   step_func.plotOn(frame, LineColor(kRed));

   // Draw all frames on a canvas
   TCanvas *c = new TCanvas("GaussExpConv", "gauss (x) exp conv", 800, 600);
   frame->GetYaxis()->SetTitleOffset(1.4);
   frame->Draw();

   c->SaveAs("gexp_conv.png");
}

// Create ROOT TH1 filled with a Gaussian distribution
TH1 *makeTH1()
{
   TH1D *hh = new TH1D("hh", "hh", 800, 0, 800);
   ifstream in("test.txt");         //207At
   int n,nx;

   while(!in.eof()){
           in >> n >> nx;
           hh->SetBinContent(n,nx);
   }
   return hh;
}