Logo ROOT  
Reference Guide
rf903_numintcache.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roofit
3 /// \notebook -js
4 /// Numeric algorithm tuning: caching of slow numeric integrals and parameterization of slow numeric integrals
5 ///
6 /// \macro_image
7 /// \macro_output
8 /// \macro_code
9 ///
10 /// \date July 2008
11 /// \author Wouter Verkerke
12 
13 #include "RooRealVar.h"
14 #include "RooDataSet.h"
15 #include "RooDataHist.h"
16 #include "RooGaussian.h"
17 #include "RooConstVar.h"
18 #include "TCanvas.h"
19 #include "TAxis.h"
20 #include "RooPlot.h"
21 #include "RooWorkspace.h"
23 #include "TFile.h"
24 #include "TH1.h"
25 
26 using namespace RooFit;
27 
28 RooWorkspace *getWorkspace(Int_t mode);
29 
30 void rf903_numintcache(Int_t mode = 0)
31 {
32  // Mode = 0 : Run plain fit (slow)
33  // Mode = 1 : Generate workspace with pre-calculated integral and store it on file (prepare for accelerated running)
34  // Mode = 2 : Run fit from previously stored workspace including cached integrals (fast, requires run in mode=1
35  // first)
36 
37  // C r e a t e , s a v e o r l o a d w o r k s p a c e w i t h p . d . f .
38  // -----------------------------------------------------------------------------------
39 
40  // Make/load workspace, exit here in mode 1
41  RooWorkspace *w1 = getWorkspace(mode);
42  if (mode == 1) {
43 
44  // Show workspace that was created
45  w1->Print();
46 
47  // Show plot of cached integral values
48  RooDataHist *hhcache = (RooDataHist *)w1->expensiveObjectCache().getObj(1);
49  if (hhcache) {
50 
51  new TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600);
52  hhcache->createHistogram("a")->Draw();
53 
54  } else {
55  Error("rf903_numintcache", "Cached histogram is not existing in workspace");
56  }
57  return;
58  }
59 
60  // U s e p . d . f . f r o m w o r k s p a c e f o r g e n e r a t i o n a n d f i t t i n g
61  // -----------------------------------------------------------------------------------
62 
63  // This is always slow (need to find maximum function value empirically in 3D space)
64  RooDataSet *d = w1->pdf("model")->generate(RooArgSet(*w1->var("x"), *w1->var("y"), *w1->var("z")), 1000);
65 
66  // This is slow in mode 0, but fast in mode 1
67  w1->pdf("model")->fitTo(*d, Verbose(kTRUE), Timer(kTRUE));
68 
69  // Projection on x (always slow as 2D integral over Y,Z at fitted value of a is not cached)
70  RooPlot *framex = w1->var("x")->frame(Title("Projection of 3D model on X"));
71  d->plotOn(framex);
72  w1->pdf("model")->plotOn(framex);
73 
74  // Draw x projection on canvas
75  auto canv = new TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600);
76  framex->Draw();
77  canv->Draw();
78 
79  // Make workspace available on command line after macro finishes
80  gDirectory->Add(w1);
81 }
82 
83 RooWorkspace *getWorkspace(Int_t mode)
84 {
85  // C r e a t e , s a v e o r l o a d w o r k s p a c e w i t h p . d . f .
86  // -----------------------------------------------------------------------------------
87  //
88  // Mode = 0 : Create workspace for plain running (no integral caching)
89  // Mode = 1 : Generate workspace with pre-calculated integral and store it on file
90  // Mode = 2 : Load previously stored workspace from file
91 
92  RooWorkspace *w(0);
93 
94  if (mode != 2) {
95 
96  // Create empty workspace workspace
97  w = new RooWorkspace("w", 1);
98 
99  // Make a difficult to normalize pdf in 3 dimensions that is integrated numerically.
100  w->factory("EXPR::model('1/((x-a)*(x-a)+0.01)+1/((y-a)*(y-a)+0.01)+1/"
101  "((z-a)*(z-a)+0.01)',x[-1,1],y[-1,1],z[-1,1],a[-5,5])");
102  }
103 
104  if (mode == 1) {
105 
106  // Instruct model to pre-calculate normalization integral that integrate at least
107  // two dimensions numerically. In this specific case the integral value for
108  // all values of parameter 'a' are stored in a histogram and available for use
109  // in subsequent fitting and plotting operations (interpolation is applied)
110 
111  // w->pdf("model")->setNormValueCaching(3) ;
112  w->pdf("model")->setStringAttribute("CACHEPARMINT", "x:y:z");
113 
114  // Evaluate pdf once to trigger filling of cache
115  RooArgSet normSet(*w->var("x"), *w->var("y"), *w->var("z"));
116  w->pdf("model")->getVal(&normSet);
117  w->writeToFile("rf903_numintcache.root");
118  }
119 
120  if (mode == 2) {
121  // Load preexisting workspace from file in mode==2
122  TFile *f = new TFile("rf903_numintcache.root");
123  w = (RooWorkspace *)f->Get("w");
124  }
125 
126  // Return created or loaded workspace
127  return w;
128 }
RooAbsRealLValue::frame
RooPlot * frame(const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Create a new RooPlot on the heap with a drawing frame initialized for this object,...
Definition: RooAbsRealLValue.cxx:199
RooWorkspace.h
RooPlot::Draw
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:691
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:91
f
#define f(i)
Definition: RSha256.hxx:104
RooGaussian.h
TCanvas.h
rf903_numintcache
Definition: rf903_numintcache.py:1
RooDataSet.h
TFile.h
RooAbsPdf::plotOn
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none(), const RooCmdArg &arg9=RooCmdArg::none(), const RooCmdArg &arg10=RooCmdArg::none()) const
Helper calling plotOn(RooPlot*, RooLinkedList&) const.
Definition: RooAbsPdf.h:121
RooWorkspace::expensiveObjectCache
RooExpensiveObjectCache & expensiveObjectCache()
Definition: RooWorkspace.h:165
RooDataHist
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:39
RooWorkspace::Print
void Print(Option_t *opts=0) const
Print contents of the workspace.
Definition: RooWorkspace.cxx:2194
RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: RooCFunction1Binding.h:29
RooFit::Verbose
RooCmdArg Verbose(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:189
RooDataHist.h
RooAbsData::createHistogram
TH1 * createHistogram(const char *name, const RooAbsRealLValue &xvar, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Calls createHistogram(const char *name, const RooAbsRealLValue& xvar, const RooLinkedList& argList) c...
Definition: RooAbsData.cxx:629
RooPlot.h
gDirectory
#define gDirectory
Definition: TDirectory.h:236
RooPlot
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:44
RooWorkspace::pdf
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
Definition: RooWorkspace.cxx:1277
RooRealVar.h
TFile
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:54
RooExpensiveObjectCache.h
RooConstVar.h
RooWorkspace
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
TCanvas
The Canvas class.
Definition: TCanvas.h:23
RooAbsPdf::generate
RooDataSet * generate(const RooArgSet &whatVars, Int_t nEvents, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none())
See RooAbsPdf::generate(const RooArgSet&,const RooCmdArg&,const RooCmdArg&,const RooCmdArg&,...
Definition: RooAbsPdf.h:58
TAxis.h
d
#define d(i)
Definition: RSha256.hxx:102
RooWorkspace::var
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found.
Definition: RooWorkspace.cxx:1295
RooDataSet
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
RooExpensiveObjectCache::getObj
const TObject * getObj(Int_t uniqueID)
Retrieve payload object of cache element with given unique ID
Definition: RooExpensiveObjectCache.cxx:154
TH1.h
RooFit::Title
RooCmdArg Title(const char *name)
Definition: RooGlobalFunc.cxx:176
RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:29
int
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:3050
RooFit::Timer
RooCmdArg Timer(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:191
Error
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:187
RooAbsPdf::fitTo
virtual RooFitResult * fitTo(RooAbsData &data, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none())
Fit PDF to given dataset.
Definition: RooAbsPdf.cxx:1264