RooAbsCachedPdf.cxx

Go to the documentation of this file.

 /***************************************************************************** 
  * Project: RooFit                                                           * 
  *                                                                           * 
  * Copyright (c) 2000-2005, Regents of the University of California          * 
  *                          and Stanford University. All rights reserved.    * 
  *                                                                           * 
  * Redistribution and use in source and binary forms,                        * 
  * with or without modification, are permitted according to the terms        * 
  * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             * 
  *****************************************************************************/ 

//////////////////////////////////////////////////////////////////////////////
// 
// BEGIN_HTML
// RooAbsCachedPdf is the abstract base class for p.d.f.s that need or
// want to cache their evaluate() output in a RooHistPdf defined in
// terms of the used observables. This base class manages the creation
// and storage of all RooHistPdf cache p.d.fs and the RooDataHists
// that define their shape. Implementations of RooAbsCachedPdf must
// define member function fillCacheObject() which serves to fill an
// already created RooDataHist with the p.d.fs function values. In
// addition the member functions actualObservables() and
// actualParameters() must be define which report what the actual
// observables to be cached are for a given set of observables passed
// by the user to getVal() and on which parameters need to be tracked
// for changes to trigger a refilling of the cache histogram.
// END_HTML
//
//
//
//

#include "Riostream.h" 
using namespace std ;

#include "RooFit.h"
#include "TString.h"
#include "RooAbsCachedPdf.h" 
#include "RooAbsReal.h" 
#include "RooMsgService.h"
#include "RooDataHist.h"
#include "RooHistPdf.h"
#include "RooGlobalFunc.h"
#include "RooRealVar.h"
#include "RooChangeTracker.h"
#include "RooExpensiveObjectCache.h"

ClassImp(RooAbsCachedPdf) 



////////////////////////////////////////////////////////////////////////////////
/// Constructor

RooAbsCachedPdf::RooAbsCachedPdf(const char *name, const char *title, Int_t ipOrder) :
  RooAbsPdf(name,title), 
  _cacheMgr(this,10),
  _ipOrder(ipOrder),
  _disableCache(kFALSE)
 { 
 } 



////////////////////////////////////////////////////////////////////////////////
/// Copy constructor

RooAbsCachedPdf::RooAbsCachedPdf(const RooAbsCachedPdf& other, const char* name) :  
   RooAbsPdf(other,name), 
   _cacheMgr(other._cacheMgr,this),
   _ipOrder(other._ipOrder),
   _disableCache(other._disableCache)
 { 
 } 



////////////////////////////////////////////////////////////////////////////////
/// Destructor

RooAbsCachedPdf::~RooAbsCachedPdf() 
{
}



////////////////////////////////////////////////////////////////////////////////
/// Implementation of getVal() overriding default implementation
/// of RooAbsPdf. Return normalized value stored in cache p.d.f
/// rather than return value of evaluate() which is undefined
/// for RooAbsCachedPdf

Double_t RooAbsCachedPdf::getValV(const RooArgSet* nset) const 
{
  if (_disableCache) {
    return RooAbsPdf::getValV(nset) ;
  }

  // Calculate current unnormalized value of object
  PdfCacheElem* cache = getCache(nset) ;

  Double_t value = cache->pdf()->getVal(nset) ;  

  _value = value ;    
  return _value ;
}



////////////////////////////////////////////////////////////////////////////////
/// Return pointer to RooHistPdf cache pdf for given choice of observables

RooAbsPdf* RooAbsCachedPdf::getCachePdf(const RooArgSet* nset) const 
{
  PdfCacheElem* cache = getCache(nset) ;

  if (cache) {
    return cache->pdf() ;
  } else {
    return 0 ;
  }
}


////////////////////////////////////////////////////////////////////////////////
/// Return pointer to RooDataHist cache histogram for given choice of observables

RooDataHist* RooAbsCachedPdf::getCacheHist(const RooArgSet* nset) const 
{
  PdfCacheElem* cache = getCache(nset) ;

  if (cache) {
    return cache->hist() ;
  } else {
    return 0 ;
  }
}


