diff --git a/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h b/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h
index 7119cc4b9a55d..7563cf881e0bc 100644
--- a/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h
+++ b/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h
@@ -728,6 +728,24 @@ inline double bernsteinIntegral(double xlo, double xhi, double xmin, double xmax
    return norm * (xmax - xmin);
 }
 
+template <class Function>
+inline double simpleNumericIntegral(double *params, double const *xlArr, Function func, double xmin, double xmax)
+{
+   const int n = 1000; // number of sampling points
+   double d = xmax - xmin;
+   double eps = d / n;
+   double out = 0.;
+   double obs[1];
+   for (int i = 0; i < n; ++i) {
+      obs[0] = xmin + eps * i;
+      double tmpA = func(params, obs, xlArr);
+      obs[0] = xmin + eps * (i + 1);
+      double tmpB = func(params, obs, xlArr);
+      out += (tmpA + tmpB) * 0.5 * eps;
+   }
+   return out;
+}
+
 } // namespace MathFuncs
 
 } // namespace Detail
diff --git a/roofit/roofitcore/src/RooRealIntegral.cxx b/roofit/roofitcore/src/RooRealIntegral.cxx
index 4e40317e9c3d7..0c7ab44b72d0e 100644
--- a/roofit/roofitcore/src/RooRealIntegral.cxx
+++ b/roofit/roofitcore/src/RooRealIntegral.cxx
@@ -1048,37 +1048,13 @@ void RooRealIntegral::translate(RooFit::CodegenContext &ctx) const
 
    auto &intVar = static_cast<RooAbsRealLValue &>(*_intList[0]);
 
-   std::string obsName = ctx.getTmpVarName();
    std::string oldIntVarResult = ctx.getResult(intVar);
    ctx.addResult(&intVar, "obs[0]");
 
    std::string funcName = ctx.buildFunction(*_function, {});
 
-   std::stringstream ss;
-
-   ss  << "double " << obsName << "[1];\n";
-
-   std::string resName = RooFit::Detail::makeValidVarName(GetName()) + "Result";
-   ctx.addResult(this, resName);
-   ctx.addToGlobalScope("double " + resName + " = 0.0;\n");
-
-   // TODO: once Clad has support for higher-order functions (follow also the
-   // Clad issue #637), we could refactor this code into an actual function
-   // instead of hardcoding it here as a string.
-   ss << "{\n"
-      << "   const int n = 1000; // number of sampling points\n"
-      << "   double d = " << intVar.getMax(intRange()) << " - " << intVar.getMin(intRange()) << ";\n"
-      << "   double eps = d / n;\n"
-      << "   for (int i = 0; i < n; ++i) {\n"
-      << "      " << obsName << "[0] = " << intVar.getMin(intRange()) << " + eps * i;\n"
-      << "      double tmpA = " << funcName << "(params, " << obsName << ", xlArr);\n"
-      << "      " << obsName << "[0] = " << intVar.getMin(intRange()) << " + eps * (i + 1);\n"
-      << "      double tmpB = " << funcName << "(params, " << obsName << ", xlArr);\n"
-      << "      " << resName << " += (tmpA + tmpB) * 0.5 * eps;\n"
-      << "   }\n"
-      << "}\n";
-
-   ctx.addToGlobalScope(ss.str());
+   ctx.addResult(this, ctx.buildCall("RooFit::Detail::MathFuncs::simpleNumericIntegral", "params", "xlArr", funcName,
+                                     intVar.getMin(intRange()), intVar.getMax(intRange())));
 
    ctx.addResult(&intVar, oldIntVarResult);
 }