[RF] Fix implementation of RooGenericPdf/RooFormulaVar::binBoundaries()

roofit/roofitcore/inc/RooHelpers.h

@@ -17,17 +17,13 @@
 #include <RooAbsArg.h>
 #include <RooAbsReal.h>
 
-#include <ROOT/RSpan.hxx>
-
 #include <sstream>
-#include <list>
 #include <vector>
 #include <string>
 #include <utility>
 
 class RooAbsPdf;
 class RooAbsData;
-class RooAbsRealLValue;
 
 namespace RooHelpers {
 
@@ -95,22 +91,6 @@ void checkRangeOfParameters(const RooAbsReal *callingClass, std::initializer_lis
 /// set all RooRealVars to constants. return true if at least one changed status
 bool setAllConstant(const RooAbsCollection &coll, bool constant = true);
 
-/// Check that `function` is constant (flat) inside each bin defined by the
-/// sorted `boundaries` when scanning the observable `obs`. Several interior
-/// points are sampled per bin and compared to the bin's first sample; if any
-/// of them deviates by more than `relTol` (relative to the value scale), the
-/// function is not flat and false is returned. The value of `obs` is restored
-/// on return.
-bool isFunctionFlatInBins(const RooAbsReal &function, RooAbsRealLValue &obs, std::span<const double> boundaries,
-                          double relTol = 1e-9);
-
-/// Return a newly allocated list with the subset of `boundaries` that lies
-/// strictly inside [`xlo`, `xhi`], with `xlo` and `xhi` added as the first and
-/// last entries. This is the form expected by RooFit's binBoundaries()
-/// interface, so the bin integrator covers exactly the integration range.
-/// The caller takes ownership of the returned list.
-std::list<double> *binBoundariesInRange(std::span<const double> boundaries, double xlo, double xhi);
-
 } // namespace RooHelpers
 
 #endif

roofit/roofitcore/res/RooFitImplHelpers.h

@@ -127,6 +127,15 @@ BinnedLOutput getBinnedL(RooAbsPdf const &pdf);
 
 void getSortedComputationGraph(RooAbsArg const &func, RooArgSet &out);
 
+/// Check that `function` is constant (flat) inside each bin defined by the
+/// sorted `boundaries` when scanning the observable `obs`. Several interior
+/// points are sampled per bin and compared to the bin's first sample; if any
+/// of them deviates by more than `relTol` (relative to the value scale), the
+/// function is not flat and false is returned. The value of `obs` is restored
+/// on return.
+bool isFunctionFlatInBins(const RooAbsReal &function, RooAbsRealLValue &obs, std::span<const double> boundaries,
+                          double relTol = 1e-9);
+
 } // namespace RooHelpers
 
 namespace RooFit::Detail {

roofit/roofitcore/src/RooFitImplHelpers.cxx

@@ -286,6 +286,42 @@ RooAbsArg *cloneTreeWithSameParametersImpl(RooAbsArg const &arg, RooArgSet const
 
 } // namespace Detail
 
+bool isFunctionFlatInBins(const RooAbsReal &function, RooAbsRealLValue &obs, std::span<const double> boundaries,
+                          double relTol)
+{
+   // Fractions of the bin width at which the function is sampled. They are kept
+   // strictly inside the bin (away from the boundaries) so that the evaluation
+   // is not affected by which side of a boundary a step function jumps.
+   const double fractions[] = {0.04, 0.27, 0.5, 0.73, 0.96};
+
+   const double savedVal = obs.getVal();
+
+   bool isFlat = true;
+   for (std::size_t i = 0; i + 1 < boundaries.size() && isFlat; ++i) {
+      const double lo = boundaries[i];
+      const double hi = boundaries[i + 1];
+      double reference = 0.0;
+      bool first = true;
+      for (double frac : fractions) {
+         obs.setVal(lo + frac * (hi - lo));
+         const double val = function.getVal();
+         if (first) {
+            reference = val;
+            first = false;
+            continue;
+         }
+         const double scale = std::max(std::abs(reference), 1e-12);
+         if (std::abs(val - reference) > relTol * scale) {
+            isFlat = false;
+            break;
+         }
+      }
+   }
+
+   obs.setVal(savedVal);
+   return isFlat;
+}
+
 } // namespace RooHelpers
 
 namespace RooFit::Detail {

roofit/roofitcore/src/RooFormulaVar.cxx

@@ -55,7 +55,7 @@
 #include "RooAbsRealLValue.h"
 #include "RooAbsBinning.h"
 #include "RooCurve.h"
-#include "RooHelpers.h"
+#include "RooFitImplHelpers.h"
 
 #ifdef ROOFIT_LEGACY_EVAL_BACKEND
 #include "RooNLLVar.h"
@@ -319,8 +319,14 @@ std::list<double>* RooFormulaVar::binBoundaries(RooAbsRealLValue& obs, double xl
    auto found = _binnings.find(_actualVars.index(obs.GetName()));
    if (found != _binnings.end()) {
       const RooAbsBinning &binning = *found->second;
-      return RooHelpers::binBoundariesInRange({binning.array(), static_cast<std::size_t>(binning.numBoundaries())}, xlo,
-                                              xhi);
+      auto hint = new std::list<double>;
+      for (int i = 0; i < binning.numBoundaries(); ++i) {
+         const double boundary = binning.array()[i];
+         if (boundary >= xlo && boundary <= xhi) {
+            hint->push_back(boundary);
+         }
+      }
+      return hint;
    }
 
   for (const auto par : _actualVars) {

roofit/roofitcore/src/RooGenericPdf.cxx

@@ -51,7 +51,7 @@ the names of the arguments are not hard coded.
 #include "RooAbsRealLValue.h"
 #include "RooAbsBinning.h"
 #include "RooCurve.h"
-#include "RooHelpers.h"
+#include "RooFitImplHelpers.h"
 
 using std::istream, std::ostream, std::endl;
 
@@ -223,8 +223,14 @@ std::list<double> *RooGenericPdf::binBoundaries(RooAbsRealLValue &obs, double xl
       return nullptr;
    }
    const RooAbsBinning &binning = *found->second;
-   return RooHelpers::binBoundariesInRange({binning.array(), static_cast<std::size_t>(binning.numBoundaries())}, xlo,
-                                           xhi);
+   auto hint = new std::list<double>;
+   for (int i = 0; i < binning.numBoundaries(); ++i) {
+      const double boundary = binning.array()[i];
+      if (boundary >= xlo && boundary <= xhi) {
+         hint->push_back(boundary);
+      }
+   }
+   return hint;
 }
 
 ////////////////////////////////////////////////////////////////////////////////

roofit/roofitcore/src/RooHelpers.cxx

@@ -29,8 +29,6 @@
 #include <ROOT/StringUtils.hxx>
 #include <TClass.h>
 
-#include <algorithm>
-#include <cmath>
 #include <unordered_map>
 
 namespace RooHelpers {
@@ -170,60 +168,4 @@ bool setAllConstant(const RooAbsCollection &coll, bool constant)
    return changed;
 }
 
-bool isFunctionFlatInBins(const RooAbsReal &function, RooAbsRealLValue &obs, std::span<const double> boundaries,
-                          double relTol)
-{
-   // Fractions of the bin width at which the function is sampled. They are kept
-   // strictly inside the bin (away from the boundaries) so that the evaluation
-   // is not affected by which side of a boundary a step function jumps.
-   const double fractions[] = {0.04, 0.27, 0.5, 0.73, 0.96};
-
-   const double savedVal = obs.getVal();
-
-   bool isFlat = true;
-   for (std::size_t i = 0; i + 1 < boundaries.size() && isFlat; ++i) {
-      const double lo = boundaries[i];
-      const double hi = boundaries[i + 1];
-      double reference = 0.0;
-      bool first = true;
-      for (double frac : fractions) {
-         obs.setVal(lo + frac * (hi - lo));
-         const double val = function.getVal();
-         if (first) {
-            reference = val;
-            first = false;
-            continue;
-         }
-         const double scale = std::max(std::abs(reference), 1e-12);
-         if (std::abs(val - reference) > relTol * scale) {
-            isFlat = false;
-            break;
-         }
-      }
-   }
-
-   obs.setVal(savedVal);
-   return isFlat;
-}
-
-std::list<double> *binBoundariesInRange(std::span<const double> boundaries, double xlo, double xhi)
-{
-   auto out = new std::list<double>;
-
-   // Small tolerance so that boundaries numerically coinciding with the range
-   // limits are not duplicated by the explicit xlo/xhi endpoints below.
-   const double delta = (xhi - xlo) * 1e-8;
-
-   for (double boundary : boundaries) {
-      if (boundary > xlo + delta && boundary < xhi - delta) {
-         out->push_back(boundary);
-      }
-   }
-
-   out->push_front(xlo);
-   out->push_back(xhi);
-
-   return out;
-}
-
 } // namespace RooHelpers

roofit/roofitcore/test/testGenericPdf.cxx

@@ -2,13 +2,15 @@
 // Authors: Stephan Hageboeck, CERN  05/2019
 //          Jonas Rembser, CERN 06/2022
 
-#include <RooRealVar.h>
-#include <RooGenericPdf.h>
-#include <RooProduct.h>
 #include <RooArgList.h>
 #include <RooBinning.h>
-#include <RooUniformBinning.h>
+#include <RooFormulaVar.h>
+#include <RooGenericPdf.h>
 #include <RooHelpers.h>
+#include <RooHistPdf.h>
+#include <RooProduct.h>
+#include <RooRealVar.h>
+#include <RooUniformBinning.h>
 #include <RooWorkspace.h>
 
 #include <gtest/gtest.h>
@@ -340,3 +342,47 @@ TEST(GenericPdf, BinnedBoundariesSurviveRename)
    EXPECT_NE(integratorName.find("RooBinIntegrator"), std::string::npos) << integratorName;
    EXPECT_DOUBLE_EQ(value, piecewiseFlatIntegral);
 }
