diff --git a/roofit/roofit/src/RooMomentMorphFuncND.cxx b/roofit/roofit/src/RooMomentMorphFuncND.cxx
index e4e2386353fa6..cbf42aecb28cf 100644
--- a/roofit/roofit/src/RooMomentMorphFuncND.cxx
+++ b/roofit/roofit/src/RooMomentMorphFuncND.cxx
@@ -850,12 +850,31 @@ void RooMomentMorphFuncND::findShape(const vector<double> &x) const
    _squareVec = output;
 
    for (int isq = 0; isq < depth; isq++) {
+      // Reset to a sentinel before searching: if no matching reference point is
+      // found below, we must not silently keep the value left over from a
+      // previous, unrelated call to findShape() (e.g. from a different
+      // hypercube visited earlier in the fit). Falling through to a stale
+      // index here would make the interpolation depend on evaluation history
+      // instead of only on the current parameter point.
+      _squareIdx[isq] = -1;
       for (int iref = 0; iref < nRef; iref++) {
          if (_squareVec[isq] == _referenceGrid._nref[iref]) {
             _squareIdx[isq] = iref;
             break;
          }
       }
+      if (_squareIdx[isq] < 0) {
+         coutE(InputArguments) << "RooMomentMorphFuncND::findShape(" << GetName()
+                               << ") ERROR: no reference pdf found for grid corner (";
+         for (unsigned int ix = 0; ix < _squareVec[isq].size(); ++ix) {
+            ccoutE(InputArguments) << (ix ? ", " : "") << _squareVec[isq][ix];
+         }
+         ccoutE(InputArguments) << ") of the hypercube enclosing the current morphing "
+                                << "parameter point. The reference grid is missing a pdf at this "
+                                << "coordinate -- check that RooMomentMorphFuncND::Grid::addPdf() was "
+                                << "called for every corner of the parameter range." << std::endl;
+         throw string("RooMomentMorphFuncND::findShape() ERROR: incomplete reference grid");
+      }
    }
 
    // std::cout << std::endl;
diff --git a/tutorials/roofit/roofit/rf617_simulation_based_inference_multidimensional.py b/tutorials/roofit/roofit/rf617_simulation_based_inference_multidimensional.py
index fc34d2242e9e9..65783f221fc0b 100644
--- a/tutorials/roofit/roofit/rf617_simulation_based_inference_multidimensional.py
+++ b/tutorials/roofit/roofit/rf617_simulation_based_inference_multidimensional.py
@@ -44,7 +44,7 @@
 def morphing(setting, n_dimensions):
     # Define binning for morphing
 
-    binning = [ROOT.RooBinning(n_bins, 0.0, n_bins - 1.0) for dim in range(n_dimensions)]
+    binning = [ROOT.RooBinning(n_bins - 1, 0.0, n_bins - 1.0) for dim in range(n_dimensions)]
     grid = ROOT.RooMomentMorphFuncND.Grid(*binning)
 
     # Set bins for each x variable
@@ -202,8 +202,7 @@ def build_ws(mu_observed):
 morphing(ROOT.RooMomentMorphFuncND.Linear, len(mu_observed))
 
 # Calculate the nll for the moprhed distribution
-# TODO: Fix RooAddPdf::fixCoefNormalization(nset) warnings with new CPU backend
-nll_morph = ws["morph"].createNLL(ws["obs_data"], EvalBackend="legacy")
+nll_morph = ws["morph"].createNLL(ws["obs_data"])
 
 # Initialize the SBI model
 model = SBI(ws, len(mu_observed))