[RF] Improve implementation of RooCurve::average()

When integrating the discretely-sampled RooCurves, the algorithm
implemented in RooCurve::average() was unnecessarily complicated.

The existing midpoints were only considered for the trapezoidal rule if
they are away from the interval limits with an arbitrary tolerance,
which seems like a premature optimization to me.

In particular, the logic was not correct if all midpoints were close to
the limits without this tolerance, resulting in issue root-project#9838.

Instead of making that case work correctly by implementing more code
branches, this commit suggests to simply don't do this tolerace check
and use all available midpoints for the trapezoidal integration rule.

A unit test with the reproducer from root-project#9838 is also implemented.

Closes root-project#9838.

roofit/roofitcore/src/RooCurve.cxx

@@ -590,31 +590,24 @@ double RooCurve::average(double xFirst, double xLast) const
   double yLast = interpolate(xLast,1e-10) ;
 
   // Find first and last mid points
-  Int_t ifirst = findPoint(xFirst,1e10) ;
-  Int_t ilast  = findPoint(xLast,1e10) ;
-  Point firstPt = getPoint(*this, ifirst);
-  Point lastPt = getPoint(*this, ilast);
+  Int_t ifirst = findPoint(xFirst, std::numeric_limits<double>::infinity());
+  Int_t ilast  = findPoint(xLast, std::numeric_limits<double>::infinity());
 
-  double tolerance=1e-3*(xLast-xFirst) ;
+  // Make sure the midpoints are actually in the interval
+  while (GetPointX(ifirst) < xFirst) {
+    ++ifirst;
+  }
+  while (GetPointX(ilast) > xLast) {
+    --ilast;
+  }
 
   // Handle trivial scenario -- no midway points, point only at or outside given range
-  if (ilast-ifirst==1 &&(firstPt.x-xFirst)<-1*tolerance && (lastPt.x-xLast)>tolerance) {
+  if (ilast < ifirst) {
     return 0.5*(yFirst+yLast) ;
   }
 
-  // If first point closest to xFirst is at xFirst or before xFirst take the next point
-  // as the first midway point
-  if ((firstPt.x-xFirst)<-1*tolerance) {
-    ifirst++ ;
-    firstPt = getPoint(*this, ifirst);
-  }
-
-  // If last point closest to yLast is at yLast or beyond yLast the previous point
-  // as the last midway point
-  if ((lastPt.x-xLast)>tolerance) {
-    ilast-- ;
-    lastPt = getPoint(*this, ilast);
-  }
+  Point firstPt = getPoint(*this, ifirst);
+  Point lastPt = getPoint(*this, ilast);
 
   // Trapezoid integration from lower edge to first midpoint
   double sum = 0.5 * (firstPt.x-xFirst)*(yFirst+firstPt.y);

roofit/roofitcore/test/CMakeLists.txt

@@ -57,12 +57,12 @@ if(clad)
     ROOT_ADD_GTEST(testRooFuncWrapper testRooFuncWrapper.cxx LIBRARIES RooFitCore RooFit HistFactory)
   endif()
 endif()
-ROOT_ADD_GTEST(testTestStatistics testTestStatistics.cxx LIBRARIES RooFitCore)
 ROOT_ADD_GTEST(testGlobalObservables testGlobalObservables.cxx LIBRARIES RooFitCore)
 ROOT_ADD_GTEST(testInterface TestStatistics/testInterface.cxx LIBRARIES RooFitCore)
 ROOT_ADD_GTEST(testLikelihoodSerial TestStatistics/testLikelihoodSerial.cxx LIBRARIES RooFitCore)
 ROOT_ADD_GTEST(testNaNPacker testNaNPacker.cxx LIBRARIES RooFitCore RooBatchCompute)
 ROOT_ADD_GTEST(testRooAbsL TestStatistics/testRooAbsL.cxx LIBRARIES RooFitCore)
+ROOT_ADD_GTEST(testRooCurve testRooCurve.cxx LIBRARIES RooFitCore)
 ROOT_ADD_GTEST(testRooHist testRooHist.cxx LIBRARIES RooFitCore)
 ROOT_ADD_GTEST(testRooHistPdf testRooHistPdf.cxx LIBRARIES RooFitCore)
 ROOT_ADD_GTEST(testRooPolyFunc testRooPolyFunc.cxx LIBRARIES Gpad RooFitCore)
@@ -72,6 +72,7 @@ ROOT_ADD_GTEST(testRooSTLRefCountList testRooSTLRefCountList.cxx LIBRARIES RooFi
 ROOT_ADD_GTEST(testRooSimultaneous testRooSimultaneous.cxx LIBRARIES RooFitCore)
 ROOT_ADD_GTEST(testRooTruthModel testRooTruthModel.cxx LIBRARIES RooFitCore RooFit)
 ROOT_ADD_GTEST(testSumW2Error testSumW2Error.cxx LIBRARIES Gpad RooFitCore)
+ROOT_ADD_GTEST(testTestStatistics testTestStatistics.cxx LIBRARIES RooFitCore)
 if (roofit_multiprocess)
   ROOT_ADD_GTEST(testTestStatisticsPlot TestStatistics/testPlot.cxx LIBRARIES RooFitMultiProcess RooFitCore RooFit
                    COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/TestStatistics/TestStatistics_ref.root)

roofit/roofitcore/test/testRooCurve.cxx

@@ -0,0 +1,48 @@
+/*
+ * Project: RooFit
+ * Authors:
+ *   Jonas Rembser, CERN 2024
+ *
+ * Copyright (c) 2024, CERN
+ *
+ * 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)
+ */
+
+#include <RooCmdArg.h>
+#include <RooGenericPdf.h>
+#include <RooHelpers.h>
+#include <RooPlot.h>
+#include <RooRealVar.h>
+
+#include "gtest_wrapper.h"
+
+// Cross-check to make sure the integration works correctly even if there is
+// only one midpoint on the RooCurve. Covers GitHub issue #9838 (the reproducer
+// in that issue was translated to this test).
+TEST(RooPlot, Average)
+{
+   // Silence the info about numeric integration because we don't care about it
+   RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING, 0u, RooFit::NumIntegration, true};
+
+   RooRealVar x("x", "x", 0, 50);
+   RooGenericPdf func("func", "Test Function", "x", x);
+
+   std::unique_ptr<RooPlot> xframe{x.frame()};
+
+   func.plotOn(xframe.get(), RooFit::Name("funcCurve"));
+
+   RooCurve *funcCurve = xframe->getCurve("funcCurve");
+
+   const double tol = 1e-10;
+
+   for (double i = 10; i < 11; i += 0.1) {
+      double avg = funcCurve->average(i, i + 0.1);
+
+      double xFirst = funcCurve->interpolate(i, tol);
+      double xLast = funcCurve->interpolate(i + 0.1, tol);
+
+      EXPECT_NEAR(avg, 0.5 * (xLast + xFirst), tol);
+   }
+}