#include "TH1.h"
#include "TF1.h"
#include "TStopwatch.h"
#include "TMinuit.h"
#include "TVirtualFitter.h"
#include "TCanvas.h"
#include "TF2.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TH2F.h"
#include "TLegend.h"
#include "math.h"
#include "TSystem.h"
#include "TStyle.h"
#include <iostream>

void Woods(Int_t & /*nPar*/, Double_t * /*grad*/ , Double_t &fval, Double_t *x, Int_t /*iflag */  )
{
  fval = 100*(x[1] - x[0]*x[0])*(x[1] - x[0]*x[0])
  		 + (x[0] - 1)*(x[0] - 1)
  		 + 90*(x[3] - x[1]*x[1])*(x[3] - x[1]*x[1])
  		 + (1 - x[2])*(1 - x[2])
  		 + 10.1*( (x[1] - 1)*(x[1] - 1) + (x[3] -1 )*(x[3] -1 ))
  		 + 19.8*(x[1]-1)(x[3]-1);
  return;
}

void DoFit(const char *fitterType="Minuit2",double par1 = 0, double par2 = 0, double par3 = 0, double par4 = 0) {

 Double_t foundMinimum[4];
 Double_t arglist[100];
 arglist[0] = 0;

  // create fitter instance using type fitter
 TVirtualFitter::SetDefaultFitter(fitterType);
 TVirtualFitter * minuit2 = TVirtualFitter::Fitter(0,2);
 arglist[0] = 1;
 // set print level
 minuit2->ExecuteCommand("SET PRINT",arglist,2);
 // set initial prameters
 minuit2->SetParameter(0, "w", par1,   0.01, 0,0);
 minuit2->SetParameter(1, "x", par2,   0.01, 0,0);
 minuit2->SetParameter(2, "y", par3,   0.01, 0,0);
 minuit2->SetParameter(3, "z", par4,   0.01, 0,0);
 // set function to minimize
 minuit2->SetFCN(Woods);

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

 // get minimum values
 foundMinimum[0] = minuit2->GetParameter(0);
 foundMinimum[1] = minuit2->GetParameter(1);
 foundMinimum[2] = minuit2->GetParameter(2);
 foundMinimum[3] = minuit2->GetParameter(3);

 std::cout << " " << fitterType << " :" << " minimum values for " << par1 << "," << par2 << "," << par3 << "," << par4 << std::endl;
 std::cout << " par0 = " << foundMinimum[0] << std::endl;
 std::cout << " par1 = " << foundMinimum[1] << std::endl;
 std::cout << " par2 = " << foundMinimum[2] << std::endl;
 std::cout << " par3 = " << foundMinimum[3] << std::endl;

 minuit2->Clear();
}

void testWoods() {
  DoFit("Minuit2",-3.0,-1.0,-3.0,-1.0);
}

#ifndef __CINT__
int main() {

  testWoods();
}
#endif