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