rf903_numintcache.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_roofit
/// \notebook -js
/// Numeric algorithm tuning: caching of slow numeric integrals and parameterization of slow numeric integrals
///
/// \macro_image
/// \macro_output
/// \macro_code
/// \author 07/2008 - Wouter Verkerke

#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooDataHist.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "RooWorkspace.h"
#include "RooExpensiveObjectCache.h"
#include "TFile.h"
#include "TH1.h"

using namespace RooFit;

RooWorkspace *getWorkspace(Int_t mode);

void rf903_numintcache(Int_t mode = 0)
{
   // Mode = 0 : Run plain fit (slow)
   // Mode = 1 : Generate workspace with pre-calculated integral and store it on file (prepare for accelerated running)
   // Mode = 2 : Run fit from previously stored workspace including cached integrals (fast, requires run in mode=1
   // first)

   // 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 .
   // -----------------------------------------------------------------------------------

   // Make/load workspace, exit here in mode 1
   RooWorkspace *w1 = getWorkspace(mode);
   if (mode == 1) {

      // Show workspace that was created
      w1->Print();

      // Show plot of cached integral values
      RooDataHist *hhcache = (RooDataHist *)w1->expensiveObjectCache().getObj(1);
      if (hhcache) {

         new TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600);
         hhcache->createHistogram("a")->Draw();

      } else {
         Error("rf903_numintcache", "Cached histogram is not existing in workspace");
      }
      return;
   }

   // 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
   // -----------------------------------------------------------------------------------

   // This is always slow (need to find maximum function value empirically in 3D space)
   RooDataSet *d = w1->pdf("model")->generate(RooArgSet(*w1->var("x"), *w1->var("y"), *w1->var("z")), 1000);

   // This is slow in mode 0, but fast in mode 1
   w1->pdf("model")->fitTo(*d, Verbose(kTRUE), Timer(kTRUE));

   // Projection on x (always slow as 2D integral over Y,Z at fitted value of a is not cached)
   RooPlot *framex = w1->var("x")->frame(Title("Projection of 3D model on X"));
   d->plotOn(framex);
   w1->pdf("model")->plotOn(framex);

   // Draw x projection on canvas
   new TCanvas("rf903_numintcache", "rf903_numintcache", 600, 600);
   framex->Draw();

   // Make workspace available on command line after macro finishes
   gDirectory->Add(w1);

   return;
}

RooWorkspace *getWorkspace(Int_t mode)
{
   // 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 .
   // -----------------------------------------------------------------------------------
   //
   // Mode = 0 : Create workspace for plain running (no integral caching)
   // Mode = 1 : Generate workspace with pre-calculated integral and store it on file
   // Mode = 2 : Load previously stored workspace from file

   RooWorkspace *w(0);

   if (mode != 2) {

      // Create empty workspace workspace
      w = new RooWorkspace("w", 1);

      // Make a difficult to normalize  p.d.f. in 3 dimensions that is integrated numerically.
      w->factory("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])");
   }

   if (mode == 1) {

      // Instruct model to pre-calculate normalization integral that integrate at least
      // two dimensions numerically. In this specific case the integral value for
      // all values of parameter 'a' are stored in a histogram and available for use
      // in subsequent fitting and plotting operations (interpolation is applied)

      // w->pdf("model")->setNormValueCaching(3) ;
      w->pdf("model")->setStringAttribute("CACHEPARMINT", "x:y:z");

      // Evaluate p.d.f. once to trigger filling of cache
      RooArgSet normSet(*w->var("x"), *w->var("y"), *w->var("z"));
      w->pdf("model")->getVal(&normSet);
      w->writeToFile("rf903_numintcache.root");
   }

   if (mode == 2) {
      // Load preexisting workspace from file in mode==2
      TFile *f = new TFile("rf903_numintcache.root");
      w = (RooWorkspace *)f->Get("w");
   }

   // Return created or loaded workspace
   return w;
}