Rf 7 0 2Efficiencyfit 2 D

Special pdf's: unbinned maximum likelihood fit of an efficiency eff(x) function to a dataset D(x,cut), cut is a category encoding a selection whose efficiency as function of x should be described by eff(x)

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 Saturday, November 28, 2020 at 11:07 AM.

In [1]:
import ROOT


flat = ROOT.kFALSE
Welcome to JupyROOT 6.23/01

Construct efficiency function e(x,y)

Declare variables x,mean, with associated name, title, value and allowed range

In [2]:
x = ROOT.RooRealVar("x", "x", -10, 10)
y = ROOT.RooRealVar("y", "y", -10, 10)
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

Efficiency function eff(x;a,b)

In [3]:
ax = ROOT.RooRealVar("ax", "ay", 0.6, 0, 1)
bx = ROOT.RooRealVar("bx", "by", 5)
cx = ROOT.RooRealVar("cx", "cy", -1, -10, 10)

ay = ROOT.RooRealVar("ay", "ay", 0.2, 0, 1)
by = ROOT.RooRealVar("by", "by", 5)
cy = ROOT.RooRealVar("cy", "cy", -1, -10, 10)

effFunc = ROOT.RooFormulaVar(
    "effFunc",
    "((1-ax)+ax*cos((x-cx)/bx))*((1-ay)+ay*cos((y-cy)/by))",
    ROOT.RooArgList(
        ax,
        bx,
        cx,
        x,
        ay,
        by,
        cy,
        y))

Acceptance state cut (1 or 0)

In [4]:
cut = ROOT.RooCategory("cut", "cutr")
cut.defineType("accept", 1)
cut.defineType("reject", 0)
Out[4]:
False

Construct conditional efficiency pdf E(cut|x,y)

Construct efficiency pdf eff(cut|x)

In [5]:
effPdf = ROOT.RooEfficiency("effPdf", "effPdf", effFunc, cut, "accept")

Generate data(x,y,cut) from a toy model

Construct global shape pdf shape(x) and product model(x,cut) = eff(cut|x)*shape(x) (These are only needed to generate some toy MC here to be used later)

In [6]:
shapePdfX = ROOT.RooPolynomial(
    "shapePdfX", "shapePdfX", x, ROOT.RooArgList(
        ROOT.RooFit.RooConst(
            0 if flat else -0.095)))
shapePdfY = ROOT.RooPolynomial(
    "shapePdfY", "shapePdfY", y, ROOT.RooArgList(
        ROOT.RooFit.RooConst(
            0 if flat else +0.095)))
shapePdf = ROOT.RooProdPdf(
    "shapePdf",
    "shapePdf",
    ROOT.RooArgList(
        shapePdfX,
        shapePdfY))
model = ROOT.RooProdPdf(
    "model",
    "model",
    ROOT.RooArgSet(shapePdf),
    ROOT.RooFit.Conditional(
        ROOT.RooArgSet(effPdf),
        ROOT.RooArgSet(cut)))

Generate some toy data from model

