#include <iostream>
#include <TH1D.h>
#include <TF1.h>
#include <TCanvas.h>
#include <TVirtualFitter.h>
#include <TList.h>
#include <TStyle.h>

TH1D *hist;
TF1 *fitFunc;
Int_t nFitPoints;

//The Fit Function
Double_t gauss(Double_t *x, Double_t *par){
  return par[0] * exp(-0.5*pow((x[0]-par[1]) / par[2],2) );
}

//The Function to be Minimized
void chi2(Int_t & /*nPar*/, Double_t * /*grad*/ , Double_t &fval, Double_t *p, Int_t /*iflag */){

  Double_t chi2(0);
  nFitPoints = 0;
  Double_t x[1];
  
  //cout <<p[0] <<" " <<p[1] <<" " <<p[2] <<endl;
  
  for (Int_t iBin=1; iBin<hist->GetNbinsX(); iBin++){
    x[0] = hist->GetBinCenter(iBin);

    chi2 += pow((hist->GetBinContent(iBin) - gauss(x,p))/(Double_t)hist->GetBinError(iBin),2);
    nFitPoints++;

    //cout <<iBin <<" " <<x[0] <<" " <<hist->GetBinContent(iBin) <<" " <<gauss(x,p) <<" " <<hist->GetBinError(iBin) <<endl;
  }

  fval = chi2;

}

//Main
void simpleSimFit(){

  gStyle->SetOptFit(1);
  
  //Create Histogram, Function, and Fill Histogram
  hist = new TH1D("hist","hist",130,-3,10);
  fitFunc = new TF1("fit",gauss,-3,10,3);
  fitFunc->SetParameters(1,4,1);
  hist->FillRandom("fit",1000);

  //Seed Parameters
  const Int_t nPars(3);
  Double_t seeds[nPars] = {50,4,1};
  Double_t lowers[nPars] = {30,3,.5};
  Double_t uppers[nPars] = {80,5,1.5};

  //Set up the fitter with seed parameters, bounds, and set the minimization function
  TVirtualFitter::SetDefaultFitter("Minuit");
  TVirtualFitter *minuit = TVirtualFitter::Fitter(0,nPars);
  for (Int_t iPar=0; iPar<nPars; iPar++){
    minuit->SetParameter(iPar,Form("par%d",iPar),seeds[iPar],0.01,lowers[iPar],uppers[iPar]);
  }
  minuit->SetFCN(chi2);

  //Set the Output Level
  Double_t arglist[100];
  arglist[0] = 0;
  minuit->ExecuteCommand("SET PRINT",arglist,2);

  //Minimize
  arglist[0] = 5000;
  arglist[1] = 0.01;
  minuit->ExecuteCommand("MIGRAD",arglist,2);

  //Collect Results
  Double_t minPars[nPars];
  Double_t minParErrors[nPars];
  Double_t chi2(0), edm(0), errdef(0);
  Int_t nvpar(0), nparx(0);
  for (Int_t iPar=0; iPar<nPars; iPar++){
    minPars[iPar] = minuit->GetParameter(iPar);
    minParErrors[iPar] = minuit->GetParError(iPar);
    cout <<iPar <<" " <<minPars[iPar] <<" " <<minParErrors[iPar] <<endl;
  }
  minuit->GetStats(chi2,edm,errdef,nvpar,nparx);

  //Pass the Results to the Function
  fitFunc->SetParameters(minPars);
  fitFunc->SetParErrors(minParErrors);
  fitFunc->SetChisquare(chi2);
  Int_t ndf = nFitPoints - nvpar;
  fitFunc->SetNDF(ndf);

  //Associate the Function with the histogram that was fit
  hist->GetListOfFunctions()->Add(fitFunc);

  //Draw
  TCanvas *canvas = new TCanvas("canvas","canvas",20,20,400,400);
  canvas->cd();
  hist->Draw("E");

}