RooIntegralMorph.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)             *
 *****************************************************************************/
 
/** \class RooIntegralMorph
    \ingroup Roofit
 
Class RooIntegralMorph is an implementation of the histogram interpolation
technique described by Alex Read in 'NIM A 425 (1999) 357-369 'Linear interpolation of histograms'
for continuous functions rather than histograms. The interpolation method, in short,
works as follows.
 
  - Given a p.d.f f1(x) with c.d.f F1(x) and p.d.f f2(x) with c.d.f F2(x)
 
  - One finds takes a value 'y' of both c.d.fs and determines the corresponding x
    values x(1,2) at which F(1,2)(x)==y.
 
  - The value of the interpolated p.d.f fbar(x) is then calculated as
    fbar(alpha*x1+(1-alpha)*x2) = f1(x1)*f2(x2) / ( alpha*f2(x2) + (1-alpha)*f1(x1) ) ;
 
From a technical point of view class RooIntegralMorph is a p.d.f that takes
two input p.d.fs f1(x,p) an f2(x,q) and an interpolation parameter to
make a p.d.f fbar(x,p,q,alpha). The shapes f1 and f2 are always taken
to be end the end-points of the parameter alpha, regardless of what
the those numeric values are.
 
Since the value of fbar(x) cannot be easily calculated for a given value
of x, class RooIntegralMorph is an implementation of RooAbsCachedPdf and
calculates the shape of the interpolated p.d.f. fbar(x) for all values
of x for a given value of alpha,p,q and caches these values in a histogram
(as implemented by RooAbsCachedPdf). The binning granularity of the cache
can be controlled by the binning named "cache" on the RooRealVar representing
the observable x. The fbar sampling algorithm first scans the range of
calculable c.d.f. values with a coarse grid and a recursive gap division
mechanism, and then polishes the c.d.f. value for each cache bin center with
Newton iterations on the monotone map X(y) = alpha*x1(y) + (1-alpha)*x2(y).
The residual inaccuracy of the cached p.d.f. is therefore dominated by the
interpolation between the cache bin centers and decreases steeply with the
number of cache bins: with O(1000) cache bins the difference between the
cached shape and the exact morphed p.d.f. is typically at the 1e-7 level
(measured as a Kolmogorov-Smirnov distance) for Gaussian input shapes.
 
Note on numeric stability of the algorithm. Since the algorithm relies
on a numeric inversion of cumulative distributions functions, some precision
may be lost at the 'edges' of the same (i.e. at regions in x where the
c.d.f. value is close to zero or one). The sampling strategy is to start
at y=0.1 (or the distance of y to 1.0) and push the y range outward by
a factor of sqrt(10) iteratively up to the point where the corresponding
x value no longer changes significantly, with a hard cutoff at y=1e-12.
For p.d.f.s with very flat tails such as Gaussians some part of the tail
may be lost due to limitations in numeric precision in the CDF inversion
step.
 
An effect related to the above limitation in numeric precision should
be anticipated when floating the alpha parameter in a fit. If a p.d.f
with such flat tails is fitted, it is likely that the dataset contains
events in the flat tail region. If the alpha parameter is varied, the
likelihood contribution from such events may exhibit discontinuities
in alpha, causing discontinuities in the summed likelihood as well
that will cause convergence problems in MINUIT. To mitigate this effect
one can use the setCacheAlpha() method to instruct RooIntegralMorph
to construct a two-dimensional cache for its output values in both
x and alpha. If linear interpolation is requested on the resulting
output histogram, the resulting interpolation of the p.d.f in the
alpha dimension will smooth out the discontinuities in the tail regions
result in a continuous likelihood distribution that can be fitted.
An added advantage of the cacheAlpha option is that if parameters
p,q of f1,f2 are fixed, the cached values in RooIntegralMorph are
valid for the entire fit session and do not need to be recalculated
for each change in alpha, which may result an considerable increase
in calculation speed.
 
**/
 
#include "Riostream.h"
 
#include "RooIntegralMorph.h"
#include "RooAbsCategory.h"
#include "RooBrentRootFinder.h"
#include "RooAbsFunc.h"
#include "RooRealVar.h"
#include "RooDataHist.h"
#include "TH1.h"
 
using std::flush, std::endl;
 