////////////////////////////////////////////////////////////////////////////////
/// Mark all bins of given cache as unitialized (value -1)

void RooAbsCachedPdf::clearCacheObject(PdfCacheElem& cache) const 
{
  cache.hist()->setAllWeights(-1) ;  
}



////////////////////////////////////////////////////////////////////////////////
/// Retrieve cache object associated with given choice of observables. If cache object
/// does not exist, create and fill and register it on the fly. If recalculate=false
/// recalculation of cache contents of existing caches that are marked dirty due to
/// dependent parameter changes is suppressed. 

RooAbsCachedPdf::PdfCacheElem* RooAbsCachedPdf::getCache(const RooArgSet* nset, Bool_t recalculate) const
{
  // Check if this configuration was created becfore
  Int_t sterileIdx(-1) ;
  PdfCacheElem* cache = (PdfCacheElem*) _cacheMgr.getObj(nset,0,&sterileIdx) ;

  // Check if we have a cache histogram in the global expensive object cache
  if (cache) {
    if (cache->paramTracker()->hasChanged(kTRUE) && (recalculate || !cache->pdf()->haveUnitNorm()) ) {
      cxcoutD(Eval) << "RooAbsCachedPdf::getCache(" << GetName() << ") cache " << cache << " pdf " 
          << cache->pdf()->GetName() << " requires recalculation as parameters changed" << endl ;     
      fillCacheObject(*cache) ;  
      cache->pdf()->setValueDirty() ;
    }
    return cache ;
  }

  // Create and fill cache
  cache = createCache(nset) ; 

  // Check if we have contents registered already in global expensive object cache 
  RooDataHist* htmp = (RooDataHist*) expensiveObjectCache().retrieveObject(cache->hist()->GetName(),RooDataHist::Class(),cache->paramTracker()->parameters()) ;

  if (htmp) {    

    cache->hist()->reset() ;
    cache->hist()->add(*htmp) ;

  } else {

    fillCacheObject(*cache) ;  

    RooDataHist* eoclone = new RooDataHist(*cache->hist()) ;
    eoclone->removeSelfFromDir() ;
    expensiveObjectCache().registerObject(GetName(),cache->hist()->GetName(),*eoclone,cache->paramTracker()->parameters()) ;
    
  } 

  
  // Store this cache configuration
  Int_t code = _cacheMgr.setObj(nset,0,((RooAbsCacheElement*)cache),0) ;

  coutI(Caching) << "RooAbsCachedPdf::getCache(" << GetName() << ") creating new cache " << cache << " with pdf "
       << cache->pdf()->GetName() << " for nset " << (nset?*nset:RooArgSet()) << " with code " << code ;
  if (htmp) {
    ccoutI(Caching) << " from preexisting content." ;
  }
  ccoutI(Caching) << endl ;
  
  return cache ;
}




////////////////////////////////////////////////////////////////////////////////
/// Constructor of cache object which owns RooDataHist cache histogram,
/// RooHistPdf pdf that represents is shape and RooChangeTracker meta
/// object that tracks changes in listed dependent parameter of cache.

