Rf 2 0 1_Composite

Addition and convolution: composite pdf with signal and background component

pdf = f_bkg * bkg(x,a0,a1) + (1-fbkg) * (f_sig1 * sig1(x,m,s1 + (1-f_sig1) * sig2(x,m,s2)))

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 Friday, May 13, 2022 at 09:18 AM.

In [1]:
import ROOT

Welcome to JupyROOT 6.27/01

Setup component pdfs

Declare observable x

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

Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters

In [3]:
mean = ROOT.RooRealVar("mean", "mean of gaussians", 5)
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)

[#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.

Build Chebychev polynomial pdf

In [4]:
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])

Method 1 - Two RooAddPdfs

Add signal components

Sum the signal components into a composite signal pdf

In [5]:
sig1frac = ROOT.RooRealVar("sig1frac", "fraction of component 1 in signal", 0.8, 0.0, 1.0)
sig = ROOT.RooAddPdf("sig", "Signal", [sig1, sig2], [sig1frac])

Add signal and background

Sum the composite signal and background

In [6]:
bkgfrac = ROOT.RooRealVar("bkgfrac", "fraction of background", 0.5, 0.0, 1.0)
model = ROOT.RooAddPdf("model", "g1+g2+a", [bkg, sig], [bkgfrac])

Sample, fit and plot model

Generate a data sample of 1000 events in x from model

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

