//  example of a simple max likelihood estimation
//  inspired in the examples from Cowan, chap.7 and 
//  with data from Vuolo, pag.168 
//  (also after Ifit.C and fitcont.C from Root tutorials)
//
//
//                  jtmn oct 2003
//

#include "TMinuit.h"
// global variables
// data from Vuolo, pg. 168
const Int_t Npoints = 11;
Double_t x[Npoints] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1};
Double_t y[Npoints] = {3.11, 2.67, 2.65, 2.43, 2.45, 2.62, 3.2, 4.05, 5.12, 
                        6.59, 9.25};
Double_t sig[Npoints] = {0.10, 0.10, 0.12, 0.12, 0.15, 0.15, 0.15, 0.20, 0.20,
                         0.20, 0.20};
// result from direct calculation from Vuolo
Double_t a1  = 3.078;
Double_t ea1 = 0.065;
Double_t a2  = 0.1011;
Double_t ea2 = 0.0019;

//_____________________________________________________________________________

Double_t yteo(Double_t x, Double_t *par)
{
  Double_t value = par[0]*TMath::Exp(-x) + par[1]*TMath::Exp(4*x);
  return value;
}

//_____________________________________________________________________________

void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
   // f is the function that MINUIT really minimizes
  Double_t chis2 = 0;
  Double_t temp;
  for (Int_t i=0; i<Npoints;i++) {
    temp = (yteo(x[i],par)-y[i])/sig[i];
    // cout << "temp, sig, i = " <<temp << " "<< sig[i] <<" "<< i << endl;
    chis2 += temp*temp;
  }
  // here we minimize, so the signal is positive 
  f = chis2;
}

//_____________________________________________________________________________

void chisq()
{
  // prepare maximization 
  TMinuit *gMinuit = new TMinuit(2);  
 //initialize TMinuit with a maximum of 2 params
   gMinuit->SetFCN(fcn);

   Double_t arglist[10];
   Int_t ierflg = 0;

   // for max likelihood = 0.5, for chisq = 1.0
   arglist[0] = 1.0;
   gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);

   // Set starting values and step sizes for parameters
   static Double_t vstart[2] = {3.,1.};
   static Double_t step[2] = {0.001 , 0.001};
   gMinuit->mnparm(0, "a1", vstart[0], step[0], 0,0,ierflg);
   gMinuit->mnparm(1, "a2", vstart[1], step[1], 0,0,ierflg);

   // Now ready for minimization step
   arglist[0] = 500;
   arglist[1] = 1.;
   gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);


   // Print results
   Double_t amin,edm,errdef;
   Int_t nvpar,nparx,icstat;
   gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
   gMinuit->mnprin(3,amin);
   
   //now make the contours (taken from fitcont.C)

   TCanvas *c2 = new TCanvas("c2","contours",10,10,600,600);
   //Get contour for parameter 0 versus parameter 1  for ERRDEF=2 
   gMinuit->SetErrorDef(2);
   TGraph *gr2 = (TGraph*)gMinuit->Contour(80,0,1);
   gr2->SetFillColor(42);
   gr2->SetTitle("Exemplo do Minuit -- chis2");
   gr2->Draw("alf");
   //Get contour for parameter 0 versus parameter 1 for ERRDEF=1  
   gMinuit->SetErrorDef(1);
   TGraph *gr1 = (TGraph*)gMinuit->Contour(80,0,1);
   gr1->SetFillColor(38);
   gr1->Draw("lf");
   // get the estimated values
   Double_t a1_est, a1_est_err;
   Double_t a2_est, a2_est_err;
   gMinuit->GetParameter(0,a1_est,a1_est_err);
   gMinuit->GetParameter(1,a2_est,a2_est_err);
   // put the square as the estimated value
   estP = new TMarker(a1_est,a2_est,8);
   estP->SetMarkerColor(kBlack);
   estP->SetMarkerSize(2);
   estP->Draw();

   // draw Vuolo's value
   // it is necessary to define arrays, even for one point
   const Int_t n = 1;
   Double_t a1v[n]={a1};
   Double_t a2v[n]={a2};
   Double_t ea1v[n]={ea1};
   Double_t ea2v[n]={ea2};
   gr = new TGraphErrors(n,a1v,a2v,ea1v,ea2v);
   gr->SetMarkerColor(kRed);
   gr->SetLineColor(kRed);
   gr->SetMarkerStyle(21);
   gr->Draw("SLP");

   gStyle->SetOptStat(0);
   // draw the data histogram with fitted curve
   TCanvas *c1 = new TCanvas("c1","dados e curva",10,10,600,600);
   TH1F *h1 = new TH1F("h1", "Dados e curva ajustada", 11,0.05,1.15);
   for (int i=0;i<Npoints;i++){
     h1->SetBinContent(i+1,y[i]);
     h1->SetBinError(i+1,sig[i]);
   }
   h1->Draw("E1P");
   TF1 *fteo= new TF1("fteo","[0]*TMath::Exp(-x)+[1]*TMath::Exp(4*x)",0,1.2);
   fteo->SetParameters(a1_est,a2_est);
   fteo->Draw("SAMEP");

}