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