rf303_conditional

'MULTIDIMENSIONAL MODELS' RooFit tutorial macro #303 Use of tailored p.d.f as conditional p.d.fs.s

pdf = gauss(x,f(y),sx | y ) with f(y) = a0 + a1*y

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 Wednesday, November 30, 2022 at 11:22 AM.

In [1]:
import ROOT


def makeFakeDataXY():
    x = ROOT.RooRealVar("x", "x", -10, 10)
    y = ROOT.RooRealVar("y", "y", -10, 10)
    coord = {x, y}

    d = ROOT.RooDataSet("d", "d", coord)

    for i in range(10000):
        tmpy = ROOT.gRandom.Gaus(0, 10)
        tmpx = ROOT.gRandom.Gaus(0.5 * tmpy, 1)
        if (abs(tmpy) < 10) and (abs(tmpx) < 10):
            x.setVal(tmpx)
            y.setVal(tmpy)
            d.add(coord)

    return d
Welcome to JupyROOT 6.27/01

Set up composed model gauss(x, m(y), s)

Create observables

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

Create function f(y) = a0 + a1*y

In [3]:
a0 = ROOT.RooRealVar("a0", "a0", -0.5, -5, 5)
a1 = ROOT.RooRealVar("a1", "a1", -0.5, -1, 1)
fy = ROOT.RooPolyVar("fy", "fy", y, [a0, a1])

Creat gauss(x,f(y),s)

In [4]:
sigma = ROOT.RooRealVar("sigma", "width of gaussian", 0.5, 0.1, 2.0)
model = ROOT.RooGaussian("model", "Gaussian with shifting mean", x, fy, sigma)

Obtain fake external experimental dataset with values for x and y

In [5]:
expDataXY = makeFakeDataXY()

Generate data from conditional p.d.f. model(x|y)

Make subset of experimental data with only y values

In [6]:
expDataY = expDataXY.reduce({y})

Generate 10000 events in x obtained from conditional model(x|y) with y values taken from experimental data

In [7]:
data = model.generate({x}, ProtoData=expDataY)
data.Print()
RooDataSet::modelData[x,y] = 6850 entries

Fit conditional p.d.f model(x|y) to data

In [8]:
model.fitTo(expDataXY, ConditionalObservables={y})
Out[8]:
<cppyy.gbl.RooFitResult object at 0x(nil)>
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
 **********
 **    1 **SET PRINT           1
 **********
 **********
 **    2 **SET NOGRAD
 **********
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 a0          -5.00000e-01  1.00000e+00   -5.00000e+00  5.00000e+00
     2 a1          -5.00000e-01  2.00000e-01   -1.00000e+00  1.00000e+00
     3 sigma        5.00000e-01  1.90000e-01    1.00000e-01  2.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=421420 FROM MIGRAD    STATUS=INITIATE       12 CALLS          13 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+00   2.02430e-01  -7.46265e+04
   2  a1          -5.00000e-01   2.00000e-01   2.35352e-01  -6.95347e+05
   3  sigma        5.00000e-01   1.90000e-01   2.52163e-01  -1.29056e+06
                               ERR DEF= 0.5
 MIGRAD MINIMIZATION HAS CONVERGED.
 MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=9659.64 FROM MIGRAD    STATUS=CONVERGED     101 CALLS         102 TOTAL
                     EDM=8.18253e-05    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  a0           9.03100e-03   1.19768e-02   1.62548e-04  -2.31829e+00
   2  a1           5.02815e-01   2.21631e-03   1.74026e-04  -2.33580e+00
   3  sigma        9.91234e-01   8.46817e-03   6.05957e-04  -4.27913e-01
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  3    ERR DEF=0.5
  1.434e-04  1.880e-07  4.948e-08 
  1.880e-07  4.912e-06  8.995e-09 
  4.948e-08  8.995e-09  7.171e-05 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2      3
        1  0.00710   1.000  0.007  0.000
        2  0.00710   0.007  1.000  0.000
        3  0.00068   0.000  0.000  1.000
 **********
 **    7 **SET ERR         0.5
 **********
 **********
 **    8 **SET PRINT           1
 **********
 **********
 **    9 **HESSE        1500
 **********
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=9659.64 FROM HESSE     STATUS=OK             16 CALLS         118 TOTAL
                     EDM=8.17764e-05    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                INTERNAL      INTERNAL  
  NO.   NAME      VALUE            ERROR       STEP SIZE       VALUE   
   1  a0           9.03100e-03   1.19768e-02   3.25095e-05   1.80620e-03
   2  a1           5.02815e-01   2.21631e-03   3.48052e-05   2.61474e+00
   3  sigma        9.91234e-01   8.46818e-03   1.21191e-04  -6.18987e-02
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  3    ERR DEF=0.5
  1.434e-04  1.880e-07  2.151e-09 
  1.880e-07  4.912e-06  2.334e-10 
  2.151e-09  2.334e-10  7.171e-05 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2      3
        1  0.00708   1.000  0.007  0.000
        2  0.00708   0.007  1.000  0.000
        3  0.00002   0.000  0.000  1.000
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization

Project conditional p.d.f on x and y dimensions

Plot x distribution of data and projection of model x = 1/Ndata sum(data(y_i)) model(x;y_i)

In [9]:
xframe = x.frame()
expDataXY.plotOn(xframe)
model.plotOn(xframe, ProjWData=expDataY)
Out[9]:
<cppyy.gbl.RooPlot object at 0x8b6c620>
[#1] INFO:Plotting -- RooAbsReal::plotOn(model) plot on x averages using data variables (y)
[#1] INFO:Plotting -- RooDataWeightedAverage::ctor(modelDataWgtAvg) constructing data weighted average of function model_Norm[x] over 6850 data points of (y) with a total weight of 6850
.........................................................................................................................................................................................................

Speed up (and approximate) projection by using binned clone of data for projection

In [10]:
binnedDataY = expDataY.binnedClone()
model.plotOn(xframe, ProjWData=binnedDataY, LineColor="c", LineStyle=":")
Out[10]:
<cppyy.gbl.RooPlot object at 0x8b6c620>
[#1] INFO:Plotting -- RooAbsReal::plotOn(model) plot on x averages using data variables (y)
[#1] INFO:Plotting -- RooDataWeightedAverage::ctor(modelDataWgtAvg) constructing data weighted average of function model_Norm[x] over 100 data points of (y) with a total weight of 6850
.........................................................................................................................................................................................................

Show effect of projection with too coarse binning

In [11]:
(expDataY.get().find("y")).setBins(5)
binnedDataY2 = expDataY.binnedClone()
model.plotOn(xframe, ProjWData=binnedDataY2, LineColor="r")
Out[11]:
<cppyy.gbl.RooPlot object at 0x8b6c620>
[#1] INFO:Plotting -- RooAbsReal::plotOn(model) plot on x averages using data variables (y)
[#1] INFO:Plotting -- RooDataWeightedAverage::ctor(modelDataWgtAvg) constructing data weighted average of function model_Norm[x] over 5 data points of (y) with a total weight of 6850
.................................................................................................................................................................................................................................................

Make canvas and draw ROOT.RooPlots

In [12]:
c = ROOT.TCanvas("rf303_conditional", "rf303_conditional", 600, 460)
ROOT.gPad.SetLeftMargin(0.15)
xframe.GetYaxis().SetTitleOffset(1.2)
xframe.Draw()

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

Draw all canvases

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