// test fitting with a user provided fcn
#include "TF1.h"
#include "TH1.h"
#include "TVirtualFitter.h"
#include "TCanvas.h"
#include "TFitResult.h"
#include "Math/DistFunc.h"
#include "TLegend.h"


int nfitpoints = 0; 

void MyChi2Function(Int_t &, Double_t *, Double_t & chi2, Double_t * p, Int_t ) { 

   // get pointer to function and histogram 
   assert(TVirtualFitter::GetFitter() );
   TF1 * func = (TF1*) (TVirtualFitter::GetFitter() )->GetUserFunc();
   TH1 * hist = (TH1*) (TVirtualFitter::GetFitter() )->GetObjectFit();
   assert(func && hist); 

   // compute chi2 using expected errors instead of observed one 
   chi2 = 0; 
   nfitpoints = 0;
   double x[1];
   int nbins = hist->GetNbinsX(); 
   for (int i = 1; i <= nbins; ++i) { 
      x[0] = hist->GetBinCenter(i); 
      double y = hist->GetBinContent(i); 
      double fval = func->EvalPar(x,p); 
      //double ey2 = std::pow(hist->GetBinError(i),2); // observed error  
      double ey2 = fval;  
         
      if (ey2 > 0)  { 
         double res = (y-fval);
         res *= res/ey2;          
         chi2 += res; 
         nfitpoints++;
      }      
   }
   return;
}

void exampleFitUser() { 

   // generate a gaussian distributed histogram 

   TH1 * h1 = new TH1D("h1","h1",100,-5,5);
   h1->FillRandom("gaus");

   // do standard chi2 fit 
   TF1 * gaus1 = new TF1("gaus1","gaus");
   gaus1->SetLineColor(kBlue); 
   TFitResultPtr r1 = h1->Fit(gaus1,"SQ");   
   r1->Print();

   std::cout << "Fit chi2 p value = " << ROOT::Math::chisquared_cdf_c(gaus1->GetChisquare(), gaus1->GetNDF() ) << std::endl;

   // create another gaussian function
   TF1 * gaus2 = new TF1("gaus2","gaus");

   // do chi2 fit using expected errors 
   TVirtualFitter::Fitter(h1)->SetFCN(MyChi2Function); 
   gaus2->SetLineColor(kRed);
   TFitResultPtr r2 = h1->Fit(gaus2,"UQS+"); 
   // need to set chi2 value and ndf in fitted function
   gaus2->SetChisquare(r2->MinFcnValue() );
   gaus2->SetNDF(nfitpoints - r2->NFreeParameters() );
   r2->Print();
   std::cout << "NDF = " << gaus2->GetNDF() << std::endl;
   std::cout << "Fit chi2 p value = " << ROOT::Math::chisquared_cdf_c(gaus2->GetChisquare(), gaus2->GetNDF() ) << std::endl;

   if (gPad) gPad->SetLogy(1);

   TLegend * l = new TLegend(0.05,0.85,0.5,0.98);
   l->AddEntry(gaus1,"chi2 fit using observed errors","l");
   l->AddEntry(gaus2,"chi2 fit using expected errors","l");
   l->Draw();

}