RooAbsCachedPdf::PdfCacheElem::PdfCacheElem(const RooAbsCachedPdf& self, const RooArgSet* nsetIn) : 
  _pdf(0), _paramTracker(0), _hist(0), _norm(0) 
{
  // Create cache object itself -- Default implementation is a RooHistPdf
  RooArgSet* nset2 = self.actualObservables(nsetIn?*nsetIn:RooArgSet()) ;

  RooArgSet orderedObs ;
  if (nset2) {
    self.preferredObservableScanOrder(*nset2,orderedObs) ;
  }

  // Create RooDataHist
  TString hname = self.GetName() ;
  hname.Append("_") ;
  hname.Append(self.inputBaseName()) ;
  hname.Append("_CACHEHIST") ;
  hname.Append(self.cacheNameSuffix(orderedObs)) ;
  hname.Append(self.histNameSuffix()) ;
  _hist = new RooDataHist(hname,hname,orderedObs,self.binningName()) ;
  _hist->removeSelfFromDir() ;

  //RooArgSet* observables= self.getObservables(orderedObs) ;
  // cout << "orderedObs = " << orderedObs << " observables = " << *observables << endl ;

  // Get set of p.d.f. observable corresponding to set of histogram observables
  RooArgSet pdfObs ;
  RooArgSet pdfFinalObs ;
  TIterator* iter = orderedObs.createIterator() ;
  RooAbsArg* harg ;
  while((harg=(RooAbsArg*)iter->Next())) {
    RooAbsArg& po = self.pdfObservable(*harg) ;
    pdfObs.add(po) ;
    if (po.isFundamental()) {
      pdfFinalObs.add(po) ;
    } else {      
      RooArgSet* tmp = po.getVariables() ;
      pdfFinalObs.add(*tmp) ;
      delete tmp ;
    }
  }
  delete iter ;

  // Create RooHistPdf
  TString pdfname = self.inputBaseName() ;
  pdfname.Append("_CACHE") ;
  pdfname.Append(self.cacheNameSuffix(pdfFinalObs)) ;
  _pdf = new RooHistPdf(pdfname,pdfname,pdfObs,orderedObs,*_hist,self.getInterpolationOrder()) ;
  if (nsetIn) {
    _nset.addClone(*nsetIn) ;
  }

  // Create pseudo-object that tracks changes in parameter values

  RooArgSet* params = self.actualParameters(pdfFinalObs) ;
  params->remove(pdfFinalObs,kTRUE,kTRUE) ;
  
  string name= Form("%s_CACHEPARAMS",_pdf->GetName()) ;
  _paramTracker = new RooChangeTracker(name.c_str(),name.c_str(),*params,kTRUE) ;
  _paramTracker->hasChanged(kTRUE) ; // clear dirty flag as cache is up-to-date upon creation

  // Introduce formal dependency of RooHistPdf on parameters so that const optimization code
  // makes the correct decisions
  _pdf->addServerList(*params) ;

  // Set initial state of cache to dirty
  _pdf->setValueDirty() ;

  //delete observables ;
  delete params ;
  delete nset2 ;

}



////////////////////////////////////////////////////////////////////////////////
/// Construct string with unique suffix for cache objects based on 
/// observable names that define cache configuration

TString RooAbsCachedPdf::cacheNameSuffix(const RooArgSet& nset) const 
{
  TString name ;
  name.Append("_Obs[") ;
  if (nset.getSize()>0) {
    TIterator* iter = nset.createIterator() ;
    RooAbsArg* arg ;
    Bool_t first(kTRUE) ;
    while((arg=(RooAbsArg*)iter->Next())) {
      if (first) {
   first=kFALSE ;
      } else {
   name.Append(",") ;
      }
      name.Append(arg->GetName()) ;
    }
    delete iter ;
  } 

  name.Append("]") ;
  const char* payloadUS = payloadUniqueSuffix() ;
  if (payloadUS) {
    name.Append(payloadUS) ;
  }
  return name ;
}



////////////////////////////////////////////////////////////////////////////////
/// Change the interpolation order that is used in RooHistPdf cache
/// representation smoothing the RooDataHist shapes.

void RooAbsCachedPdf::setInterpolationOrder(Int_t order) 
{
  _ipOrder = order ;

  Int_t i ;
  for (i=0 ; i<_cacheMgr.cacheSize() ; i++) {
    PdfCacheElem* cache = (PdfCacheElem*) _cacheMgr.getObjByIndex(i) ;
    if (cache) {
      cache->pdf()->setInterpolationOrder(order) ;
    }
  }
}



////////////////////////////////////////////////////////////////////////////////
/// Returns all RooAbsArg objects contained in the cache element

RooArgList RooAbsCachedPdf::PdfCacheElem::containedArgs(Action) 
{
  RooArgList ret(*_pdf) ;
  ret.add(*_paramTracker) ;
  if (_norm) ret.add(*_norm) ;
  return ret ;
}



////////////////////////////////////////////////////////////////////////////////
/// Cache element destructor

RooAbsCachedPdf::PdfCacheElem::~PdfCacheElem() 
{
  if (_norm) {
    delete _norm ;
  }
  if (_pdf) {
    delete _pdf ;
  }
  if (_paramTracker) {
    delete _paramTracker ;
  }
  if (_hist) {
    delete _hist ;
  }
}



