/*---------------------------------------------
Example of fit with MINUIT
 
 Work only in compiled mode:
  .L minuitFit.C+;
  minuitFit();
-------------------------------------------- */
#include "TH1.h"
#include "TF1.h"
#include "TMinuit.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TLegend.h"


//------------------------------------------------------------------------------
// Global variables

const Double_t trueTau = 100.;
const Int_t numberOfPoints=500; // statistics size
const Bool_t fitWithTruncature = true;
const Double_t measLimits[2] = {200., 600.}; // measurement range
const Double_t fitLimits[2] = {0., 1010.};
const Int_t nBins = 100;
const Int_t nExperiences = 1000; // for toy MC

//------------------------------------------------------------------------------
//Define here base function for fit and drawing
Double_t pdf( Double_t x, Double_t tau ) {

  if( fitWithTruncature ) {
    return exp(-x/tau)/tau/(exp(-measLimits[0]/tau)-exp(-measLimits[1]/tau));
  }
  else {
    return exp(-x/tau)/tau;    
  }
}

//------------------------------------------------------------------------------
// Wrapper to pdf using 
// First parameter is used as normalisation
Double_t pdfWrap( Double_t *x, Double_t *par ) {
  
  return par[0] * pdf( x[0], par[1]);

}

//------------------------------------------------------------------------------
// Function used by MINUIT
void fcn(Int_t &numberOfParameters, Double_t *inputGradient, Double_t &value, Double_t *parameters, Int_t flag){

  static Double_t *x = new Double_t[numberOfPoints]; // declared as static to keep the same value in the different call to the function
  static TH1F *histoToFit = new TH1F("histoToFit", "Original data", nBins, fitLimits[0], fitLimits[1]);
  histoToFit->SetMarkerStyle(20);
  histoToFit->SetMarkerSize(0.5);
  
  TF1 *fToGetRandom = new TF1("getRandom", pdfWrap, measLimits[0], measLimits[1], 2);
  fToGetRandom->SetNpx(1000);

  // init
  // Here fill the array of measurements randomly from a function
  if( flag==1 ) {
    //printf( "FCN: (flag=%1d) Initializing data to fit\n\n", flag );
    fToGetRandom->SetParameters( 1., parameters[0]);
    histoToFit->Reset();
    for( Int_t i=0; i<numberOfPoints; i++ ) {
      x[i] = fToGetRandom->GetRandom(  );
      histoToFit->Fill( x[i] );
    }
  }

  // initialize gradient
  if( flag==2 ) {
    printf( "FCN: (flag=%1d) nop (to compute first input gradient)\n", flag );
  }

    
  // Compute here the -log( Likelyhood ) = value
  // loop on all points
  Double_t prob;
  value = 0.0;
  for( Int_t i=0; i<numberOfPoints; i++ ) {
    prob = pdf( x[i], parameters[0] );
    value = value - log( prob );
    //printf( "point %d, x=%.1f, param (%.1f, %.1e), prob=%e, value=%e\n", i, x[i], parameters[0], parameters[1], prob, value);
  } // end loop on all points
 

  // end
  if( flag==3 ) {
    printf( "FCN: (flag=%1d) ending fit\n", flag );
    TF1 *hFunction = new TF1("fitFunction", pdfWrap, fitLimits[0], fitLimits[1], 2);
    hFunction->SetParameters( numberOfPoints*(fitLimits[1]-fitLimits[0])/nBins, parameters[0] );
    histoToFit->Draw("PE");
    hFunction->Draw("same");
  }

}


