#include "TMinuit.h"
#include "TRandom.h"
#include "TH1F.h"
#include "TMath.h"
#include <iostream>
#include <stdlib.h>
#include "TGraphErrors.h"

const double TWOPI=TMath::TwoPi();
const double RANDM=RAND_MAX;

// Sample definition
const int Nmax=1000;
float M[Nmax];    //array of gaussian random deviates

// Range of the sample
float Mmin = 20.;
float Mmax = 70.;

// Summary in a Histogram
TH1F *hd = new TH1F("hd","The Data",100,20.,70.);

// Likelihood distribution
TH1F *ld = new TH1F("ld","Likelihood Distribution",100,0.,10.);


//---------------------------------------------------
double gauss(double x, double mean, double sigma, double integral)
//---------------------------------------------------
{
  return (integral/sqrt(TWOPI)/sigma)*exp(-(x-mean)*(x-mean)/2/sigma/sigma);
}

//---------------------------------------------------
double Likelihood(double mu, double s)
//---------------------------------------------------
{
  double proba=0.;
  for (int i=0;i!=Nmax;i++)
    {      
      proba = proba-2.*log(gauss(M[i],mu,s,1.));
      // proba = proba*gauss(M[i],mu,s);
      //      cout << M[i] << "  " <<  gauss(M[i],mu,s) << endl;
    }
  //  cout << proba << endl;
  
  return proba;
}


//---------------------------------------------------
void fitFCN(Int_t &npar, Double_t *gin, Double_t &fcn, Double_t *x, Int_t iflag)
//---------------------------------------------------
{
  fcn=Likelihood(x[0],x[1]);
}


//----------------------------------------------------------------------
void Fit()
//----------------------------------------------------------------------
{
  float mu = 45.;
  float s = 6.2;
  for (int i=0;i!=Nmax;i++)
    {  
      // The trial and rejection method : for sampling a gausian random deviate
      float rndG1 = 1.;
      float rndG2 = 1.;
      double m,fm = 0;
      while (fm<rndG1) {
	rndG1=rand()/RANDM;
	rndG2=rand()/RANDM;
	m = rndG2*(Mmax-Mmin)+Mmin;      
	fm = gauss(m,mu,s,1.);
      }
      M[i]=m;
      hd->Fill(m);
      // When understood, can use directly this 
      // (where an optimised trial and rejection method is coded for you) 
      // M[i] = gRandom->Gaus(mu,s);
    }
  
  const int nPlotMax=100;
  float p[nPlotMax];
  float L[nPlotMax];
  float ep[nPlotMax];
  float eL[nPlotMax];
  double Lmin=9999999.;
  for (int iPlot=0;iPlot!=nPlotMax;iPlot++)
    {
      p[iPlot]=44.8+iPlot*0.4/nPlotMax;
      L[iPlot]=Likelihood(p[iPlot],s);
      ep[iPlot]=0.0001;
      eL[iPlot]=0.0001;
      if (L[iPlot] < Lmin) Lmin=L[iPlot];
      cout << iPlot << "  " << p[iPlot] << "  " << L[iPlot] << endl;
    }
    for (int iPlot=0;iPlot!=nPlotMax;iPlot++)
    {
      L[iPlot]=L[iPlot]-Lmin;
    }

  int ierflg = 0;
  double arglist[10], edm, errdef, dum1,q1, dum2,q2;
  
  TMinuit *gMinuit = new TMinuit(2);
  gMinuit->Command("SET PRI 0");
  gMinuit->Command("SET NOW");
  gMinuit->Command("SET ERR 1");
  gMinuit->SetFCN(fitFCN);
  // Set starting values and step sizes for parameters
  gMinuit->mnparm(0, "Mean", 45., 0.01, 30., 60., ierflg);
  gMinuit->mnparm(1, "Sigma", 1., 0.01, 0., 10., ierflg);
  // Now ready for minimization step
  gMinuit->Command("MIGRAD 500 1");
  gMinuit->GetParameter(0, q1, dum1);  
  gMinuit->GetParameter(1, q2, dum2);
  cout << "Mass " << q1 << " +/- " << dum1 << endl; 
  cout <<  "Width " << q2 << " +/- " << dum2 << endl;

  // Plot the result 
  gROOT->SetStyle("Plain");
  gStyle->SetCanvasColor(0);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetPadColor(0);
  gStyle->SetLineWidth(1.);
  gStyle->SetCanvasDefH(596);
  gStyle->SetCanvasDefW(596);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetFrameBorderMode(0);
  gStyle->SetTitleBorderSize(0);
  gStyle->SetStatBorderSize(0);
  gStyle->SetPadTickX(1);
  gStyle->SetPadTickY(1);
  gStyle->SetPalette(1);
  gStyle->SetOptStat(0);
  
  c0 = new TCanvas("c0", "", 0, 0, 600, 450);
  c0->SetLeftMargin(0.12);
  c0->SetRightMargin(0.12);
  c0->SetBottomMargin(0.2);
  
  // Here perform a simple binned fit
  hd->SetMarkerStyle(20);
  hd->SetMarkerColor(2);
  hd->SetMarkerSize(0.8);
  hd->Fit("gaus");  
  hd->Draw("epsame");

  c1 = new TCanvas("c1", "", 50, 50, 650, 500);
  c1->SetLeftMargin(0.12);
  c1->SetRightMargin(0.12);
  c1->SetBottomMargin(0.2);
  
  TMultiGraph *mg = new TMultiGraph();
  TGraphErrors *LikelihoodGraph = new TGraphErrors(nPlotMax,p,L,ep,eL);
  LikelihoodGraph->SetMarkerStyle(20);
  LikelihoodGraph->SetMarkerColor(2);
  LikelihoodGraph->SetMarkerSize(1.);
  mg->Add(LikelihoodGraph);
  mg->SetMinimum(0.);
  mg->SetMaximum(100.);
  mg->Draw("AP");
cout<<"finished";
}