ROOT   Reference Guide
Searching...
No Matches
rf211_paramconv.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_roofit
3## \notebook
4## 'ADDITION AND CONVOLUTION' RooFit tutorial macro #211
5## Working a with a p.d.f. with a convolution operator in terms
6## of a parameter
7##
8## (require ROOT to be compiled with --enable-fftw3)
9##
10## \macro_image
11## \macro_code
12## \macro_output
13##
14## \date February 2018
15## \authors Clemens Lange, Wouter Verkerke (C version)
16
17import ROOT
18
19# Set up component pdfs
20# ---------------------------------------
21
22# Gaussian g(x ; mean,sigma)
23x = ROOT.RooRealVar("x", "x", -10, 10)
24mean = ROOT.RooRealVar("mean", "mean", -3, 3)
25sigma = ROOT.RooRealVar("sigma", "sigma", 0.5, 0.1, 10)
26modelx = ROOT.RooGaussian("gx", "gx", x, mean, sigma)
27
28# Block function in mean
29a = ROOT.RooRealVar("a", "a", 2, 1, 10)
30model_mean = ROOT.RooGenericPdf("model_mean", "abs(mean)<a", [mean, a])
31
32# Convolution in mean model = g(x,mean,sigma) (x) block(mean)
33x.setBins(1000, "cache")
34mean.setBins(50, "cache")
35model = ROOT.RooFFTConvPdf("model", "model", mean, modelx, model_mean)
36
37# Configure convolution to construct a 2-D cache in (x,mean)
38# rather than a 1-d cache in mean that needs to be recalculated
39# for each value of x
40model.setCacheObservables({x})
41model.setBufferFraction(1.0)
42
43# Integrate model over projModel = Int model dmean
44projModel = model.createProjection({mean})
45
46# Generate 1000 toy events
47d = projModel.generateBinned({x}, 1000)
48
49# Fit p.d.f. to toy data
50projModel.fitTo(d, Verbose=True, PrintLevel=-1)
51
52# Plot data and fitted p.d.f.
53frame = x.frame(Bins=25)
54d.plotOn(frame)
55projModel.plotOn(frame)
56
57# Make 2d histogram of model(x;mean)
58hh = model.createHistogram(
59 "hh",
60 x,
61 Binning=50,
62 YVar=dict(var=mean, Binning=50),
63 ConditionalObservables={mean},
64)
65hh.SetTitle("histogram of model(x|mean)")
66hh.SetLineColor(ROOT.kBlue)
67
68# Draw frame on canvas
69c = ROOT.TCanvas("rf211_paramconv", "rf211_paramconv", 800, 400)
70c.Divide(2)
71c.cd(1)