//----------- The code starts here -------------
#include "TH1.h"
#include "TF1.h"
#include <iostream>
#include <stdlib.h>
#include "TSystem.h"
#include "TMath.h"
#include "TRandom.h"
#include "TMinuit.h"
using namespace std;

Bool_t fit = true;
Int_t ncall = 0;


void fastfit(TH1F *h, Double_t *pars)
{
        if (fit)
   {
      h->Fit("gaus");
      //h->Fit("gaus","0L","",0.,50.);
      TF1 *f=h->GetFunction("gaus");
      if (f == 0) cout << "Error: no function gaus exist\n";
      else for (Int_t i = 0; i < 3; i++) pars[i] = f->GetParameter(i);
   }
   else
   {
      pars[0] = h->GetMaximum();
      pars[1] = h->GetMean();
      pars[2] = h->GetRMS();
   }
   return;
}

void getpar(Double_t *x, Float_t &mass, Float_t &width)
{
   Double_t par[6];
//
   TH1F h("peak","peak",50,0.,50.);
   TRandom myrandom(654321);
   for (Int_t i=0; i<10000;i++) // Fill the histogram with a Gaussian 
   {
     Float_t M = myrandom.Gaus(25.*x[0],5.0+10.*(x[1]-1.1)*(x[1]-1.1));
     h.Fill(M);
   }
//
   fastfit(&h,par);
   mass = par[1]; 
   width = par[2];
   return;
}

void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
   Float_t mass, width;
   getpar(par, mass, width);
   ncall++;
//
   f = 1000.*width + (mass - 22.5)*(mass - 22.5);
   if (ncall < 10) cout << "iflag, mass, width, f = " << iflag << ", "<< mass << ", " << width << ", " << f << endl;
//
   return;
}
void fittest(Int_t flag = 0)
{
   Double_t x[2], err[2]; //, tmp;
   //char name[10];
//
   ncall = 0;
   fit=true;
   if (flag != 0) fit = false;
//
   TMinuit *minuit = new TMinuit();
   minuit->DefineParameter(0, "Mass Scale",  1., 0.05, 0.8, 1.2);
   minuit->DefineParameter(1, "Width Scale", 1., 0.05, 0.8, 1.2);
// 
   minuit->SetFCN(fcn);
   Double_t arglist[10];
   Int_t ierflg = 0;
//
   arglist[0] = 1;
   minuit->mnexcm("SET ERR", arglist ,1, ierflg);
//
   arglist[0] = 5000;
   arglist[1] = 1.;
   if (ierflg == 0) minuit->mnexcm("MINIMIZE", arglist ,1, ierflg);
//
   if (ierflg != 0) cout << "Error in mnexcm: " << ierflg << endl;
//
   for (Int_t i = 0; i < 2; i++)
   {
      minuit->GetParameter(i, x[i], err[i]);
      cout << "Parameter " << i << " = " << x[i] << " +- " << err[i] <<endl;
   }
}