diff --git a/roofit/histfactory/inc/RooStats/HistFactory/ParamHistFunc.h b/roofit/histfactory/inc/RooStats/HistFactory/ParamHistFunc.h
index 07dbbfd48a1d9..7efad55d99417 100644
--- a/roofit/histfactory/inc/RooStats/HistFactory/ParamHistFunc.h
+++ b/roofit/histfactory/inc/RooStats/HistFactory/ParamHistFunc.h
@@ -48,7 +48,7 @@ class ParamHistFunc : public RooAbsReal {
 
   RooDataHist const& dataHist() const { return _dataSet; }
 
-  double binVolume() const { return _dataSet.binVolume(); }
+  double binVolume(std::size_t i) const { return _dataSet.binVolume(i); }
 
   bool forceAnalyticalInt(const RooAbsArg&) const override { return true ; }
 
diff --git a/roofit/histfactory/src/RooBarlowBeestonLL.cxx b/roofit/histfactory/src/RooBarlowBeestonLL.cxx
index d5bdb47cb6e4a..94251d16ef722 100644
--- a/roofit/histfactory/src/RooBarlowBeestonLL.cxx
+++ b/roofit/histfactory/src/RooBarlowBeestonLL.cxx
@@ -192,7 +192,7 @@ void RooStats::HistFactory::RooBarlowBeestonLL::initializeBarlowCache() {
       cache.bin_center = bin_center;
       cache.observables = obsSet;
 
-      cache.binVolume = param_func->binVolume();
+      cache.binVolume = param_func->binVolume(bin_index);
 
       // Delete this part, simply a comment
       RooArgList obs_list( *(cache.bin_center) );
diff --git a/roofit/roofit/src/RooHistConstraint.cxx b/roofit/roofit/src/RooHistConstraint.cxx
index fcff8e9307f4e..80cfc3f7113d9 100644
--- a/roofit/roofit/src/RooHistConstraint.cxx
+++ b/roofit/roofit/src/RooHistConstraint.cxx
@@ -60,11 +60,11 @@ RooHistConstraint::RooHistConstraint(const char *name, const char *title,
 
     // Now populate nominal with parameters
     for (int i=0 ; i<phf->_dh.numEntries() ; i++) {
-      phf->_dh.get(i) ;
-      if (phf->_dh.weight()<threshold && phf->_dh.weight() != 0.) {
+      const double wi = phf->_dh.weight(i) ;
+      if (wi<threshold && wi != 0.) {
         const char* vname = Form("%s_nominal_bin_%i",GetName(),i) ;
         auto var = std::make_unique<RooRealVar>(vname,vname,0,1.E30);
-        var->setVal(phf->_dh.weight()) ;
+        var->setVal(wi) ;
         var->setConstant(true);
 
         auto gamma = static_cast<RooRealVar*>(phf->_p.at(i));
diff --git a/roofit/roofit/src/RooIntegralMorph.cxx b/roofit/roofit/src/RooIntegralMorph.cxx
index 9e2e1cc6cfdcb..1bc8519f7adee 100644
--- a/roofit/roofit/src/RooIntegralMorph.cxx
+++ b/roofit/roofit/src/RooIntegralMorph.cxx
@@ -392,8 +392,8 @@ void RooIntegralMorph::MorphCacheElem::calculate(TIterator* dIter)
   // Zero output histogram below lowest calculable X value
   for (int i=0; i<_yatXmin ; i++) {
     dIter->Next() ;
-    //_hist->get(i) ;
-    hist()->set(0) ;
+    const std::size_t binIdx = hist()->getIndex(*hist()->get(), /*fast=*/true);
+    hist()->set(binIdx, 0, -1) ;
   }
 
   double x1 = _x->getMin("cache");
@@ -420,14 +420,16 @@ void RooIntegralMorph::MorphCacheElem::calculate(TIterator* dIter)
     double fbarX = f1x1*f2x2 / ( _alpha->getVal()*f2x2 + (1-_alpha->getVal())*f1x1 ) ;
 
     dIter->Next() ;
-    //_hist->get(i) ;
-    hist()->set(fbarX) ;
+    {
+      const std::size_t binIdx = hist()->getIndex(*hist()->get(), /*fast=*/true);
+      hist()->set(binIdx, fbarX, -1) ;
+    }
   }
   // Zero output histogram above highest calculable X value
   for (int i=_yatXmax+1 ; i<nbins ; i++) {
     dIter->Next() ;
-    //_hist->get(i) ;
-    hist()->set(0) ;
+    const std::size_t binIdx = hist()->getIndex(*hist()->get(), /*fast=*/true);
+    hist()->set(binIdx, 0, -1) ;
   }
 
   pdf()->setUnitNorm(true) ;
diff --git a/roofit/roofit/src/RooLagrangianMorphFunc.cxx b/roofit/roofit/src/RooLagrangianMorphFunc.cxx
index c380ccb006221..94473baec0a1f 100644
--- a/roofit/roofit/src/RooLagrangianMorphFunc.cxx
+++ b/roofit/roofit/src/RooLagrangianMorphFunc.cxx
@@ -2572,12 +2572,12 @@ TH1 *RooLagrangianMorphFunc::createTH1(const std::string &name, bool correlateEr
             continue;
          }
          const RooDataHist &dhist = hf->dataHist();
-         dhist.get(i);
          RooAbsReal *formula = dynamic_cast<RooAbsReal *>(prod->components().find(Form("w_%s", prod->GetName())));
          double weight = formula->getVal();
-         unc2 += dhist.weightSquared() * weight * weight;
-         unc += sqrt(dhist.weightSquared()) * weight;
-         val += dhist.weight() * weight;
+         const double w2 = dhist.weightSquared(i);
+         unc2 += w2 * weight * weight;
+         unc += sqrt(w2) * weight;
+         val += dhist.weight(i) * weight;
       }
       hist->SetBinContent(i + 1, val);
       hist->SetBinError(i + 1, correlateErrors ? unc : sqrt(unc2));
@@ -2827,8 +2827,7 @@ double RooLagrangianMorphFunc::expectedUncertainty() const
       if (hf) {
          const RooDataHist &hist = hf->dataHist();
          for (Int_t j = 0; j < observable->getBins(); ++j) {
-            hist.get(j);
-            newunc2 += hist.weightSquared();
+            newunc2 += hist.weightSquared(j);
          }
       } else if (rv) {
          newunc2 = pow(rv->getError(), 2);
diff --git a/roofit/roofit/src/RooParamHistFunc.cxx b/roofit/roofit/src/RooParamHistFunc.cxx
index 2890c083c2c45..3f37cd548ece1 100644
--- a/roofit/roofit/src/RooParamHistFunc.cxx
+++ b/roofit/roofit/src/RooParamHistFunc.cxx
@@ -44,11 +44,11 @@ RooParamHistFunc::RooParamHistFunc(const char *name, const char *title, RooDataH
    // Now populate p with parameters
    RooArgSet allVars;
    for (Int_t i = 0; i < _dh.numEntries(); i++) {
-      _dh.get(i);
+      const double wi = _dh.weight(i);
       const char *vname = Form("%s_gamma_bin_%i", GetName(), i);
       RooRealVar *var = new RooRealVar(vname, vname, 0, 1000);
-      var->setVal(_relParam ? 1 : _dh.weight());
-      var->setError(_relParam ? 1 / sqrt(_dh.weight()) : sqrt(_dh.weight()));
+      var->setVal(_relParam ? 1 : wi);
+      var->setError(_relParam ? 1 / sqrt(wi) : sqrt(wi));
       var->setConstant(true);
       allVars.add(*var);
       _p.add(*var);
@@ -94,8 +94,7 @@ void RooParamHistFunc::setActual(Int_t ibin, double newVal)
 
 double RooParamHistFunc::getNominal(Int_t ibin) const
 {
-  _dh.get(ibin) ;
-  return _dh.weight() ;
+  return _dh.weight(ibin) ;
 }
 
 ////////////////////////////////////////////////////////////////////////////////
diff --git a/roofit/roofitcore/inc/RooEvaluatorWrapper.h b/roofit/roofitcore/inc/RooEvaluatorWrapper.h
index dea80a4179602..ce394fc629152 100644
--- a/roofit/roofitcore/inc/RooEvaluatorWrapper.h
+++ b/roofit/roofitcore/inc/RooEvaluatorWrapper.h
@@ -21,8 +21,10 @@
 #include <RooRealProxy.h>
 #include <RooSetProxy.h>
 
+#include <memory>
 #include <stack>
 
+class ChangeOperModeRAII;
 class RooAbsArg;
 class RooAbsCategory;
 class RooAbsPdf;
@@ -71,9 +73,10 @@ class RooEvaluatorWrapper final : public RooAbsReal {
    void generateHessian();
 
    void setUseGeneratedFunctionCode(bool);
-
    void writeDebugMacro(std::string const &) const;
 
+   std::unique_ptr<ChangeOperModeRAII> setOperModes(RooAbsArg::OperMode opMode);
+
 protected:
    double evaluate() const override;
 
diff --git a/roofit/roofitcore/inc/RooFit/Evaluator.h b/roofit/roofitcore/inc/RooFit/Evaluator.h
index 281e4ab7213f3..f50424c634abf 100644
--- a/roofit/roofitcore/inc/RooFit/Evaluator.h
+++ b/roofit/roofitcore/inc/RooFit/Evaluator.h
@@ -20,7 +20,6 @@
 #include <RConfig.h>
 
 #include <memory>
-#include <stack>
 
 class ChangeOperModeRAII;
 class RooAbsArg;
@@ -45,6 +44,8 @@ class Evaluator {
 
    void setOffsetMode(RooFit::EvalContext::OffsetMode);
 
+   std::unique_ptr<ChangeOperModeRAII> setOperModes(RooAbsArg::OperMode opMode);
+
 private:
    void processVariable(NodeInfo &nodeInfo);
    void processCategory(NodeInfo &nodeInfo);
@@ -66,7 +67,7 @@ class Evaluator {
    RooFit::EvalContext _evalContextCUDA;
    std::vector<NodeInfo> _nodes;                             // the ordered computation graph
    std::unordered_map<TNamed const *, NodeInfo *> _nodesMap; // for quick lookup of nodes
-   std::stack<std::unique_ptr<ChangeOperModeRAII>> _changeOperModeRAIIs;
+   std::unique_ptr<ChangeOperModeRAII> _operModeChanges;
 };
 
 } // end namespace RooFit
diff --git a/roofit/roofitcore/inc/RooMinimizer.h b/roofit/roofitcore/inc/RooMinimizer.h
index fcf2dfe960a64..95ff438588f3e 100644
--- a/roofit/roofitcore/inc/RooMinimizer.h
+++ b/roofit/roofitcore/inc/RooMinimizer.h
@@ -234,6 +234,7 @@ class RooMinimizer : public TObject {
    void fillCorrMatrix(RooFitResult &fitRes);
    void updateErrors();
 
+   RooAbsReal &_function;
    ROOT::Fit::FitConfig _config;                      ///< fitter configuration (options and parameter settings)
    std::unique_ptr<FitResult> _result;                ///<! pointer to the object containing the result of the fit
    std::unique_ptr<ROOT::Math::Minimizer> _minimizer; ///<! pointer to used minimizer
diff --git a/roofit/roofitcore/res/RooFitImplHelpers.h b/roofit/roofitcore/res/RooFitImplHelpers.h
index 6d04134bc85ec..d14bc905cd21e 100644
--- a/roofit/roofitcore/res/RooFitImplHelpers.h
+++ b/roofit/roofitcore/res/RooFitImplHelpers.h
@@ -46,25 +46,51 @@ class DisableCachingRAII {
    bool _oldState;
 };
 
-/// Struct to temporarily change the operation mode of a RooAbsArg until it
-/// goes out of scope.
+/// Scope guard that temporarily changes the operation mode of one or more
+/// RooAbsArg instances. Each call to change() records the arg's current
+/// operMode before flipping it to the requested mode (non-recursively, i.e.
+/// value clients are not touched). Destruction (or an explicit clear())
+/// restores every recorded mode in LIFO order.
+///
+/// The class is movable but not copyable, so it can be returned from
+/// functions that build up a batch of changes to hand to the caller.
 class ChangeOperModeRAII {
 public:
-   ChangeOperModeRAII(RooAbsArg *arg, RooAbsArg::OperMode opMode) : _arg{arg}, _oldOpMode(arg->operMode())
+   ChangeOperModeRAII() = default;
+
+   /// Convenience ctor: behaves like a scope guard for a single arg.
+   ChangeOperModeRAII(RooAbsArg *arg, RooAbsArg::OperMode opMode) { change(arg, opMode); }
+
+   ~ChangeOperModeRAII() { clear(); }
+
+   ChangeOperModeRAII(ChangeOperModeRAII &&) = default;
+   ChangeOperModeRAII &operator=(ChangeOperModeRAII &&) = default;
+   ChangeOperModeRAII(ChangeOperModeRAII const &) = delete;
+   ChangeOperModeRAII &operator=(ChangeOperModeRAII const &) = delete;
+
+   /// Record arg's current operMode and flip it to opMode. If the current
+   /// mode already equals opMode, this is a no-op (nothing to restore).
+   void change(RooAbsArg *arg, RooAbsArg::OperMode opMode)
    {
-      arg->setOperMode(opMode, /*recurse=*/false);
+      if (opMode == arg->operMode())
+         return;
+      _entries.emplace_back(arg, arg->operMode());
+      arg->setOperMode(opMode, /*recurseADirty=*/false);
    }
 
-   ChangeOperModeRAII(ChangeOperModeRAII const &other) = delete;
-   ChangeOperModeRAII &operator=(ChangeOperModeRAII const &other) = delete;
-   ChangeOperModeRAII(ChangeOperModeRAII &&other) = delete;
-   ChangeOperModeRAII &operator=(ChangeOperModeRAII &&other) = delete;
+   /// Restore every recorded change right away, emptying this guard.
+   void clear()
+   {
+      for (auto it = _entries.rbegin(); it != _entries.rend(); ++it) {
+         it->first->setOperMode(it->second, /*recurseADirty=*/false);
+      }
+      _entries.clear();
+   }
 
-   ~ChangeOperModeRAII() { _arg->setOperMode(_oldOpMode, /*recurse=*/false); }
+   bool empty() const { return _entries.empty(); }
 
 private:
-   RooAbsArg *_arg = nullptr;
-   RooAbsArg::OperMode _oldOpMode;
+   std::vector<std::pair<RooAbsArg *, RooAbsArg::OperMode>> _entries;
 };
 
 namespace RooHelpers {
diff --git a/roofit/roofitcore/src/RooAbsPdf.cxx b/roofit/roofitcore/src/RooAbsPdf.cxx
index 7aad2f55d3b9f..296649b6aae4e 100644
--- a/roofit/roofitcore/src/RooAbsPdf.cxx
+++ b/roofit/roofitcore/src/RooAbsPdf.cxx
@@ -1703,26 +1703,27 @@ RooFit::OwningPtr<RooDataHist> RooAbsPdf::generateBinned(const RooArgSet &whatVa
   Int_t histOutSum(0) ;
   for (int i=0 ; i<hist->numEntries() ; i++) {
     hist->get(i) ;
+    const double wi = hist->weight(i) ;
     if (expectedData) {
 
       // Expected data, multiply p.d.f by nEvents
-      double w=hist->weight()*nEvents ;
+      double w=wi*nEvents ;
       hist->set(i, w, sqrt(w));
 
     } else if (extended) {
 
       // Extended mode, set contents to Poisson(pdf*nEvents)
-      double w = RooRandom::randomGenerator()->Poisson(hist->weight()*nEvents) ;
-      hist->set(w,sqrt(w)) ;
+      double w = RooRandom::randomGenerator()->Poisson(wi*nEvents) ;
+      hist->set(i, w, sqrt(w)) ;
 
     } else {
 
       // Regular mode, fill array of weights with Poisson(pdf*nEvents), but to not fill
       // histogram yet.
-      if (hist->weight()>histMax) {
-        histMax = hist->weight() ;
+      if (wi>histMax) {
+        histMax = wi ;
       }
-      histOut[i] = RooRandom::randomGenerator()->Poisson(hist->weight()*nEvents) ;
+      histOut[i] = RooRandom::randomGenerator()->Poisson(wi*nEvents) ;
       histOutSum += histOut[i] ;
     }
   }
@@ -1745,7 +1746,7 @@ RooFit::OwningPtr<RooDataHist> RooAbsPdf::generateBinned(const RooArgSet &whatVa
       hist->get(ibinRand) ;
       double ranY = RooRandom::randomGenerator()->Uniform(histMax) ;
 
-      if (ranY<hist->weight()) {
+      if (ranY<hist->weight(ibinRand)) {
         if (wgt==1) {
           histOut[ibinRand]++ ;
         } else {
@@ -1768,8 +1769,7 @@ RooFit::OwningPtr<RooDataHist> RooAbsPdf::generateBinned(const RooArgSet &whatVa
 
     // Transfer working array to histogram
     for (int i=0 ; i<hist->numEntries() ; i++) {
-      hist->get(i) ;
-      hist->set(histOut[i],sqrt(1.0*histOut[i])) ;
+      hist->set(i, histOut[i], sqrt(1.0*histOut[i])) ;
     }
 
   } else if (expectedData) {
@@ -1779,8 +1779,8 @@ RooFit::OwningPtr<RooDataHist> RooAbsPdf::generateBinned(const RooArgSet &whatVa
     // bin average and bin integral in sampling bins
     double corr = nEvents/hist->sumEntries() ;
     for (int i=0 ; i<hist->numEntries() ; i++) {
-      hist->get(i) ;
-      hist->set(hist->weight()*corr,sqrt(hist->weight()*corr)) ;
+      const double wnew = hist->weight(i)*corr ;
+      hist->set(i, wnew, sqrt(wnew)) ;
     }
 
   }
diff --git a/roofit/roofitcore/src/RooAbsReal.cxx b/roofit/roofitcore/src/RooAbsReal.cxx
index b6e5a1f901749..da039f60240cd 100644
--- a/roofit/roofitcore/src/RooAbsReal.cxx
+++ b/roofit/roofitcore/src/RooAbsReal.cxx
@@ -1151,7 +1151,7 @@ RooDataHist* RooAbsReal::fillDataHist(RooDataHist *hist, const RooArgSet* normSe
     const RooArgSet* obs = hist->get(i) ;
     double binVal = theClone->getVal(normSet?normSet:obs)*scaleFactor ;
     if (correctForBinSize) {
-      binVal*= hist->binVolume() ;
+      binVal*= hist->binVolume(i) ;
     }
     hist->set(i, binVal, 0.);
   }
diff --git a/roofit/roofitcore/src/RooBinnedGenContext.cxx b/roofit/roofitcore/src/RooBinnedGenContext.cxx
index dbcc5652e1bf3..e844a59dc6148 100644
--- a/roofit/roofitcore/src/RooBinnedGenContext.cxx
+++ b/roofit/roofitcore/src/RooBinnedGenContext.cxx
@@ -141,27 +141,28 @@ RooDataSet *RooBinnedGenContext::generate(double nEvt, bool /*skipInit*/, bool e
   double histMax(-1) ;
   Int_t histOutSum(0) ;
   for (int i=0 ; i<_hist->numEntries() ; i++) {
-    _hist->get(i) ;
+    const RooArgSet* coords = _hist->get(i) ;
+    const double wi = _hist->weight(i) ;
     if (_expectedData) {
 
       // Expected data, multiply p.d.f by nEvents
-      double w=_hist->weight()*nEvents ;
-      wudata->add(*_hist->get(),w) ;
+      double w=wi*nEvents ;
+      wudata->add(*coords,w) ;
 
     } else if (extended) {
 
       // Extended mode, set contents to Poisson(pdf*nEvents)
-      double w = RooRandom::randomGenerator()->Poisson(_hist->weight()*nEvents) ;
-      wudata->add(*_hist->get(),w) ;
+      double w = RooRandom::randomGenerator()->Poisson(wi*nEvents) ;
+      wudata->add(*coords,w) ;
 
     } else {
 
       // Regular mode, fill array of weights with Poisson(pdf*nEvents), but to not fill
       // histogram yet.
-      if (_hist->weight()>histMax) {
-   histMax = _hist->weight() ;
+      if (wi>histMax) {
+   histMax = wi ;
       }
-      histOut[i] = RooRandom::randomGenerator()->Poisson(_hist->weight()*nEvents) ;
+      histOut[i] = RooRandom::randomGenerator()->Poisson(wi*nEvents) ;
       histOutSum += histOut[i] ;
     }
   }
@@ -182,7 +183,7 @@ RooDataSet *RooBinnedGenContext::generate(double nEvt, bool /*skipInit*/, bool e
       _hist->get(ibinRand) ;
       double ranY = RooRandom::randomGenerator()->Uniform(histMax) ;
 
-      if (ranY<_hist->weight()) {
+      if (ranY<_hist->weight(ibinRand)) {
    if (wgt==1) {
      histOut[ibinRand]++ ;
    } else {
diff --git a/roofit/roofitcore/src/RooChi2Var.cxx b/roofit/roofitcore/src/RooChi2Var.cxx
index b6df3e207403b..89020e32ece56 100644
--- a/roofit/roofitcore/src/RooChi2Var.cxx
+++ b/roofit/roofitcore/src/RooChi2Var.cxx
@@ -107,9 +107,9 @@ double RooChi2Var::evaluatePartition(std::size_t firstEvent, std::size_t lastEve
     }
     if (!doSelect) continue ;
 
-    const double nData = hdata->weight() ;
+    const double nData = hdata->weight(i) ;
 
-    const double nPdf = _funcClone->getVal(_normSet) * normFactor * hdata->binVolume() ;
+    const double nPdf = _funcClone->getVal(_normSet) * normFactor * hdata->binVolume(i) ;
 
     const double eExt = nPdf-nData ;
 
diff --git a/roofit/roofitcore/src/RooDataHist.cxx b/roofit/roofitcore/src/RooDataHist.cxx
index 5d6be15717fd8..c54b845d9fd88 100644
--- a/roofit/roofitcore/src/RooDataHist.cxx
+++ b/roofit/roofitcore/src/RooDataHist.cxx
@@ -660,7 +660,7 @@ void RooDataHist::importDHistSet(const RooArgList & /*vars*/, RooCategory &index
       // Transfer contents
       for (Int_t i = 0; i < dhist->numEntries(); i++) {
          _vars.assign(*dhist->get(i));
-         add(_vars, dhist->weight() * initWgt, pow(dhist->weightError(SumW2), 2));
+         add(_vars, dhist->weight(i) * initWgt, pow(dhist->weightError(SumW2), 2));
       }
    }
 }
@@ -955,7 +955,7 @@ std::unique_ptr<RooAbsData> RooDataHist::reduceEng(const RooArgSet& varSubset, c
 
     if (!cloneVar || cloneVar->getVal()) {
       weightError(lo,hi,SumW2) ;
-      rdh->add(*row,weight(),lo*lo) ;
+      rdh->add(*row,weight(i),lo*lo) ;
     }
   }
 
@@ -1572,7 +1572,7 @@ void RooDataHist::weightError(double& lo, double& hi, ErrorType etype) const
     throw std::invalid_argument("RooDataHist::weightError(" + std::string(GetName()) + ") error type Expected not allowed here");
     break ;
 
-  case Poisson:
+  case Poisson: {
     if (_errLo && _errLo[_curIndex] >= 0.0) {
       // Weight is preset or precalculated
       lo = _errLo[_curIndex];
@@ -1586,12 +1586,14 @@ void RooDataHist::weightError(double& lo, double& hi, ErrorType etype) const
     // Calculate poisson errors
     double ym;
     double yp;
-    RooHistError::instance().getPoissonInterval(Int_t(weight()+0.5),ym,yp,1) ;
-    _errLo[_curIndex] = weight()-ym;
-    _errHi[_curIndex] = yp-weight();
+    const double w = weight(_curIndex);
+    RooHistError::instance().getPoissonInterval(Int_t(w+0.5),ym,yp,1) ;
+    _errLo[_curIndex] = w-ym;
+    _errHi[_curIndex] = yp-w;
     lo = _errLo[_curIndex];
     hi = _errHi[_curIndex];
     return ;
+  }
 
   case SumW2:
     lo = std::sqrt(weightSquared(_curIndex));
@@ -2348,7 +2350,7 @@ void RooDataHist::printContents(std::ostream& os) const
 
         double lo, hi;
         weightError(lo, hi, RooAbsData::SumW2);
-        os << ", weight=" << weight() << " +/- [" << lo << "," << hi << "]"
+        os << ", weight=" << weight(i) << " +/- [" << lo << "," << hi << "]"
            << std::endl;
     }
 }
diff --git a/roofit/roofitcore/src/RooEvaluatorWrapper.cxx b/roofit/roofitcore/src/RooEvaluatorWrapper.cxx
index 9f4d3be01b532..c0a2a66c0b1e5 100644
--- a/roofit/roofitcore/src/RooEvaluatorWrapper.cxx
+++ b/roofit/roofitcore/src/RooEvaluatorWrapper.cxx
@@ -741,6 +741,11 @@ void RooEvaluatorWrapper::writeDebugMacro(std::string const &filename) const
       return _funcWrapper->writeDebugMacro(filename);
 }
 
+std::unique_ptr<ChangeOperModeRAII> RooEvaluatorWrapper::setOperModes(RooAbsArg::OperMode opMode)
+{
+   return _evaluator->setOperModes(opMode);
+}
+
 } // namespace RooFit::Experimental
 
 /// \endcond
diff --git a/roofit/roofitcore/src/RooFFTConvPdf.cxx b/roofit/roofitcore/src/RooFFTConvPdf.cxx
index 15bfde3f3883a..790ea9aef1f5b 100644
--- a/roofit/roofitcore/src/RooFFTConvPdf.cxx
+++ b/roofit/roofitcore/src/RooFFTConvPdf.cxx
@@ -657,10 +657,11 @@ void RooFFTConvPdf::fillCacheSlice(FFTCacheElem& aux, const RooArgSet& slicePos)
     while (j>=N2) j-= N2 ;
 
     iter->Next() ;
+    const std::size_t binIdx = cacheHist.getIndex(*cacheHist.get(), /*fast=*/true);
 #ifndef ROOFIT_MATH_FFTW3
-    cacheHist.set(output[j]);
+    cacheHist.set(binIdx, output[j], -1.);
 #else
-    cacheHist.set(aux.fftc2r->GetPointReal(j));
+    cacheHist.set(binIdx, aux.fftc2r->GetPointReal(j), -1.);
 #endif
   }
 }
diff --git a/roofit/roofitcore/src/RooFit/Evaluator.cxx b/roofit/roofitcore/src/RooFit/Evaluator.cxx
index 0ad61d2774a26..b17fc94cbe598 100644
--- a/roofit/roofitcore/src/RooFit/Evaluator.cxx
+++ b/roofit/roofitcore/src/RooFit/Evaluator.cxx
@@ -46,6 +46,7 @@ RooAbsPdf::fitTo() is called and gets destroyed when the fitting ends.
 #include <iomanip>
 #include <numeric>
 #include <thread>
+#include <unordered_set>
 
 namespace RooFit {
 
@@ -325,16 +326,16 @@ void Evaluator::updateOutputSizes()
    for (auto &info : _nodes) {
       info.outputSize = outputSizeMap.at(info.absArg);
       info.isDirty = true;
-
-      // In principle we don't need dirty flag propagation because the driver
-      // takes care of deciding which node needs to be re-evaluated. However,
-      // disabling it also for scalar mode results in very long fitting times
-      // for specific models (test 14 in stressRooFit), which still needs to be
-      // understood. TODO.
-      if (!info.isScalar()) {
+      // We don't need dirty flag propagation because the evaluator takes care
+      // of deciding what needs to be re-evaluated. We can disable the regular
+      // dirty state propagation. However, fundamental variables like
+      // RooRealVars and RooCategories are usually shared with other
+      // computation graphs outside the evaluator, so they can't be mutated.
+      // See also the code of the RooMinimizer, which ensures that dirty state
+      // propagation is temporarily disabled during minimization to really
+      // eliminate any overhead from the dirty flag propagation.
+      if (!info.absArg->isFundamental()) {
          setOperMode(info.absArg, RooAbsArg::ADirty);
-      } else {
-         setOperMode(info.absArg, info.originalOperMode);
       }
    }
 
@@ -627,9 +628,52 @@ void Evaluator::markGPUNodes()
 /// Evaluator gets deleted.
 void Evaluator::setOperMode(RooAbsArg *arg, RooAbsArg::OperMode opMode)
 {
-   if (opMode != arg->operMode()) {
-      _changeOperModeRAIIs.emplace(std::make_unique<ChangeOperModeRAII>(arg, opMode));
+   if (!_operModeChanges)
+      _operModeChanges = std::make_unique<ChangeOperModeRAII>();
+   _operModeChanges->change(arg, opMode);
+}
+
+// Change the operation modes of all RooAbsArgs in the computation graph.
+// The changes are reset when the returned RAII object goes out of scope.
+//
+// We also walk transitively through value clients of the nodes to cover any
+// node that RooAbsReal::doEval (the fallback scalar implementation) might
+// inadvertently propagate the ADirty mode to via its recursive restore: that
+// helper sets servers temporarily to AClean and then calls
+// setOperMode(oldOperMode) to restore, which recurses to value clients when
+// oldOperMode is ADirty. If we did not protect those clients here, any node
+// outside the computation graph that shares a fundamental (e.g. a parameter
+// like a RooRealVar) would be left permanently in ADirty after the first
+// minimization, dramatically slowing down later scalar evaluations (for
+// example on pdfs held by the legacy test statistics' internal cache).
+std::unique_ptr<ChangeOperModeRAII> Evaluator::setOperModes(RooAbsArg::OperMode opMode)
+{
+   auto out = std::make_unique<ChangeOperModeRAII>();
+   std::unordered_set<RooAbsArg *> visited;
+
+   std::vector<RooAbsArg *> queue;
+   queue.reserve(_nodes.size());
+   for (auto &info : _nodes) {
+      queue.push_back(info.absArg);
+   }
+
+   while (!queue.empty()) {
+      RooAbsArg *node = queue.back();
+      queue.pop_back();
+      if (!visited.insert(node).second)
+         continue;
+
+      out->change(node, opMode);
+
+      // Only follow value-client links: that is exactly the propagation path
+      // used by RooAbsArg::setOperMode with mode==ADirty.
+      if (opMode == RooAbsArg::ADirty) {
+         for (auto *client : node->valueClients()) {
+            queue.push_back(client);
+         }
+      }
    }
+   return out;
 }
 
 void Evaluator::print(std::ostream &os)
diff --git a/roofit/roofitcore/src/RooHistFunc.cxx b/roofit/roofitcore/src/RooHistFunc.cxx
index 54f9bf9e472cd..043df4cd492e4 100644
--- a/roofit/roofitcore/src/RooHistFunc.cxx
+++ b/roofit/roofitcore/src/RooHistFunc.cxx
@@ -247,8 +247,7 @@ double RooHistFunc::maxVal(Int_t code) const
 
   double max(-1) ;
   for (Int_t i=0 ; i<_dataHist->numEntries() ; i++) {
-    _dataHist->get(i) ;
-    double wgt = _dataHist->weight() ;
+    double wgt = _dataHist->weight(i) ;
     if (wgt>max) max=wgt ;
   }
 
@@ -487,9 +486,7 @@ bool RooHistFunc::areIdentical(const RooDataHist& dh1, const RooDataHist& dh2)
   if (std::abs(dh1.sumEntries()-dh2.sumEntries())>1e-8) return false ;
   if (dh1.numEntries() != dh2.numEntries()) return false ;
   for (int i=0 ; i < dh1.numEntries() ; i++) {
-    dh1.get(i) ;
-    dh2.get(i) ;
-    if (std::abs(dh1.weight()-dh2.weight())>1e-8) return false ;
+    if (std::abs(dh1.weight(i)-dh2.weight(i))>1e-8) return false ;
   }
   using RooHelpers::getColonSeparatedNameString;
   if (getColonSeparatedNameString(*dh1.get()) != getColonSeparatedNameString(*dh2.get())) return false ;
diff --git a/roofit/roofitcore/src/RooHistPdf.cxx b/roofit/roofitcore/src/RooHistPdf.cxx
index dd3555827e854..4879f782596b3 100644
--- a/roofit/roofitcore/src/RooHistPdf.cxx
+++ b/roofit/roofitcore/src/RooHistPdf.cxx
@@ -530,8 +530,7 @@ double RooHistPdf::maxVal(Int_t code) const
 
   double max(-1) ;
   for (Int_t i=0 ; i<_dataHist->numEntries() ; i++) {
-    _dataHist->get(i) ;
-    double wgt = _dataHist->weight() ;
+    double wgt = _dataHist->weight(i) ;
     if (wgt>max) max=wgt ;
   }
 
@@ -548,9 +547,7 @@ bool RooHistPdf::areIdentical(const RooDataHist& dh1, const RooDataHist& dh2)
   if (std::abs(dh1.sumEntries()-dh2.sumEntries())>1e-8) return false ;
   if (dh1.numEntries() != dh2.numEntries()) return false ;
   for (int i=0 ; i < dh1.numEntries() ; i++) {
-    dh1.get(i) ;
-    dh2.get(i) ;
-    if (std::abs(dh1.weight()-dh2.weight())>1e-8) return false ;
+    if (std::abs(dh1.weight(i)-dh2.weight(i))>1e-8) return false ;
   }
   return true ;
 }
diff --git a/roofit/roofitcore/src/RooMinimizer.cxx b/roofit/roofitcore/src/RooMinimizer.cxx
index 718a13880ac2b..af81c50102902 100644
--- a/roofit/roofitcore/src/RooMinimizer.cxx
+++ b/roofit/roofitcore/src/RooMinimizer.cxx
@@ -40,27 +40,30 @@ automatic PDF optimization.
 #include "RooMinimizer.h"
 
 #include "RooAbsMinimizerFcn.h"
-#include "RooArgSet.h"
-#include "RooArgList.h"
 #include "RooAbsReal.h"
+#include "RooArgList.h"
+#include "RooArgSet.h"
+#include "RooCategory.h"
 #include "RooDataSet.h"
-#include "RooRealVar.h"
-#include "RooSentinel.h"
+#include "RooEvaluatorWrapper.h"
+#include "RooFit/TestStatistics/RooAbsL.h"
+#include "RooFit/TestStatistics/RooRealL.h"
+#include "RooFitResult.h"
+#include "RooHelpers.h"
+#include "RooMinimizerFcn.h"
 #include "RooMsgService.h"
-#include "RooCategory.h"
 #include "RooMultiPdf.h"
 #include "RooPlot.h"
-#include "RooHelpers.h"
-#include "RooMinimizerFcn.h"
-#include "RooFitResult.h"
-#include "RooFit/TestStatistics/RooAbsL.h"
-#include "RooFit/TestStatistics/RooRealL.h"
+#include "RooRealVar.h"
+#include "RooSentinel.h"
 #ifdef ROOFIT_MULTIPROCESS
 #include "TestStatistics/MinuitFcnGrad.h"
 #include "RooFit/MultiProcess/Config.h"
 #include "RooFit/MultiProcess/ProcessTimer.h"
 #endif
 
+#include "RooFitImplHelpers.h"
+
 #include <Fit/BasicFCN.h>
 #include <Math/Minimizer.h>
 #include <TClass.h>
@@ -120,6 +123,22 @@ void reorderCombinations(std::vector<std::vector<int>> &combos, const std::vecto
    }
 }
 
+// The RooEvaluatorWrapper uses its own logic to decide what needs to be
+// re-evaluated. We can therefore disable the regular dirty state propagation
+// temporarily during minimization. However, some RooAbsArgs shared with other
+// regular RooFit computation graphs outside the minimized likelihood, so we
+// have to make sure that the operation mode is reset after the minimization.
+//
+// This should be called before running any routine via the _minimizer data
+// member. The RAII object should only be destructed after the routine is done.
+std::unique_ptr<ChangeOperModeRAII> setOperModesDirty(RooAbsReal &function)
+{
+   if (auto *wrapper = dynamic_cast<RooFit::Experimental::RooEvaluatorWrapper *>(&function)) {
+      return wrapper->setOperModes(RooAbsArg::ADirty);
+   }
+   return {};
+}
+
 } // namespace
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -135,7 +154,7 @@ void reorderCombinations(std::vector<std::vector<int>> &combos, const std::vecto
 /// value of the input function.
 
 /// Constructor that accepts all configuration in struct with RooAbsReal likelihood
-RooMinimizer::RooMinimizer(RooAbsReal &function, Config const &cfg) : _cfg(cfg)
+RooMinimizer::RooMinimizer(RooAbsReal &function, Config const &cfg) : _function{function}, _cfg(cfg)
 {
    initMinimizerFirstPart();
    auto nll_real = dynamic_cast<RooFit::TestStatistics::RooRealL *>(&function);
@@ -692,6 +711,7 @@ RooPlot *RooMinimizer::contour(RooRealVar &var1, RooRealVar &var2, double n1, do
    n[4] = n5;
    n[5] = n6;
 
+   auto operModeRAII = setOperModesDirty(_function);
    for (int ic = 0; ic < 6; ic++) {
       if (n[ic] > 0) {
 
@@ -906,6 +926,8 @@ bool RooMinimizer::fitFCN()
    // fit a user provided FCN function
    // create fit parameter settings
 
+   auto operModeRAII = setOperModesDirty(_function);
+
    // Check number of parameters
    unsigned int npar = getNPar();
    if (npar == 0) {
@@ -1045,6 +1067,8 @@ bool RooMinimizer::calculateHessErrors()
    // compute the Hesse errors according to configuration
    // set in the parameters and append value in fit result
 
+   auto operModeRAII = setOperModesDirty(_function);
+
    // update  minimizer (recreate if not done or if name has changed
    if (!updateMinimizerOptions()) {
       coutE(Minimization) << "RooMinimizer::calculateHessErrors() Error re-initializing the minimizer" << std::endl;
@@ -1079,6 +1103,8 @@ bool RooMinimizer::calculateMinosErrors()
    // (in DoMinimization) aftewr minimizing if the
    //  FitConfig::MinosErrors() flag is set
 
+   auto operModeRAII = setOperModesDirty(_function);
+
    // update  minimizer (but cannot re-create in this case). Must use an existing one
    if (!updateMinimizerOptions(false)) {
       coutE(Minimization) << "RooMinimizer::calculateHessErrors() Error re-initializing the minimizer" << std::endl;
diff --git a/roofit/roofitcore/src/RooRealIntegral.cxx b/roofit/roofitcore/src/RooRealIntegral.cxx
index c580bb0b2c10e..a1aa4bc459e54 100644
--- a/roofit/roofitcore/src/RooRealIntegral.cxx
+++ b/roofit/roofitcore/src/RooRealIntegral.cxx
@@ -841,12 +841,12 @@ double RooRealIntegral::evaluate() const
         // from caching subgraph results (e.g. for nested numeric integrals).
         RooArgList serverList;
         _function->treeNodeServerList(&serverList, nullptr, true, true, false, true);
-        std::list<ChangeOperModeRAII> operModeRAII;
+        ChangeOperModeRAII operModeRAII;
 
         for (auto *arg : serverList) {
            arg->syncCache();
            if (arg->operMode() == RooAbsArg::AClean) {
-              operModeRAII.emplace_back(arg, RooAbsArg::Auto);
+              operModeRAII.change(arg, RooAbsArg::Auto);
            }
         }
 
diff --git a/roofit/roofitcore/test/testRooDataHist.cxx b/roofit/roofitcore/test/testRooDataHist.cxx
index e8d6084e75541..d2b78a81f94f3 100644
--- a/roofit/roofitcore/test/testRooDataHist.cxx
+++ b/roofit/roofitcore/test/testRooDataHist.cxx
@@ -96,7 +96,7 @@ TEST(RooDataHist, UnWeightedEntries)
    EXPECT_EQ(static_cast<RooRealVar *>(coordsAtZero.find(x))->getVal(),
              static_cast<RooRealVar *>(coordsAtPoint9.find(x))->getVal());
 
-   const double weight = dataHist.weight();
+   const double weight = dataHist.weight(dataHist.getIndex(coordinates));
    EXPECT_EQ(weight, targetBinContent);
 
    EXPECT_NEAR(dataHist.weightError(RooAbsData::Poisson), sqrt(targetBinContent),
@@ -124,7 +124,7 @@ TEST(RooDataHist, UnWeightedEntries)
 
    RooArgSet coordsAt10;
    dataHist.get(10)->snapshot(coordsAt10);
-   const double weightBin10 = dataHist.weight();
+   const double weightBin10 = dataHist.weight(10);
 
    EXPECT_NEAR(static_cast<RooRealVar *>(coordsAt10.find(x))->getVal(), 0.5, 1.E-1);
    EXPECT_EQ(weight, weightBin10);
@@ -149,7 +149,7 @@ TEST(RooDataHist, WeightedEntries)
 
    x.setVal(0.);
    dataHist.get(coordinates);
-   const double weight = dataHist.weight();
+   const double weight = dataHist.weight(dataHist.getIndex(coordinates));
    ASSERT_EQ(weight, targetBinContent);
 
    const double targetError = sqrt(10 * 4.);
@@ -178,7 +178,7 @@ TEST(RooDataHist, WeightedEntries)
 
    RooArgSet coordsAt10;
    dataHist.get(10)->snapshot(coordsAt10);
-   const double weightBin10 = dataHist.weight();
+   const double weightBin10 = dataHist.weight(10);
 
    EXPECT_NEAR(static_cast<RooRealVar *>(coordsAt10.find(x))->getVal(), 0.5, 1.E-1);
    EXPECT_EQ(weight, weightBin10);
@@ -217,8 +217,7 @@ TEST_P(RooDataHistIO, ReadLegacy)
    EXPECT_EQ(dataHist.sumEntries(), 20 * targetBinContent);
 
    static_cast<RooRealVar *>(legacyVals.find("x"))->setVal(0.);
-   dataHist.get(legacyVals); // trigger side effect for weight below.
-   const double weight = dataHist.weight();
+   const double weight = dataHist.weight(dataHist.getIndex(legacyVals));
    ASSERT_EQ(weight, targetBinContent);
 
    const double targetError = sqrt(10 * 4.);
@@ -247,7 +246,7 @@ TEST_P(RooDataHistIO, ReadLegacy)
 
    RooArgSet coordsAt10;
    dataHist.get(10)->snapshot(coordsAt10);
-   const double weightBin10 = dataHist.weight();
+   const double weightBin10 = dataHist.weight(10);
 
    EXPECT_NEAR(static_cast<RooRealVar *>(coordsAt10.find("x"))->getVal(), 0.5, 1.E-1);
    EXPECT_EQ(weight, weightBin10);
@@ -804,11 +803,11 @@ TEST(RooDataHist, SplitDataHistWithSumW2)
 
    data1.get(0);
    data2.get(0);
-   EXPECT_DOUBLE_EQ(data2.weightSquared(), data1.weightSquared());
+   EXPECT_DOUBLE_EQ(data2.weightSquared(), data1.weightSquared(0));
    EXPECT_DOUBLE_EQ(data2.weightError(), data1.weightError());
 
    data1.get(1);
    data2.get(1);
-   EXPECT_DOUBLE_EQ(data2.weightSquared(), data1.weightSquared());
+   EXPECT_DOUBLE_EQ(data2.weightSquared(), data1.weightSquared(1));
    EXPECT_DOUBLE_EQ(data2.weightError(), data1.weightError());
 }
diff --git a/roofit/roofitcore/test/testRooDataSet.cxx b/roofit/roofitcore/test/testRooDataSet.cxx
index c3d72473f4deb..7e1e45443feaa 100644
--- a/roofit/roofitcore/test/testRooDataSet.cxx
+++ b/roofit/roofitcore/test/testRooDataSet.cxx
@@ -374,10 +374,9 @@ TEST(RooDataSet, ImportDataHist)
    RooDataSet ds{"ds", "ds", x, RooFit::Import(dh)};
 
    for (int i = 0; i < x.numBins(); ++i) {
-      dh.get(i);
       ds.get(i);
-      EXPECT_FLOAT_EQ(ds.weight(), dh.weight()) << "weight() is off in bin " << i;
-      EXPECT_FLOAT_EQ(ds.weightSquared(), dh.weightSquared()) << "weightSquared() is off in bin " << i;
+      EXPECT_FLOAT_EQ(ds.weight(), dh.weight(i)) << "weight() is off in bin " << i;
+      EXPECT_FLOAT_EQ(ds.weightSquared(), dh.weightSquared(i)) << "weightSquared() is off in bin " << i;
    }
 }
 
diff --git a/roofit/roostats/inc/RooStats/BayesianCalculator.h b/roofit/roostats/inc/RooStats/BayesianCalculator.h
index f0800cb866abf..5529a456d362a 100644
--- a/roofit/roostats/inc/RooStats/BayesianCalculator.h
+++ b/roofit/roostats/inc/RooStats/BayesianCalculator.h
@@ -151,8 +151,6 @@ namespace RooStats {
 
       void ComputeIntervalFromCdf(double c1, double c2) const;
 
-      void ComputeIntervalUsingRooFit(double c1, double c2) const;
-
       void ComputeShortestInterval() const;
 
    private:
diff --git a/roofit/roostats/src/BayesianCalculator.cxx b/roofit/roostats/src/BayesianCalculator.cxx
index d836e39948c75..3d5d2e513231d 100644
--- a/roofit/roostats/src/BayesianCalculator.cxx
+++ b/roofit/roostats/src/BayesianCalculator.cxx
@@ -54,7 +54,6 @@ credible interval from the given function.
 #include "RooAbsReal.h"
 #include "RooRealVar.h"
 #include "RooArgSet.h"
-#include "RooBrentRootFinder.h"
 #include "RooFormulaVar.h"
 #include "RooGenericPdf.h"
 #include "RooPlot.h"
@@ -167,7 +166,6 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction {
       fFunctor(nll, bindParams, RooArgList() ),              // functor
       fPriorFunc(nullptr),
       fLikelihood(fFunctor, nullptr, nllMinimum),         // integral of exp(-nll) function
-      fIntegrator(ROOT::Math::IntegratorMultiDim::GetType(integType) ),  // integrator
       fXmin(bindParams.size()),               // vector of parameters (min values)
       fXmax(bindParams.size())
    {
@@ -176,7 +174,15 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction {
          fLikelihood.SetPrior(fPriorFunc.get() );
       }
 
-      fIntegrator.SetFunction(fLikelihood, bindParams.size() );
+      // ROOT::Math::IntegratorMultiDim does not support dim = 1, so fall back to
+      // the 1D integrator when there are no nuisance parameters.
+      if (bindParams.size() == 1) {
+         fIntegratorOneDim = std::make_unique<ROOT::Math::Integrator>(ROOT::Math::IntegratorOneDim::GetType(integType));
+         fIntegratorOneDim->SetFunction(fLikelihood);
+      } else {
+         fIntegrator = std::make_unique<ROOT::Math::IntegratorMultiDim>(ROOT::Math::IntegratorMultiDim::GetType(integType));
+         fIntegrator->SetFunction(fLikelihood, bindParams.size());
+      }
 
       ooccoutD(nullptr,NumericIntegration) << "PosteriorCdfFunction::Compute integral of posterior in nuisance and poi. "
                                            << " nllMinimum is " << nllMinimum << std::endl;
@@ -191,7 +197,10 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction {
                                               << " in interval [ " <<  fXmin[i] << " , " << fXmax[i] << " ] " << std::endl;
       }
 
-      fIntegrator.Options().Print(ooccoutD(nullptr,NumericIntegration));
+      if (fIntegrator)
+         fIntegrator->Options().Print(ooccoutD(nullptr,NumericIntegration));
+      else
+         fIntegratorOneDim->Options().Print(ooccoutD(nullptr,NumericIntegration));
 
       // store max POI value because it will be changed when evaluating the function
       fMaxPOI = fXmax[0];
@@ -214,7 +223,6 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction {
       //fPriorFunc(std::shared_ptr<RooFunctor>((RooFunctor*)0)),
       fPriorFunc(rhs.fPriorFunc),
       fLikelihood(fFunctor, fPriorFunc.get(), rhs.fLikelihood.fOffset),
-      fIntegrator(ROOT::Math::IntegratorMultiDim::GetType( rhs.fIntegrator.Name().c_str() ) ),  // integrator
       fXmin( rhs.fXmin),
       fXmax( rhs.fXmax),
       fNorm( rhs.fNorm),
@@ -226,7 +234,13 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction {
       fError(rhs.fError),
       fNormCdfValues(rhs.fNormCdfValues)
    {
-      fIntegrator.SetFunction(fLikelihood, fXmin.size() );
+      if (rhs.fIntegratorOneDim) {
+         fIntegratorOneDim = std::make_unique<ROOT::Math::Integrator>(ROOT::Math::IntegratorOneDim::GetType(rhs.fIntegratorOneDim->Name().c_str()));
+         fIntegratorOneDim->SetFunction(fLikelihood);
+      } else {
+         fIntegrator = std::make_unique<ROOT::Math::IntegratorMultiDim>(ROOT::Math::IntegratorMultiDim::GetType(rhs.fIntegrator->Name().c_str()));
+         fIntegrator->SetFunction(fLikelihood, fXmin.size());
+      }
       // need special treatment for the smart pointer
       // if (rhs.fPriorFunc.get() ) {
       //    fPriorFunc = std::shared_ptr<RooFunctor>(new RooFunctor(*(rhs.fPriorFunc) ) );
@@ -277,8 +291,15 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction {
 
       fFunctor.binding().resetNumCall();  // reset number of calls for debug
 
-      double cdf = fIntegrator.Integral(&fXmin[0],&fXmax[0]);
-      double error = fIntegrator.Error();
+      double cdf;
+      double error;
+      if (fIntegratorOneDim) {
+         cdf = fIntegratorOneDim->Integral(fXmin[0], fXmax[0]);
+         error = fIntegratorOneDim->Error();
+      } else {
+         cdf = fIntegrator->Integral(&fXmin[0], &fXmax[0]);
+         error = fIntegrator->Error();
+      }
       double normcdf =  cdf/fNorm;  // normalize the cdf
 
       ooccoutD(nullptr,NumericIntegration) << "PosteriorCdfFunction: poi = [" << fXmin[0] << " , "
@@ -324,7 +345,8 @@ class PosteriorCdfFunction : public ROOT::Math::IGenFunction {
    mutable RooFunctor fFunctor;                   // functor binding nll
    mutable std::shared_ptr<RooFunctor> fPriorFunc;  // functor binding the prior
    LikelihoodFunction fLikelihood;               // likelihood function
-   mutable ROOT::Math::IntegratorMultiDim  fIntegrator; // integrator  (mutable because Integral() is not const
+   mutable std::unique_ptr<ROOT::Math::IntegratorMultiDim> fIntegrator; // multi-dim integrator (for >=2 dims)
+   mutable std::unique_ptr<ROOT::Math::Integrator> fIntegratorOneDim;   // 1D integrator (for POI-only case)
    mutable std::vector<double> fXmin;    // min value of parameters (poi+nuis) -
    mutable std::vector<double> fXmax;   // max value of parameters (poi+nuis) - max poi changes so it is mutable
    double fNorm = 1.0;                  // normalization value (computed in constructor)
@@ -1127,14 +1149,7 @@ SimpleInterval* BayesianCalculator::GetInterval() const
       }
 
       else {
-         // use integration method if there are nuisance parameters
-         if (!fNuisanceParameters.empty()) {
-            ComputeIntervalFromCdf(lowerCutOff, upperCutOff);
-         }
-         else {
-            // case of no nuisance - just use createCdf from roofit
-            ComputeIntervalUsingRooFit(lowerCutOff, upperCutOff);
-         }
+         ComputeIntervalFromCdf(lowerCutOff, upperCutOff);
          // case cdf failed (scan then the posterior)
          if (!fValidInterval) {
             fNScanBins = 100;
@@ -1187,55 +1202,6 @@ double BayesianCalculator::GetMode() const {
    //return  fApproxPosterior->GetMaximumX();
 }
 
-////////////////////////////////////////////////////////////////////////////////
-/// internal function compute the interval using RooFit
-
-void BayesianCalculator::ComputeIntervalUsingRooFit(double lowerCutOff, double upperCutOff ) const {
-
-   coutI(Eval) <<  "BayesianCalculator: Compute interval using RooFit:  posteriorPdf + createCdf + RooBrentRootFinder " << std::endl;
-
-   RooRealVar* poi = dynamic_cast<RooRealVar*>(fPOI.first() );
-   assert(poi);
-
-   fValidInterval = false;
-   if (!fPosteriorPdf) fPosteriorPdf = (RooAbsPdf*) GetPosteriorPdf();
-   if (!fPosteriorPdf) return;
-
-   std::unique_ptr<RooAbsReal> cdf{fPosteriorPdf->createCdf(fPOI,RooFit::ScanNoCdf())};
-   if (!cdf) return;
-
-   std::unique_ptr<RooAbsFunc> cdf_bind{cdf->bindVars(fPOI,&fPOI)};
-   if (!cdf_bind) return;
-
-   RooBrentRootFinder brf(*cdf_bind);
-   brf.setTol(fBrfPrecision); // set the brf precision
-
-   double tmpVal = poi->getVal();  // patch used because findRoot changes the value of poi
-
-   bool ret = true;
-   if (lowerCutOff > 0) {
-      double y = lowerCutOff;
-      ret &= brf.findRoot(fLower,poi->getMin(),poi->getMax(),y);
-   }
-   else
-      fLower = poi->getMin();
-
-   if (upperCutOff < 1.0) {
-      double y=upperCutOff;
-      ret &= brf.findRoot(fUpper,poi->getMin(),poi->getMax(),y);
-   }
-   else
-      fUpper = poi->getMax();
-   if (!ret) {
-      coutE(Eval) << "BayesianCalculator::GetInterval "
-                  << "Error returned from Root finder, estimated interval is not fully correct" << std::endl;
-   } else {
-      fValidInterval = true;
-   }
-
-   poi->setVal(tmpVal); // patch: restore the original value of poi
-}
-
 ////////////////////////////////////////////////////////////////////////////////
 /// internal function compute the interval using Cdf integration
 
diff --git a/roofit/roostats/src/MCMCInterval.cxx b/roofit/roostats/src/MCMCInterval.cxx
index acfe5261503a9..3f246ccd79912 100644
--- a/roofit/roostats/src/MCMCInterval.cxx
+++ b/roofit/roostats/src/MCMCInterval.cxx
@@ -111,10 +111,8 @@ MCMCInterval::~MCMCInterval() = default;
 struct CompareDataHistBins {
    CompareDataHistBins(RooDataHist* hist) : fDataHist(hist) {}
    bool operator() (Int_t bin1 , Int_t bin2) {
-      fDataHist->get(bin1);
-      double n1 = fDataHist->weight();
-      fDataHist->get(bin2);
-      double n2 = fDataHist->weight();
+      double n1 = fDataHist->weight(bin1);
+      double n2 = fDataHist->weight(bin2);
       return (n1 < n2);
    }
    RooDataHist* fDataHist;
@@ -185,8 +183,7 @@ bool MCMCInterval::IsInInterval(const RooArgSet& point) const
             // >= cutoff
             Int_t bin;
             bin = fDataHist->getIndex(point);
-            fDataHist->get(bin);
-            return (fDataHist->weight() >= (double)fHistCutoff);
+            return (fDataHist->weight(bin) >= (double)fHistCutoff);
          }
       }
    } else if (fIntervalType == kTailFraction) {
@@ -839,8 +836,7 @@ void MCMCInterval::DetermineByDataHist()
    double content;
    Int_t i;
    for (i = numBins - 1; i >= 0; i--) {
-      fDataHist->get(bins[i]);
-      content = fDataHist->weight();
+      content = fDataHist->weight(bins[i]);
       if ((sum + content) / nEntries >= fConfidenceLevel) {
          fHistCutoff = content;
          if (fIsHistStrict) {
@@ -858,8 +854,7 @@ void MCMCInterval::DetermineByDataHist()
    if (fIsHistStrict) {
       // keep going to find the sum
       for ( ; i >= 0; i--) {
-         fDataHist->get(bins[i]);
-         content = fDataHist->weight();
+         content = fDataHist->weight(bins[i]);
          if (content == fHistCutoff) {
             sum += content;
          } else {
@@ -869,8 +864,7 @@ void MCMCInterval::DetermineByDataHist()
    } else {
       // backtrack to find the cutoff and sum
       for ( ; i < numBins; i++) {
-         fDataHist->get(bins[i]);
-         content = fDataHist->weight();
+         content = fDataHist->weight(bins[i]);
          if (content > fHistCutoff) {
             fHistCutoff = content;
             break;
@@ -1070,9 +1064,8 @@ double MCMCInterval::LowerLimitByDataHist(RooRealVar& param)
          double lowerLimit = param.getMax();
          double val;
          for (Int_t i = 0; i < numBins; i++) {
-            fDataHist->get(i);
-            if (fDataHist->weight() >= fHistCutoff) {
-               val = fDataHist->get()->getRealValue(param.GetName());
+            if (fDataHist->weight(i) >= fHistCutoff) {
+               val = fDataHist->get(i)->getRealValue(param.GetName());
                if (val < lowerLimit)
                   lowerLimit = val;
             }
@@ -1147,9 +1140,8 @@ double MCMCInterval::UpperLimitByDataHist(RooRealVar& param)
          double upperLimit = param.getMin();
          double val;
          for (Int_t i = 0; i < numBins; i++) {
-            fDataHist->get(i);
-            if (fDataHist->weight() >= fHistCutoff) {
-               val = fDataHist->get()->getRealValue(param.GetName());
+            if (fDataHist->weight(i) >= fHistCutoff) {
+               val = fDataHist->get(i)->getRealValue(param.GetName());
                if (val > upperLimit)
                   upperLimit = val;
             }
@@ -1187,9 +1179,8 @@ double MCMCInterval::LowerLimitByKeys(RooRealVar& param)
          double lowerLimit = param.getMax();
          double val;
          for (Int_t i = 0; i < numBins; i++) {
-            fKeysDataHist->get(i);
-            if (fKeysDataHist->weight() >= fKeysCutoff) {
-               val = fKeysDataHist->get()->getRealValue(param.GetName());
+            if (fKeysDataHist->weight(i) >= fKeysCutoff) {
+               val = fKeysDataHist->get(i)->getRealValue(param.GetName());
                if (val < lowerLimit)
                   lowerLimit = val;
             }
@@ -1227,9 +1218,8 @@ double MCMCInterval::UpperLimitByKeys(RooRealVar& param)
          double upperLimit = param.getMin();
          double val;
          for (Int_t i = 0; i < numBins; i++) {
-            fKeysDataHist->get(i);
-            if (fKeysDataHist->weight() >= fKeysCutoff) {
-               val = fKeysDataHist->get()->getRealValue(param.GetName());
+            if (fKeysDataHist->weight(i) >= fKeysCutoff) {
+               val = fKeysDataHist->get(i)->getRealValue(param.GetName());
                if (val > upperLimit)
                   upperLimit = val;
             }
@@ -1264,8 +1254,7 @@ double MCMCInterval::GetKeysMax()
    double max = 0;
    double w;
    for (Int_t i = 0; i < numBins; i++) {
-      fKeysDataHist->get(i);
-      w = fKeysDataHist->weight();
+      w = fKeysDataHist->weight(i);
       if (w > max)
          max = w;
    }
diff --git a/roofit/roostats/test/CMakeLists.txt b/roofit/roostats/test/CMakeLists.txt
index 75eff086a88d6..7770dc075e93e 100644
--- a/roofit/roostats/test/CMakeLists.txt
+++ b/roofit/roostats/test/CMakeLists.txt
@@ -1,4 +1,5 @@
 ROOT_ADD_GTEST(testAsymptoticCalculator testAsymptoticCalculator.cxx LIBRARIES RooStats)
+ROOT_ADD_GTEST(testBayesianCalculator testBayesianCalculator.cxx LIBRARIES RooStats)
 ROOT_ADD_GTEST(testHypoTestInvResult testHypoTestInvResult.cxx
   LIBRARIES RooStats
   COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/testHypoTestInvResult_1.root)
diff --git a/roofit/roostats/test/testBayesianCalculator.cxx b/roofit/roostats/test/testBayesianCalculator.cxx
new file mode 100644
index 0000000000000..164ed8a1357ae
--- /dev/null
+++ b/roofit/roostats/test/testBayesianCalculator.cxx
@@ -0,0 +1,65 @@
+#include <RooAddPdf.h>
+#include <RooArgSet.h>
+#include <RooDataSet.h>
+#include <RooExponential.h>
+#include <RooGaussian.h>
+#include <RooRandom.h>
+#include <RooRealVar.h>
+#include <RooStats/BayesianCalculator.h>
+#include <RooStats/SimpleInterval.h>
+#include <RooUniform.h>
+
+#include <gtest/gtest.h>
+
+/// Regression test for https://github.com/root-project/root/issues/17567
+
+TEST(BayesianCalculator, UpperLimitNoNuisance_Issue17567)
+{
+   RooRandom::randomGenerator()->SetSeed(12345);
+
+   // Background-only generation model.
+   RooRealVar mass("mass", "mass", 6500, 7500);
+   RooRealVar slopeGen("slopeGen", "slopeGen", -8e-4);
+   RooExponential expoGen("expoGen", "expoGen", mass, slopeGen);
+   std::unique_ptr<RooDataSet> data{expoGen.generate(mass, 1000)};
+
+   // Fit model: Gaussian signal on top of exponential background.
+   // All parameters except the POI ngauss are fixed so there are no nuisances.
+   RooRealVar mean("mean", "mean", 7024);
+   RooRealVar sigma("sigma", "sigma", 44);
+   RooGaussian gauss("gauss", "gauss", mass, mean, sigma);
+   RooRealVar slope("slope", "slope", -8e-4);
+   RooExponential expo("expo", "expo", mass, slope);
+   RooRealVar ngauss("ngauss", "ngauss", 10, 0, 50);
+   RooRealVar nexpo("nexpo", "nexpo", 1e3);
+   RooAddPdf pdftot("pdftot", "pdftot", RooArgList(gauss, expo), RooArgList(ngauss, nexpo));
+
+   RooUniform prior("prior", "prior", ngauss);
+
+   // Run with the default interval method (numerical integration instead of
+   // brute-force scan).
+   RooStats::BayesianCalculator bc(*data, pdftot, RooArgSet(ngauss), prior);
+   bc.SetConfidenceLevel(0.95);
+   bc.SetLeftSideTailFraction(0.);
+   std::unique_ptr<RooStats::SimpleInterval> bcInterval{bc.GetInterval()};
+   ASSERT_NE(bcInterval, nullptr);
+   const double limitDefault = bcInterval->UpperLimit();
+   EXPECT_EQ(bcInterval->LowerLimit(), 0.0);
+
+   // Reference: explicit scan with many bins gives the same answer up to
+   // numerical integration and scan discretization errors.
+   RooStats::BayesianCalculator bcScan(*data, pdftot, RooArgSet(ngauss), prior);
+   bcScan.SetConfidenceLevel(0.95);
+   bcScan.SetLeftSideTailFraction(0.);
+   bcScan.SetScanOfPosterior(200);
+   std::unique_ptr<RooStats::SimpleInterval> scanInterval{bcScan.GetInterval()};
+   ASSERT_NE(scanInterval, nullptr);
+   const double limitScan = scanInterval->UpperLimit();
+
+   // Both limits should sit well inside the POI range (if the old bug were
+   // present the limit would be at or near poi->getMax() = 50, or pulled down
+   // by a broken CDF root) and must agree within 1%.
+   EXPECT_GT(limitDefault, 20.0);
+   EXPECT_LT(limitDefault, 35.0);
+   EXPECT_NEAR(limitDefault, limitScan, 0.01 * limitScan);
+}
diff --git a/roofit/xroofit/src/xRooNode.cxx b/roofit/xroofit/src/xRooNode.cxx
index 2bdeebf94c663..ec8d5034a8aa8 100644
--- a/roofit/xroofit/src/xRooNode.cxx
+++ b/roofit/xroofit/src/xRooNode.cxx
@@ -2136,7 +2136,7 @@ xRooNode xRooNode::Add(const xRooNode &child, Option_t *opt)
             // need to divide by bin widths first
             for (int i = 0; i < _f->dataHist().numEntries(); i++) {
                auto bin_pars = _f->dataHist().get(i);
-               _f->dataHist().set(*bin_pars, _f->dataHist().weight() / _f->dataHist().binVolume(*bin_pars));
+               _f->dataHist().set(*bin_pars, _f->dataHist().weight(i) / _f->dataHist().binVolume(*bin_pars));
             }
             _f->setAttribute("autodensity", false);
             _f->setValueDirty();
@@ -3329,7 +3329,7 @@ xRooNode xRooNode::Multiply(const xRooNode &child, Option_t *opt)
             // need to divide by bin widths first
             for (int i = 0; i < _f->dataHist().numEntries(); i++) {
                auto bin_pars = _f->dataHist().get(i);
-               _f->dataHist().set(*bin_pars, _f->dataHist().weight() / _f->dataHist().binVolume(*bin_pars));
+               _f->dataHist().set(*bin_pars, _f->dataHist().weight(i) / _f->dataHist().binVolume(*bin_pars));
             }
             _f->setValueDirty();
 
@@ -4893,7 +4893,7 @@ bool xRooNode::SetBinError(int bin, double value)
          TString origName = (f->getStringAttribute("origName")) ? f->getStringAttribute("origName") : GetName();
          rrv->setStringAttribute(Form("sumw2_%s", origName.Data()), TString::Format("%f", pow(value, 2)));
          auto bin_pars = f->dataHist().get(bin - 1);
-         auto _binContent = f->dataHist().weight();
+         auto _binContent = f->dataHist().weight(bin - 1);
          if (f->getAttribute("density")) {
             _binContent *= f->dataHist().binVolume(*bin_pars);
          }