rf903_numintcache.C File Reference

Detailed Description

'NUMERIC ALGORITHM TUNING' RooFit tutorial macro #903

Caching of slow numeric integrals and parameterization of slow numeric integrals

Processing /mnt/vdb/lsf/workspace/root-makedoc-v608/rootspi/rdoc/src/v6-08-00-patches/tutorials/roofit/rf903_numintcache.C...

RooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby 
                Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
                All rights reserved, please read http://roofit.sourceforge.net/license.txt

[#1] INFO:ObjectHandling -- RooWorkspace::exportToCint(w) INFO: references to all objects in this workspace will be created in CINT in 'namespace w'
[#1] INFO:NumericIntegration -- RooRealIntegral::init(model_Int[x,y,z]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(x,y,z)
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
[#0] WARNING:Minization -- RooMinimizerFcn::synchronize: WARNING: no initial error estimate available for a: using 1
 **********
 **    1 **SET PRINT           1
 **********
 **********
 **    2 **SET NOGRAD
 **********
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 a            0.00000e+00  1.00000e+00   -5.00000e+00  5.00000e+00
 **********
 **    3 **SET ERR         0.5
 **********
 **********
 **    4 **SET PRINT           1
 **********
 **********
 **    5 **SET STR           1
 **********
 NOW USING STRATEGY  1: TRY TO BALANCE SPEED AGAINST RELIABILITY
 **********
 **    6 **MIGRAD         500           1
 **********
 FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.

prevFCN = 1659.930708   START MIGRAD MINIMIZATION.  STRATEGY  1.  CONVERGENCE WHEN EDM .LT. 1.00e-03
a=0.02003, [#1] INFO:NumericIntegration -- RooRealIntegral::init(model_Int[x,y,z]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(x,y,z)

prevFCN = 1668.090831  a=-0.02003, 
prevFCN = 1666.341378  a=0.002003, 
prevFCN = 1660.094465  a=-0.002003, 
prevFCN = 1659.913535   FCN=1659.93 FROM MIGRAD    STATUS=INITIATE        4 CALLS           5 TOTAL
                     EDM= unknown      STRATEGY= 1      NO ERROR MATRIX       
  EXT PARAMETER               CURRENT GUESS       STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  a            0.00000e+00   1.00000e+00   2.01358e-01   2.25805e+02
                               ERR DEF= 0.5
a=-0.001236, 
prevFCN = 1659.902781  a=-0.001036, 
prevFCN = 1659.903514  a=-0.001437, 
prevFCN = 1659.903515  a=-0.001089, 
prevFCN = 1659.903177  a=-0.001383, 
prevFCN = 1659.903178   MIGRAD MINIMIZATION HAS CONVERGED.
 MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.
a=-0.001236, 
prevFCN = 1659.902781  a=-0.001089, 
prevFCN = 1659.903177  a=-0.001383, 
prevFCN = 1659.903178  a=-0.001207, 
prevFCN = 1659.902797  a=-0.001266, 
prevFCN = 1659.902797   COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=1659.9 FROM MIGRAD    STATUS=CONVERGED      14 CALLS          15 TOTAL
                     EDM=2.26676e-10    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                   STEP         FIRST   
  NO.   NAME      VALUE            ERROR          SIZE      DERIVATIVE 
   1  a           -1.23627e-03   5.23021e-03   2.94345e-05  -1.43931e-02
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  1    ERR DEF=0.5
  2.736e-05 
[#1] INFO:Minization -- Command timer: Real time 0:00:13, CP time 13.780
[#1] INFO:Minization -- Session timer: Real time 0:00:13, CP time 13.780
a=-0.001236,  **********
 **    7 **SET ERR         0.5
 **********
 **********
 **    8 **SET PRINT           1
 **********
 **********
 **    9 **HESSE         500
 **********

prevFCN = 1659.902781  a=-0.001207, 
prevFCN = 1659.902797  a=-0.001266, 
prevFCN = 1659.902797  a=-0.00123, 
prevFCN = 1659.902782  a=-0.001242, 
prevFCN = 1659.902782   COVARIANCE MATRIX CALCULATED SUCCESSFULLY
 FCN=1659.9 FROM HESSE     STATUS=OK              5 CALLS          20 TOTAL
                     EDM=2.26062e-10    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                                INTERNAL      INTERNAL  
  NO.   NAME      VALUE            ERROR       STEP SIZE       VALUE   
   1  a           -1.23627e-03   5.23023e-03   5.88690e-06  -2.47254e-04
                               ERR DEF= 0.5
 EXTERNAL ERROR MATRIX.    NDIM=  25    NPAR=  1    ERR DEF=0.5
  2.736e-05 
[#1] INFO:Minization -- Command timer: Real time 0:00:04, CP time 4.850
[#1] INFO:Minization -- Session timer: Real time 0:00:18, CP time 18.630, 2 slices
a=-0.001236, [#1] INFO:Minization -- RooMinimizer::optimizeConst: deactivating const optimization
[#0] WARNING:NumericIntegration -- RooAdaptiveIntegratorND::dtor(model) WARNING: Number of suppressed warningings about integral evaluations where target precision was not reached is 20
[#1] INFO:Plotting -- RooAbsReal::plotOn(model) plot on x integrates over variables (y,z)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(model_Int[y,z]_Norm[x,y,z]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(y,z)

#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 ;
}

Author

07/2008 - Wouter Verkerke

Definition in file rf903_numintcache.C.