[RF] Support analytical integration of RooStepFunction

roofit/roofit/inc/RooStepFunction.h

@@ -38,6 +38,9 @@ class RooStepFunction : public RooAbsReal {
 
   std::list<double>* plotSamplingHint(RooAbsRealLValue& obs, double xlo, double xhi) const override ;
 
+  int getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=nullptr) const override;
+  double analyticalIntegral(int code, const char* rangeName=nullptr) const override;
+
  protected:
 
   double evaluate() const override;
@@ -48,6 +51,8 @@ class RooStepFunction : public RooAbsReal {
   RooListProxy _coefList ;
   RooListProxy _boundaryList ;
   bool       _interpolate = false;
+  mutable std::vector<double> _coefCache; //!
+  mutable std::vector<double> _boundaryCache; //!
 
   ClassDefOverride(RooStepFunction,1) //  Step Function
 };

roofit/roofit/src/RooStepFunction.cxx

@@ -32,12 +32,26 @@ but a not-normalized function (RooAbsReal)
 
 #include <RooArgList.h>
 #include <RooCurve.h>
-#include <RooMsgService.h>
+#include <RooFit/Detail/MathFuncs.h>
 #include <RooMath.h>
-#include <RooRealVar.h>
+#include <RooMsgService.h>
 
 ClassImp(RooStepFunction);
 
+namespace {
+
+// Get the values for each element in a RooListProxy, using the correct normalization set.
+std::span<double const> values(RooListProxy const &listProxy, std::vector<double> &valueCache)
+{
+   valueCache.resize(listProxy.size());
+   for (std::size_t i = 0; i < listProxy.size(); ++i) {
+      valueCache[i] = static_cast<RooAbsReal &>(listProxy[i]).getVal(listProxy.nset());
+   }
+   return valueCache;
+}
+
+} // namespace
+
 ////////////////////////////////////////////////////////////////////////////////
 /// Constructor
 
@@ -77,54 +91,51 @@ RooStepFunction::RooStepFunction(const RooStepFunction& other, const char* name)
 
 double RooStepFunction::evaluate() const
 {
-  std::vector<double> b(_boundaryList.size()) ;
-  std::vector<double> c(_coefList.size()+3) ;
-  Int_t nb(0) ;
-  for (auto * boundary : static_range_cast<RooAbsReal*>(_boundaryList)) {
-    b[nb++] = boundary->getVal() ;
-  }
+  std::span<const double> b = values(_boundaryList, _boundaryCache);
+  int nb = b.size();
 
   // Return zero if outside any boundaries
   if ((_x<b[0]) || (_x>b[nb-1])) return 0 ;
 
   if (!_interpolate) {
 
     // No interpolation -- Return values bin-by-bin
-    for (Int_t i=0;i<nb-1;i++){
+    for (int i=0;i<nb-1;i++){
       if (_x>b[i]&&_x<=b[i+1]) {
    return (static_cast<RooAbsReal*>(_coefList.at(i)))->getVal() ;
       }
     }
     return 0 ;
 
-  } else {
+  }
+
+    std::vector<double> c(_coefList.size()+3) ;
 
     // Interpolation
 
     // Make array of (b[0],bin centers,b[last])
     c[0] = b[0] ; c[nb] = b[nb-1] ;
-    for (Int_t i=0 ; i<nb-1 ; i++) {
+    for (int i=0 ; i<nb-1 ; i++) {
       c[i+1] = (b[i]+b[i+1])/2 ;
     }
 
     // Make array of (0,coefficient values,0)
-    Int_t nc(0) ;
+    int nc(0) ;
     std::vector<double> y(_coefList.size()+3) ;
     y[nc++] = 0 ;
     for(auto * coef : static_range_cast<RooAbsReal*>(_coefList)) {
       y[nc++] = coef->getVal() ;
     }
     y[nc++] = 0 ;
 
-    for (Int_t i=0;i<nc-1;i++){
+    for (int i=0;i<nc-1;i++){
       if (_x>c[i]&&_x<=c[i+1]) {
    double xx[2] ; xx[0]=c[i] ; xx[1]=c[i+1] ;
    double yy[2] ; yy[0]=y[i] ; yy[1]=y[i+1] ;
    return RooMath::interpolate(xx,yy,2,_x) ;
       }
     }
     return 0;
-  }
 }
 
 
@@ -134,14 +145,23 @@ std::list<double> *RooStepFunction::plotSamplingHint(RooAbsRealLValue &obs, doub
       return nullptr;
    }
 
-   // Retrieve position of all bin boundaries
-   std::vector<double> boundaries;
-   boundaries.reserve(_boundaryList.size());
-   for (auto *boundary : static_range_cast<RooAbsReal *>(_boundaryList)) {
-      boundaries.push_back(boundary->getVal());
-   }
-
    // Use the helper function from RooCurve to make sure to get sampling hints
    // that work with the RooFitPlotting.
-   return RooCurve::plotSamplingHintForBinBoundaries(boundaries, xlo, xhi);
+   return RooCurve::plotSamplingHintForBinBoundaries(values(_boundaryList, _boundaryCache), xlo, xhi);
+}
+
+int RooStepFunction::getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char * /*rangeName*/) const
+{
+   // We only support analytical integration if we integrate over x and there is
+   // no interpolation.
+   if (!_interpolate && matchArgs(allVars, analVars, _x))
+      return 1;
+   return 0;
+}
+
+double RooStepFunction::analyticalIntegral(int /*code*/, const char *rangeName) const
+{
+   return RooFit::Detail::MathFuncs::stepFunctionIntegral(_x.min(rangeName), _x.max(rangeName), _coefList.size(),
+                                                          values(_boundaryList, _boundaryCache).data(),
+                                                          values(_coefList, _coefCache).data());
 }

