rf608_fitresultaspdf.C File Reference

Detailed Description

'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.

Processing /mnt/build/workspace/root-makedoc-v612/rootspi/rdoc/src/v6-12-00-patches/tutorials/roofit/rf608_fitresultaspdf.C...

RooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby 
                Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
                All rights reserved, please read http://roofit.sourceforge.net/license.txt

[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#1] INFO:Minization --  The following expressions will be evaluated in cache-and-track mode: (g1,g2)
 **********
 **    1 **SET PRINT           1
 **********
 **********
 **    2 **SET NOGRAD
 **********
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 frac         5.00000e-01  1.00000e-01    0.00000e+00  1.00000e+00
     2 mean         0.00000e+00  2.00000e-01   -1.00000e+00  1.00000e+00
     3 sigma_g2     4.00000e+00  2.00000e-01    3.00000e+00  5.00000e+00
 **********
 **    3 **SET ERR         0.5
 **********
 **********
 **    4 **SET PRINT           1
 **********
 **********
 **    5 **SET STR           1
 **********
 NOW USING STRATEGY  1: TRY TO BALANCE SPEED AGAINST RELIABILITY
 **********
 **    6 **MIGRAD        1500           1
 **********
 FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
 START MIGRAD MINIMIZATION.  STRATEGY  1.  CONVERGENCE WHEN EDM .LT. 1.00e-03
 FCN=2523.5 FROM MIGRAD    STATUS=INITIATE       10 CALLS          11 TOTAL
                     EDM= unknown      STRATEGY= 1      NO ERROR MATRIX       
  EXT PARAMETER               CURRENT GUESS       STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  frac         5.00000e-01   1.00000e-01   2.01358e-01  -1.31418e+01
   2  mean         0.00000e+00   2.00000e-01   2.01358e-01   4.68351e+00
   3  sigma_g2     4.00000e+00   2.00000e-01   2.01358e-01   5.45887e+00
                               ERR DEF= 0.5
 MIGRAD MINIMIZATION HAS CONVERGED.
 MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=2522.85 FROM MIGRAD    STATUS=CONVERGED      50 CALLS          51 TOTAL
                     EDM=4.43971e-05    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  frac         5.43820e-01   5.60284e-02   2.70800e-03   8.34315e-02
   2  mean        -4.18828e-02   9.05559e-02   3.14649e-03  -5.37978e-03
   3  sigma_g2     4.01280e+00   2.08953e-01   4.84455e-03  -4.10552e-02
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  3    ERR DEF=0.5
  3.153e-03 -3.621e-04  8.647e-03 
 -3.621e-04  8.223e-03 -1.611e-03 
  8.647e-03 -1.611e-03  4.431e-02 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2      3
        1  0.73165   1.000 -0.071  0.732
        2  0.08551  -0.071  1.000 -0.084
        3  0.73231   0.732 -0.084  1.000
 **********
 **    7 **SET ERR         0.5
 **********
 **********
 **    8 **SET PRINT           1
 **********
 **********
 **    9 **HESSE        1500
 **********
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=2522.85 FROM HESSE     STATUS=OK             16 CALLS          67 TOTAL
                     EDM=4.43683e-05    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                INTERNAL      INTERNAL  
  NO.   NAME      VALUE            ERROR       STEP SIZE       VALUE   
   1  frac         5.43820e-01   5.60802e-02   5.41601e-04   8.77525e-02
   2  mean        -4.18828e-02   9.05561e-02   6.29298e-04  -4.18951e-02
   3  sigma_g2     4.01280e+00   2.09149e-01   9.68910e-04   1.27960e-02
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  3    ERR DEF=0.5
  3.158e-03 -3.614e-04  8.670e-03 
 -3.614e-04  8.223e-03 -1.615e-03 
  8.670e-03 -1.615e-03  4.440e-02 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2      3
        1  0.73224   1.000 -0.071  0.732
        2  0.08555  -0.071  1.000 -0.085
        3  0.73291   0.732 -0.085  1.000
[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization

#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'
   new TCanvas("rf608_fitresultaspdf_1","rf608_fitresultaspdf_1",600,600) ;
   hh_3d->Draw("iso") ; 

   // 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") ;

}

Author

07/2008 - Wouter Verkerke

Definition in file rf608_fitresultaspdf.C.