[#0] WARNING:Eval -- Evaluating RooAddPdf without a defined normalization set. This can lead to ambiguos coefficients definition and incorrect results. Use RooAddPdf::fixCoefNormalization(nset) to provide a normalization set for defining uniquely RooAddPdf coefficients!

Fit model to data

In [8]:
model.fitTo(data)
Out[8]:
<cppyy.gbl.RooFitResult object at 0x(nil)>
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization --  The following expressions have been identified as constant and will be precalculated and cached: (sig1,sig2)
[#1] INFO:Minimization --  The following expressions will be evaluated in cache-and-track mode: (bkg)
 **********
 **    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 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        2000           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=1926.2 FROM MIGRAD    STATUS=INITIATE       31 CALLS          32 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   2.31850e+00
   2  a1           7.50936e-02   1.00000e-01   5.53213e-01  -4.13475e-03
   3  bkgfrac      5.00000e-01   1.00000e-01   0.00000e+00   2.67741e+01
   4  sig1frac     8.00000e-01   1.00000e-01   0.00000e+00   3.58316e-01
                               ERR DEF= 0.5
 MIGRAD MINIMIZATION HAS CONVERGED.
 MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=1924.17 FROM MIGRAD    STATUS=CONVERGED     119 CALLS         120 TOTAL
                     EDM=2.77609e-07    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  a0           5.07642e-01   7.96101e-02   4.71991e-03  -9.40696e-04
   2  a1           2.66515e-01   1.34621e-01   5.72412e-03  -6.75205e-05
   3  bkgfrac      4.27793e-01   3.57760e-02   1.24946e-03  -1.00877e-02
   4  sig1frac     6.41234e-01   9.72166e-02   4.23189e-03   1.73342e-03
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  4    ERR DEF=0.5
  6.392e-03  2.434e-03 -5.449e-04 -1.352e-03 
  2.434e-03  1.871e-02 -3.804e-03 -8.768e-03 
 -5.449e-04 -3.804e-03  1.282e-03  2.494e-03 
 -1.352e-03 -8.768e-03  2.494e-03  9.583e-03 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2      3      4
        1  0.22587   1.000  0.223 -0.190 -0.173
        2  0.79303   0.223  1.000 -0.777 -0.655
        3  0.82183  -0.190 -0.777  1.000  0.712
        4  0.73007  -0.173 -0.655  0.712  1.000
 **********
 **    7 **SET ERR         0.5
 **********
 **********
 **    8 **SET PRINT           1
 **********
 **********
 **    9 **HESSE        2000
 **********
 COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=1924.17 FROM HESSE     STATUS=OK             23 CALLS         143 TOTAL
                     EDM=2.75649e-07    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                INTERNAL      INTERNAL  
  NO.   NAME      VALUE            ERROR       STEP SIZE       VALUE   
   1  a0           5.07642e-01   7.95738e-02   1.88796e-04   1.52838e-02
   2  a1           2.66515e-01   1.33780e-01   2.28965e-04  -4.85861e-01
   3  bkgfrac      4.27793e-01   3.55511e-02   2.49893e-04  -1.44921e-01
   4  sig1frac     6.41234e-01   9.67457e-02   8.46378e-04   2.86365e-01
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  4    ERR DEF=0.5
  6.386e-03  2.397e-03 -5.349e-04 -1.328e-03 
  2.397e-03  1.847e-02 -3.740e-03 -8.613e-03 
 -5.349e-04 -3.740e-03  1.266e-03  2.455e-03 
 -1.328e-03 -8.613e-03  2.455e-03  9.489e-03 
 PARAMETER  CORRELATION COEFFICIENTS  
       NO.  GLOBAL      1      2      3      4
        1  0.22394   1.000  0.221 -0.188 -0.171
        2  0.78997   0.221  1.000 -0.773 -0.651
        3  0.81931  -0.188 -0.773  1.000  0.708
        4  0.72690  -0.171 -0.651  0.708  1.000
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization

Plot data and PDF overlaid

In [9]:
xframe = x.frame(Title="Example of composite pdf=(sig1+sig2)+bkg")
data.plotOn(xframe)
model.plotOn(xframe)
Out[9]:
<cppyy.gbl.RooPlot object at 0x27a12a0>

Overlay the background component of model with a dashed line

In [10]:
model.plotOn(xframe, Components={bkg}, LineStyle="--")
Out[10]:
<cppyy.gbl.RooPlot object at 0x27a12a0>
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (bkg)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()

Overlay the background+sig2 components of model with a dotted line

In [11]:
model.plotOn(xframe, Components={bkg, sig2}, LineStyle=":")
Out[11]:
<cppyy.gbl.RooPlot object at 0x27a12a0>
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (bkg,sig2)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: (sig)

Print structure of composite pdf

In [12]:
model.Print("t")

0x8d5b0e0 RooAddPdf::model = 0.885987 [Auto,Dirty] 
  0x8ce9c60/V- RooChebychev::bkg = 0.733485 [Auto,Dirty] 
    0x88b50b0/V- RooRealVar::x = 5
    0x8976450/V- RooRealVar::a0 = 0.507642 +/- 0.0795738
    0x89975f0/V- RooRealVar::a1 = 0.266515 +/- 0.13378
  0x8d59540/V- RooRealVar::bkgfrac = 0.427793 +/- 0.0355511
  0x8d18140/V- RooAddPdf::sig = 1 [Auto,Dirty] 
    0x8a025c0/V- RooGaussian::sig1 = 1 [Auto,Dirty] 
      0x88b50b0/V- RooRealVar::x = 5
      0x88ac5c0/V- RooRealVar::mean = 5
      0x7802f30/V- RooRealVar::sigma1 = 0.5
    0x8c69910/V- RooRealVar::sig1frac = 0.641234 +/- 0.0967457
    0x89c4720/V- RooGaussian::sig2 = 1 [Auto,Dirty] 
      0x88b50b0/V- RooRealVar::x = 5
      0x88ac5c0/V- RooRealVar::mean = 5
      0x877f230/V- RooRealVar::sigma2 = 1

Method 2 - One RooAddPdf with recursive fractions

Construct sum of models on one go using recursive fraction interpretations

model2 = bkg + (sig1 + sig2)

In [13]:
model2 = ROOT.RooAddPdf("model", "g1+g2+a", [bkg, sig1, sig2], [bkgfrac, sig1frac], True)

NB: Each coefficient is interpreted as the fraction of the left-hand component of the i-th recursive sum, i.e.

sum4 = A + ( B + ( C + D) with fraction fA, and fC expands to

sum4 = fAA + (1-fA)(fBB + (1-fB)(fCC + (1-fC)D))

Plot recursive addition model

In [14]:
model2.plotOn(xframe, LineColor="r", LineStyle="--")
model2.plotOn(xframe, Components={bkg, sig2}, LineColor="r", LineStyle="--")
model2.Print("t")

[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) directly selected PDF components: (bkg,sig2)
[#1] INFO:Plotting -- RooAbsPdf::plotOn(model) indirectly selected PDF components: ()
0x9548fc0 RooAddPdf::model = 0.885987 [Auto,Dirty] 
  0x8ce9c60/V- RooChebychev::bkg = 0.733485 [Auto,Dirty] 
    0x88b50b0/V- RooRealVar::x = 5
    0x8976450/V- RooRealVar::a0 = 0.507642 +/- 0.0795738
    0x89975f0/V- RooRealVar::a1 = 0.266515 +/- 0.13378
  0x8d59540/V- RooRealVar::bkgfrac = 0.427793 +/- 0.0355511
  0x8a025c0/V- RooGaussian::sig1 = 1 [Auto,Dirty] 
    0x88b50b0/V- RooRealVar::x = 5
    0x88ac5c0/V- RooRealVar::mean = 5
    0x7802f30/V- RooRealVar::sigma1 = 0.5
  0x9360ea0/V- RooRecursiveFraction::model_recursive_fraction_sig1 = 0.366918 [Auto,Clean] 
    0x8c69910/V- RooRealVar::sig1frac = 0.641234 +/- 0.0967457
    0x8d59540/V- RooRealVar::bkgfrac = 0.427793 +/- 0.0355511
  0x89c4720/V- RooGaussian::sig2 = 1 [Auto,Dirty] 
    0x88b50b0/V- RooRealVar::x = 5
    0x88ac5c0/V- RooRealVar::mean = 5
    0x877f230/V- RooRealVar::sigma2 = 1
  0x9368f50/V- RooRecursiveFraction::model_recursive_fraction_sig2 = 0.205289 [Auto,Clean] 
    0x830ca70/V- RooConstVar::1 = 1
    0x8c69910/V- RooRealVar::sig1frac = 0.641234 +/- 0.0967457
    0x8d59540/V- RooRealVar::bkgfrac = 0.427793 +/- 0.0355511

Draw the frame on the canvas

In [15]:
c = ROOT.TCanvas("rf201_composite", "rf201_composite", 600, 600)
ROOT.gPad.SetLeftMargin(0.15)
xframe.GetYaxis().SetTitleOffset(1.4)
xframe.Draw()

c.SaveAs("rf201_composite.png")

Info in <TCanvas::Print>: png file rf201_composite.png has been created

Draw all canvases

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