Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
35 if mode == 1:
36 # Instruct model to precalculate normalization integral that integrate at least
37 # two dimensions numerically. In self specific case the integral value for
38 # all values of parameter 'a' are stored in a histogram and available for use
39 # in subsequent fitting and plotting operations (interpolation is
40 # applied)
41
42 # w.pdf("model").setNormValueCaching(3)
43 model = w["model"]
44 model.setStringAttribute("CACHEPARMINT", "x:y:z")
45
46 # Evaluate pdf once to trigger filling of cache
47 normSet = {w["x"], w["y"], w["z"]}
48 model.getVal(normSet)
49 w.writeToFile("rf903_numintcache.root")
50
51 if mode == 2:
52 # Load preexisting workspace from file in mode==2
53 f = ROOT.TFile("rf903_numintcache.root")
54 w = f.Get("w")
55
56 # Return created or loaded workspace
57 return w
58
59
60mode = 0
61# Mode = 0 : Run plain fit (slow)
62# Mode = 1 : Generate workspace with precalculated integral and store it on file (prepare for accelerated running)
63# Mode = 2 : Run fit from previously stored workspace including cached
64# integrals (fast, run in mode=1 first)
65
66# Create, save or load workspace with pdf
67# -----------------------------------------------------------------------------------
68
69# Make/load workspace, here in mode 1
70w = getWorkspace(mode)
71if mode == 1:
72 # Show workspace that was created
73 w.Print()
74
75 # Show plot of cached integral values
76 hhcache = w.expensiveObjectCache().getObj(1)
77 if hhcache:
78 ROOT.TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600)
79 hhcache.createHistogram("a").Draw()
80 else:
81 ROOT.RooFit.Error("rf903_numintcache", "Cached histogram is not existing in workspace")
82 sys.exit()
83
84# Use pdf from workspace for generation and fitting
85# -----------------------------------------------------------------------------------
86
87# ROOT.This is always slow (need to find maximum function value
88# empirically in 3D space)
89model = w["model"]
90d = model.generate({w["x"], w["y"], w["z"]}, 1000)
91
92# ROOT.This is slow in mode 0, fast in mode 1
93model.fitTo(d, Verbose=True, Timer=True)
94
95# Projection on x (always slow as 2D integral over Y, at fitted value of a
96# is not cached)
97framex = w["x"].frame(Title="Projection of 3D model on X")
98d.plotOn(framex)
99model.plotOn(framex)
100
101# Draw x projection on canvas
102c = ROOT.TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600)
103framex.Draw()
104
105c.SaveAs("rf903_numintcache.png")
106
107# Make workspace available on command line after macro finishes
108ROOT.gDirectory.Add(w)
th1 Draw()