#include <VerticalInterpPdf.h>

Classes
class	CacheElem
Public Member Functions
Double_t	analyticalIntegralWN (Int_t code, const RooArgSet normSet, const char rangeName=0) const
virtual Bool_t	checkObservables (const RooArgSet *nset) const
virtual TObject *	clone (const char *newname) const
const RooArgList &	coefList () const
Double_t	evaluate () const
virtual Bool_t	forceAnalyticalInt (const RooAbsArg &) const
const RooArgList &	funcList () const
Int_t	getAnalyticalIntegralWN (RooArgSet &allVars, RooArgSet &numVars, const RooArgSet normSet, const char rangeName=0) const
	VerticalInterpPdf (const char name, const char title, const RooArgList &funcList, const RooArgList &coefList, Double_t quadraticRegion=0., Int_t quadraticAlgo=0)
	VerticalInterpPdf (const VerticalInterpPdf &other, const char *name=0)
	VerticalInterpPdf ()
virtual	~VerticalInterpPdf ()
Protected Member Functions
Double_t	interpolate (Double_t coeff, Double_t central, RooAbsReal fUp, RooAbsReal fDown) const
	Iterator over coefficient list.
Protected Attributes
TIterator *	_coefIter
	Iterator over FUNC list.
RooListProxy	_coefList
TIterator *	_funcIter
RooListProxy	_funcList
RooObjCacheManager	_normIntMgr
Int_t	_quadraticAlgo
Double_t	_quadraticRegion

Detailed Description

Vertical interpolation between multiple histograms (or pdfs). Based on RooRealSumPdf

Definition at line 12 of file VerticalInterpPdf.h.

Constructor & Destructor Documentation

VerticalInterpPdf::VerticalInterpPdf ( )

Referenced by clone().

VerticalInterpPdf::VerticalInterpPdf	(	const char *	name,
		const char *	title,
		const RooArgList &	funcList,
		const RooArgList &	coefList,
		Double_t	quadraticRegion = `0.`,
		Int_t	quadraticAlgo = `0`
	)

Definition at line 33 of file VerticalInterpPdf.cc.

References _coefIter, _coefList, _funcIter, _funcList, and _quadraticAlgo.

                                                                                                                                                                                   :
  RooAbsPdf(name,title),
  _normIntMgr(this,10),
  _funcList("!funcList","List of functions",this),
  _coefList("!coefList","List of coefficients",this),
  _quadraticRegion(quadraticRegion),
  _quadraticAlgo(quadraticAlgo)
{ 

  if (inFuncList.getSize()!=2*inCoefList.getSize()+1) {
    coutE(InputArguments) << "VerticalInterpPdf::VerticalInterpPdf(" << GetName() 
                          << ") number of pdfs and coefficients inconsistent, must have Nfunc=1+2*Ncoef" << endl ;
    assert(0);
  }

  TIterator* funcIter = inFuncList.createIterator() ;
  RooAbsArg* func;
  while((func = (RooAbsArg*)funcIter->Next())) {
    if (!dynamic_cast<RooAbsReal*>(func)) {
      coutE(InputArguments) << "ERROR: VerticalInterpPdf::VerticalInterpPdf(" << GetName() << ") function  " << func->GetName() << " is not of type RooAbsReal" << endl;
      assert(0);
    }
    _funcList.add(*func) ;
  }
  delete funcIter;

  TIterator* coefIter = inCoefList.createIterator() ;
  RooAbsArg* coef;
  while((coef = (RooAbsArg*)coefIter->Next())) {
    if (!dynamic_cast<RooAbsReal*>(coef)) {
      coutE(InputArguments) << "ERROR: VerticalInterpPdf::VerticalInterpPdf(" << GetName() << ") coefficient " << coef->GetName() << " is not of type RooAbsReal" << endl;
      assert(0);
    }
    _coefList.add(*coef) ;    
  }
  delete coefIter;

  _funcIter  = _funcList.createIterator() ;
  _coefIter = _coefList.createIterator() ;

  if (_quadraticAlgo == -1) { 
    // multiplicative morphing: no way to do analytical integrals.
    _forceNumInt = kTRUE; 
  } else if (_quadraticAlgo >= 100) {
      _quadraticAlgo -= 100;
      _forceNumInt = kTRUE; 
  }
}

