Logo ROOT  
Reference Guide
rf903_numintcache.py
Go to the documentation of this file.
1## \file
2## \ingroup tutorial_roofit
3## \notebook
4## Numeric algorithm tuning: caching of slow numeric integrals and parameterizations of slow numeric integrals
5##
6## \macro_code
7##
8## \date February 2018
9## \authors Clemens Lange, Wouter Verkerke (C++ version)
10
11import sys
12import ROOT
13
14
15def getWorkspace(mode):
16 # Create, save or load workspace with pdf
17 # -----------------------------------------------------------------------------------
18 #
19 # Mode = 0 : Create workspace for plain running (no integral caching)
20 # Mode = 1 : Generate workspace with precalculated integral and store it on file
21 # Mode = 2 : Load previously stored workspace from file
22
23 w = ROOT.RooWorkspace()
24
25 if mode != 2:
26 # Create empty workspace workspace
27 w = ROOT.RooWorkspace("w", 1)
28
29 # Make a difficult to normalize pdf in 3 dimensions that is
30 # integrated numerically.
31 w.factory(
32 "EXPR::model('1/((x-a)*(x-a)+0.01)+1/((y-a)*(y-a)+0.01)+1/((z-a)*(z-a)+0.01)',x[-1,1],y[-1,1],z[-1,1],a[-5,5])")
33
34 if mode == 1:
35 # Instruct model to precalculate normalization integral that integrate at least
36 # two dimensions numerically. In self specific case the integral value for
37 # all values of parameter 'a' are stored in a histogram and available for use
38 # in subsequent fitting and plotting operations (interpolation is
39 # applied)
40
41 # w.pdf("model").setNormValueCaching(3)
42 w.pdf("model").setStringAttribute("CACHEPARMINT", "x:y:z")
43
44 # Evaluate pdf once to trigger filling of cache
45 normSet = ROOT.RooArgSet(w.var("x"), w.var("y"), w.var("z"))
46 w.pdf("model").getVal(normSet)
47 w.writeToFile("rf903_numintcache.root")
48
49 if (mode == 2):
50 # Load preexisting workspace from file in mode==2
51 f = ROOT.TFile("rf903_numintcache.root")
52 w = f.Get("w")
53
54 # Return created or loaded workspace
55 return w
56
57
58mode = 0
59# Mode = 0 : Run plain fit (slow)
60# Mode = 1 : Generate workspace with precalculated integral and store it on file (prepare for accelerated running)
61# Mode = 2 : Run fit from previously stored workspace including cached
62# integrals (fast, run in mode=1 first)
63
64# Create, save or load workspace with pdf
65# -----------------------------------------------------------------------------------
66
67# Make/load workspace, here in mode 1
68w = getWorkspace(mode)
69if mode == 1:
70 # Show workspace that was created
71 w.Print()
72
73 # Show plot of cached integral values
74 hhcache = w.expensiveObjectCache().getObj(1)
75 if (hhcache):
76 ROOT.TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600)
77 hhcache.createHistogram("a").Draw()
78 else:
79 ROOT.RooFit.Error("rf903_numintcache",
80 "Cached histogram is not existing in workspace")
81 sys.exit()
82
83# Use pdf from workspace for generation and fitting
84# -----------------------------------------------------------------------------------
85
86# ROOT.This is always slow (need to find maximum function value
87# empirically in 3D space)
88d = w.pdf("model").generate(
89 ROOT.RooArgSet(
90 w.var("x"),
91 w.var("y"),
92 w.var("z")),
93 1000)
94
95# ROOT.This is slow in mode 0, fast in mode 1
96w.pdf("model").fitTo(
97 d, ROOT.RooFit.Verbose(
98 ROOT.kTRUE), ROOT.RooFit.Timer(
99 ROOT.kTRUE))
100
101# Projection on x (always slow as 2D integral over Y, at fitted value of a
102# is not cached)
103framex = w.var("x").frame(ROOT.RooFit.Title("Projection of 3D model on X"))
104d.plotOn(framex)
105w.pdf("model").plotOn(framex)
106
107# Draw x projection on canvas
108c = ROOT.TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600)
109framex.Draw()
110
111c.SaveAs("rf903_numintcache.png")
112
113# Make workspace available on command line after macro finishes
114ROOT.gDirectory.Add(w)
th1 Draw()