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 
11 import sys
12 import ROOT
13 
14 
15 def 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 
58 mode = 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
68 w = getWorkspace(mode)
69 if 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)
88 d = 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
96 w.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)
103 framex = w.var("x").frame(ROOT.RooFit.Title("Projection of 3D model on X"))
104 d.plotOn(framex)
105 w.pdf("model").plotOn(framex)
106 
107 # Draw x projection on canvas
108 c = ROOT.TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600)
109 framex.Draw()
110 
111 c.SaveAs("rf903_numintcache.png")
112 
113 # Make workspace available on command line after macro finishes
114 ROOT.gDirectory.Add(w)
Draw
th1 Draw()