VerticalInterpPdf::VerticalInterpPdf	(	const VerticalInterpPdf &	other,
		const char *	name = `0`
	)

Definition at line 86 of file VerticalInterpPdf.cc.

References _coefIter, _coefList, _funcIter, and _funcList.

                                                                                     :
  RooAbsPdf(other,name),
  _normIntMgr(other._normIntMgr,this),
  _funcList("!funcList",this,other._funcList),
  _coefList("!coefList",this,other._coefList),
  _quadraticRegion(other._quadraticRegion),
  _quadraticAlgo(other._quadraticAlgo)
{
  // Copy constructor

  _funcIter  = _funcList.createIterator() ;
  _coefIter = _coefList.createIterator() ;
}

VerticalInterpPdf::~VerticalInterpPdf ( ) [virtual]

Definition at line 103 of file VerticalInterpPdf.cc.

References _coefIter, and _funcIter.

{
  // Destructor
  delete _funcIter ;
  delete _coefIter ;
}

Member Function Documentation

Double_t VerticalInterpPdf::analyticalIntegralWN	(	Int_t	code,
		const RooArgSet *	normSet,
		const char *	rangeName = `0`
	)		const

Definition at line 238 of file VerticalInterpPdf.cc.

References _coefIter, VerticalInterpPdf::CacheElem::_funcIntList, VerticalInterpPdf::CacheElem::_funcNormList, _normIntMgr, interpolate(), and relativeConstraints::value.

{
  // Implement analytical integrations by deferring integration of component
  // functions to integrators of components

  // Handle trivial passthrough scenario
  if (code==0) return getVal(normSet2) ;

  RooAbsReal *coef;
  Double_t value = 0;

  // WVE needs adaptation for rangeName feature
  CacheElem* cache = (CacheElem*) _normIntMgr.getObjByIndex(code-1) ;

  TIterator* funcIntIter = cache->_funcIntList.createIterator() ;
  RooAbsReal *funcInt = (RooAbsReal *) funcIntIter->Next();
  Double_t central = funcInt->getVal();
  value += central;

  _coefIter->Reset() ;
  while((coef=(RooAbsReal*)_coefIter->Next())) {
    Double_t coefVal = coef->getVal(normSet2) ;
    RooAbsReal * funcIntUp = (RooAbsReal*)funcIntIter->Next() ;
    RooAbsReal * funcIntDn = (RooAbsReal*)funcIntIter->Next() ;
    value += interpolate(coefVal, central, funcIntUp, funcIntDn);
  }
  
  delete funcIntIter ;
  
  Double_t normVal(1) ;
  if (normSet2) {
    normVal = 0 ;

    TIterator* funcNormIter = cache->_funcNormList.createIterator() ;

    RooAbsReal* funcNorm = (RooAbsReal*) funcNormIter->Next();
    central = funcNorm->getVal(normSet2) ;

    _coefIter->Reset() ;
    while((coef=(RooAbsReal*)_coefIter->Next())) {
      RooAbsReal *funcNormUp = (RooAbsReal*)funcNormIter->Next() ;
      RooAbsReal *funcNormDn = (RooAbsReal*)funcNormIter->Next() ;
      Double_t coefVal = coef->getVal(normSet2) ;
      normVal += interpolate(coefVal, central, funcNormUp, funcNormDn);
    }
    
    delete funcNormIter ;      
  }

  return ( value > 0 ? value : 1E-9 ) / normVal;
}

Bool_t VerticalInterpPdf::checkObservables ( const RooArgSet * nset ) const [virtual]

Definition at line 150 of file VerticalInterpPdf.cc.

References _coefIter, _funcIter, and _quadraticAlgo.