+
+// Ensure the implementation of RooGenericPdf::binBoundaries() (and
+// RooFormulaVar::binBoundaries()) is consistent with equivalent RooHistPdf.
+TEST(GenericPdf, BinnedBoundariesConsistentWithHistPdf)
+{
+   int nBins = 5;
+
+   RooRealVar x{"x", "x", 0, 0, 5};
+   x.setBins(nBins);
+
+   RooDataHist dh{"dh", "", x};
+
+   RooHistPdf hpdf{"hpdf", "", x, dh};
+   RooGenericPdf gpdf{"gpdf", "x[0] - x[0] + 1", {x}}; // uniform dummy
+   gpdf.setBinning(x, x.getBinning());
+
+   RooFormulaVar fvar{"fvar", "x[0] - x[0] + 1", {x}}; // uniform dummy
+   fvar.setBinning(x, x.getBinning());
+
+   // intentionally beyond the bin boundaries
+   double lo = -10;
+   double hi = 10;
+
+   std::unique_ptr<std::list<double>> boundsHistPdf{hpdf.binBoundaries(x, lo, hi)};
+   std::unique_ptr<std::list<double>> boundsGenericPdf{gpdf.binBoundaries(x, lo, hi)};
+   std::unique_ptr<std::list<double>> boundsFormulaVar{fvar.binBoundaries(x, lo, hi)};
+
+   EXPECT_EQ(boundsHistPdf->size(), nBins + 1);
+   EXPECT_EQ(boundsGenericPdf->size(), nBins + 1);
+   EXPECT_EQ(boundsFormulaVar->size(), nBins + 1);
+
+   auto iterHistPdf = boundsHistPdf->begin();
+   auto iterGenericPdf = boundsGenericPdf->begin();
+   auto iterFormulaVar = boundsFormulaVar->begin();
+
+   for (int i = 0; i < nBins + 1; ++i) {
+      EXPECT_EQ(*iterHistPdf, static_cast<double>(i));
+      EXPECT_EQ(*iterGenericPdf, static_cast<double>(i));
+      EXPECT_EQ(*iterFormulaVar, static_cast<double>(i));
+      iterHistPdf++;
+      iterGenericPdf++;
+      iterFormulaVar++;
+   }
+}