rf608_fitresultaspdf

Likelihood and minimization: representing the parabolic approximation of the fit as a multi-variate Gaussian on the parameters of the fitted pdf

Author: Clemens Lange, Wouter Verkerke (C++ version)
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Sunday, November 27, 2022 at 11:08 AM.

In [1]:
import ROOT
Welcome to JupyROOT 6.27/01

Create model and dataset

Observable

In [2]:
x = ROOT.RooRealVar("x", "x", -20, 20)

Model (intentional strong correlations)

In [3]:
mean = ROOT.RooRealVar("mean", "mean of g1 and g2", 0, -1, 1)
sigma_g1 = ROOT.RooRealVar("sigma_g1", "width of g1", 2)
g1 = ROOT.RooGaussian("g1", "g1", x, mean, sigma_g1)

sigma_g2 = ROOT.RooRealVar("sigma_g2", "width of g2", 4, 3.0, 5.0)
g2 = ROOT.RooGaussian("g2", "g2", x, mean, sigma_g2)

frac = ROOT.RooRealVar("frac", "frac", 0.5, 0.0, 1.0)
model = ROOT.RooAddPdf("model", "model", [g1, g2], [frac])
[#0] WARNING:InputArguments -- The parameter 'sigma_g1' with range [-1e+30, 1e+30] of the RooGaussian 'g1' exceeds the safe range of (0, inf). Advise to limit its range.

Generate 1000 events

In [4]:
data = model.generate({x}, 1000)

Fit model to data

In [5]:
r = model.fitTo(data, Save=True)
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization --  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:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization

Create MV Gaussian pdf of fitted parameters

In [6]:
parabPdf = r.createHessePdf({frac, mean, sigma_g2})

Some exercises with the parameter pdf

Generate 100K points in the parameter space, from the MVGaussian pdf

In [7]:
d = parabPdf.generate({mean, sigma_g2, frac}, 100000)

Sample a 3-D histogram of the pdf to be visualized as an error ellipsoid using the GLISO draw option

In [8]:
hh_3d = parabPdf.createHistogram("mean,sigma_g2,frac", 25, 25, 25)
hh_3d.SetFillColor(ROOT.kBlue)

Project 3D parameter pdf down to 3 permutations of two-dimensional pdfs The integrations corresponding to these projections are performed analytically by the MV Gaussian pdf

In [9]:
pdf_sigmag2_frac = parabPdf.createProjection({mean})
pdf_mean_frac = parabPdf.createProjection({sigma_g2})
pdf_mean_sigmag2 = parabPdf.createProjection({frac})

Make 2D plots of the 3 two-dimensional pdf projections

In [10]:
hh_sigmag2_frac = pdf_sigmag2_frac.createHistogram("sigma_g2,frac", 50, 50)
hh_mean_frac = pdf_mean_frac.createHistogram("mean,frac", 50, 50)
hh_mean_sigmag2 = pdf_mean_sigmag2.createHistogram("mean,sigma_g2", 50, 50)
hh_mean_frac.SetLineColor(ROOT.kBlue)
hh_sigmag2_frac.SetLineColor(ROOT.kBlue)
hh_mean_sigmag2.SetLineColor(ROOT.kBlue)

Draw the 'sigar'

In [11]:
ROOT.gStyle.SetCanvasPreferGL(True)
ROOT.gStyle.SetPalette(1)
c1 = ROOT.TCanvas("rf608_fitresultaspdf_1", "rf608_fitresultaspdf_1", 600, 600)
hh_3d.Draw("gliso")

c1.SaveAs("rf608_fitresultaspdf_1.png")
Info in <TCanvas::Print>: png file rf608_fitresultaspdf_1.png has been created

Draw the 2D projections of the 3D pdf

In [12]:
c2 = ROOT.TCanvas("rf608_fitresultaspdf_2", "rf608_fitresultaspdf_2", 900, 600)
c2.Divide(3, 2)
c2.cd(1)
ROOT.gPad.SetLeftMargin(0.15)
hh_mean_sigmag2.GetZaxis().SetTitleOffset(1.4)
hh_mean_sigmag2.Draw("surf3")
c2.cd(2)
ROOT.gPad.SetLeftMargin(0.15)
hh_sigmag2_frac.GetZaxis().SetTitleOffset(1.4)
hh_sigmag2_frac.Draw("surf3")
c2.cd(3)
ROOT.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 pdf

In [13]:
tmp1 = d.createHistogram(mean, sigma_g2, 50, 50)
tmp2 = d.createHistogram(sigma_g2, frac, 50, 50)
tmp3 = d.createHistogram(mean, frac, 50, 50)

c2.cd(4)
ROOT.gPad.SetLeftMargin(0.15)
tmp1.GetZaxis().SetTitleOffset(1.4)
tmp1.Draw("lego3")
c2.cd(5)
ROOT.gPad.SetLeftMargin(0.15)
tmp2.GetZaxis().SetTitleOffset(1.4)
tmp2.Draw("lego3")
c2.cd(6)
ROOT.gPad.SetLeftMargin(0.15)
tmp3.GetZaxis().SetTitleOffset(1.4)
tmp3.Draw("lego3")

c2.SaveAs("rf608_fitresultaspdf_2.png")
Info in <TCanvas::Print>: png file rf608_fitresultaspdf_2.png has been created

Draw all canvases

In [14]:
from ROOT import gROOT 
gROOT.GetListOfCanvases().Draw()