//------------------------------------------------------------------------------
// Actual macro to make the fit
void minuitFit(){

  Double_t fitParameters[3], paramValues[3], paramErrors[3];
  Int_t error = 0;
  Double_t argument;
  TGraph *gcontour, *gscan, *g;

  // MINUIT settings
  TMinuit *minimization = new TMinuit(3);
  minimization->SetFCN( fcn );
  minimization->DefineParameter( 0, "#tau", trueTau, .5, 0., 20*trueTau);
  
  // Generation step
  argument = 1;
  minimization->mnexcm("CALL FCN", &argument, 1, error);

  // Just scan the function to minimize
  //  range is +/- 5 x expected uncertainty
  printf( "\n************\n Scanning\n ****************\n\n" );
  TCanvas * clike = new TCanvas("cmin","Likelihood", 500, 500);
  clike->SetGrid(1,1);
  Double_t arg[4] = {1, 20, trueTau/10, trueTau*10};
  minimization->mnexcm("SCAN", arg, 4, error);
  gscan = (TGraph*)gMinuit->GetPlot();
  gscan->SetTitle("Basic scan");
  TF1 *flikelihood;
  Double_t expectedUncertainty;
  if( fitWithTruncature ) {
    expectedUncertainty = trueTau / sqrt( numberOfPoints) / sqrt(1- pow( (measLimits[1]-measLimits[0])/trueTau/2 / sinh((measLimits[1]-measLimits[0])/trueTau/2) ,2) );
   flikelihood= new TF1("flikelihood","[0]*(log(x)+[1]/x+log(exp(-[2]/x)-exp(-[3]/x))-log([1])-1-log(exp(-[2]/[1])-exp(-[3]/[1])))",trueTau-5*expectedUncertainty, trueTau+5*expectedUncertainty);
    flikelihood->SetParameters(numberOfPoints,trueTau,fitLimits[0],fitLimits[1]);
  }
  else {
    expectedUncertainty = trueTau / sqrt( numberOfPoints);
    flikelihood= new TF1("flikelihood","[0]*(log(x)+[1]/x-log([1])-1)",trueTau-5*expectedUncertainty, trueTau+5*expectedUncertainty);    
    flikelihood->SetParameters(numberOfPoints,trueTau);
  }
  printf("\n Expected uncertainty = %.3f\n\n", expectedUncertainty);
  flikelihood->SetTitle("Analytical likelihood");
  flikelihood->Draw();
  TF1 *fparabolic = new TF1("fparabolic","[0]*(x/[1]-1)*(x/[1]-1)/2",trueTau-5*expectedUncertainty, trueTau+5*expectedUncertainty);
  fparabolic->SetTitle("Parabolic likelihood");
  fparabolic->SetParameters(numberOfPoints,trueTau);
  fparabolic->SetLineStyle(2);
  fparabolic->Draw("same");
  gscan->Draw("PL");
  TLegend * leg = (TLegend*)clike->BuildLegend();
  leg->SetFillColor(0);
  

  // one Minimization for display
  printf( "\n************\n Starting minimization\n ****************\n\n" );
  minimization->DefineParameter( 0, "#tau", 4*trueTau, 1., 0., 20*trueTau);
  argument = 0;
  minimization->mnexcm("MIGRAD", &argument, 0, error);

  minimization->GetParameter( 0, paramValues[0], paramErrors[0] );
  printf( "\n FIT RESULT : tau=%5.2f+/-%5.3f\n\n", paramValues[0], paramErrors[0] );

  // A second step could be used to improve errors
  argument = 0;
  minimization->mnexcm("MINOS", &argument, 0, error);
  minimization->GetParameter( 0, paramValues[0], paramErrors[0] );
  printf( "\n AFTER MINOS : tau=%5.2f+/-%5.3f\n\n", paramValues[0], paramErrors[0] );

  // Does not work...
//  Double_t xpt[100], ypt[100];
//  minimization->mnplot( xpt, ypt, "x", 100, 10, 1);
//  cmin->cd(2);
//  g = (TGraph*)gMinuit->GetPlot();
//  g->Draw("AP");
  
  // Useless with only one parameter
//  minimization->SetErrorDef(2);
//  gcontour = (TGraph*)minimization->Contour(50,0,1);
//  gcontour->SetTitle("Contour with ErrDef=2");
//  cmin->cd(3);
//  gcontour->Draw("ALP");
  

  TCanvas * cfit = new TCanvas("cfit","Fit", 500, 500);
  cfit->SetLogy();
  argument = 3;
  minimization->mnexcm("CALL FCN", &argument, 1, error);
  
  
  // Simulate many experiments = toy MC
  
  TH1F *hfitresult = new TH1F( "hfitresult", "Distribution of fit results", 50, trueTau-20*expectedUncertainty, trueTau+5*expectedUncertainty);
  hfitresult->SetFillColor(29);
  TH1F *herrresult = new TH1F( "herrresult", "Distribution of uncertainty results", 50, 0, 5*expectedUncertainty);
  herrresult->SetFillColor(29);
  
  argument = -1; // make minuit silent
  minimization->mnexcm("SET PRINTOUT ", &argument, 1, error);
  for ( Int_t iExp=0; iExp<nExperiences; iExp++) {
    minimization->DefineParameter( 0, "#tau", trueTau, .5, 0., 20*trueTau);
    argument = 1;
    minimization->mnexcm("CALL FCN", &argument, 1, error);
    minimization->DefineParameter( 0, "#tau", 4*trueTau, 1., 0., 20*trueTau);
    argument = 0;
    minimization->mnexcm("MIGRAD", &argument, 0, error);
    minimization->GetParameter( 0, paramValues[0], paramErrors[0] );
    hfitresult->Fill( paramValues[0]);
    herrresult->Fill( paramErrors[0]);
  }
  
  TCanvas * cmc = new TCanvas("cmc","MC distribution", 500, 500);
  cmc->Divide(2,1);
  cmc->cd(1);
  hfitresult->Draw();
  hfitresult->Fit("gaus");
  printf( "\n\n Estimator average = %.2f, RMS = %.2f\n", hfitresult->GetFunction("gaus")->GetParameter(1), hfitresult->GetFunction("gaus")->GetParameter(2));
  cmc->cd(2);
  herrresult->Draw();
  herrresult->Fit("gaus");
  printf( "\n\n Uncertainty average = %.2f, RMS = %.2f\n", herrresult->GetFunction("gaus")->GetParameter(1), herrresult->GetFunction("gaus")->GetParameter(2));


  delete minimization;

}