#include "TH1F.h"
#include "TF1.h"
#include "TVirtualFitter.h"
#include "TGraph2D.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TH2F.h"
#include "TStyle.h"

void LikelihoodCurve(TH1F* h, TF1* f)
{

   // we first do a 1D likelihood curve for the lifetime
   // the output will be a TGraph
   Float_t x[40];
   Float_t y[40];
   
   Float_t p0 = f->GetParameter(1); // this is the best fit.
   Float_t step = f->GetParError(1); // this is the estimate of sigma.
   Float_t a0 = f->GetParameter(0); // this is the best fit.
   Float_t astep = f->GetParError(0); // this is the estimate of sigma.

   step /= 4.;
   astep /= 4.;
   Float_t p = p0-(20*step);
   Float_t a = a0-(20*astep);
   Double_t amin,edm,errdef;
   Int_t nvpar,nparx;
   TVirtualFitter *fitter = TVirtualFitter::Fitter(h);
   fitter->GetStats(amin,edm,errdef,nvpar,nparx);
   Float_t minimum = amin;
   for(int i=0;i<40;++i,p+=step) {
     f->FixParameter(1,p);
     h->Fit(f->GetName(),"LQ0","",500,3500);
     x[i] = p;
     TVirtualFitter *fitter = TVirtualFitter::Fitter(h);
     fitter->GetStats(amin,edm,errdef,nvpar,nparx);
     y[i] = amin-minimum;
   }
   new TCanvas;
   TGraph* g1 = new TGraph(40,x,y);
   g1->Draw("al");
   
   // we do also a 2D likelihood including the amplitude as a nuisance parameter 
   // to visualize the impact of correlation
   TGraph2D *g2 = new TGraph2D();
   TH2F* h2 = new TH2F("likelihood","likelihood contour",40,p0-(20*step),p0+(20*step),40,a0-(20*astep),a0+(20*astep));
   p = p0-(20*step);
   for(int i=0;i<40;++i,p+=step) {
     a = a0-(20*astep);
     for(int j=0;j<40;++j,a+=astep) {
       f->FixParameter(0,a);
       f->FixParameter(1,p);
       h->Fit(f->GetName(),"LQ0","",500,3500);
       TVirtualFitter *fitter = TVirtualFitter::Fitter(h);
       fitter->GetStats(amin,edm,errdef,nvpar,nparx);
       g2->SetPoint(i+j*40,a,p,amin-minimum);
       h2->SetBinContent(i,j,amin-minimum);
     }
   }
   TCanvas* c3 = new TCanvas;
   c3->Divide(2,2);
   c3->cd(1);
   g2->Draw("surf1");
   c3->cd(2);
   gStyle->SetPalette(1);
   g2->Draw("cont4z");
   c3->cd(3);
   Double_t level[20];
   for(int i=1;i<=20;++i) level[i-1]=i*i/2.;
   h2->SetContour(20,level);
   h2->Draw("cont4z");
}