{
  // Check if FUNC is valid for given normalization set.
  // Coeffient and FUNC must be non-overlapping, but func-coefficient 
  // pairs may overlap each other
  //
  // In the present implementation, coefficients may not be observables or derive
  // from observables
  if (_quadraticAlgo == -1) return false; // multiplicative morphing. we don't care.
  

  _coefIter->Reset() ;
  RooAbsReal* coef ;
  while((coef=(RooAbsReal*)_coefIter->Next())) {
    if (coef->dependsOn(*nset)) {
      coutE(InputArguments) << "RooRealPdf::checkObservables(" << GetName() << "): ERROR coefficient " << coef->GetName() 
                            << " depends on one or more of the following observables" ; nset->Print("1") ;
      return true;
    }
  }

  _funcIter->Reset() ;
  RooAbsReal* func ;
  while((func = (RooAbsReal*)_funcIter->Next())) { 
    if (func->observableOverlaps(nset,*coef)) {
      coutE(InputArguments) << "VerticalInterpPdf::checkObservables(" << GetName() << "): ERROR: coefficient " << coef->GetName() 
                            << " and FUNC " << func->GetName() << " have one or more observables in common" << endl ;
      return true;
    }
  }
  
  return false;
}

virtual TObject* VerticalInterpPdf::clone ( const char * newname ) const [inline, virtual]

Definition at line 18 of file VerticalInterpPdf.h.

References VerticalInterpPdf().

{ return new VerticalInterpPdf(*this,newname) ; }

const RooArgList& VerticalInterpPdf::coefList ( ) const [inline]

Definition at line 29 of file VerticalInterpPdf.h.

References _coefList.

{ return _coefList ; }

Double_t VerticalInterpPdf::evaluate ( ) const

Definition at line 111 of file VerticalInterpPdf.cc.

References _coefIter, _funcIter, _quadraticAlgo, interpolate(), and relativeConstraints::value.

{
  // Calculate the current value
  Double_t value(0) ;

  // Do running sum of coef/func pairs, calculate lastCoef.
  _funcIter->Reset() ;
  _coefIter->Reset() ;
  RooAbsReal* coef ;
  RooAbsReal* func = (RooAbsReal*)_funcIter->Next();

  Double_t central = func->getVal();
  value = central;

  if (_quadraticAlgo >= 0) {
      // additive interpolation
      while((coef=(RooAbsReal*)_coefIter->Next())) {
          Double_t coefVal = coef->getVal() ;
          RooAbsReal* funcUp = (RooAbsReal*)_funcIter->Next() ;
          RooAbsReal* funcDn = (RooAbsReal*)_funcIter->Next() ;
          value += interpolate(coefVal, central, funcUp, funcDn);
      }
  } else {
      // multiplicative interpolation
      while((coef=(RooAbsReal*)_coefIter->Next())) {
          Double_t coefVal = coef->getVal() ;
          RooAbsReal* funcUp = (RooAbsReal*)_funcIter->Next() ;
          RooAbsReal* funcDn = (RooAbsReal*)_funcIter->Next() ;
          value *= interpolate(coefVal, central, funcUp, funcDn);
      }
  }
   
  return value > 0 ? value : 1E-9 ;
}

virtual Bool_t VerticalInterpPdf::forceAnalyticalInt ( const RooAbsArg & ) const [inline, virtual]

Definition at line 24 of file VerticalInterpPdf.h.

{ return kTRUE ; }

const RooArgList& VerticalInterpPdf::funcList ( ) const [inline]

Definition at line 28 of file VerticalInterpPdf.h.

References _funcList.

{ return _funcList ; }

Int_t VerticalInterpPdf::getAnalyticalIntegralWN	(	RooArgSet &	allVars,
		RooArgSet &	numVars,
		const RooArgSet *	normSet,
		const char *	rangeName = `0`
	)		const

Definition at line 188 of file VerticalInterpPdf.cc.

References VerticalInterpPdf::CacheElem::_funcIntList, _funcIter, VerticalInterpPdf::CacheElem::_funcNormList, and _normIntMgr.