////////////////////////////////////////////////////////////////////////////////
/// Constructor with observables x, pdf shapes pdf1 and pdf2 which represent
/// the shapes at the end points of the interpolation parameter alpha
/// If doCacheAlpha is true, a two-dimensional cache is constructed in
/// both alpha and x
 
RooIntegralMorph::RooIntegralMorph(const char *name, const char *title, RooAbsReal &_pdf1, RooAbsReal &_pdf2,
                                   RooAbsReal &_x, RooAbsReal &_alpha, bool doCacheAlpha)
   : RooAbsCachedPdf(name, title, 2),
     pdf1("pdf1", "pdf1", this, _pdf1),
     pdf2("pdf2", "pdf2", this, _pdf2),
     x("x", "x", this, _x),
     alpha("alpha", "alpha", this, _alpha),
     _cacheAlpha(doCacheAlpha)
{
}
 
////////////////////////////////////////////////////////////////////////////////
/// Copy constructor
 
RooIntegralMorph::RooIntegralMorph(const RooIntegralMorph &other, const char *name)
   : RooAbsCachedPdf(other, name),
     pdf1("pdf1", this, other.pdf1),
     pdf2("pdf2", this, other.pdf2),
     x("x", this, other.x),
     alpha("alpha", this, other.alpha),
     _cacheAlpha(other._cacheAlpha)
{
}
 
////////////////////////////////////////////////////////////////////////////////
/// Observable to be cached for given choice of normalization.
/// Returns the 'x' observable unless doCacheAlpha is set in which
/// case a set with both x and alpha
 
RooFit::OwningPtr<RooArgSet> RooIntegralMorph::actualObservables(const RooArgSet & /*nset*/) const
{
   RooArgSet *obs = new RooArgSet;
   if (_cacheAlpha) {
      obs->add(alpha.arg());
   }
   obs->add(x.arg());
   return RooFit::OwningPtr<RooArgSet>{obs};
}
 
////////////////////////////////////////////////////////////////////////////////
/// Parameters of the cache. Returns parameters of both pdf1 and pdf2
/// and parameter cache, in case doCacheAlpha is not set.
 
