Logo ROOT  
Reference Guide
rf903_numintcache.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roofit
3 /// \notebook -js
4 ///
5 /// Numeric algorithm tuning: caching of slow numeric integrals and parameterization of slow numeric integrals
6 ///
7 /// \macro_image
8 /// \macro_output
9 /// \macro_code
10 ///
11 /// \date 07/2008
12 /// \author Wouter Verkerke
13 
14 #include "RooRealVar.h"
15 #include "RooDataSet.h"
16 #include "RooDataHist.h"
17 #include "RooGaussian.h"
18 #include "RooConstVar.h"
19 #include "TCanvas.h"
20 #include "TAxis.h"
21 #include "RooPlot.h"
22 #include "RooWorkspace.h"
24 #include "TFile.h"
25 #include "TH1.h"
26 
27 using namespace RooFit;
28 
29 RooWorkspace *getWorkspace(Int_t mode);
30 
31 void rf903_numintcache(Int_t mode = 0)
32 {
33  // Mode = 0 : Run plain fit (slow)
34  // Mode = 1 : Generate workspace with pre-calculated integral and store it on file (prepare for accelerated running)
35  // Mode = 2 : Run fit from previously stored workspace including cached integrals (fast, requires run in mode=1
36  // first)
37 
38  // 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 .
39  // -----------------------------------------------------------------------------------
40 
41  // Make/load workspace, exit here in mode 1
42  RooWorkspace *w1 = getWorkspace(mode);
43  if (mode == 1) {
44 
45  // Show workspace that was created
46  w1->Print();
47 
48  // Show plot of cached integral values
49  RooDataHist *hhcache = (RooDataHist *)w1->expensiveObjectCache().getObj(1);
50  if (hhcache) {
51 
52  new TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600);
53  hhcache->createHistogram("a")->Draw();
54 
55  } else {
56  Error("rf903_numintcache", "Cached histogram is not existing in workspace");
57  }
58  return;
59  }
60 
61  // 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
62  // -----------------------------------------------------------------------------------
63 
64  // This is always slow (need to find maximum function value empirically in 3D space)
65  RooDataSet *d = w1->pdf("model")->generate(RooArgSet(*w1->var("x"), *w1->var("y"), *w1->var("z")), 1000);
66 
67  // This is slow in mode 0, but fast in mode 1
68  w1->pdf("model")->fitTo(*d, Verbose(kTRUE), Timer(kTRUE));
69 
70  // Projection on x (always slow as 2D integral over Y,Z at fitted value of a is not cached)
71  RooPlot *framex = w1->var("x")->frame(Title("Projection of 3D model on X"));
72  d->plotOn(framex);
73  w1->pdf("model")->plotOn(framex);
74 
75  // Draw x projection on canvas
76  new TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600);
77  framex->Draw();
78 
79  // Make workspace available on command line after macro finishes
80  gDirectory->Add(w1);
81 
82  return;
83 }
84 
85 RooWorkspace *getWorkspace(Int_t mode)
86 {
87  // 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 .
88  // -----------------------------------------------------------------------------------
89  //
90  // Mode = 0 : Create workspace for plain running (no integral caching)
91  // Mode = 1 : Generate workspace with pre-calculated integral and store it on file
92  // Mode = 2 : Load previously stored workspace from file
93 
94  RooWorkspace *w(0);
95 
96  if (mode != 2) {
97 
98  // Create empty workspace workspace
99  w = new RooWorkspace("w", 1);
100 
101  // Make a difficult to normalize p.d.f. in 3 dimensions that is integrated numerically.
102  w->factory("EXPR::model('1/((x-a)*(x-a)+0.01)+1/((y-a)*(y-a)+0.01)+1/"
103  "((z-a)*(z-a)+0.01)',x[-1,1],y[-1,1],z[-1,1],a[-5,5])");
104  }
105 
106  if (mode == 1) {
107 
108  // Instruct model to pre-calculate normalization integral that integrate at least
109  // two dimensions numerically. In this specific case the integral value for
110  // all values of parameter 'a' are stored in a histogram and available for use
111  // in subsequent fitting and plotting operations (interpolation is applied)
112 
113  // w->pdf("model")->setNormValueCaching(3) ;
114  w->pdf("model")->setStringAttribute("CACHEPARMINT", "x:y:z");
115 
116  // Evaluate p.d.f. once to trigger filling of cache
117  RooArgSet normSet(*w->var("x"), *w->var("y"), *w->var("z"));
118  w->pdf("model")->getVal(&normSet);
119  w->writeToFile("rf903_numintcache.root");
120  }
121 
122  if (mode == 2) {
123  // Load preexisting workspace from file in mode==2
124  TFile *f = new TFile("rf903_numintcache.root");
125  w = (RooWorkspace *)f->Get("w");
126  }
127 
128  // Return created or loaded workspace
129  return w;
130 }
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:122
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:118
RooWorkspace::expensiveObjectCache
RooExpensiveObjectCache & expensiveObjectCache()
Definition: RooWorkspace.h:165
RooDataHist
Definition: RooDataHist.h:39
RooWorkspace::Print
void Print(Option_t *opts=0) const
Print contents of the workspace.
Definition: RooWorkspace.cxx:2194
RooFit
Definition: RooCFunction1Binding.h:29
RooFit::Verbose
RooCmdArg Verbose(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:186
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
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
Definition: TFile.h:54
RooExpensiveObjectCache.h
RooConstVar.h
RooWorkspace
Definition: RooWorkspace.h:43
TCanvas
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:55
TAxis.h
d
#define d(i)
Definition: RSha256.hxx:120
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
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:173
RooArgSet
Definition: RooArgSet.h:28
int
TH1::Draw
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2997
RooFit::Timer
RooCmdArg Timer(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:188
Error
void Error(const char *location, const char *msgfmt,...)
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:1261