#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooAddPdf.h"
#include "TCanvas.h"
#include "RooPlot.h"

#ifndef RooStats_ConfInterval
#include "RooStats/ConfInterval.h"
#endif

using namespace RooFit ;
using namespace RooStats;

int gammatest(){
    const double Emin=0.5;
    const double Emax=20;
    const double nbins=100;
    const int nphot=10000;
    const int nrep=1000;
    double lEmin=log10(Emin);
    double lEmax=log10(Emax);

    const double gamma0=-3;
    
    RooRealVar a0("a0","a0",0,nphot*10);
    RooRealVar lE("lE","log_{10}(E/TeV)",lEmin,lEmax);
    RooRealVar gamma("gamma","#Gamma",-10,10);

    RooBinning linbins(nbins,lEmin,lEmax);
    lE.setBinning(linbins);

    RooGenericPdf bkg("bkg","pow( pow(10,@0) ,@1) *(pow(10,@0))",RooArgList(lE,gamma) );
    RooExtendPdf modelEx("modelEx","modelEx",bkg,a0);
    
    TH1D gammahist("gammahist","#Gamma reco- unbinned in blue, binned in red",100,-3.2,-2.8);
    TH1D gammaBhist("gammaBhist","#Gamma reco, binned",                       100,-3.2,-2.8);
    
    double xmin = linbins.binLow( linbins.binNumber(lEmin+0.1) );
    double xmax = linbins.binHigh( linbins.binNumber(lEmax-0.1) );

    for (int i=0;i<nrep;i++){ //write the 
        gamma=gamma0;
        RooDataSet *dat = modelEx.generate(lE,nphot);
        modelEx.fitTo(*dat,Range(xmin,xmax));
        gammahist.Fill(gamma.getValV());

        gamma=gamma0;
        RooDataHist *datH = modelEx.generateBinned(lE,nphot);
        modelEx.fitTo(*datH,Range(xmin,xmax));
        gammaBhist.Fill(gamma.getValV());

        delete dat;
        delete datH;
    }

    TCanvas cgamma("cgamma","cgamma",600,600);
    gammahist.Draw();
    gammaBhist.SetLineColor(kRed);
    gammaBhist.Draw("SAME");
    cgamma.Print("cgammaPlot.pdf");



    return 0;
}