RooFit::OwningPtr<RooArgSet> RooIntegralMorph::actualParameters(const RooArgSet & /*nset*/) const
{
   std::unique_ptr<RooArgSet> par1{pdf1->getParameters(static_cast<RooArgSet *>(nullptr))};
   RooArgSet par2;
   pdf2->getParameters(nullptr, par2);
   par1->add(par2, true);
   par1->remove(x.arg(), true, true);
   if (!_cacheAlpha) {
      par1->add(alpha.arg());
   }
   return RooFit::makeOwningPtr(std::move(par1));
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return base name component for cache components in this case
/// a string encoding the names of both end point p.d.f.s
 
const char *RooIntegralMorph::inputBaseName() const
{
   static TString name;
 
   name = pdf1.arg().GetName();
   name.Append("_MORPH_");
   name.Append(pdf2.arg().GetName());
   return name.Data();
}
 
////////////////////////////////////////////////////////////////////////////////
/// Fill the cache with the interpolated shape.
 
void RooIntegralMorph::fillCacheObject(PdfCacheElem &cache) const
{
   MorphCacheElem &mcache = static_cast<MorphCacheElem &>(cache);
 
   // If cacheAlpha is true employ slice iterator here to fill all slices
 
   if (!_cacheAlpha) {
 
      std::unique_ptr<TIterator> dIter{cache.hist()->sliceIterator(const_cast<RooAbsReal &>(x.arg()), RooArgSet())};
      mcache.calculate(dIter.get());
 
   } else {
      std::unique_ptr<TIterator> slIter{
         cache.hist()->sliceIterator(const_cast<RooAbsReal &>(alpha.arg()), RooArgSet())};
 
      double alphaSave = alpha;
      RooArgSet alphaSet(alpha.arg());
      coutP(Eval) << "RooIntegralMorph::fillCacheObject(" << GetName() << ") filling multi-dimensional cache";
      while (slIter->Next()) {
         alphaSet.assign(*cache.hist()->get());
         std::unique_ptr<TIterator> dIter{
            cache.hist()->sliceIterator(const_cast<RooAbsReal &>(x.arg()), RooArgSet(alpha.arg()))};
         mcache.calculate(dIter.get());
         ccoutP(Eval) << "." << flush;
      }
      ccoutP(Eval) << std::endl;
 
      const_cast<RooIntegralMorph *>(this)->alpha = alphaSave;
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// Create and return a derived MorphCacheElem.
 
RooAbsCachedPdf::PdfCacheElem *RooIntegralMorph::createCache(const RooArgSet *nset) const
{
   return new MorphCacheElem(const_cast<RooIntegralMorph &>(*this), nset);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return all RooAbsArg components contained in this cache
 
RooArgList RooIntegralMorph::MorphCacheElem::containedArgs(Action action)
{
   RooArgList ret;
   ret.add(PdfCacheElem::containedArgs(action));
   ret.add(*_self);
   ret.add(*_pdf1);
   ret.add(*_pdf2);
   ret.add(*_x);
   ret.add(*_alpha);
   ret.add(*_c1);
   ret.add(*_c2);
 
   return ret;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Construct of cache element, copy relevant input from RooIntegralMorph,
/// create the cdfs from the input p.d.fs and instantiate the root finders
/// on the cdfs to perform the inversion.
 
RooIntegralMorph::MorphCacheElem::MorphCacheElem(RooIntegralMorph &self, const RooArgSet *nsetIn)
   : PdfCacheElem(self, nsetIn),
     _self(&self),
     _pdf1(static_cast<RooAbsPdf *>(self.pdf1.absArg())),
     _pdf2(static_cast<RooAbsPdf *>(self.pdf2.absArg())),
     _x(static_cast<RooRealVar *>(self.x.absArg())),
     _alpha(static_cast<RooAbsReal *>(self.alpha.absArg())),
     _yatXmin(0),
     _yatXmax(0),
     _ccounter(0),
     _ycutoff(1e-12)
{
   // Mark in base class that normalization of cached pdf is invariant under pdf parameters
 
   _nset = std::make_unique<RooArgSet>(*_x);
 
   _c1 = std::unique_ptr<RooAbsReal>{_pdf1->createCdf(*_x)};
   _c2 = std::unique_ptr<RooAbsReal>{_pdf2->createCdf(*_x)};
   _cb1 = std::unique_ptr<RooAbsFunc>{_c1->bindVars(*_x, _nset.get())};
   _cb2 = std::unique_ptr<RooAbsFunc>{_c2->bindVars(*_x, _nset.get())};
 
   _rf1 = std::make_unique<RooBrentRootFinder>(*_cb1);
   _rf2 = std::make_unique<RooBrentRootFinder>(*_cb2);
 
   _rf1->setTol(1e-12);
   _rf2->setTol(1e-12);
 
   // _yatX = 0 ;
   // _calcX = 0 ;
 
   // Must do this here too: fillCache() may not be called if cache contents is retrieved from EOcache
   pdf()->setUnitNorm(true);
}
 
////////////////////////////////////////////////////////////////////////////////
/// Destructor
 
RooIntegralMorph::MorphCacheElem::~MorphCacheElem() {}
 
////////////////////////////////////////////////////////////////////////////////
/// Calculate the x value of the output p.d.f at the given cdf value y.
/// The ok boolean is filled with the success status of the operation.
 
double RooIntegralMorph::MorphCacheElem::calcX(double y, bool &ok)
{
   if (y < 0 || y > 1) {
      oocoutW(_self, Eval)
         << "RooIntegralMorph::MorphCacheElem::calcX() WARNING: requested root finding for unphysical CDF value " << y
         << std::endl;
   }
   double x1;
   double x2;
 
   double xmax = _x->getMax("cache");
   double xmin = _x->getMin("cache");
 
   ok = true;
   ok &= _rf1->findRoot(x1, xmin, xmax, y);
   ok &= _rf2->findRoot(x2, xmin, xmax, y);
   if (!ok)
      return 0;
   _ccounter++;
 
   return _alpha->getVal() * x1 + (1 - _alpha->getVal()) * x2;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Return the bin number enclosing the given x value
 
Int_t RooIntegralMorph::MorphCacheElem::binX(double X)
{
   double xmax = _x->getMax("cache");
   double xmin = _x->getMin("cache");
   return (Int_t)(_x->numBins("cache") * (X - xmin) / (xmax - xmin));
}
 
////////////////////////////////////////////////////////////////////////////////
/// Calculate shape of p.d.f for x,alpha values
/// defined by dIter iterator over cache histogram
 
void RooIntegralMorph::MorphCacheElem::calculate(TIterator *dIter)
{
   double xsave = _self->x;
 
   // if (!_yatX) {
   //   _yatX = new double[_x->numBins("cache")+1] ;
   //   _calcX = new double[_x->numBins("cache")+1] ;
   // }
 
   _yatX.resize(_x->numBins("cache") + 1);
   _calcX.resize(_x->numBins("cache") + 1);
 
   _ccounter = 0;
 
   // Get number of bins from PdfCacheElem histogram
   Int_t nbins = _x->numBins("cache");
   if (nbins < 2) {
      oocoutE(_self, Eval) << "RooIntegralMorph::MorphCacheElem::calculate(" << _self->GetName()
                           << ") ERROR: observable " << _x->GetName()
                           << " has an empty binning for the cache histogram."
                           << " Define one with RooRealVar::setBins(nbins, \"cache\")." << std::endl;
      return;
   }
 
   // Initialize yatX array to 'un-calculated values (-1)'
   for (int i = 0; i < nbins; i++) {
      _yatX[i] = -1;
      _calcX[i] = 0;
   }
 
   // Find low and high point
   findRange();
 
   // Perform initial scan of 100 points
   for (int i = 0; i < 10; i++) {
 
      // Take a point in y
      double offset = _yatX[_yatXmin];
      double delta = (_yatX[_yatXmax] - _yatX[_yatXmin]) / 10;
      double y = offset + i * delta;
 
      // Calculate corresponding X
      bool ok;
      double X = calcX(y, ok);
      if (ok) {
         Int_t iX = binX(X);
         _yatX[iX] = y;
         _calcX[iX] = X;
      }
   }
 
   // Now take an iteration filling the 'gaps'
   Int_t igapLow = _yatXmin + 1;
   while (true) {
      // Find next gap
      Int_t igapHigh = igapLow + 1;
      while (igapHigh < (_yatXmax) && _yatX[igapHigh] < 0)
         igapHigh++;
 
      // Fill the gap (iteratively and/or using interpolation)
      fillGap(igapLow - 1, igapHigh);
 
      // Terminate after processing of last gap
      if (igapHigh >= _yatXmax - 1)
         break;
      igapLow = igapHigh + 1;
   }
 
   // Make one more iteration to recalculate Y value at bin centers
   double xmax = _x->getMax("cache");
   double xmin = _x->getMin("cache");
   double binw = (xmax - xmin) / _x->numBins("cache");
   for (int i = _yatXmin + 1; i < _yatXmax - 1; i++) {
 
      // Calculate additional offset to apply if bin ixlo does not have X value calculated at bin center
      double xBinC = xmin + (i + 0.5) * binw;
      double xOffset = xBinC - _calcX[i];
      if (std::abs(xOffset / binw) > 1e-3) {
         double slope = (_yatX[i + 1] - _yatX[i - 1]) / (_calcX[i + 1] - _calcX[i - 1]);
         double newY = _yatX[i] + slope * xOffset;
         // cout << "bin " << i << " needs to be re-centered " << xOffset/binw << " slope = " << slope << " origY = " <<
         // _yatX[i] << " newY = " << newY << std::endl ;
         _yatX[i] = newY;
      }
   }
 
   // Zero output histogram below lowest calculable X value
   for (int i = 0; i < _yatXmin; i++) {
      dIter->Next();
      const std::size_t binIdx = hist()->getIndex(*hist()->get(), /*fast=*/true);
      hist()->set(binIdx, 0, -1);
   }
 
   double xMax = _x->getMax("cache");
   const double aval = _alpha->getVal();
 
   // Lower bounds for the root finding in the loop below, exploiting the fact
   // that the cumulative distribution functions increase monotonically: as y
   // increases from bin to bin, the x values found in the previous bin are
   // valid lower bounds for the current bin.
   double x1lo = _x->getMin("cache");
   double x2lo = _x->getMin("cache");
 
   // Transfer calculated values to histogram
   for (int i = _yatXmin; i <= _yatXmax; i++) {
 
      double y = _yatX[i];
      const double xBinC = xmin + (i + 0.5) * binw;
 
      double x1 = x1lo;
      double x2 = x2lo;
      double f1x1 = 0;
      double f2x2 = 0;
 
      // The y values obtained from the recursive gap filling are only
      // approximate solutions of X(y) == xBinC. Polish them with Newton
      // iterations on the monotone map X(y) = alpha*x1(y) + (1-alpha)*x2(y),
      // whose derivative is dX/dy = alpha/f1(x1) + (1-alpha)/f2(x2), so that
      // the stored p.d.f. value corresponds to the bin center to full precision.
      for (int iter = 0; iter < 10; iter++) {
         bool ok = _rf1->findRoot(x1, x1lo, xMax, y);
         ok &= _rf2->findRoot(x2, x2lo, xMax, y);
         _x->setVal(x1);
         f1x1 = _pdf1->getVal(_nset.get());
         _x->setVal(x2);
         f2x2 = _pdf2->getVal(_nset.get());
         if (!ok || f1x1 <= 0 || f2x2 <= 0)
            break;
         const double X = aval * x1 + (1 - aval) * x2;
         if (std::abs(X - xBinC) < 1e-12 * (xMax - xmin))
            break;
         const double dXdy = aval / f1x1 + (1 - aval) / f2x2;
         const double yNew = y + (xBinC - X) / dXdy;
         if (!(yNew > 0.) || !(yNew < 1.) || yNew == y)
            break;
         y = yNew;
      }
      _yatX[i] = y;
 
      double fbarX = f1x1 * f2x2 / (aval * f2x2 + (1 - aval) * f1x1);
 
      dIter->Next();
      {
         const std::size_t binIdx = hist()->getIndex(*hist()->get(), /*fast=*/true);
         hist()->set(binIdx, fbarX, -1);
      }
 
      x1lo = x1;
      x2lo = x2;
   }
   // Zero output histogram above highest calculable X value
   for (int i = _yatXmax + 1; i < nbins; i++) {
      dIter->Next();
      const std::size_t binIdx = hist()->getIndex(*hist()->get(), /*fast=*/true);
      hist()->set(binIdx, 0, -1);
   }
 
   pdf()->setUnitNorm(true);
   _self->x = xsave;
 
   oocxcoutD(_self, Eval) << "RooIntegralMorph::MorphCacheElem::calculate(" << _self->GetName()
                          << ") calculation required " << _ccounter << " samplings of cdfs" << std::endl;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Fill all empty histogram bins between bins ixlo and ixhi. The value of 'splitPoint'
/// defines the split point for the recursive division strategy to fill the gaps
/// If the midpoint value of y is very close to the midpoint in x, use interpolation
/// to fill the gaps, otherwise the intervals again.
 
void RooIntegralMorph::MorphCacheElem::fillGap(Int_t ixlo, Int_t ixhi, double splitPoint)
{
   // CONVENTION: _yatX[ixlo] is filled, _yatX[ixhi] is filled, elements in between are empty
   //   std::cout << "fillGap: gap from _yatX[" << ixlo << "]=" << _yatX[ixlo] << " to _yatX[" << ixhi << "]=" <<
   //   _yatX[ixhi] << ", size = " << ixhi-ixlo << std::endl ;
 
   if (_yatX[ixlo] < 0) {
      oocoutE(_self, Eval) << "RooIntegralMorph::MorphCacheElme::fillGap(" << _self->GetName() << "): ERROR in fillgap "
                           << ixlo << " = " << ixhi << " splitPoint= " << splitPoint << " _yatX[ixlo] = " << _yatX[ixlo]
                           << std::endl;
   }
   if (_yatX[ixhi] < 0) {
      oocoutE(_self, Eval) << "RooIntegralMorph::MorphCacheElme::fillGap(" << _self->GetName() << "): ERROR in fillgap "
                           << ixlo << " = " << ixhi << " splitPoint " << splitPoint << " _yatX[ixhi] = " << _yatX[ixhi]
                           << std::endl;
   }
 
   // Determine where half-way Y value lands
   double ymid = _yatX[ixlo] * splitPoint + _yatX[ixhi] * (1 - splitPoint);
   bool ok;
   double Xmid = calcX(ymid, ok);
   if (!ok) {
      oocoutW(_self, Eval) << "RooIntegralMorph::MorphCacheElem::fillGap(" << _self->GetName()
                           << ") unable to calculate midpoint in gap [" << ixlo << "," << ixhi
                           << "], resorting to interpolation" << std::endl;
      interpolateGap(ixlo, ixhi);
   }
 
   Int_t iX = binX(Xmid);
   double cq = (Xmid - _calcX[ixlo]) / (_calcX[ixhi] - _calcX[ixlo]) - 0.5;
 
   // Store midway point
   _yatX[iX] = ymid;
   _calcX[iX] = Xmid;
 
   // Policy: If centration quality is better than 1% OR better than 1/10 of a bin, fill interval with linear
   // interpolation
   if (std::abs(cq) < 0.01 || std::abs(cq * (ixhi - ixlo)) < 0.1 || ymid < _ycutoff) {
 
      // Fill remaining gaps on either side with linear interpolation
      if (iX - ixlo > 1) {
         interpolateGap(ixlo, iX);
      }
      if (ixhi - iX > 1) {
         interpolateGap(iX, ixhi);
      }
 
   } else {
 
      if (iX == ixlo) {
 
         if (splitPoint < 0.95) {
            // Midway value lands on lowest bin, retry split with higher split point
            double newSplit = splitPoint + 0.5 * (1 - splitPoint);
            fillGap(ixlo, ixhi, newSplit);
         } else {
            // Give up and resort to interpolation
            interpolateGap(ixlo, ixhi);
         }
 
      } else if (iX == ixhi) {
 
         // Midway value lands on highest bin, retry split with lower split point
         if (splitPoint > 0.05) {
            double newSplit = splitPoint / 2;
            fillGap(ixlo, ixhi, newSplit);
         } else {
            // Give up and resort to interpolation
            interpolateGap(ixlo, ixhi);
         }
 
      } else {
 
         // Midway point reasonable, iterate on interval on both sides
         if (iX - ixlo > 1) {
            fillGap(ixlo, iX);
         }
         if (ixhi - iX > 1) {
            fillGap(iX, ixhi);
         }
      }
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// Fill empty histogram bins between ixlo and ixhi with values obtained
/// from linear interpolation of ixlo,ixhi elements.
 
void RooIntegralMorph::MorphCacheElem::interpolateGap(Int_t ixlo, Int_t ixhi)
{
   // cout << "filling gap with linear interpolation ixlo=" << ixlo << " ixhi=" << ixhi << std::endl ;
 
   double xmax = _x->getMax("cache");
   double xmin = _x->getMin("cache");
   double binw = (xmax - xmin) / _x->numBins("cache");
 
   // Calculate deltaY in terms of actual X difference calculate, not based on nominal bin width
   double deltaY = (_yatX[ixhi] - _yatX[ixlo]) / ((_calcX[ixhi] - _calcX[ixlo]) / binw);
 
   // Calculate additional offset to apply if bin ixlo does not have X value calculated at bin center
   double xBinC = xmin + (ixlo + 0.5) * binw;
   double xOffset = xBinC - _calcX[ixlo];
 
   for (int j = ixlo + 1; j < ixhi; j++) {
      _yatX[j] = _yatX[ixlo] + (xOffset / binw + (j - ixlo)) * deltaY;
      _calcX[j] = xmin + (j + 0.5) * binw;
   }
}
 
////////////////////////////////////////////////////////////////////////////////
/// Determine which range of y values can be mapped to x values
/// from the numeric inversion of the input c.d.fs.
/// Start with a y range of [0.1-0.9] and push boundaries
/// outward with a factor of 1/sqrt(10). Stop iteration if
/// inverted x values no longer change
 
void RooIntegralMorph::MorphCacheElem::findRange()
{
   double xmin = _x->getMin("cache");
   double xmax = _x->getMax("cache");
   Int_t nbins = _x->numBins("cache");
 
   double x1;
   double x2;
   bool ok = true;
   double ymin = 0.1;
   double yminSave(-1);
   double Xsave(-1);
   double Xlast = xmax;
 
   // Find lowest Y value that can be measured
   // Start at 0.1 and iteratively lower limit by sqrt(10)
   while (true) {
      ok &= _rf1->findRoot(x1, xmin, xmax, ymin);
      ok &= _rf2->findRoot(x2, xmin, xmax, ymin);
      oocxcoutD(_self, Eval) << "RooIntegralMorph::MorphCacheElem::findRange(" << _self->GetName()
                             << ") findMin: x1 = " << x1 << " x2 = " << x2 << " ok = " << (ok ? "T" : "F") << std::endl;
 
      // Terminate in case of non-convergence
      if (!ok)
         break;
 
      // Terminate if value of X no longer moves by >0.1 bin size
      double X = _alpha->getVal() * x1 + (1 - _alpha->getVal()) * x2;
      if (std::abs(X - Xlast) / (xmax - xmin) < 0.0001) {
         break;
      }
      Xlast = X;
 
      // Store new Y value
      _yatXmin = (Int_t)(nbins * (X - xmin) / (xmax - xmin));
      _yatX[_yatXmin] = ymin;
      _calcX[_yatXmin] = X;
      yminSave = ymin;
      Xsave = X;
 
      // Reduce ymin by half an order of magnitude
      ymin /= sqrt(10.);
 
      // Emergency break
      if (ymin < _ycutoff)
         break;
   }
   _yatX[_yatXmin] = yminSave;
   _calcX[_yatXmin] = Xsave;
 
   // Find highest Y value that can be measured
   // Start at 1 - 0.1 and iteratively lower delta by sqrt(10)
   ok = true;
   double deltaymax = 0.1;
   double deltaymaxSave(-1);
   Xlast = xmin;
   while (true) {
      ok &= _rf1->findRoot(x1, xmin, xmax, 1 - deltaymax);
      ok &= _rf2->findRoot(x2, xmin, xmax, 1 - deltaymax);
 
      oocxcoutD(_self, Eval) << "RooIntegralMorph::MorphCacheElem::findRange(" << _self->GetName()
                             << ") findMax: x1 = " << x1 << " x2 = " << x2 << " ok = " << (ok ? "T" : "F") << std::endl;
 
      // Terminate in case of non-convergence
      if (!ok)
         break;
 
      // Terminate if value of X no longer moves by >0.1 bin size
      double X = _alpha->getVal() * x1 + (1 - _alpha->getVal()) * x2;
      if (std::abs(X - Xlast) / (xmax - xmin) < 0.0001) {
         break;
      }
      Xlast = X;
 
      // Store new Y value
      _yatXmax = (Int_t)(nbins * (X - xmin) / (xmax - xmin));
      _yatX[_yatXmax] = 1 - deltaymax;
      _calcX[_yatXmax] = X;
      deltaymaxSave = deltaymax;
      Xsave = X;
 
      // Reduce ymin by half an order of magnitude
      deltaymax /= sqrt(10.);
 
      // Emergency break
      if (deltaymax < _ycutoff)
         break;
   }
 
   _yatX[_yatXmax] = 1 - deltaymaxSave;
   _calcX[_yatXmax] = Xsave;
 
   // Initialize values out of range to 'out-of-range' (-2)
   for (int i = 0; i < _yatXmin; i++)
      _yatX[i] = -2;
   for (int i = _yatXmax + 1; i < nbins; i++)
      _yatX[i] = -2;
   oocxcoutD(_self, Eval) << "RooIntegralMorph::findRange(" << _self->GetName() << "): ymin = " << _yatX[_yatXmin]
                          << " ymax = " << _yatX[_yatXmax] << std::endl;
   oocxcoutD(_self, Eval) << "RooIntegralMorph::findRange(" << _self->GetName() << "): xmin = " << _calcX[_yatXmin]
                          << " xmax = " << _calcX[_yatXmax] << std::endl;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Dummy
 
double RooIntegralMorph::evaluate() const
{
   return 0;
}
 
////////////////////////////////////////////////////////////////////////////////
/// Indicate to the RooAbsCachedPdf base class that for the filling of the
/// cache the traversal of the x should be in the innermost loop, to minimize
/// recalculation of the one-dimensional internal cache for a fixed value of alpha
 
void RooIntegralMorph::preferredObservableScanOrder(const RooArgSet &obs, RooArgSet &orderedObs) const
{
   // Put x last to minimize cache faulting
   orderedObs.removeAll();
 
   orderedObs.add(obs);
   RooAbsArg *obsX = obs.find(x.arg().GetName());
   if (obsX) {
      orderedObs.remove(*obsX);
      orderedObs.add(*obsX);
   }
}