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
26using namespace RooFit;
27
28RooWorkspace *getWorkspace(Int_t mode);
29
30void 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
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
83RooWorkspace *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}
#define d(i)
Definition: RSha256.hxx:102
#define f(i)
Definition: RSha256.hxx:104
int Int_t
Definition: RtypesCore.h:45
const Bool_t kTRUE
Definition: RtypesCore.h:91
#define gDirectory
Definition: TDirectory.h:290
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:187
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
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
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:1410
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
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,...
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:29
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:37
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:33
const TObject * getObj(Int_t uniqueID)
Retrieve payload object of cache element with given unique ID
A RooPlot is a plot frame and a container for graphics objects within that frame.
Definition: RooPlot.h:44
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:691
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
void Print(Option_t *opts=0) const
Print contents of the workspace.
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found.
RooExpensiveObjectCache & expensiveObjectCache()
Definition: RooWorkspace.h:165
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
The Canvas class.
Definition: TCanvas.h:23
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:54
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:3073
RooCmdArg Verbose(Bool_t flag=kTRUE)
RooCmdArg Timer(Bool_t flag=kTRUE)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
const char * Title
Definition: TXMLSetup.cxx:68