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 Sunday, February 05, 2023 at 11:14 AM.
import ROOT
Welcome to JupyROOT 6.29/01
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)
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
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])
Add signal components
Sum the signal components into a composite signal pdf
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])
Generate a data sample of 1000 events in x from model
data = model.generate({x}, 1000)
Fit model to data
model.fitTo(data)
<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
xframe = x.frame(Title="Example of composite pdf=(sig1+sig2)+bkg")
data.plotOn(xframe)
model.plotOn(xframe)
<cppyy.gbl.RooPlot object at 0xef5d7c0>
Overlay the background component of model with a dashed line
model.plotOn(xframe, Components={bkg}, LineStyle="--")
<cppyy.gbl.RooPlot object at 0xef5d7c0>
[#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
model.plotOn(xframe, Components={bkg, sig2}, LineStyle=":")
<cppyy.gbl.RooPlot object at 0xef5d7c0>
[#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
model.Print("t")
0xb648d50 RooAddPdf::model = 0.885987/1 [Auto,Clean] 0xa9c3f70/V- RooChebychev::bkg = 0.733485 [Auto,Dirty] 0xd2790c0/V- RooRealVar::x = 5 0xdd8e4a0/V- RooRealVar::a0 = 0.507642 +/- 0.0795738 0xddc1250/V- RooRealVar::a1 = 0.266515 +/- 0.13378 0x96328c0/V- RooRealVar::bkgfrac = 0.427793 +/- 0.0355511 0x962df70/V- RooAddPdf::sig = 1/1 [Auto,Clean] 0xdd934a0/V- RooGaussian::sig1 = 1 [Auto,Dirty] 0xd2790c0/V- RooRealVar::x = 5 0xda376e0/V- RooRealVar::mean = 5 0x1019740/V- RooRealVar::sigma1 = 0.5 0xe100e20/V- RooRealVar::sig1frac = 0.641234 +/- 0.0967457 0xdddef00/V- RooGaussian::sig2 = 1 [Auto,Dirty] 0xd2790c0/V- RooRealVar::x = 5 0xda376e0/V- RooRealVar::mean = 5 0x10112f0/V- RooRealVar::sigma2 = 1
Construct sum of models on one go using recursive fraction interpretations
model2 = bkg + (sig1 + sig2)
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))
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: () 0xefd7d30 RooAddPdf::model = 0.885987/1 [Auto,Clean] 0xa9c3f70/V- RooChebychev::bkg = 0.733485 [Auto,Dirty] 0xd2790c0/V- RooRealVar::x = 5 0xdd8e4a0/V- RooRealVar::a0 = 0.507642 +/- 0.0795738 0xddc1250/V- RooRealVar::a1 = 0.266515 +/- 0.13378 0x96328c0/V- RooRealVar::bkgfrac = 0.427793 +/- 0.0355511 0xdd934a0/V- RooGaussian::sig1 = 1 [Auto,Dirty] 0xd2790c0/V- RooRealVar::x = 5 0xda376e0/V- RooRealVar::mean = 5 0x1019740/V- RooRealVar::sigma1 = 0.5 0xeffbc80/V- RooRecursiveFraction::model_recursive_fraction_sig1_2 = 0.366918 [Auto,Clean] 0xe100e20/V- RooRealVar::sig1frac = 0.641234 +/- 0.0967457 0x96328c0/V- RooRealVar::bkgfrac = 0.427793 +/- 0.0355511 0xdddef00/V- RooGaussian::sig2 = 1 [Auto,Dirty] 0xd2790c0/V- RooRealVar::x = 5 0xda376e0/V- RooRealVar::mean = 5 0x10112f0/V- RooRealVar::sigma2 = 1 0xf04d530/V- RooRecursiveFraction::model_recursive_fraction_sig2_3 = 0.205289 [Auto,Clean] 0xee26590/V- RooConstVar::1 = 1 0xe100e20/V- RooRealVar::sig1frac = 0.641234 +/- 0.0967457 0x96328c0/V- RooRealVar::bkgfrac = 0.427793 +/- 0.0355511
Draw the frame on the canvas
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
from ROOT import gROOT
gROOT.GetListOfCanvases().Draw()