'LIKELIHOOD AND MINIMIZATION' RooFit tutorial macro #602
Setting up a chi^2 fit to a binned dataset
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, April 17, 2024 at 11:19 AM.
from __future__ import print_function
import ROOT
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 [-inf, inf] of the RooGaussian 'sig1' exceeds the safe range of (0, inf). Advise to limit its range. [#0] WARNING:InputArguments -- The parameter 'sigma2' with range [-inf, inf] of the RooGaussian 'sig2' exceeds the safe range of (0, inf). Advise to limit its range.
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])
d = model.generate({x}, 10000)
dh = d.binnedClone()
Construct a chi^2 of the data and the model. When a p.d.f. is used in a chi^2 fit, probability density scaled by the number of events in the dataset to obtain the fit function If model is an extended p.d.f, expected number events is used instead of the observed number of events.
ll = ROOT.RooLinkedList()
model.chi2FitTo(dh, ll)
<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) Minuit2Minimizer: Minimize with max-calls 2000 convergence for edm < 1 strategy 1 Minuit2Minimizer : Valid minimum - status = 0 FVAL = 104.639633447510988 Edm = 0.000778057047730882148 Nfcn = 70 a0 = 0.501526 +/- 0.0229096 (limited) a1 = 0.158456 +/- 0.0368354 (limited) bkgfrac = 0.506609 +/- 0.011349 (limited) sig1frac = 0.815448 +/- 0.0373695 (limited) [#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
Info in <Minuit2>: MnSeedGenerator Computing seed using NumericalGradient calculator Info in <Minuit2>: MnSeedGenerator Initial state: FCN = 106.1754789 Edm = 1.260070585 NCalls = 17 Info in <Minuit2>: MnSeedGenerator Initial state Minimum value : 106.1754789 Edm : 1.260070585 Internal parameters: [ 0 -0.6435011088 0 0.6435011088] Internal gradient : [ -9.229193845 35.69317399 24.70460566 9.713828671] Internal covariance matrix: [[ 0.0021104285 0 0 0] [ 0 0.0034850852 0 0] [ 0 0 0.0001669053 0] [ 0 0 0 0.0033769849]]] Info in <Minuit2>: VariableMetricBuilder Start iterating until Edm is < 0.002 with call limit = 2000 Info in <Minuit2>: VariableMetricBuilder 0 - FCN = 106.1754789 Edm = 1.260070585 NCalls = 17 Info in <Minuit2>: VariableMetricBuilder 1 - FCN = 105.1128131 Edm = 0.08526485556 NCalls = 27 Info in <Minuit2>: VariableMetricBuilder 2 - FCN = 104.86813 Edm = 0.1538210632 NCalls = 37 Info in <Minuit2>: VariableMetricBuilder 3 - FCN = 104.6396334 Edm = 0.0006219311146 NCalls = 47 Info in <Minuit2>: VariableMetricBuilder After Hessian Info in <Minuit2>: VariableMetricBuilder 4 - FCN = 104.6396334 Edm = 0.0007780570477 NCalls = 70 Info in <Minuit2>: Minuit2Minimizer::Hesse Using max-calls 2000 Info in <Minuit2>: Minuit2Minimizer::Hesse Hesse is valid - matrix is accurate
NB: It is also possible to fit a ROOT.RooAbsReal function to a ROOT.RooDataHist using chi2FitTo().
Note that entries with zero bins are not allowed for a proper chi^2 calculation and will give error messages
dsmall = d.reduce(ROOT.RooFit.EventRange(1, 100))
dhsmall = dsmall.binnedClone()
chi2_lowstat = model.createChi2(dhsmall)
print(chi2_lowstat.getVal())
90.86495991394004
Draw all canvases
from ROOT import gROOT
gROOT.GetListOfCanvases().Draw()