rf510_wsnamedsets

'ORGANIZATION AND SIMULTANEOUS FITS' RooFit tutorial macro #510

Working with named parameter sets and parameter snapshots in workspaces

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:07 AM.

In [1]:
import ROOT


def fillWorkspace(w):
    # Create model
    # -----------------------

    # Declare observable x
    x = ROOT.RooRealVar("x", "x", 0, 10)

    # Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and
    # their parameters
    mean = ROOT.RooRealVar("mean", "mean of gaussians", 5, 0, 10)
    sigma1 = ROOT.RooRealVar("sigma1", "width of gaussians", 0.5)
    sigma2 = ROOT.RooRealVar("sigma2", "width of gaussians", 1)

    sig1 = ROOT.RooGaussian("sig1", "Signal component 1", x, mean, sigma1)
    sig2 = ROOT.RooGaussian("sig2", "Signal component 2", x, mean, sigma2)

    # Build Chebychev polynomial p.d.f.
    a0 = ROOT.RooRealVar("a0", "a0", 0.5, 0.0, 1.0)
    a1 = ROOT.RooRealVar("a1", "a1", -0.2, 0.0, 1.0)
    bkg = ROOT.RooChebychev("bkg", "Background", x, [a0, a1])

    # Sum the signal components into a composite signal p.d.f.
    sig1frac = ROOT.RooRealVar("sig1frac", "fraction of component 1 in signal", 0.8, 0.0, 1.0)
    sig = ROOT.RooAddPdf("sig", "Signal", [sig1, sig2], [sig1frac])

    # Sum the composite signal and background
    bkgfrac = ROOT.RooRealVar("bkgfrac", "fraction of background", 0.5, 0.0, 1.0)
    model = ROOT.RooAddPdf("model", "g1+g2+a", [bkg, sig], [bkgfrac])

    # Import model into p.d.f.
    w.Import(model)

    # Encode definition of parameters in workspace
    # ---------------------------------------------------------------------------------------

    # Define named sets "parameters" and "observables", list which variables should be considered
    # parameters and observables by the users convention
    #
    # Variables appearing in sets _must_ live in the workspace already, the autoImport flag
    # of defineSet must be set to import them on the fly. Named sets contain only references
    # to the original variables, the value of observables in named sets already
    # reflect their 'current' value
    params = model.getParameters({x})
    w.defineSet("parameters", params)
    w.defineSet("observables", {x})

    # Encode reference value for parameters in workspace
    # ---------------------------------------------------------------------------------------------------

    # Define a parameter 'snapshot' in the p.d.f.
    # Unlike a named set, parameter snapshot stores an independent set of values for
    # a given set of variables in the workspace. The values can be stored and reloaded
    # into the workspace variable objects using the loadSnapshot() and saveSnapshot()
    # methods. A snapshot saves the value of each variable, errors that are stored
    # with it as well as the 'Constant' flag that is used in fits to determine if a
    # parameter is kept fixed or not.

    # Do a dummy fit to a (supposedly) reference dataset here and store the results
    # of that fit into a snapshot
    refData = model.generate({x}, 10000)
    model.fitTo(refData, PrintLevel=-1)

    # The kTRUE flag imports the values of the objects in (*params) into the workspace
    # If not set, present values of the workspace parameters objects are stored
    w.saveSnapshot("reference_fit", params, True)

    # Make another fit with the signal componentforced to zero
    # and save those parameters too

    bkgfrac.setVal(1)
    bkgfrac.setConstant(True)
    bkgfrac.removeError()
    model.fitTo(refData, PrintLevel=-1)

    w.saveSnapshot("reference_fit_bkgonly", params, True)
Welcome to JupyROOT 6.27/01

Create model and dataset