{
  // Advertise that all integrals can be handled internally.

  // Handle trivial no-integration scenario
  if (allVars.getSize()==0) return 0 ;
  if (_forceNumInt) return 0 ;

  // Select subset of allVars that are actual dependents
  analVars.add(allVars) ;
  RooArgSet* normSet = normSet2 ? getObservables(normSet2) : 0 ;


  // Check if this configuration was created before
  Int_t sterileIdx(-1) ;
  CacheElem* cache = (CacheElem*) _normIntMgr.getObj(normSet,&analVars,&sterileIdx,(const char *)0);
  if (cache) {
    return _normIntMgr.lastIndex()+1 ;
  }
  
  // Create new cache element
  cache = new CacheElem ;

  // Make list of function projection and normalization integrals 
  _funcIter->Reset() ;
  RooAbsReal *func ;
  while((func=(RooAbsReal*)_funcIter->Next())) {
    RooAbsReal* funcInt = func->createIntegral(analVars) ;
    cache->_funcIntList.addOwned(*funcInt) ;
    if (normSet && normSet->getSize()>0) {
      RooAbsReal* funcNorm = func->createIntegral(*normSet) ;
      cache->_funcNormList.addOwned(*funcNorm) ;
    }
  }

  // Store cache element
  Int_t code = _normIntMgr.setObj(normSet,&analVars,(RooAbsCacheElement*)cache,0) ;

  if (normSet) {
    delete normSet ;
  }

  return code+1 ; 
}

Double_t VerticalInterpPdf::interpolate	(	Double_t	coeff,
		Double_t	central,
		RooAbsReal *	fUp,
		RooAbsReal *	fDown
	)		const `[protected]`

Iterator over coefficient list.

Definition at line 290 of file VerticalInterpPdf.cc.

References _quadraticAlgo, _quadraticRegion, and funct::pow().

Referenced by analyticalIntegralWN(), and evaluate().

{
    if (_quadraticAlgo == -1) {
        Double_t kappa = (coeff > 0 ? fUp->getVal()/central : central/fDn->getVal());
        return pow(kappa,fabs(coeff));
    }

    if (fabs(coeff) >= _quadraticRegion) {
        return coeff * (coeff > 0 ? fUp->getVal() - central : central - fDn->getVal());
    } else {
        // quadratic interpolation coefficients between the three
        if (_quadraticAlgo != 1) {
            // quadratic interpolation null in zero and continuos at boundaries, but not differentiable at boundaries
            // conditions:
            //   c_up (+_quadraticRegion) = +_quadraticRegion
            //   c_cen(+_quadraticRegion) = -_quadraticRegion
            //   c_dn (+_quadraticRegion) = 0
            //   c_up (-_quadraticRegion) = 0 
            //   c_cen(-_quadraticRegion) = -_quadraticRegion
            //   c_dn (-_quadraticRegion) = +_quadraticRegion
            //   c_up(0) = c_dn(0) = c_cen(0) = 0
            Double_t c_up  = + coeff * (_quadraticRegion + coeff) / (2 * _quadraticRegion);
            Double_t c_dn  = - coeff * (_quadraticRegion - coeff) / (2 * _quadraticRegion);
            Double_t c_cen = - coeff * coeff / _quadraticRegion;
            return c_up * fUp->getVal() + c_dn * fDn->getVal() + c_cen * central;
        } else { //if (_quadraticAlgo == 1) { 
            // quadratic interpolation that is everywhere differentiable, but it's not null in zero
            // conditions on the function
            //   c_up (+_quadraticRegion) = +_quadraticRegion
            //   c_cen(+_quadraticRegion) = -_quadraticRegion
            //   c_dn (+_quadraticRegion) = 0
            //   c_up (-_quadraticRegion) = 0 
            //   c_cen(-_quadraticRegion) = -_quadraticRegion
            //   c_dn (-_quadraticRegion) = +_quadraticRegion
            // conditions on the derivatives
            //   c_up '(+_quadraticRegion) = +1
            //   c_cen'(+_quadraticRegion) = -1
            //   c_dn '(+_quadraticRegion) = 0
            //   c_up '(-_quadraticRegion) = 0
            //   c_cen'(-_quadraticRegion) = +1
            //   c_dn '(-_quadraticRegion) = -1
            Double_t c_up  = (_quadraticRegion + coeff) * (_quadraticRegion + coeff) / (4 * _quadraticRegion);
            Double_t c_dn  = (_quadraticRegion - coeff) * (_quadraticRegion - coeff) / (4 * _quadraticRegion);
            Double_t c_cen = - c_up - c_dn;
            return c_up * fUp->getVal() + c_dn * fDn->getVal() + c_cen * central;
        } 
    }

}