diff --git a/interpreter/cling/tools/plugins/clad/CMakeLists.txt b/interpreter/cling/tools/plugins/clad/CMakeLists.txt
index ab93259e12fa1..0298aad36a1f2 100644
--- a/interpreter/cling/tools/plugins/clad/CMakeLists.txt
+++ b/interpreter/cling/tools/plugins/clad/CMakeLists.txt
@@ -71,7 +71,7 @@ if (DEFINED CLAD_SOURCE_DIR)
   list(APPEND _clad_extra_settings SOURCE_DIR ${CLAD_SOURCE_DIR})
 else()
   list(APPEND _clad_extra_settings GIT_REPOSITORY https://github.com/vgvassilev/clad.git)
-  list(APPEND _clad_extra_settings GIT_TAG v2.4)
+  list(APPEND _clad_extra_settings GIT_TAG master)
 endif()
 
 # list(APPEND _clad_patches_list "patch1.patch" "patch2.patch")
@@ -86,7 +86,7 @@ set(_clad_patch_command
 ExternalProject_Add(
   clad
   UPDATE_COMMAND ""
-  PATCH_COMMAND ${_clad_patch_command}
+  # PATCH_COMMAND ${_clad_patch_command}
   CMAKE_ARGS -G ${CMAKE_GENERATOR}
              -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}
              -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
diff --git a/math/mathcore/inc/Math/FitMethodFunction.h b/math/mathcore/inc/Math/FitMethodFunction.h
index 63eb814ac42aa..cb6aa7526d169 100644
--- a/math/mathcore/inc/Math/FitMethodFunction.h
+++ b/math/mathcore/inc/Math/FitMethodFunction.h
@@ -73,12 +73,13 @@ class BasicFitMethodFunction : public FunctionType {
    virtual double DataElement(const double *x, unsigned int i, double *g = nullptr, double *h = nullptr, bool fullHessian = false) const = 0;
 
    // flag to indicate if full Hessian computation is supported
-   virtual bool HasHessian() const { return false;}
+   bool HasHessian() const override { return false; }
 
    /**
     * Computes the full Hessian. Return false if Hessian is not supported
     */
-   virtual bool Hessian(const double * x, double * hess) const {
+   bool Hessian(const double *x, double *hess) const override
+   {
       //return full Hessian of  the objective function which is Sum(F(i))
       unsigned int np = NPoints();
       unsigned int ndim = NDim();
diff --git a/math/mathcore/inc/Math/Functor.h b/math/mathcore/inc/Math/Functor.h
index 8e6f848eec2c2..e7f0857ca3719 100644
--- a/math/mathcore/inc/Math/Functor.h
+++ b/math/mathcore/inc/Math/Functor.h
@@ -23,6 +23,7 @@
 #include <algorithm>
 #include <memory>
 #include <functional>
+#include <stdexcept>
 #include <type_traits>
 #include <vector>
 
@@ -202,10 +203,13 @@ class GradFunctor : public IGradientFunctionMultiDim  {
     * @param f  : function object computing the function value
     * @param dim : number of function dimension
     * @param g   : function object computing the function gradient
+    * @param h   : optional function object computing the function Hessian,
+    *              filling its second argument as a row-major dim * dim array
     */
-   GradFunctor(std::function<double(double const *)> const&f, unsigned int dim,
-               std::function<void(double const *, double *)> const& g)
-      : fDim{dim}, fFunc{f}, fGradFunc{g}
+   GradFunctor(std::function<double(double const *)> const &f, unsigned int dim,
+               std::function<void(double const *, double *)> const &g,
+               std::function<void(double const *, double *)> const &h = nullptr)
+      : fDim{dim}, fFunc{f}, fGradFunc{g}, fHessFunc{h}
    {}
 
    // Clone of the function handler (use copy-ctor).
@@ -224,6 +228,17 @@ class GradFunctor : public IGradientFunctionMultiDim  {
       fGradFunc(x, g);
    }
 
+   bool HasHessian() const override { return bool(fHessFunc); }
+
+   bool Hessian(const double *x, double *hess) const override
+   {
+      if (!fHessFunc) {
+         return false;
+      }
+      fHessFunc(x, hess);
+      return true;
+   }
+
 private :
 
    inline double DoEval (const double * x) const override {
@@ -246,6 +261,7 @@ private :
    std::function<double(const double *)> fFunc;
    std::function<double(double const *, unsigned int)> fDerivFunc;
    std::function<void(const double *, double*)> fGradFunc;
+   std::function<void(const double *, double *)> fHessFunc;
 };
 
 
@@ -285,29 +301,42 @@ class GradFunctor1D : public IGradientFunctionOneDim  {
        the function evaluation and the other for the derivative.
        The member functions must take a double as argument and return a double
     */
-   template <class PtrObj, typename MemFn, typename GradMemFn>
-   GradFunctor1D(const PtrObj& p, MemFn memFn, GradMemFn gradFn)
+   template <class PtrObj, typename MemFn, typename GradMemFn,
+             std::enable_if_t<std::is_member_pointer<MemFn>::value, bool> = true>
+   GradFunctor1D(const PtrObj &p, MemFn memFn, GradMemFn gradFn)
       : fFunc{std::bind(memFn, p, std::placeholders::_1)}, fDerivFunc{std::bind(gradFn, p, std::placeholders::_1)}
    {}
 
-
    /// Specialized constructor from 2 function objects implementing double
    /// operator()(double x). The first one for the function evaluation and the
-   /// second one implementing the function derivative.
-   GradFunctor1D(std::function<double(double)> const& f, std::function<double(double)> const& g)
-      : fFunc{f}, fDerivFunc{g}
+   /// second one implementing the function derivative. Optionally, a third
+   /// function object implementing the second derivative can be passed.
+   GradFunctor1D(std::function<double(double)> const &f, std::function<double(double)> const &g,
+                 std::function<double(double)> const &g2 = nullptr)
+      : fFunc{f}, fDerivFunc{g}, fSecondDerivFunc{g2}
    {}
 
    // clone of the function handler (use copy-ctor)
    GradFunctor1D * Clone() const override { return new GradFunctor1D(*this); }
 
+   bool HasHessian() const override { return bool(fSecondDerivFunc); }
+
 private :
 
    inline double DoEval (double x) const override { return fFunc(x); }
    inline double DoDerivative (double x) const override { return fDerivFunc(x); }
+   inline double DoSecondDerivative(double x) const override
+   {
+      if (!fSecondDerivFunc) {
+         throw std::runtime_error(
+            "The second derivative evaluation is not implemented for this function (see HasHessian())");
+      }
+      return fSecondDerivFunc(x);
+   }
 
    std::function<double(double)> fFunc;
    std::function<double(double)> fDerivFunc;
+   std::function<double(double)> fSecondDerivFunc;
 };
 
 
diff --git a/math/mathcore/inc/Math/IFunction.h b/math/mathcore/inc/Math/IFunction.h
index bc0b029d4bceb..a24a970c130d7 100644
--- a/math/mathcore/inc/Math/IFunction.h
+++ b/math/mathcore/inc/Math/IFunction.h
@@ -33,6 +33,7 @@
 
 #include "Math/IFunctionfwd.h"
 
+#include <stdexcept>
 
 namespace ROOT {
    namespace Math {
@@ -111,6 +112,15 @@ namespace ROOT {
             Gradient(x, df);
          }
 
+         // Indicate whether this class supports second derivative (Hessian) calculations,
+         // i.e., if the Hessian() method is implemented.
+         virtual bool HasHessian() const { return false; }
+
+         /// Evaluate all second derivatives (the Hessian matrix) at a point x,
+         /// filling hess as a row-major NDim() * NDim() array.
+         /// Return false if the second derivatives are not implemented (see HasHessian()).
+         virtual bool Hessian(const T * /*x*/, T * /*hess*/) const { return false; }
+
          /// Return the partial derivative with respect to the passed coordinate.
          T Derivative(const T *x, unsigned int icoord = 0) const { return DoDerivative(x, icoord); }
 
@@ -178,6 +188,10 @@ namespace ROOT {
          // if it inherits from ROOT::Math::IGradientFunctionOneDim.
          virtual bool HasGradient() const { return false; }
 
+         // Indicate whether this class supports second derivative calculations,
+         // i.e., if the SecondDerivative() method is implemented.
+         virtual bool HasHessian() const { return false; }
+
          /// Return the derivative of the function at a point x
          /// Use the private method DoDerivative
          double Derivative(double x) const { return DoDerivative(x); }
@@ -188,6 +202,21 @@ namespace ROOT {
          /// Compatibility method with multi-dimensional interface for Gradient.
          void Gradient(const double *x, double *g) const { g[0] = DoDerivative(*x); }
 
+         /// Return the second derivative of the function at a point x.
+         /// Check HasHessian() to see whether it is implemented.
+         double SecondDerivative(double x) const { return DoSecondDerivative(x); }
+
+         /// Compatibility method with multi-dimensional interface for the Hessian.
+         /// Return false if the second derivative is not implemented (see HasHessian()).
+         bool Hessian(const double *x, double *hess) const
+         {
+            if (!HasHessian()) {
+               return false;
+            }
+            hess[0] = DoSecondDerivative(*x);
+            return true;
+         }
+
          /// Optimized method to evaluate at the same time the function value and derivative at a point x.
          /// Often both value and derivatives are needed and it is often more efficient to compute them at the same time.
          /// Derived class should implement this method if performances play an important role and if it is faster to
@@ -208,6 +237,13 @@ namespace ROOT {
 
          /// Function to evaluate the derivative with respect each coordinate. To be implemented by the derived class.
          virtual double  DoDerivative(double) const { return 0.; }
+
+         /// Function to evaluate the second derivative. To be implemented by the derived class (see HasHessian()).
+         virtual double DoSecondDerivative(double) const
+         {
+            throw std::runtime_error(
+               "The second derivative evaluation is not implemented for this function (see HasHessian())");
+         }
       };
 
 
diff --git a/roofit/codegen/src/CodegenImpl.cxx b/roofit/codegen/src/CodegenImpl.cxx
index 750f6912bc380..574fec2480a0c 100644
--- a/roofit/codegen/src/CodegenImpl.cxx
+++ b/roofit/codegen/src/CodegenImpl.cxx
@@ -483,6 +483,67 @@ void functorCodegenImpl(RooArg_t &arg, RooArgList const &variables, CodegenConte
               "}\n"
               "} // namespace clad::custom_derivatives\n";
 
+      // If the functor also provides second derivatives, emit the
+      // forward-mode derivative and its pullback, which Clad needs to
+      // generate second derivatives of code calling the wrapper (it
+      // implements Hessians as forward-mode derivatives that are then
+      // differentiated in reverse mode).
+      if (arg.function()->HasHessian()) {
+         code += "#include <plugins/include/clad/Differentiator/BuiltinDerivatives.h>\n"
+                 "namespace clad::custom_derivatives {\n\n"
+                 "clad::ValueAndPushforward<double, double> " +
+                 wrapperName +
+                 "_pushforward(double const *x, double const *d_x) {\n"
+                 "   double grad[" +
+                 nStr +
+                 "]{};\n"
+                 "   " +
+                 funcAddrCasted +
+                 "->Gradient(x, grad);\n"
+                 "   double dot = 0.;\n"
+                 "   for (int i = 0; i < " +
+                 nStr +
+                 "; ++i) {\n"
+                 "      dot += grad[i] * d_x[i];\n"
+                 "   }\n"
+                 "   return {" +
+                 funcAddrCasted +
+                 "->operator()(x), dot};\n"
+                 "}\n\n"
+                 "void " +
+                 wrapperName +
+                 "_pushforward_pullback(double const *x, double const *d_x, clad::ValueAndPushforward<double, "
+                 "double> d_y, double *d_out_x, double *d_out_d_x) {\n"
+                 "   double grad[" +
+                 nStr +
+                 "]{};\n"
+                 "   " +
+                 funcAddrCasted +
+                 "->Gradient(x, grad);\n"
+                 "   double hess[" +
+                 nStr + " * " + nStr +
+                 "]{};\n"
+                 "   " +
+                 funcAddrCasted +
+                 "->Hessian(x, hess);\n"
+                 "   for (int i = 0; i < " +
+                 nStr +
+                 "; ++i) {\n"
+                 "      double hdot = 0.;\n"
+                 "      for (int j = 0; j < " +
+                 nStr +
+                 "; ++j) {\n"
+                 "         hdot += hess[i * " +
+                 nStr +
+                 " + j] * d_x[j];\n"
+                 "      }\n"
+                 "      d_out_x[i] += d_y.value * grad[i] + d_y.pushforward * hdot;\n"
+                 "      d_out_d_x[i] += d_y.pushforward * grad[i];\n"
+                 "   }\n"
+                 "}\n"
+                 "} // namespace clad::custom_derivatives\n";
+      }
+
       gInterpreter->Declare(code.c_str());
    }
 
@@ -590,7 +651,7 @@ void codegenImpl(RooLandau &arg, CodegenContext &ctx)
 
 void codegenImpl(RooLognormal &arg, CodegenContext &ctx)
 {
-   std::string funcName = arg.useStandardParametrization() ? "logNormalEvaluateStandard" : "logNormal";
+   std::string funcName = arg.useStandardParametrization() ? "logNormalStandard" : "logNormal";
    ctx.addResult(&arg, ctx.buildCall(mathFunc(funcName), arg.getX(), arg.getShapeK(), arg.getMedian()));
 }
 
diff --git a/roofit/histfactory/test/testHistFactory.cxx b/roofit/histfactory/test/testHistFactory.cxx
index 29acd2c4ded64..58c1725efa9e3 100644
--- a/roofit/histfactory/test/testHistFactory.cxx
+++ b/roofit/histfactory/test/testHistFactory.cxx
@@ -12,17 +12,19 @@
 #include <RooFitHS3/JSONIO.h>
 #include <RooFitHS3/RooJSONFactoryWSTool.h>
 
-#include <RooFit/Detail/NormalizationHelpers.h>
-#include <RooDataHist.h>
-#include <RooWorkspace.h>
 #include <RooArgSet.h>
-#include <RooSimultaneous.h>
-#include <RooRealSumPdf.h>
-#include <RooRealVar.h>
-#include <RooHelpers.h>
+#include <RooDataHist.h>
+#include <RooEvaluatorWrapper.h>
+#include <RooFit/Detail/NormalizationHelpers.h>
+#include <RooFit/Evaluator.h>
 #include <RooFitResult.h>
+#include <RooHelpers.h>
+#include <RooMinimizer.h>
 #include <RooPlot.h>
-#include <RooFit/Evaluator.h>
+#include <RooRealSumPdf.h>
+#include <RooRealVar.h>
+#include <RooSimultaneous.h>
+#include <RooWorkspace.h>
 
 #include <TROOT.h>
 #include <TFile.h>
@@ -610,8 +612,21 @@ TEST_P(HFFixtureFit, Fit)
    }
 
    using namespace RooFit;
-   std::unique_ptr<RooFitResult> fitResult{simPdf->fitTo(
-      *data, evalBackend, GlobalObservables(*mc->GetGlobalObservables()), Save(), PrintLevel(verbose ? 1 : -1))};
+   std::unique_ptr<RooAbsReal> nll{
+      simPdf->createNLL(*data, evalBackend, GlobalObservables(*mc->GetGlobalObservables()))};
+   RooMinimizer::Config cfg;
+   if (evalBackend == RooFit::EvalBackend::Codegen()) {
+      // Make sure we use both analytical gradient and Hessian
+      static_cast<RooFit::Experimental::RooEvaluatorWrapper &>(*nll).generateGradient();
+      static_cast<RooFit::Experimental::RooEvaluatorWrapper &>(*nll).generateHessian();
+      cfg.useGradient = true;
+      cfg.useHessian = true;
+   }
+   RooMinimizer minim{*nll, cfg};
+   minim.setPrintLevel(verbose ? 1 : -1);
+   minim.minimize("Minuit2", "Migrad");
+   minim.hesse();
+   std::unique_ptr<RooFitResult> fitResult{minim.save()};
    ASSERT_NE(fitResult, nullptr);
    if (verbose)
       fitResult->Print("v");
diff --git a/roofit/roofit/src/RooONNXFunc.cxx b/roofit/roofit/src/RooONNXFunc.cxx
index b3f9e86506948..ac8a7d54b6e04 100644
--- a/roofit/roofit/src/RooONNXFunc.cxx
+++ b/roofit/roofit/src/RooONNXFunc.cxx
@@ -192,6 +192,7 @@ struct RooONNXFunc::RuntimeCache {
 
    RooFit::Detail::AnyWithVoidPtr _session;
    RooFit::Detail::AnyWithVoidPtr _d_session;
+   RooFit::Detail::AnyWithVoidPtr _zero_session; ///< Zero-weight session, the tangent for forward-mode derivatives.
    Func _func;
 };
 
@@ -295,20 +296,26 @@ std::string _RooONNXFunc_onnxToCppWithSofie(std::uint8_t const *onnxBytes, std::
       gInterpreter->ProcessLine(("std::size(" + namespaceName + "::inputTensorDims);").c_str()));
 
    // Per-input-tensor parameter / argument lists used by the JIT'd code below.
-   std::string innerParams;       // "float const *input0, float const *input1, ..."
-   std::string innerArgs;         // "input0, input1, ..."
-   std::string outerDoubleParams; // "double const *input0, double const *input1, ..."
-   std::string cladInputs;        // "input0, input1, ..."  (for clad::gradient param spec)
+   std::string innerParams;        // "float const *input0, float const *input1, ..."
+   std::string innerArgs;          // "input0, input1, ..."
+   std::string innerTangentParams; // "float const *d_input0, float const *d_input1, ..."
+   std::string innerTangentArgs;   // "d_input0, d_input1, ..."
+   std::string outerDoubleParams;  // "double const *input0, double const *input1, ..."
+   std::string cladInputs;         // "input0, input1, ..."  (for clad::gradient param spec)
    for (std::size_t i = 0; i < nInputTensors; ++i) {
       std::string istr = std::to_string(i);
       if (i > 0) {
          innerParams += ", ";
          innerArgs += ", ";
+         innerTangentParams += ", ";
+         innerTangentArgs += ", ";
          outerDoubleParams += ", ";
          cladInputs += ", ";
       }
       innerParams += "float const *input" + istr;
       innerArgs += "input" + istr;
+      innerTangentParams += "float const *d_input" + istr;
+      innerTangentArgs += "d_input" + istr;
       outerDoubleParams += "double const *input" + istr;
       cladInputs += "input" + istr;
    }
@@ -363,6 +370,47 @@ std::string _RooONNXFunc_onnxToCppWithSofie(std::uint8_t const *onnxBytes, std::
 
    gInterpreter->ProcessLine(("clad::gradient(" + namespaceName + "::roo_wrapper, \"" + cladInputs + "\");").c_str());
 
+   // Also generate roo_inner_wrapper_pushforward: the forward-mode derivative
+   // of the inner wrapper, which takes the session and input tangents as
+   // function arguments (requesting the derivative of the outer wrapper
+   // generates it, the same trick as for the pullback above). It is the basis
+   // of the exact second derivatives below.
+   gInterpreter->ProcessLine(("clad::differentiate(" + namespaceName + "::roo_wrapper, \"input0[0]\");").c_str());
+
+   // The session tangent for the pushforward: the weights are constants, so
+   // their tangent is zero (a default-constructed Session holds the actual
+   // weights). Intermediate tensors inside this session are used as tangent
+   // scratch space during the forward pass.
+   _runtime->_zero_session.emplace(sessionName);
+   auto ptrZeroSession = toPtrString(_runtime->_zero_session.ptr, sessionName);
+   gInterpreter->ProcessLine((ptrZeroSession + "->SetWeightsToZero();").c_str());
+
+   // Directional-derivative wrapper: reverse differentiation of it w.r.t.
+   // both the inputs and the tangent direction gives the exact
+   // Hessian-vector product H . d_input (in the input adjoints) and the
+   // gradient (in the tangent adjoints). Both must be requested as active:
+   // for a non-varied argument clad would look up the custom helper pullbacks
+   // with reduced signatures, not find them, and silently fall back to
+   // differentiating the BLAS calls, which yields zero.
+   {
+      std::ostringstream ss;
+      ss << "namespace " << namespaceName << " {\n\n"
+         << "float roo_dir(Session const &session, Session const &zeroSession, " << innerParams << ", "
+         << innerTangentParams << ") {\n"
+         << "   return roo_inner_wrapper_pushforward(session, " << innerArgs << ", zeroSession, " << innerTangentArgs
+         << ").pushforward;\n"
+         << "}\n\n"
+         << "float roo_dir_outer(Session const &session, Session const &zeroSession, " << innerParams << ", "
+         << innerTangentParams << ") {\n"
+         << "   return roo_dir(session, zeroSession, " << innerArgs << ", " << innerTangentArgs << ");\n"
+         << "}\n\n"
+         << "} // namespace " << namespaceName << "\n";
+      gInterpreter->Declare(ss.str().c_str());
+   }
+   gInterpreter->ProcessLine(
+      ("clad::gradient(" + namespaceName + "::roo_dir_outer, \"" + cladInputs + ", " + innerTangentArgs + "\");")
+         .c_str());
+
    // The codegen call site (CodegenImpl::codegenImpl(RooONNXFunc)) passes one
    // double-array argument per input tensor. Emit roo_outer_wrapper and the matching
    // custom-derivative pullback with the corresponding number of parameters.
@@ -416,6 +464,108 @@ std::string _RooONNXFunc_onnxToCppWithSofie(std::uint8_t const *onnxBytes, std::
             << "       d_input" << i << "[i] += d_inputFlt" << i << "[i];\n"
             << "    }\n";
       }
+      ss << "}\n\n";
+
+      // Custom derivatives for Clad Hessians, which are implemented in clad
+      // as forward-mode derivatives that are then differentiated in reverse
+      // mode: the forward pass needs a "pushforward" for the opaque
+      // roo_outer_wrapper call, and the reverse pass over it needs the
+      // matching "pushforward_pullback". Both are exact: they are built on
+      // the clad-generated roo_inner_wrapper_pushforward and roo_dir_pullback
+      // (see above).
+
+      // Comma-separated parameter helper for the emitted code.
+      std::string tangentDoubleParams; // "double const *d_input0, ..."
+      for (std::size_t i = 0; i < nInputTensors; ++i) {
+         std::string istr = std::to_string(i);
+         if (i > 0) {
+            tangentDoubleParams += ", ";
+         }
+         tangentDoubleParams += "double const *d_input" + istr;
+      }
+
+      // Emits float conversion buffers for the inputs and the tangents.
+      auto emitFloatBuffers = [&]() {
+         for (std::size_t i = 0; i < nInputTensors; ++i) {
+            ss << "    float inputFlt" << i << "[inputTensorDims[" << i << "].total_size()];\n"
+               << "    float dInputFlt" << i << "[inputTensorDims[" << i << "].total_size()];\n"
+               << "    for (::std::size_t i = 0; i < ::std::size(inputFlt" << i << "); ++i) {\n"
+               << "       inputFlt" << i << "[i] = input" << i << "[i];\n"
+               << "       dInputFlt" << i << "[i] = d_input" << i << "[i];\n"
+               << "    }\n";
+         }
+      };
+
+      // The pushforward evaluates the function value together with the exact
+      // directional derivative grad . d_input, in one forward pass.
+      ss << "clad::ValueAndPushforward<double, double> roo_outer_wrapper_pushforward(" << outerDoubleParams << ", "
+         << tangentDoubleParams << ") {\n"
+         << "    using namespace ::" << namespaceName << ";\n"
+         << "    auto &session = *" << ptrSession << ";\n"
+         << "    auto &zeroSession = *" << ptrZeroSession << ";\n";
+      emitFloatBuffers();
+      ss << "    auto vp = roo_inner_wrapper_pushforward(session, ";
+      for (std::size_t i = 0; i < nInputTensors; ++i) {
+         ss << "inputFlt" << i << ", ";
+      }
+      ss << "zeroSession";
+      for (std::size_t i = 0; i < nInputTensors; ++i) {
+         ss << ", dInputFlt" << i;
+      }
+      ss << ");\n"
+         << "    return {vp.value, vp.pushforward};\n"
+         << "}\n\n";
+
+      // The pullback of the pushforward needs second derivatives only in the
+      // form of the Hessian-vector product H . d_input, which roo_dir_pullback
+      // evaluates exactly (together with the gradient, as the adjoint of the
+      // tangent direction).
+      ss << "void roo_outer_wrapper_pushforward_pullback(" << outerDoubleParams << ", " << tangentDoubleParams
+         << ", clad::ValueAndPushforward<double, double> d_y";
+      for (std::size_t i = 0; i < nInputTensors; ++i) {
+         ss << ", double *d_out_input" << i;
+      }
+      for (std::size_t i = 0; i < nInputTensors; ++i) {
+         ss << ", double *d_out_d_input" << i;
+      }
+      ss << ") {\n"
+         << "    using namespace ::" << namespaceName << ";\n"
+         << "    auto &session = *" << ptrSession << ";\n"
+         << "    auto &zeroSession = *" << ptrZeroSession << ";\n"
+         << "    // Fresh adjoint sessions for every call: the clad-generated\n"
+         << "    // second-order pullback does not restore the adjoint state inside\n"
+         << "    // the sessions to zero (clad bug, see the \"clad referenced\n"
+         << "    // '_tracker...' before its declaration\" warnings it prints during\n"
+         << "    // generation), so reusing them would leak state between calls.\n"
+         << "    Session dSession1;\n"
+         << "    Session dSession2;\n";
+      emitFloatBuffers();
+      for (std::size_t i = 0; i < nInputTensors; ++i) {
+         ss << "    float hvpFlt" << i << "[inputTensorDims[" << i << "].total_size()] = {};\n"
+            << "    float gradFlt" << i << "[inputTensorDims[" << i << "].total_size()] = {};\n";
+      }
+      ss << "    roo_dir_pullback(session, zeroSession, ";
+      for (std::size_t i = 0; i < nInputTensors; ++i) {
+         ss << "inputFlt" << i << ", ";
+      }
+      for (std::size_t i = 0; i < nInputTensors; ++i) {
+         ss << "dInputFlt" << i << ", ";
+      }
+      ss << "1.F, &dSession1, &dSession2";
+      for (std::size_t i = 0; i < nInputTensors; ++i) {
+         ss << ", hvpFlt" << i;
+      }
+      for (std::size_t i = 0; i < nInputTensors; ++i) {
+         ss << ", gradFlt" << i;
+      }
+      ss << ");\n";
+      for (std::size_t i = 0; i < nInputTensors; ++i) {
+         ss << "    for (::std::size_t i = 0; i < ::std::size(inputFlt" << i << "); ++i) {\n"
+            << "       d_out_input" << i << "[i] += d_y.value * gradFlt" << i << "[i] + d_y.pushforward * hvpFlt" << i
+            << "[i];\n"
+            << "       d_out_d_input" << i << "[i] += d_y.pushforward * gradFlt" << i << "[i];\n"
+            << "    }\n";
+      }
       ss << "}\n\n"
          << "} // namespace " << namespaceName << "\n"
          << "} // namespace clad::custom_derivatives\n";
diff --git a/roofit/roofit/test/create_onnx_model.py b/roofit/roofit/test/create_onnx_model.py
index b5568dc0fbce6..90b38d729bba2 100644
--- a/roofit/roofit/test/create_onnx_model.py
+++ b/roofit/roofit/test/create_onnx_model.py
@@ -1,19 +1,18 @@
 import numpy as np
 import torch
 import torch.nn as nn
-import torch.optim as optim
 
 
 class SmallMLP(nn.Module):
     """Single-input MLP regression model."""
 
-    def __init__(self, in_features=10, hidden=32, out_features=1):
+    def __init__(self, in_features=10, hidden=32, out_features=1, activation=nn.ReLU):
         super().__init__()
         self.net = nn.Sequential(
             nn.Linear(in_features, hidden),
-            nn.ReLU(),
+            activation(),
             nn.Linear(hidden, hidden),
-            nn.ReLU(),
+            activation(),
             nn.Linear(hidden, out_features),
         )
 
@@ -24,13 +23,13 @@ def forward(self, x):
 class TwoInputMLP(nn.Module):
     """Dual-input MLP: concatenates two input tensors then runs an MLP."""
 
-    def __init__(self, in_features_a=10, in_features_b=5, hidden=32, out_features=1):
+    def __init__(self, in_features_a=10, in_features_b=5, hidden=32, out_features=1, activation=nn.ReLU):
         super().__init__()
         self.net = nn.Sequential(
             nn.Linear(in_features_a + in_features_b, hidden),
-            nn.ReLU(),
+            activation(),
             nn.Linear(hidden, hidden),
-            nn.ReLU(),
+            activation(),
             nn.Linear(hidden, out_features),
         )
 
@@ -68,7 +67,7 @@ def export_model(
     return model
 
 
-def run_inference_and_save(model, example_inputs, name):
+def run_inference_and_save(model, example_inputs, name, with_hessian=False):
     """Run forward+backward with fixed inputs and save prediction + gradients."""
     model = model.cpu()
     inputs = [t.clone().detach().requires_grad_(True) for t in example_inputs]
@@ -83,6 +82,19 @@ def run_inference_and_save(model, example_inputs, name):
     for i, x in enumerate(inputs):
         np.savetxt(f"{name}_grad_{i}.txt", x.grad.detach().numpy())
 
+    if with_hessian:
+        # Full Hessian over the concatenation of all (flattened) input
+        # tensors, matching the parameter ordering used on the RooFit side.
+        blocks = torch.autograd.functional.hessian(
+            lambda *args: model(*args).squeeze(), tuple(t.clone().detach() for t in example_inputs)
+        )
+        n_per_input = [t.numel() for t in example_inputs]
+        hess = np.block([
+            [blocks[i][j].reshape(n_per_input[i], n_per_input[j]).numpy() for j in range(len(example_inputs))]
+            for i in range(len(example_inputs))
+        ])
+        np.savetxt(f"{name}_hessian.txt", hess)
+
 
 def main():
     # ---- Model 1: single input ----
@@ -114,6 +126,39 @@ def main():
         name="regression_mlp_two_input",
     )
 
+    # ---- Models 3 and 4: like models 1 and 2, but with a smooth activation
+    # function so that the Hessian with respect to the inputs is non-trivial
+    # (a ReLU MLP is piecewise-linear in its inputs, so its input Hessian
+    # vanishes almost everywhere). Used for the Hessian tests.
+    model3 = SmallMLP(in_features=10, hidden=32, out_features=1, activation=nn.Tanh)
+    model3 = export_model(
+        model3,
+        input_shapes=[(10,)],
+        onnx_path="regression_mlp_tanh.onnx",
+    )
+    run_inference_and_save(
+        model3,
+        example_inputs=[torch.tensor([[0.1] * 10])],
+        name="regression_mlp_tanh",
+        with_hessian=True,
+    )
+
+    model4 = TwoInputMLP(in_features_a=10, in_features_b=5, hidden=32, out_features=1, activation=nn.Tanh)
+    model4 = export_model(
+        model4,
+        input_shapes=[(10,), (5,)],
+        onnx_path="regression_mlp_two_input_tanh.onnx",
+    )
+    run_inference_and_save(
+        model4,
+        example_inputs=[
+            torch.tensor([[0.1] * 10]),
+            torch.tensor([[0.2] * 5]),
+        ],
+        name="regression_mlp_two_input_tanh",
+        with_hessian=True,
+    )
+
 
 if __name__ == "__main__":
     main()
diff --git a/roofit/roofit/test/testRooONNXFunc.cxx b/roofit/roofit/test/testRooONNXFunc.cxx
index 06b2638a84ffe..f2f2d5975261e 100644
--- a/roofit/roofit/test/testRooONNXFunc.cxx
+++ b/roofit/roofit/test/testRooONNXFunc.cxx
@@ -2,6 +2,7 @@
 // Authors: Jonas Rembser, CERN 2026
 
 #include <RooONNXFunc.h>
+#include <RooProduct.h>
 #include <RooRealVar.h>
 #include <RooDataSet.h>
 #include <RooEvaluatorWrapper.h>
@@ -197,4 +198,119 @@ TEST(RooONNXFunc, Basic_CodegenAD_2Tensors)
       EXPECT_NEAR(output_vec[10 + i], refGrad1[i], 1e-5);
    }
 }
+
+/// Validate the Clad Hessian of a RooONNXFunc against a PyTorch reference.
+/// The Hessian tests use models with tanh activations: a ReLU MLP is
+/// piecewise-linear in its inputs, so its input Hessian vanishes almost
+/// everywhere and would not validate anything.
+TEST(RooONNXFunc, CodegenHessian)
+{
+   RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING, 0u, RooFit::Fitting, true};
+
+   double refPred = readDoublesFromFile("regression_mlp_tanh_pred.txt")[0];
+   std::vector<double> refGrad = readDoublesFromFile("regression_mlp_tanh_grad_0.txt");
+   std::vector<double> refHess = readDoublesFromFile("regression_mlp_tanh_hessian.txt");
+
+   const std::size_t n = 10;
+
+   RooArgList args;
+   fillArgs(args, n);
+
+   RooONNXFunc roo_func{"func", "", {args}, "regression_mlp_tanh.onnx"};
+
+   RooDataSet data("data", "data", {});
+
+   RooFit::Experimental::RooEvaluatorWrapper roo_final{roo_func, &data, false, "", nullptr, false};
+
+   EXPECT_NEAR(roo_final.getVal(), refPred, 1e-5);
+
+   roo_final.generateGradient();
+
+   std::vector<double> grad(n);
+   roo_final.gradient(grad.data());
+   for (std::size_t i = 0; i < n; ++i) {
+      EXPECT_NEAR(grad[i], refGrad[i], 1e-5);
+   }
+
+   roo_final.generateHessian();
+
+   // The second derivatives are exact Hessian-vector products (reverse-mode
+   // derivatives of the forward-mode derivative of the generated code), so
+   // they merit the same tolerance as the gradient checks.
+   std::vector<double> hess(n * n);
+   roo_final.hessian(hess.data());
+   for (std::size_t i = 0; i < n * n; ++i) {
+      EXPECT_NEAR(hess[i], refHess[i], 1e-5);
+   }
+}
+
+/// Test the Clad Hessian of a RooONNXFunc with two input tensors, including
+/// the Hessian blocks that mix the two tensors.
+TEST(RooONNXFunc, CodegenHessian_2Tensors)
+{
+   RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING, 0u, RooFit::Fitting, true};
+
+   double refPred = readDoublesFromFile("regression_mlp_two_input_tanh_pred.txt")[0];
+   std::vector<double> refHess = readDoublesFromFile("regression_mlp_two_input_tanh_hessian.txt");
+
+   RooArgList args0;
+   fillArgs(args0, 10, 0.1, "a");
+   RooArgList args1;
+   fillArgs(args1, 5, 0.2, "b");
+
+   RooONNXFunc roo_func{"func", "", {args0, args1}, "regression_mlp_two_input_tanh.onnx"};
+
+   RooDataSet data("data", "data", {});
+
+   RooFit::Experimental::RooEvaluatorWrapper roo_final{roo_func, &data, false, "", nullptr, false};
+
+   EXPECT_NEAR(roo_final.getVal(), refPred, 1e-5);
+
+   roo_final.generateHessian();
+
+   const std::size_t nTotal = 10 + 5;
+   std::vector<double> hess(nTotal * nTotal);
+   roo_final.hessian(hess.data());
+   for (std::size_t i = 0; i < nTotal * nTotal; ++i) {
+      EXPECT_NEAR(hess[i], refHess[i], 1e-5);
+   }
+}
+
+/// Test the Clad Hessian of a compound expression: the square of a
+/// RooONNXFunc. When the RooONNXFunc is the top-level function, the adjoint
+/// of the primal value in the emitted pushforward pullback stays zero; the
+/// product rule in H(f^2) = 2 * (f * H + grad * grad^T) exercises it.
+TEST(RooONNXFunc, CodegenHessianCompound)
+{
+   RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING, 0u, RooFit::Fitting, true};
+
+   double refPred = readDoublesFromFile("regression_mlp_tanh_pred.txt")[0];
+   std::vector<double> refGrad = readDoublesFromFile("regression_mlp_tanh_grad_0.txt");
+   std::vector<double> refHess = readDoublesFromFile("regression_mlp_tanh_hessian.txt");
+
+   const std::size_t n = 10;
+
+   RooArgList args;
+   fillArgs(args, n);
+
+   RooONNXFunc roo_func{"func", "", {args}, "regression_mlp_tanh.onnx"};
+   RooProduct square{"square", "", {roo_func, roo_func}};
+
+   RooDataSet data("data", "data", {});
+
+   RooFit::Experimental::RooEvaluatorWrapper roo_final{square, &data, false, "", nullptr, false};
+
+   EXPECT_NEAR(roo_final.getVal(), refPred * refPred, 1e-5);
+
+   roo_final.generateHessian();
+
+   std::vector<double> hess(n * n);
+   roo_final.hessian(hess.data());
+   for (std::size_t i = 0; i < n; ++i) {
+      for (std::size_t j = 0; j < n; ++j) {
+         const double ref = 2. * (refPred * refHess[i * n + j] + refGrad[i] * refGrad[j]);
+         EXPECT_NEAR(hess[i * n + j], ref, 1e-5);
+      }
+   }
+}
 #endif
diff --git a/roofit/roofitcore/src/RooFactoryWSTool.cxx b/roofit/roofitcore/src/RooFactoryWSTool.cxx
index adcd24acf0889..ef36fb8b71188 100644
--- a/roofit/roofitcore/src/RooFactoryWSTool.cxx
+++ b/roofit/roofitcore/src/RooFactoryWSTool.cxx
@@ -1873,7 +1873,14 @@ const char* RooFactoryWSTool::asSTRING(const char* arg)
 
 Int_t RooFactoryWSTool::asINT(const char* arg)
 {
-  return atoi(arg) ;
+   // This is also the conversion used for bool constructor arguments, where
+   // atoi() would silently turn the spelled-out literals into zero.
+   const std::string s{arg};
+   if (s == "true")
+      return 1;
+   if (s == "false")
+      return 0;
+   return atoi(arg);
 }
 
 
diff --git a/roofit/roofitcore/test/testRooFuncWrapper.cxx b/roofit/roofitcore/test/testRooFuncWrapper.cxx
index 9a5250e8b25b7..91095a201fa4f 100644
--- a/roofit/roofitcore/test/testRooFuncWrapper.cxx
+++ b/roofit/roofitcore/test/testRooFuncWrapper.cxx
@@ -78,7 +78,9 @@ void randomizeParameters(const RooArgSet &parameters)
       double mul = rng.Uniform(lowerBound, upperBound);
 
       auto par = dynamic_cast<RooAbsRealLValue *>(param);
-      if (!par)
+      // Constant parameters must not be varied. Note that they can also have
+      // an infinite range, so the randomization logic below would break.
+      if (!par || par->isConstant())
          continue;
       double val = par->getVal();
       val = val + mul * (mul > 0 ? (par->getMax() - val) : (val - par->getMin()));
@@ -97,12 +99,15 @@ class FactoryTestParams {
 public:
    FactoryTestParams() = default;
    FactoryTestParams(std::string const &name, WorkspaceSetupFunc setupWorkspace, CreateNLLFunc createNLL,
-                     double fitResultTolerance, bool randomizeParameters)
+                     double fitResultTolerance, bool randomizeParameters, double hesseTolerance = -1.,
+                     bool enableHessian = true)
       : _name{name},
         _setupWorkspace{setupWorkspace},
         _createNLL{createNLL},
         _fitResultTolerance{fitResultTolerance},
-        _randomizeParameters{randomizeParameters}
+        _hesseTolerance{hesseTolerance},
+        _randomizeParameters{randomizeParameters},
+        _enableHessian{enableHessian}
    {
    }
 
@@ -110,7 +115,15 @@ class FactoryTestParams {
    WorkspaceSetupFunc _setupWorkspace;
    CreateNLLFunc _createNLL;
    double _fitResultTolerance = 1e-4;
+   // Tolerance for comparing parameter errors, which are less precise than
+   // the values because Minuits numeric Hesse is only accurate to about the
+   // percent level. If negative, _fitResultTolerance is used.
+   double _hesseTolerance = -1.;
    bool _randomizeParameters = true;
+   // Whether to test the minimization also with the analytic Hessian from
+   // Clad. Disabled for models where the Hessian is known to come out wrong
+   // or can't be generated at all (see comments at the test definitions).
+   bool _enableHessian = true;
 };
 
 class FactoryTest : public testing::TestWithParam<FactoryTestParams> {
@@ -129,14 +142,16 @@ class FactoryTest : public testing::TestWithParam<FactoryTestParams> {
    std::unique_ptr<RooHelpers::LocalChangeMsgLevel> _changeMsgLvl;
 };
 
-std::unique_ptr<RooFitResult> runMinimizer(RooAbsReal &absReal, bool useGradient = true)
+std::unique_ptr<RooFitResult> runMinimizer(RooAbsReal &absReal, bool useAD = true, bool useHessian = true)
 {
    RooMinimizer::Config cfg;
-   cfg.useGradient = useGradient;
+   cfg.useGradient = useAD;
+   cfg.useHessian = useAD && useHessian;
    RooMinimizer m{absReal, cfg};
    m.setPrintLevel(-1);
    m.setStrategy(0);
    m.minimize("Minuit2");
+   m.hesse(); // to use the analytical Hessian for error estimation
    return std::unique_ptr<RooFitResult>{m.save()};
 }
 
@@ -164,8 +179,13 @@ TEST_P(FactoryTest, NLLFit)
    // We want to use the generated code also for the nominal likelihood. Like
    // this, we make sure to validate also the NLL values of the generated code.
    static_cast<RooFit::Experimental::RooEvaluatorWrapper &>(*nllFunc).setUseGeneratedFunctionCode(true);
+   // Let's also validate the Hessian
+   if (_params._enableHessian) {
+      static_cast<RooFit::Experimental::RooEvaluatorWrapper &>(*nllFunc).generateHessian();
+   }
 
    double tol = _params._fitResultTolerance;
+   double hesseTol = _params._hesseTolerance < 0. ? tol : _params._hesseTolerance;
 
    EXPECT_NEAR(nllRef->getVal(observables), nllFunc->getVal(), tol);
 
@@ -202,7 +222,7 @@ TEST_P(FactoryTest, NLLFit)
    paramsRefNll.assign(parametersOrig);
 
    // Minimize the RooFuncWrapper Implementation with AD
-   auto resultAd = runMinimizer(*nllFunc);
+   auto resultAd = runMinimizer(*nllFunc, true, _params._enableHessian);
    paramsRefNll.assign(parametersOrig);
 
    // Minimize the reference NLL
@@ -214,8 +234,8 @@ TEST_P(FactoryTest, NLLFit)
    // because for very small correlations it's usually not the same within the
    // relative tolerance because you would compare two small values that are
    // only different from zero because of noise.
-   EXPECT_TRUE(result->isIdenticalNoCov(*resultRef, tol, tol));
-   EXPECT_TRUE(resultAd->isIdenticalNoCov(*resultRef, tol, tol));
+   EXPECT_TRUE(result->isIdenticalNoCov(*resultRef, tol, hesseTol));
+   EXPECT_TRUE(resultAd->isIdenticalNoCov(*resultRef, tol, hesseTol));
 }
 
 /// Initial minimization that was not based on any other tutorial/test.
@@ -223,7 +243,10 @@ FactoryTestParams param1{"Gaussian",
                          [](RooWorkspace &ws) {
                             constexpr double inf = std::numeric_limits<double>::infinity();
                             ws.import(RooRealVar{"mu", "mu", -inf, inf});
-                            ws.factory("sum::mu_shifted(mu, shift[1.0, -10, 10])");
+                            // The shift parameter needs to be constant, because otherwise the
+                            // model would only depend on the sum "mu + shift": the Hessian at the
+                            // minimum is singular and the parameter errors are meaningless.
+                            ws.factory("sum::mu_shifted(mu, shift[1.0])");
                             ws.factory("prod::sigma_scaled(sigma[3.0, 0.01, 10], 1.5)");
                             ws.factory("Gaussian::model(x[0, -10, 10], mu_shifted, sigma_scaled)");
 
@@ -234,7 +257,8 @@ FactoryTestParams param1{"Gaussian",
                             return std::unique_ptr<RooAbsReal>{pdf.createNLL(data, backend)};
                          },
                          1e-4,
-                         /*randomizeParameters=*/false};
+                         /*randomizeParameters=*/false,
+                         /*hesseTolerance=*/1e-3};
 
 /// Test based on the rf301 tutorial.
 FactoryTestParams param2{"PolyVar",
@@ -269,7 +293,8 @@ FactoryTestParams param4{"ConstraintSum",
                                pdf.createNLL(data, ExternalConstraints(*ws.pdf("fconstext")), backend)};
                          },
                          1e-4,
-                         /*randomizeParameters=*/true};
+                         /*randomizeParameters=*/true,
+                         /*hesseTolerance=*/1e-3};
 
 namespace {
 
@@ -341,12 +366,14 @@ void getSimPdfModel(RooWorkspace &ws)
 } // namespace
 
 /// Test based on the simultaneous fit shown in CHEP'23 results
-FactoryTestParams param5{"SimPdf", getSimPdfModel,
+FactoryTestParams param5{"SimPdf",
+                         getSimPdfModel,
                          [](RooAbsPdf &pdf, RooAbsData &data, RooWorkspace &, RooFit::EvalBackend backend) {
                             return std::unique_ptr<RooAbsReal>{pdf.createNLL(data, backend)};
                          },
                          5e-3,
-                         /*randomizeParameters=*/true};
+                         /*randomizeParameters=*/true,
+                         /*hesseTolerance=*/5e-2};
 
 FactoryTestParams param6{"GaussianExtended",
                          [](RooWorkspace &ws) {
@@ -388,12 +415,14 @@ void getDataHistModel(RooWorkspace &ws)
 } // namespace
 
 /// Test based on rf706 tutorial
-FactoryTestParams param7{"HistPdf", getDataHistModel,
+FactoryTestParams param7{"HistPdf",
+                         getDataHistModel,
                          [](RooAbsPdf &pdf, RooAbsData &data, RooWorkspace &, RooFit::EvalBackend backend) {
                             return std::unique_ptr<RooAbsReal>{pdf.createNLL(data, backend)};
                          },
                          1e-4,
-                         /*randomizeParameters=*/true};
+                         /*randomizeParameters=*/true,
+                         /*hesseTolerance=*/2e-3};
 
 FactoryTestParams param8{"Lognormal",
                          [](RooWorkspace &ws) {
@@ -404,7 +433,8 @@ FactoryTestParams param8{"Lognormal",
                             return std::unique_ptr<RooAbsReal>{pdf.createNLL(data, backend)};
                          },
                          1e-4,
-                         /*randomizeParameters=*/true};
+                         /*randomizeParameters=*/true,
+                         /*hesseTolerance=*/1e-2};
 
 FactoryTestParams param8p1{"LognormalStandard",
                            [](RooWorkspace &ws) {
@@ -416,7 +446,8 @@ FactoryTestParams param8p1{"LognormalStandard",
                               return std::unique_ptr<RooAbsReal>{pdf.createNLL(data, backend)};
                            },
                            3e-4,
-                           /*randomizeParameters=*/true};
+                           /*randomizeParameters=*/true,
+                           /*hesseTolerance=*/3e-2};
 
 FactoryTestParams param9{"Poisson",
                          [](RooWorkspace &ws) {
@@ -427,9 +458,10 @@ FactoryTestParams param9{"Poisson",
                             return std::unique_ptr<RooAbsReal>{pdf.createNLL(data, backend)};
                          },
                          1e-4,
-                         /*randomizeParameters=*/true};
+                         /*randomizeParameters=*/true,
+                         /*hesseTolerance=*/1e-2};
 
-// A RooPoisson where x is not rounded, like it is used in HistFactory
+// A RooPoisson where x is not rounded, like it is used in HistFactory.
 FactoryTestParams param10{"PoissonNoRounding",
                           [](RooWorkspace &ws) {
                              constexpr double inf = std::numeric_limits<double>::infinity();
@@ -444,7 +476,8 @@ FactoryTestParams param10{"PoissonNoRounding",
                              return std::unique_ptr<RooAbsReal>{pdf.createNLL(data, backend)};
                           },
                           1e-4,
-                          /*randomizeParameters=*/true};
+                          /*randomizeParameters=*/true,
+                          /*hesseTolerance=*/1e-2};
 
 FactoryTestParams param11{"ClassFactory1D",
                           [](RooWorkspace &ws) {
@@ -503,10 +536,14 @@ void func_1_grad(double const * /*x*/, double *grad_out)
    grad_out[0] = 1;
    grad_out[1] = -1;
 }
+void func_1_hess(double const * /*x*/, double *hess_out)
+{
+   std::fill(hess_out, hess_out + 4, 0.0);
+}
 
 auto &functor_1()
 {
-   static ROOT::Math::GradFunctor functor{func_1, 2, func_1_grad};
+   static ROOT::Math::GradFunctor functor{func_1, 2, func_1_grad, func_1_hess};
    return functor;
 }
 
@@ -530,9 +567,35 @@ void func_gaussian_grad(double const *x, double *grad_out)
    grad_out[2] = f * arg * arg * inv_sig2 / sig;
 }
 
+void func_gaussian_hess(double const *x, double *hess_out)
+{
+   const double arg = x[0] - x[1];
+   const double sig = x[2];
+
+   const double invSig2 = 1.0 / (sig * sig);
+   const double t = arg * arg * invSig2;
+   const double f = std::exp(-0.5 * t);
+
+   const double h00 = f * (t - 1.0) * invSig2;
+   const double h02 = f * arg * (2.0 - t) * invSig2 / sig;
+   const double h22 = f * t * (t - 3.0) * invSig2;
+
+   // The function is invariant under exchanging x[0] - x[1] with x[1] - x[0],
+   // so the Hessian rows/columns for x[0] and x[1] only differ in sign.
+   hess_out[0] = h00;
+   hess_out[1] = -h00;
+   hess_out[2] = h02;
+   hess_out[3] = -h00;
+   hess_out[4] = h00;
+   hess_out[5] = -h02;
+   hess_out[6] = h02;
+   hess_out[7] = -h02;
+   hess_out[8] = h22;
+}
+
 auto &functor_gaussian()
 {
-   static ROOT::Math::GradFunctor functor{func_gaussian, 3, func_gaussian_grad};
+   static ROOT::Math::GradFunctor functor{func_gaussian, 3, func_gaussian_grad, func_gaussian_hess};
    return functor;
 }
 
@@ -544,10 +607,14 @@ double func_1_1D_diff(double /*x*/)
 {
    return 1.;
 }
+double func_1_1D_diff2(double /*x*/)
+{
+   return 0.;
+}
 
 auto &functor_1_1D()
 {
-   static ROOT::Math::GradFunctor1D functor{func_1_1D, func_1_1D_diff};
+   static ROOT::Math::GradFunctor1D functor{func_1_1D, func_1_1D_diff, func_1_1D_diff2};
    return functor;
 }
 
@@ -556,12 +623,17 @@ auto &functor_1_1D()
 FactoryTestParams param13{"RooFunctor",
                           [](RooWorkspace &ws) {
                              RooRealVar x("x", "", 0.0, -4, 4);
-                             RooRealVar mu("mu", "", 0.0, -4, 4);
+                             // The initial mu value is away from zero such that relative
+                             // comparisons of the fitted values are meaningful.
+                             RooRealVar mu("mu", "", 2.0, -4, 4);
                              RooRealVar shift("shift", "", 1.0, -4, 4);
                              shift.setConstant(true);
                              RooFunctorBinding mu_shifted("mu_shifted", "", functor_1(), {mu, shift});
                              RooFunctor1DBinding mu_shifted_1D("mu_shifted_1D", "", functor_1_1D(), {mu});
-                             RooRealVar sigma("sigma", "", 4., 0.01, 10.);
+                             // The initial sigma value is small enough for the Gaussian to be
+                             // well-contained in the observable range, such that the generated
+                             // data constrains all parameters and the fits are well-defined.
+                             RooRealVar sigma("sigma", "", 1.5, 0.01, 10.);
 
                              RooFunctorPdfBinding gauss_1("model_1", "", functor_gaussian(), {x, mu_shifted, sigma});
                              RooFunctorPdfBinding gauss_2("model_2", "", functor_gaussian(), {x, mu_shifted_1D, sigma});
@@ -572,15 +644,27 @@ FactoryTestParams param13{"RooFunctor",
 
                              ws.import(model);
                              ws.defineSet("observables", vars);
+
+                             // The model is statistically very flat around the minimum, so with
+                             // only the default 100 events, the fitted values with and without
+                             // AD can differ at the percent level just because the minimizers
+                             // stop at slightly different points in the flat valley. Generate
+                             // more events to make the minimum well-defined.
+                             std::unique_ptr<RooDataSet> data{ws.pdf("model")->generate(*ws.set("observables"), 1000)};
+                             std::unique_ptr<RooAbsData> binned{data->binnedClone()};
+                             binned->SetName("data");
+                             ws.import(*binned);
                           },
                           [](RooAbsPdf &pdf, RooAbsData &data, RooWorkspace &, RooFit::EvalBackend backend) {
                              return std::unique_ptr<RooAbsReal>{pdf.createNLL(data, backend)};
                           },
                           1e-3,
-                          /*randomizeParameters=*/true};
+                          /*randomizeParameters=*/true,
+                          /*hesseTolerance=*/1e-2};
 
 FactoryTestParams makeTestParams(const char *name, std::vector<std::string> const &expressions,
-                                 double fitResultTolerance, bool randomizeParameters = true)
+                                 double fitResultTolerance, bool randomizeParameters = true,
+                                 double hesseTolerance = -1., bool enableHessian = true, std::size_t nEvents = 0)
 {
    return {name,
            [=](RooWorkspace &ws) {
@@ -588,12 +672,23 @@ FactoryTestParams makeTestParams(const char *name, std::vector<std::string> cons
                  ws.factory(expr.c_str());
               }
               ws.defineSet("observables", "x");
+              // Provide a custom dataset if the default number of events that
+              // the test generates otherwise is not right for this model.
+              if (nEvents > 0) {
+                 std::unique_ptr<RooDataSet> data{ws.pdf("model")->generate(*ws.set("observables"), nEvents)};
+                 std::unique_ptr<RooAbsData> binned{data->binnedClone()};
+                 binned->SetName("data");
+                 ws.import(*binned);
+              }
            },
            [](RooAbsPdf &pdf, RooAbsData &data, RooWorkspace &, RooFit::EvalBackend backend) {
               using namespace RooFit;
               return std::unique_ptr<RooAbsReal>{pdf.createNLL(data, backend)};
            },
-           fitResultTolerance, randomizeParameters};
+           fitResultTolerance,
+           randomizeParameters,
+           hesseTolerance,
+           enableHessian};
 }
 
 auto testValues = testing::Values(
@@ -601,12 +696,15 @@ auto testValues = testing::Values(
    makeTestParams(
       "BifurGauss",
       {"x[0, -10, 10]", "mu[0, -10, 10]", "BifurGauss::model(x, mu, sigmaL[3.0, 0.01, 10], sigmaR[2.0, 0.01, 10])"},
-      1e-4, false),
+      1e-4, false, /*hesseTolerance=*/2e-2),
+   // The shift parameter needs to be constant, because otherwise the model
+   // would only depend on the sum "mu + shift": the Hessian at the minimum is
+   // singular and the parameter errors are meaningless.
    makeTestParams("RooFormulaVar",
-                  {"expr::mu_shifted('mu+shift',{mu[0, -10, 10], shift[1.0, -10, 10]})",
+                  {"expr::mu_shifted('mu+shift',{mu[0, -10, 10], shift[1.0]})",
                    "expr::sigma_scaled('sigma*1.5',{sigma[3.0, 0.01, 10]})",
                    "Gaussian::model(x[0, -10, 10], mu_shifted, sigma_scaled)"},
-                  1e-4, false),
+                  1e-4, false, /*hesseTolerance=*/5e-3),
    // Test for the uniform pdf. Since it doesn't depend on any parameters, we need
    // to add it to some other model like a Gaussian to get a meaningful fit model.
    makeTestParams("Uniform",
@@ -620,22 +718,39 @@ auto testValues = testing::Values(
                    "Gaussian::sig4(x, 6.0, sigma4[2.0, .01, 10])",
                    "RecursiveFraction::recfrac({a1[0.25, 0.0, 1.0], a2[0.25, 0.0, 1.0]})",
                    "SUM::model(a1 * sig1, a2 * sig2, recfrac * sig3, sig4)"},
-                  5e-3, true),
+                  5e-3, true, /*hesseTolerance=*/3e-2),
    makeTestParams("RooCBShape",
                   {"x[0., -200., 200.]", "x0[100., -200., 200.]",
                    "CBShape::model(x, x0, sigma[2., 1.E-1, 100.], alpha[1., 1.E-1, 100.], n[1., 1.E-1, 100.])"},
                   6e-3, true),
-   makeTestParams("RooBernstein",
-                  {"Bernstein::model(x[0., 100.], {c0[0.3, 0., 10.], c1[0.7, 0., 10.], c2[0.2, 0., 10.]})"}, 6e-3,
-                  true),
+   // The first Bernstein coefficient needs to be constant, because the pdf is
+   // invariant under a common rescaling of all coefficients: the Hessian at
+   // the minimum is singular and the parameter errors are meaningless.
+   makeTestParams("RooBernstein", {"Bernstein::model(x[0., 100.], {c0[0.3], c1[0.7, 0., 10.], c2[0.2, 0., 10.]})"},
+                  6e-3, true,
+                  /*hesseTolerance=*/2e-2),
    // We're testing several Landau configurations, because the underlying
    // ROOT::Math::landau_cdf is defined piecewise. Like this, we're covering
    // all possible code paths in the pullback.
-   makeTestParams("RooLandau1", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[1., 0.01, 50.])"}, 7e-3, false),
-   makeTestParams("RooLandau2", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[2.1, 0.01, 50.])"}, 7e-3, false),
-   makeTestParams("RooLandau3", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[10., 0.01, 50.])"}, 7e-3, false),
-   makeTestParams("RooLandau4", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[0.3, 0.01, 50.])"}, 7e-3, false),
-   makeTestParams("RooLandau5", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[0.07, 0.01, 50.])"}, 7e-3, false),
+   makeTestParams("RooLandau1", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[1., 0.01, 50.])"}, 7e-3, false,
+                  /*hesseTolerance=*/1e-2),
+   makeTestParams("RooLandau2", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[2.1, 0.01, 50.])"}, 7e-3, false,
+                  /*hesseTolerance=*/1e-2),
+   // With sl = 10, the Landau is nearly uniform over the observable window,
+   // so with only 100 events the fit is unidentifiable: the minimizers stop
+   // at arbitrary points in an almost flat valley, and the fitted values
+   // can't be meaningfully compared. Generate more events to make the
+   // minimum well-defined.
+   // The Hesse tolerance is wider than for the other Landau tests: with sl =
+   // 10, the model is poorly conditioned and Minuits numeric Hesse in the
+   // reference fit is less accurate (the Clad Hessian was validated against
+   // finite differences of the exact gradient).
+   makeTestParams("RooLandau3", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[10., 0.01, 50.])"}, 7e-3, false,
+                  /*hesseTolerance=*/3e-2, /*enableHessian=*/true, /*nEvents=*/10000),
+   makeTestParams("RooLandau4", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[0.3, 0.01, 50.])"}, 7e-3, false,
+                  /*hesseTolerance=*/1e-2),
+   makeTestParams("RooLandau5", {"Landau::model(x[5., 0., 30.], ml[6., 1., 30.], sl[0.07, 0.01, 50.])"}, 7e-3, false,
+                  /*hesseTolerance=*/1e-2),
    makeTestParams(
       "RooRealSumPdf1",
       {"Gaussian::gx(x[-10,10],m[0],1.0)", "Chebychev::ch(x,{0.1,0.2,-0.3})", "RealSumPdf::model({gx, ch}, {f[0,1]})"},