// root.exe TestPLLFit.C+

#include "TF1.h"
#include "TH1.h"
#include "TList.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TMath.h"
#include "Fit/Chi2FCN.h"
#include "Fit/PoissonLikelihoodFCN.h"
#include "HFitInterface.h"
#include "Fit/BinData.h"
#include "Fit/Fitter.h"
#include "Math/WrappedMultiTF1.h"
#include "TMatrixDSym.h"

#include <iostream>
using std::cout;

/******************************************************************************/
Double_t Gaus(Double_t * x, Double_t * par)
{
 return par[0]*TMath::Gaus(x[0], par[1], par[2], kTRUE);
}// end of function
/******************************************************************************/
void TestPLLFit(void) 
{ 

 // ====================================================
 // Creating toy fit input
 //
 const Double_t xMin = -100.0, xMax = +100.0;
 TH1D * h1 = new TH1D("h1", "", 200, xMin, xMax);

 const Int_t nFitFunPar = 3;
 TF1 * f1 = new TF1("f1", Gaus, xMin, xMax, nFitFunPar);
 f1->SetParameters(1.0, 2.0, 20.0);

 const Int_t nEvents = 100;
 //const Int_t nEvents = 100*100;
 h1->FillRandom("f1", nEvents);



 // ====================================================
 // Performing the fit using ROOT::Fit::PoissonLLFunction
 //
 ROOT::Fit::Fitter fitter;
 const int nPar = f1->GetNpar();
 fitter.Config().SetParamsSettings(nPar, f1->GetParameters());
 fitter.Config().ParSettings(0).SetLimits(0.0, 20.0*nEvents);
 fitter.Config().ParSettings(1).SetLimits(xMin, xMax);
 fitter.Config().ParSettings(2).SetLimits(xMin, xMax);
 fitter.Config().MinimizerOptions().SetPrintLevel(2);
 fitter.Config().SetMinimizer("Minuit2", "Migrad");

 TF1 * fitFun = new TF1("fitFun", Gaus, xMin, xMax, nFitFunPar);
 ROOT::Math::WrappedMultiTF1 * wf = new ROOT::Math::WrappedMultiTF1(*fitFun, 1);
 ROOT::Fit::DataOptions opt;
 opt.fUseEmpty = true;
 //opt.fBinVolume = true;
 //opt.fIntegral = true;
 //opt.fExpErrors = true; // probably, only in new versions of ROOT...
 ROOT::Fit::DataRange range;
 range.SetRange(xMin, xMax);
 ROOT::Fit::BinData * data = new ROOT::Fit::BinData(opt, range);
 ROOT::Fit::FillData(*data, h1);

 const bool chi2Fit = false;
 ROOT::Math::FitMethodFunction * localFitFun = new ROOT::Fit::PoissonLLFunction(*data, *wf);
 //ROOT::Math::FitMethodFunction * localFitFun = new ROOT::Fit::PoissonLLFunction(*data, *wf, 0, true);
 //ROOT::Math::FitMethodFunction * localFitFun = new ROOT::Fit::PoissonLLFunction(*data, *wf, 0, false);
 //ROOT::Math::FitMethodFunction * localFitFun = new ROOT::Fit::PoissonLLFunction(*data, *wf, 1, true);
 //ROOT::Math::FitMethodFunction * localFitFun = new ROOT::Fit::PoissonLLFunction(*data, *wf, 1, false);

 fitter.FitFCN(fitter.Config().NPar(), *localFitFun, 0, data->Size(), chi2Fit);
 const ROOT::Fit::FitResult result = fitter.Result();
 result.Print(cout);



 // ====================================================
 // Creating a canvas to draw the fit results
 //

 gStyle->SetOptFit(1111);
 gStyle->SetStatW(0.15);
 gStyle->SetStatH(0.1);
 gStyle->SetStatY(0.93);

 TCanvas * c1 = new TCanvas("Simfit", "Fit of a histogram with low statistics", 10,10,700,700);
 c1->Divide(1,2);



 // ====================================================
 // Drawing the fit using ROOT::Fit::PoissonLLFunction
 //

 c1->cd(1);
 f1->SetParameters(result.GetParams());
 f1->SetChisquare( result.Chi2() );
 f1->SetNDF( result.Ndf() );
 TMatrixDSym covM(nPar); 
 result.GetCovarianceMatrix(covM);
 covM.Print();
 for(Int_t n=0; n<nPar; n++) f1->SetParError(n, sqrt(covM(n,n)));
 f1->SetLineColor(kBlue);
 h1->GetListOfFunctions()->Add(f1);
 h1->Draw();

 const Double_t fitFunIntegral1 = f1->Integral(xMin, xMax);
 cout<<"fitFunIntegral1 = "<<fitFunIntegral1<<" (ratio="<<fitFunIntegral1/nEvents<<")\n";



 // ====================================================
 // Performing and drawing the fit using TH1::Fit
 //

 TH1D * h2 = new TH1D(*h1);
 TF1 * f2 = new TF1(*f1);

 c1->cd(2);
 h2->Fit(f2, "L");
 const Double_t fitFunIntegral2 = f2->Integral(xMin, xMax);
 cout<<"fitFunIntegral2 = "<<fitFunIntegral2<<" (ratio="<<fitFunIntegral2/nEvents<<")\n";

}
/******************************************************************************/