////////////////////////////////////////////////////////////////////////////////
/// Print contents of cache when printing self as part of object tree

void RooAbsCachedPdf::PdfCacheElem::printCompactTreeHook(ostream& os, const char* indent, Int_t curElem, Int_t maxElem) 
{
  if (curElem==0) {
    os << indent << "--- RooAbsCachedPdf begin cache ---" << endl ;
  }

  TString indent2(indent) ;
  os << Form("[%d] Configuration for observables ",curElem) << _nset << endl ;
  indent2 += Form("[%d] ",curElem) ;
  _pdf->printCompactTree(os,indent2) ;
  if (_norm) {
    os << Form("[%d] Norm ",curElem) ;
    _norm->printStream(os,kName|kArgs,kSingleLine) ;
  }
  
  if (curElem==maxElem) {
    os << indent << "--- RooAbsCachedPdf end cache --- " << endl ;
  }
}



////////////////////////////////////////////////////////////////////////////////
/// Force RooRealIntegral to offer all our actual observable for internal
/// integration

Bool_t RooAbsCachedPdf::forceAnalyticalInt(const RooAbsArg& dep) const 
{
  RooArgSet* actObs = actualObservables(dep) ;
  Bool_t ret = (actObs->getSize()>0) ;
  delete actObs ;
  return ret ;
}



////////////////////////////////////////////////////////////////////////////////
/// Advertises internal (analytical) integration capabilities. Call
/// is forwarded to RooHistPdf cache p.d.f of cache that is used for
/// given choice of observables

Int_t RooAbsCachedPdf::getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& analVars, const RooArgSet* normSet, const char* rangeName) const 
{
  if (allVars.getSize()==0) {
    return 0 ;
  }

  PdfCacheElem* cache = getCache(normSet?normSet:&allVars) ;
  Int_t code = cache->pdf()->getAnalyticalIntegralWN(allVars,analVars,normSet,rangeName) ;

  if (code==0) {
    return 0 ;
  }

  RooArgSet* all = new RooArgSet ;
  RooArgSet* ana = new RooArgSet ;
  RooArgSet* nrm = new RooArgSet ;
  all->addClone(allVars) ;
  ana->addClone(analVars) ;
  if (normSet) {
    nrm->addClone(*normSet) ;
  }
  std::vector<Int_t> codeList(2);
  codeList[0] = code ;
  codeList[1] = cache->pdf()->haveUnitNorm() ? 1 : 0 ;
  Int_t masterCode = _anaReg.store(codeList,all,ana,nrm)+1 ; // takes ownership of all sets

  
  // Mark all observables as internally integrated 
  if (cache->pdf()->haveUnitNorm()) {
    analVars.add(allVars,kTRUE) ;
  }

  return masterCode ;
}



////////////////////////////////////////////////////////////////////////////////
/// Implements internal (analytical) integration capabilities. Call
/// is forwarded to RooHistPdf cache p.d.f of cache that is used for
/// given choice of observables

Double_t RooAbsCachedPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
{  
  if (code==0) {
    return getVal(normSet) ; 
  }  

  RooArgSet *allVars(0),*anaVars(0),*normSet2(0),*dummy(0) ;
  const std::vector<Int_t> codeList = _anaReg.retrieve(code-1,allVars,anaVars,normSet2,dummy) ;
  
  PdfCacheElem* cache = getCache(normSet2?normSet2:anaVars,kFALSE) ;
  Double_t ret = cache->pdf()->analyticalIntegralWN(codeList[0],normSet,rangeName) ;

  if (codeList[1]>0) {
    RooArgSet factObs(*allVars) ;
    factObs.remove(*anaVars,kTRUE,kTRUE) ;
    TIterator* iter = factObs.createIterator() ;
    RooAbsLValue* arg ;
    while((arg=dynamic_cast<RooAbsLValue*>(iter->Next()))) {
      ret *= arg->volume(rangeName) ;
    }
    delete iter ;
  }
  
  return ret ;
}