//////////////////////////////////////////////////////////////////////////
//
// 'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #608
// 
// Representing the parabolic approximation of the fit as
// a multi-variate Gaussian on the parameters of the fitted p.d.f.
//
//
// 07/2008 - Wouter Verkerke 
// 
/////////////////////////////////////////////////////////////////////////

#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooAddPdf.h"
#include "RooChebychev.h"
#include "RooFitResult.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "TFile.h"
#include "TStyle.h"
#include "TH2.h"
#include "TH3.h"

using namespace RooFit ;


void rf608_fitresultaspdf()
{
  // C r e a t e   m o d e l   a n d   d a t a s e t 
  // -----------------------------------------------

  // Observable
  RooRealVar x("x","x",-20,20) ;

  // Model (intentional strong correlations)
  RooRealVar mean("mean","mean of g1 and g2",0,-1,1) ;
  RooRealVar sigma_g1("sigma_g1","width of g1",2) ; 
  RooGaussian g1("g1","g1",x,mean,sigma_g1) ;

  RooRealVar sigma_g2("sigma_g2","width of g2",4,3.0,5.0) ;
  RooGaussian g2("g2","g2",x,mean,sigma_g2) ;

  RooRealVar frac("frac","frac",0.5,0.0,1.0) ;
  RooAddPdf model("model","model",RooArgList(g1,g2),frac) ;

  // Generate 1000 events
  RooDataSet* data = model.generate(x,1000) ;


  // F i t   m o d e l   t o   d a t a 
  // ----------------------------------

  RooFitResult* r = model.fitTo(*data,Save()) ;


  // C r e a t e M V   G a u s s i a n   p d f   o f   f i t t e d    p a r a m e t e r s
  // ------------------------------------------------------------------------------------

  RooAbsPdf* parabPdf = r->createHessePdf(RooArgSet(frac,mean,sigma_g2)) ;


  // S o m e   e x e c e r c i s e s   w i t h   t h e   p a r a m e t e r   p d f 
  // -----------------------------------------------------------------------------

  // Generate 100K points in the parameter space, sampled from the MVGaussian p.d.f.
  RooDataSet* d = parabPdf->generate(RooArgSet(mean,sigma_g2,frac),100000) ;


  // Sample a 3-D histogram of the p.d.f. to be visualized as an error ellipsoid using the GLISO draw option
  TH3* hh_3d = (TH3*) parabPdf->createHistogram("mean,sigma_g2,frac",25,25,25) ;
  hh_3d->SetFillColor(kBlue) ;  


  // Project 3D parameter p.d.f. down to 3 permutations of two-dimensional p.d.f.s 
  // The integrations corresponding to these projections are performed analytically
  // by the MV Gaussian p.d.f.
  RooAbsPdf* pdf_sigmag2_frac = parabPdf->createProjection(mean) ;
  RooAbsPdf* pdf_mean_frac    = parabPdf->createProjection(sigma_g2) ;
  RooAbsPdf* pdf_mean_sigmag2 = parabPdf->createProjection(frac) ;


  // Make 2D plots of the 3 two-dimensional p.d.f. projections
  TH2* hh_sigmag2_frac = (TH2*) pdf_sigmag2_frac->createHistogram("sigma_g2,frac",50,50) ;
  TH2* hh_mean_frac    = (TH2*) pdf_mean_frac->createHistogram("mean,frac",50,50) ;
  TH2* hh_mean_sigmag2 = (TH2*) pdf_mean_sigmag2->createHistogram("mean,sigma_g2",50,50) ;
  hh_mean_frac->SetLineColor(kBlue) ;
  hh_sigmag2_frac->SetLineColor(kBlue) ;
  hh_mean_sigmag2->SetLineColor(kBlue) ;


  // Draw the 'sigar'
  gStyle->SetCanvasPreferGL(true);
  gStyle->SetPalette(1) ;
  new TCanvas("rf608_fitresultaspdf_1","rf608_fitresultaspdf_1",600,600) ;
  hh_3d->Draw("gliso") ; 

  // Draw the 2D projections of the 3D p.d.f.
  TCanvas* c2 = new TCanvas("rf608_fitresultaspdf_2","rf608_fitresultaspdf_2",900,600) ;
  c2->Divide(3,2) ;
  c2->cd(1) ; gPad->SetLeftMargin(0.15) ; hh_mean_sigmag2->GetZaxis()->SetTitleOffset(1.4) ; hh_mean_sigmag2->Draw("surf3") ; 
  c2->cd(2) ; gPad->SetLeftMargin(0.15) ; hh_sigmag2_frac->GetZaxis()->SetTitleOffset(1.4) ; hh_sigmag2_frac->Draw("surf3") ; 
  c2->cd(3) ; gPad->SetLeftMargin(0.15) ; hh_mean_frac->GetZaxis()->SetTitleOffset(1.4) ; hh_mean_frac->Draw("surf3") ; 

  // Draw the distributions of parameter points sampled from the p.d.f.
  TH1* tmp1 = d->createHistogram("mean,sigma_g2",50,50) ;
  TH1* tmp2 = d->createHistogram("sigma_g2,frac",50,50) ;
  TH1* tmp3 = d->createHistogram("mean,frac",50,50) ;

  c2->cd(4) ; gPad->SetLeftMargin(0.15) ; tmp1->GetZaxis()->SetTitleOffset(1.4) ; tmp1->Draw("lego3") ;
  c2->cd(5) ; gPad->SetLeftMargin(0.15) ; tmp2->GetZaxis()->SetTitleOffset(1.4) ; tmp2->Draw("lego3") ;
  c2->cd(6) ; gPad->SetLeftMargin(0.15) ; tmp3->GetZaxis()->SetTitleOffset(1.4) ; tmp3->Draw("lego3") ;

}