In [2]:
w = ROOT.RooWorkspace("w")
fillWorkspace(w)
[#0] WARNING:InputArguments -- The parameter 'sigma1' with range [-1e+30, 1e+30] of the RooGaussian 'sig1' exceeds the safe range of (0, inf). Advise to limit its range.
[#0] WARNING:InputArguments -- The parameter 'sigma2' with range [-1e+30, 1e+30] of the RooGaussian 'sig2' exceeds the safe range of (0, inf). Advise to limit its range.
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooAddPdf::model
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooChebychev::bkg
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooRealVar::x
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooRealVar::a0
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooRealVar::a1
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooRealVar::bkgfrac
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooAddPdf::sig
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooGaussian::sig1
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooRealVar::mean
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooRealVar::sigma1
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooRealVar::sig1frac
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooGaussian::sig2
[#1] INFO:ObjectHandling -- RooWorkspace::import(w) importing RooRealVar::sigma2
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization --  The following expressions will be evaluated in cache-and-track mode: (bkg,sig1,sig2)
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization --  The following expressions will be evaluated in cache-and-track mode: (bkg,sig1,sig2)
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization

Exploit convention encoded in named set "parameters" and "observables" to use workspace contents w/o need for introspected

In [3]:
model = w["model"]

Generate data from p.d.f. in given observables

In [4]:
data = model.generate(w.set("observables"), 1000)

Fit model to data

In [5]:
model.fitTo(data)
Out[5]:
<cppyy.gbl.RooFitResult object at 0x(nil)>
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization --  The following expressions will be evaluated in cache-and-track mode: (bkg,sig1,sig2)
 **********
 **    1 **SET PRINT           1
 **********
 **********
 **    2 **SET NOGRAD
 **********
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 a0           5.00000e-01  1.00000e-01    0.00000e+00  1.00000e+00
     2 a1           0.00000e+00  1.00000e-01    0.00000e+00  1.00000e+00
 MINUIT WARNING IN PARAM DEF
 ============== STARTING VALUE IS AT LIMIT.
 MINUIT WARNING IN PARAMETR
 ============== VARIABLE2 IS AT ITS LOWER ALLOWED LIMIT.
 MINUIT WARNING IN PARAMETR
 ============== VARIABLE2 BROUGHT BACK INSIDE LIMITS.
     3 bkgfrac      5.00000e-01  1.00000e-01    0.00000e+00  1.00000e+00
     4 mean         5.00000e+00  1.00000e+00    0.00000e+00  1.00000e+01
     5 sig1frac     8.00000e-01  1.00000e-01    0.00000e+00  1.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        2500           1
 **********
 FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
 MINUIT WARNING IN MIGrad    
 ============== VARIABLE2 IS AT ITS LOWER ALLOWED LIMIT.
 START MIGRAD MINIMIZATION.  STRATEGY  1.  CONVERGENCE WHEN EDM .LT. 1.00e-03
 FCN=1921.72 FROM MIGRAD    STATUS=INITIATE       37 CALLS          38 TOTAL
                     EDM= unknown      STRATEGY= 1      NO ERROR MATRIX       
  EXT PARAMETER               CURRENT GUESS       STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  a0           5.00000e-01   1.00000e-01   0.00000e+00   1.05962e+01
   2  a1           1.20069e-01   1.00000e-01   7.05742e-01   5.42891e-01
   3  bkgfrac      5.00000e-01   1.00000e-01   0.00000e+00   3.05354e+01
   4  mean         5.00000e+00   1.00000e+00   0.00000e+00  -7.51085e+01
   5  sig1frac     8.00000e-01   1.00000e-01   0.00000e+00  -6.73860e+00
                               ERR DEF= 0.5
 MIGRAD MINIMIZATION HAS CONVERGED.
 MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=1919.24 FROM MIGRAD    STATUS=CONVERGED     130 CALLS         131 TOTAL
                     EDM=3.30126e-05    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  a0           3.74916e-01   7.69275e-02   4.69528e-03   1.66844e-02
   2  a1           8.17702e-02   1.21821e-01   8.86749e-03  -6.33902e-04
   3  bkgfrac      4.91367e-01   3.69596e-02   1.23375e-03   3.69153e-02
   4  mean         5.01943e+00   2.99569e-02   1.80527e-04   6.61892e-01
   5  sig1frac     8.84317e-01   1.16911e-01   6.63479e-03  -1.33559e-02
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  5    ERR DEF=0.5
  5.968e-03  2.115e-03 -4.931e-04 -2.130e-04 -1.575e-03 
  2.115e-03  1.594e-02 -3.476e-03 -5.180e-05 -1.046e-02 
 -4.931e-04 -3.476e-03  1.369e-03 -1.294e-05  3.488e-03 
 -2.130e-04 -5.180e-05 -1.294e-05  8.974e-04 -4.427e-05 
 -1.575e-03 -1.046e-02  3.488e-03 -4.427e-05  1.432e-02 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2      3      4      5
        1  0.23649   1.000  0.217 -0.173 -0.092 -0.170
        2  0.76815   0.217  1.000 -0.744 -0.014 -0.692
        3  0.83467  -0.173 -0.744  1.000 -0.012  0.788
        4  0.09976  -0.092 -0.014 -0.012  1.000 -0.012
        5  0.80365  -0.170 -0.692  0.788 -0.012  1.000
 **********
 **    7 **SET ERR         0.5
 **********
 **********
 **    8 **SET PRINT           1
 **********
 **********
 **    9 **HESSE        2500
 **********
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=1919.24 FROM HESSE     STATUS=OK             31 CALLS         162 TOTAL
                     EDM=3.28745e-05    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                INTERNAL      INTERNAL  
  NO.   NAME      VALUE            ERROR       STEP SIZE       VALUE   
   1  a0           3.74916e-01   7.68702e-02   9.39056e-04  -2.52855e-01
   2  a1           8.17702e-02   1.20796e-01   3.54700e-04  -9.90791e-01
   3  bkgfrac      4.91367e-01   3.66875e-02   2.46750e-04  -1.72673e-02
   4  mean         5.01943e+00   2.99563e-02   3.61054e-05   3.88610e-03
   5  sig1frac     8.84317e-01   1.16255e-01   1.32696e-03   8.76705e-01
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  5    ERR DEF=0.5
  5.959e-03  2.063e-03 -4.793e-04 -2.126e-04 -1.535e-03 
  2.063e-03  1.565e-02 -3.395e-03 -4.999e-05 -1.022e-02 
 -4.793e-04 -3.395e-03  1.348e-03 -1.340e-05  3.430e-03 
 -2.126e-04 -4.999e-05 -1.340e-05  8.974e-04 -4.452e-05 
 -1.535e-03 -1.022e-02  3.430e-03 -4.452e-05  1.416e-02 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2      3      4      5
        1  0.23348   1.000  0.214 -0.169 -0.092 -0.167
        2  0.76323   0.214  1.000 -0.739 -0.013 -0.687
        3  0.83196  -0.169 -0.739  1.000 -0.012  0.785
        4  0.09958  -0.092 -0.013 -0.012  1.000 -0.012
        5  0.80103  -0.167 -0.687  0.785 -0.012  1.000
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization

Plot fitted model and data on frame of first (only) observable

In [6]:
frame = (w.set("observables").first()).frame()
data.plotOn(frame)
model.plotOn(frame)
Out[6]:
<cppyy.gbl.RooPlot object at 0x9fe1810>

Overlay plot with model with reference parameters as stored in snapshots

In [7]:
w.loadSnapshot("reference_fit")
model.plotOn(frame, LineColor="r")
w.loadSnapshot("reference_fit_bkgonly")
model.plotOn(frame, LineColor="r", LineStyle="--")
Out[7]:
<cppyy.gbl.RooPlot object at 0x9fe1810>

Draw the frame on the canvas

In [8]:
c = ROOT.TCanvas("rf510_wsnamedsets", "rf503_wsnamedsets", 600, 600)
ROOT.gPad.SetLeftMargin(0.15)
frame.GetYaxis().SetTitleOffset(1.4)
frame.Draw()

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

Print workspace contents

In [9]:
w.Print()
RooWorkspace(w) w contents

variables
---------
(a0,a1,bkgfrac,mean,sig1frac,sigma1,sigma2,x)

p.d.f.s
-------
RooChebychev::bkg[ x=x coefList=(a0,a1) ] = 1
RooAddPdf::model[ bkgfrac * bkg + [%] * sig ] = 1/1
RooAddPdf::sig[ sig1frac * sig1 + [%] * sig2 ] = 0.0131991/1
RooGaussian::sig1[ x=x mean=mean sigma=sigma1 ] = 1.11947e-05
RooGaussian::sig2[ x=x mean=mean sigma=sigma2 ] = 0.0578433

parameter snapshots
-------------------
reference_fit = (a0=0.484427 +/- 0.0229437,a1=1.85852e-06 +/- 0.615506,bkgfrac=0.489779 +/- 0.00797979,mean=5.01099 +/- 0.0103497,sigma1=0.5[C],sig1frac=0.771962 +/- 0.0274623,sigma2=1[C])
reference_fit_bkgonly = (a0=0.469238 +/- 0.0224872,a1=0 +/- 0.000154899,bkgfrac=1[C],mean=2.61253 +/- 2.39847,sigma1=0.5[C],sig1frac=0.771962 +/- 0.353053,sigma2=1[C])

named sets
----------
observables:(x)
parameters:(a0,a1,bkgfrac,mean,sig1frac,sigma1,sigma2)

Workspace will remain in memory after macro finishes

In [10]:
ROOT.gDirectory.Add(w)

Draw all canvases

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