#include "TH2D.h"
#include "TF2.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TRandom3.h"
#include "TVirtualFitter.h"
#include "TList.h"

#include <iostream>

double gFun(double *x, double *par) {
   double z1 = double((x[0]-par[1])/par[2]);
   return par[0]*exp(-0.5*(z1*z1));
}

// data need to be globals to be visible by fcn
//TRandom3 rndm;
TH1D *h1;
Int_t npfits;

void myFcn(Int_t & /*nPar*/, Double_t * /*grad*/ , Double_t &fval, Double_t *p, Int_t /*iflag */  )
{
   TAxis *xaxis1  = h1->GetXaxis();
   TAxis *yaxis1  = h1->GetYaxis();

   int nbinX1 = h1->GetNbinsX();
   //int nbinY1 = h1->GetNbinsY();

   double chi2 = 0;
   double x[1];
   double tmp;
   npfits = 0;
   for (int ix = 1; ix <= nbinX1; ++ix) {
     x[0] = xaxis1->GetBinCenter(ix);
     if ( h1->GetBinError(ix) > 0 ) {
         tmp = (h1->GetBinContent(ix) - gFun(x,p))/h1->GetBinError(ix);
         chi2 += tmp*tmp;
         npfits++;
     }
   }


   fval = chi2;
}

int fit2dHist(int option=1) {

  h1 = new TH1D("h1","",20,-5.0,5.0);

   double iniParams[3] = { 100, 6.0, 2.0};
   // create fit function                                                       
   TF1 * func = new TF1("func",gFun,-30.0,30.0,3);
   func->SetParameters(iniParams);

   // fill Histos
   UInt_t mySeed = 3;
   TRandom3 *myRandGen = new TRandom3(mySeed);
   Double_t myRandVal;
   int nToPlot = 10000;
   for (int myInt = 0; myInt < nToPlot; myInt++){
     double xFlat = myRandGen->Uniform(-3.0,3.0);
     double xGaus = myRandGen->Gaus(0.0,1.0);
     //hFlat->Fill(xFlat);
     h1->Fill(xGaus);
   }



   // fill data structure for fit (coordinates + values + errors)
   //std::cout << "Do global fit" << std::endl;
   // fit now all the function together

   //The default minimizer is Minuit, you can also try Minuit2

   if (option%10 == 2) TVirtualFitter::SetDefaultFitter("Minuit2");
   else                TVirtualFitter::SetDefaultFitter("Minuit");
   TVirtualFitter * minuit = TVirtualFitter::Fitter(0,3);

   for (int i = 0; i < 3; ++i) {
     minuit->SetParameter(i, func->GetParName(i), func->GetParameter(i), 0.01, 0,0);
   }
   minuit->SetFCN(myFcn);

   double arglist[100];
   arglist[0] = 0;
   //arglist[1] = 1.1;
   // set print level
   minuit->ExecuteCommand("SET PRINT",arglist,2);

   // minimize (DO FIT)
   arglist[0] = 5000; // number of function calls
   arglist[1] = 0.01; // tolerance
   minuit->ExecuteCommand("MIGRAD",arglist,2);


   //get result
   double minParams[3];
   double parErrors[3];
   for (int i = 0; i < 3; ++i) {
     minParams[i] = minuit->GetParameter(i);
     parErrors[i] = minuit->GetParError(i);
   }
   double chi2, edm, errdef;
   int nvpar, nparx;
   //minuit->GetStats(chi2,edm,errdef,nvpar,nparx);


   func->SetParameters(minParams);
   func->SetParErrors(parErrors);
   func->SetChisquare(chi2);
   int ndf = npfits-nvpar;
   func->SetNDF(ndf);


   // add to list of functions
   h1->GetListOfFunctions()->Add(func);
   
   
   // Create a new canvas.
   TCanvas * c1 = new TCanvas("c1","Two HIstogram Fit example",100,10,900,800);
   h1->Draw();
   
   return 0;
}