roofit/roofit/test/CMakeLists.txt

@@ -19,6 +19,7 @@ ROOT_ADD_GTEST(testRooKeysPdf testRooKeysPdf.cxx LIBRARIES Gpad RooFit)
 ROOT_ADD_GTEST(testRooLandau testRooLandau.cxx LIBRARIES RooFit)
 ROOT_ADD_GTEST(testRooParamHistFunc testRooParamHistFunc.cxx LIBRARIES Gpad RooFit)
 ROOT_ADD_GTEST(testRooPoisson testRooPoisson.cxx LIBRARIES RooFit)
+ROOT_ADD_GTEST(testRooStepFunc testRooStepFunc.cxx LIBRARIES RooFit)
 
 if(mathmore)
   ROOT_EXECUTABLE(testRooFit testRooFit.cxx LIBRARIES RooFit MathMore)

roofit/roofit/test/testRooStepFunc.cxx

@@ -0,0 +1,35 @@
+// Tests for the RooStepFunc
+// Authors: Jonas Rembser, CERN 2025
+
+#include <RooStepFunction.h>
+#include <RooWrapperPdf.h>
+#include <RooRealVar.h>
+#include <RooRealIntegral.h>
+
+#include <gtest/gtest.h>
+
+// Test that RooStepFunction can be wrapped into a pdf object, where it is
+// normalized without any analytical integrals.
+TEST(RooStepFunc, InPdfWrapper)
+{
+   RooRealVar x{"x", "x", 5.0, 0.0, 10.0};
+   RooRealVar a{"a", "a", 3.0, 0.0, 10.0};
+   RooRealVar b{"b", "b", 7.0, 0.0, 10.0};
+
+   RooStepFunction stepFunc{"step", "", x, RooArgList{1.0}, {a, b}};
+   RooWrapperPdf stepPdf{"pdf", "", stepFunc};
+
+   RooArgSet normSet{x};
+
+   RooArgSet intSet{x};
+   RooArgSet numSet{};
+
+   RooRealIntegral integ{"integ", "", stepPdf, x};
+
+   // Make sure that the normalization integral for the pdf is analytical
+   EXPECT_EQ(integ.anaIntVars().size(), 1);
+   EXPECT_DOUBLE_EQ(integ.getVal(), 4.0);
+
+   // Check that the pdf in correctly normalized
+   EXPECT_DOUBLE_EQ(stepPdf.getVal(normSet), 0.25);
+}

roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h

@@ -765,6 +765,20 @@ inline double multiVarGaussian(int n, const double *x, const double *mu, const d
    return std::exp(-0.5 * result);
 }
 
+// Integral of a step function defined by `nBins` intervals, where the
+// intervals have values `coefs` and the boundary on the interval `iBin` is
+// given by `[boundaries[i], boundaries[i+1])`.
+inline double stepFunctionIntegral(double xmin, double xmax, std::size_t nBins, double const *boundaries, double const *coefs)
+{
+   double out = 0.0;
+   for (std::size_t i = 0; i < nBins; ++i) {
+      double a = boundaries[i];
+      double b = boundaries[i + 1];
+      out += coefs[i] * std::max(0.0, std::min(b, xmax) - std::max(a, xmin));
+   }
+   return out;
+}
+
 } // namespace MathFuncs
 } // namespace Detail
 } // namespace RooFit

roofit/roofitcore/inc/RooWrapperPdf.h

@@ -53,8 +53,11 @@ class RooWrapperPdf final : public RooAbsPdf {
   bool selfNormalized() const override { return _selfNormalized; }
 
   // Analytical Integration handling
-  bool forceAnalyticalInt(const RooAbsArg& dep) const override {
-    return _func->forceAnalyticalInt(dep);
+  bool forceAnalyticalInt(const RooAbsArg& /*dep*/) const override {
+     // Just like with other wrapper classes like RooExtendPdf, we can safely
+     // use the analytical integration capabilities of the wrapped object,
+     // becuase we don't do any no-linear transformation.
+     return true;
   }
   Int_t getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& analVars, const RooArgSet* normSet,
       const char* rangeName=nullptr) const override {