In [7]:
data = model.generate(ROOT.RooArgSet(x, y, cut), 10000)
[#0] WARNING:Generation -- RooAcceptReject::ctor(effPdf_Int[]_Norm[cut]) WARNING: performing accept/reject sampling on a p.d.f in 2 dimensions without prior knowledge on maximum value of p.d.f. Determining maximum value by taking 200000 trial samples. If p.d.f contains sharp peaks smaller than average distance between trial sampling points these may be missed and p.d.f. may be sampled incorrectly.
[#0] WARNING:Generation -- RooAcceptReject::ctor(effPdf_Int[]_Norm[cut]): WARNING: 200000 trial samples requested by p.d.f for 2-dimensional accept/reject sampling, this may take some time

Fit conditional efficiency pdf to data

Fit conditional efficiency pdf to data

In [8]:
effPdf.fitTo(data, ROOT.RooFit.ConditionalObservables(ROOT.RooArgSet(x, y)))
Out[8]:
<cppyy.gbl.RooFitResult object at 0x(nil)>
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
 **********
 **    1 **SET PRINT           1
 **********
 **********
 **    2 **SET NOGRAD
 **********
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 ax           6.00000e-01  1.00000e-01    0.00000e+00  1.00000e+00
     2 ay           2.00000e-01  1.00000e-01    0.00000e+00  1.00000e+00
     3 cx          -1.00000e+00  2.00000e+00   -1.00000e+01  1.00000e+01
     4 cy          -1.00000e+00  2.00000e+00   -1.00000e+01  1.00000e+01
 **********
 **    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        2000           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=5447.45 FROM MIGRAD    STATUS=INITIATE       16 CALLS          17 TOTAL
                     EDM= unknown      STRATEGY= 1      NO ERROR MATRIX       
  EXT PARAMETER               CURRENT GUESS       STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  ax           6.00000e-01   1.00000e-01   2.05758e-01   2.07611e+01
   2  ay           2.00000e-01   1.00000e-01   2.57889e-01   6.35766e+01
   3  cx          -1.00000e+00   2.00000e+00   2.02430e-01   3.98906e+02
   4  cy          -1.00000e+00   2.00000e+00   2.02430e-01  -1.37299e+02
                               ERR DEF= 0.5
 MIGRAD MINIMIZATION HAS CONVERGED.
 MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=5441.98 FROM MIGRAD    STATUS=CONVERGED      79 CALLS          80 TOTAL
                     EDM=2.10777e-05    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  ax           6.10442e-01   9.85749e-03   8.90018e-04   1.06907e-01
   2  ay           2.06325e-01   1.06819e-02   9.76478e-04   2.34590e-01
   3  cx          -1.13975e+00   5.97322e-02   2.77681e-04   6.58278e-02
   4  cy          -5.30427e-01   2.15761e-01   8.34515e-04  -2.03235e-01
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  4    ERR DEF=0.5
  9.718e-05 -2.419e-05 -2.246e-04  2.779e-05 
 -2.419e-05  1.141e-04 -2.514e-05  1.447e-03 
 -2.246e-04 -2.514e-05  3.568e-03  1.894e-05 
  2.779e-05  1.447e-03  1.894e-05  4.656e-02 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2      3      4
        1  0.50176   1.000 -0.230 -0.381  0.013
        2  0.68607  -0.230  1.000 -0.039  0.628
        3  0.42040  -0.381 -0.039  1.000  0.001
        4  0.65573   0.013  0.628  0.001  1.000
 **********
 **    7 **SET ERR         0.5
 **********
 **********
 **    8 **SET PRINT           1
 **********
 **********
 **    9 **HESSE        2000
 **********
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=5441.98 FROM HESSE     STATUS=OK             23 CALLS         103 TOTAL
                     EDM=2.1074e-05    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                INTERNAL      INTERNAL  
  NO.   NAME      VALUE            ERROR       STEP SIZE       VALUE   
   1  ax           6.10442e-01   9.85365e-03   1.78004e-04   2.22720e-01
   2  ay           2.06325e-01   1.06842e-02   1.95296e-04  -6.27781e-01
   3  cx          -1.13975e+00   5.97056e-02   5.55362e-05  -1.14223e-01
   4  cy          -5.30427e-01   2.15797e-01   1.66903e-04  -5.30676e-02
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  4    ERR DEF=0.5
  9.711e-05 -2.421e-05 -2.238e-04  2.765e-05 
 -2.421e-05  1.142e-04 -2.530e-05  1.448e-03 
 -2.238e-04 -2.530e-05  3.565e-03  1.965e-05 
  2.765e-05  1.448e-03  1.965e-05  4.658e-02 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2      3      4
        1  0.50117   1.000 -0.230 -0.380  0.013
        2  0.68624  -0.230  1.000 -0.040  0.628
        3  0.41953  -0.380 -0.040  1.000  0.002
        4  0.65587   0.013  0.628  0.002  1.000
[#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization

Plot fitted, data efficiency

Make 2D histograms of all data, data and efficiency function

In [9]:
hh_data_all = ROOT.RooAbsData.createHistogram(
    data, "hh_data_all", x, ROOT.RooFit.Binning(8), ROOT.RooFit.YVar(
        y, ROOT.RooFit.Binning(8)))
hh_data_sel = ROOT.RooAbsData.createHistogram(
    data, "hh_data_sel", x, ROOT.RooFit.Binning(8), ROOT.RooFit.YVar(
        y, ROOT.RooFit.Binning(8)), ROOT.RooFit.Cut("cut==cut::accept"))
hh_eff = effFunc.createHistogram(
    "hh_eff", x, ROOT.RooFit.Binning(
        50), ROOT.RooFit.YVar(y, ROOT.RooFit.Binning(50)))
[#0] WARNING:InputArguments -- RooAbsReal::createHistogram(effFunc) WARNING extended mode requested for a non-pdf object, ignored

Some adjustsment for good visualization

In [10]:
hh_data_all.SetMinimum(0)
hh_data_sel.SetMinimum(0)
hh_eff.SetMinimum(0)
hh_eff.SetLineColor(ROOT.kBlue)

Draw all frames on a canvas

In [11]:
ca = ROOT.TCanvas("rf702_efficiency_2D", "rf702_efficiency_2D", 1200, 400)
ca.Divide(3)
ca.cd(1)
ROOT.gPad.SetLeftMargin(0.15)
hh_data_all.GetZaxis().SetTitleOffset(1.8)
hh_data_all.Draw("lego")
ca.cd(2)
ROOT.gPad.SetLeftMargin(0.15)
hh_data_sel.GetZaxis().SetTitleOffset(1.8)
hh_data_sel.Draw("lego")
ca.cd(3)
ROOT.gPad.SetLeftMargin(0.15)
hh_eff.GetZaxis().SetTitleOffset(1.8)
hh_eff.Draw("surf")

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

Draw all canvases

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