From 7c74defce6d5e4c43403a2aded265b809648cff5 Mon Sep 17 00:00:00 2001
From: Jonas Rembser <jonas.rembser@cern.ch>
Date: Tue, 2 Dec 2025 21:33:50 +0100
Subject: [PATCH 1/6] [RF] More verbose errors in `RooPlot::findObject`

It greatly helps debugging if the error message tells you which objects
are actually existing.
---
 roofit/roofitcore/src/RooPlot.cxx | 32 ++++++++++++++++++-------------
 1 file changed, 19 insertions(+), 13 deletions(-)

diff --git a/roofit/roofitcore/src/RooPlot.cxx b/roofit/roofitcore/src/RooPlot.cxx
index 89cf8818476eb..a3a833b2e7cb0 100644
--- a/roofit/roofitcore/src/RooPlot.cxx
+++ b/roofit/roofitcore/src/RooPlot.cxx
@@ -898,22 +898,28 @@ bool RooPlot::drawAfter(const char *after, const char *target)
 /// methods to change the drawing style attributes of a contained
 /// object directly.
 
-TObject *RooPlot::findObject(const char *name, const TClass* tClass) const
+TObject *RooPlot::findObject(const char *name, const TClass *tClass) const
 {
-  TObject *ret = nullptr;
+   TObject *ret = nullptr;
 
-  for(auto const& item : _items) {
-    TObject &obj = *item.first;
-    if ((!name || name[0] == '\0' || !TString(name).CompareTo(obj.GetName()))
-        && (!tClass || (obj.IsA()==tClass))) {
-      ret = &obj ;
-    }
-  }
+   for (auto const &item : _items) {
+      TObject &obj = *item.first;
+      if ((!name || name[0] == '\0' || !TString(name).CompareTo(obj.GetName())) && (!tClass || (obj.IsA() == tClass))) {
+         ret = &obj;
+      }
+   }
 
-  if (ret == nullptr) {
-    coutE(InputArguments) << "RooPlot::findObject(" << GetName() << ") cannot find object " << (name?name:"<last>") << std::endl ;
-  }
-  return ret ;
+   if (ret == nullptr) {
+      std::stringstream error;
+      error << "RooPlot::findObject(" << GetName() << ") cannot find object " << (name ? name : "<last>") << "\n"
+            << "Available objects are:\n";
+      for (auto const &item : _items) {
+         TObject &obj = *item.first;
+         error << "  - " << obj.IsA()->GetName() << " \"" << obj.GetName() << "\"\n";
+      }
+      coutE(InputArguments) << error.str();
+   }
+   return ret;
 }
 
 

From 99cebb731f791ea23a818afa659ce4fcf9f53c57 Mon Sep 17 00:00:00 2001
From: Jonas Rembser <jonas.rembser@cern.ch>
Date: Tue, 2 Dec 2025 21:34:51 +0100
Subject: [PATCH 2/6] [RF] Avoid confusing `plotOn(Slice(sample),
 ProjWData(same))` pattern

In the documentation of the RooSimultaneous, it says about the
`ProjWData()` argument:

> For observables present in given dataset projection of PDF is achieved
  by constructing an average over all observable values in given set.

And about `Slice()`, it says:

> Override default projection behaviour by omitting the specified
  category observable from the projection, i.e., by not integrating over
  all states of this category.

Starting from this explanation, it is highly unintuitive that one should
do things like `simPdf.plotOn(Slice(sample, "state"), ProjWData(same))` when
plotting the pdf of state `"state"` from a RooSimultaneous.

I think that instead, we should promote easy patterns in the tutorials.

That is, if you want to plot a slice pdf from a RooSimultaneous, you are
retreiving it with `getPdf()` and then plot it. The result is the same.
---
 roofit/roofitcore/test/stressRooFit_tests.h   | 25 +++++++------
 .../roofit/roofit/rf501_simultaneouspdf.C     | 27 +++++++-------
 .../roofit/roofit/rf501_simultaneouspdf.py    | 35 +++++++++----------
 3 files changed, 42 insertions(+), 45 deletions(-)

diff --git a/roofit/roofitcore/test/stressRooFit_tests.h b/roofit/roofitcore/test/stressRooFit_tests.h
index 652414c05acd3..6cb169eedf035 100644
--- a/roofit/roofitcore/test/stressRooFit_tests.h
+++ b/roofit/roofitcore/test/stressRooFit_tests.h
@@ -2972,22 +2972,21 @@ class TestBasic501 : public RooUnitTest {
       RooPlot *frame1 = x.frame(Bins(30), Title("Physics sample"));
 
       // Plot all data tagged as physics sample
-      combData.plotOn(frame1, Cut("sample==sample::physics"));
+      std::unique_ptr<RooAbsData> slicedData1{combData.reduce(Cut("sample==sample::physics"))};
+      slicedData1->SetName("combData_Cut[sample==sample::physics]"); // to be consistent with reference file
+      slicedData1->plotOn(frame1);
 
       // Plot "physics" slice of simultaneous pdf.
-      // NBL You _must_ project the sample index category with data using ProjWData
-      // as a RooSimultaneous makes no prediction on the shape in the index category
-      // and can thus not be integrated
-      simPdf.plotOn(frame1, Slice(sample, "physics"), ProjWData(sample, combData));
-      simPdf.plotOn(frame1, Slice(sample, "physics"), Components("px"), ProjWData(sample, combData),
-                    LineStyle(kDashed));
-
-      // The same plot for the control sample slice
+      simPdf.getPdf("physics")->plotOn(frame1);
+      simPdf.getPdf("physics")->plotOn(frame1, Components("px"), LineStyle(kDashed));
+
+      // The same plot for the control sample slice.
       RooPlot *frame2 = x.frame(Bins(30), Title("Control sample"));
-      combData.plotOn(frame2, Cut("sample==sample::control"));
-      simPdf.plotOn(frame2, Slice(sample, "control"), ProjWData(sample, combData));
-      simPdf.plotOn(frame2, Slice(sample, "control"), Components("px_ctl"), ProjWData(sample, combData),
-                    LineStyle(kDashed));
+      std::unique_ptr<RooAbsData> slicedData2{combData.reduce(Cut("sample==sample::control"))};
+      slicedData2->SetName("combData_Cut[sample==sample::control]"); // to be consistent with reference file
+      slicedData2->plotOn(frame2);
+      simPdf.getPdf("control")->plotOn(frame2);
+      simPdf.getPdf("control")->plotOn(frame2, Components("px_ctl"), LineStyle(kDashed));
 
       regPlot(frame1, "rf501_plot1");
       regPlot(frame2, "rf501_plot2");
diff --git a/tutorials/roofit/roofit/rf501_simultaneouspdf.C b/tutorials/roofit/roofit/rf501_simultaneouspdf.C
index a658fc66ea844..fd666c6b88039 100644
--- a/tutorials/roofit/roofit/rf501_simultaneouspdf.C
+++ b/tutorials/roofit/roofit/rf501_simultaneouspdf.C
@@ -102,29 +102,28 @@ void rf501_simultaneouspdf()
    RooPlot *frame1 = x.frame(Title("Physics sample"));
 
    // Plot all data tagged as physics sample
-   combData.plotOn(frame1, Cut("sample==sample::physics"));
+   std::unique_ptr<RooAbsData> slicedData1{combData.reduce(Cut("sample==sample::physics"))};
+   slicedData1->plotOn(frame1);
 
    // Plot "physics" slice of simultaneous pdf.
+   simPdf.getPdf("physics")->plotOn(frame1);
+   simPdf.getPdf("physics")->plotOn(frame1, Components("px"), LineStyle(kDashed));
+
+   // The same plot for the control sample slice. We do this with a different
+   // approach, using the data slice for the projection. This approach is more
+   // general, because you can plot sums of slices by using logical or in the
+   // Cut() command.
    // NBL You _must_ project the sample index category with data using ProjWData
    // as a RooSimultaneous makes no prediction on the shape in the index category
    // and can thus not be integrated.
    // In other words: Since the PDF doesn't know the number of events in the different
    // category states, it doesn't know how much of each component it has to project out.
    // This information is read from the data.
-   simPdf.plotOn(frame1, Slice(sample, "physics"), ProjWData(sample, combData));
-   simPdf.plotOn(frame1, Slice(sample, "physics"), Components("px"), ProjWData(sample, combData), LineStyle(kDashed));
-
-   // The same plot for the control sample slice. We do this with a different
-   // approach this time, for illustration purposes. Here, we are slicing the
-   // dataset and then use the data slice for the projection, because then the
-   // RooFit::Slice() becomes unnecessary. This approach is more general,
-   // because you can plot sums of slices by using logical or in the Cut()
-   // command.
    RooPlot *frame2 = x.frame(Bins(30), Title("Control sample"));
-   std::unique_ptr<RooAbsData> slicedData{combData.reduce(Cut("sample==sample::control"))};
-   slicedData->plotOn(frame2);
-   simPdf.plotOn(frame2, ProjWData(sample, *slicedData));
-   simPdf.plotOn(frame2, Components("px_ctl"), ProjWData(sample, *slicedData), LineStyle(kDashed));
+   std::unique_ptr<RooAbsData> slicedData2{combData.reduce(Cut("sample==sample::control"))};
+   slicedData2->plotOn(frame2);
+   simPdf.plotOn(frame2, ProjWData(sample, *slicedData2));
+   simPdf.plotOn(frame2, Components("px_ctl"), ProjWData(sample, *slicedData2), LineStyle(kDashed));
 
    // The same plot for all the phase space. Here, we can just use the original
    // combined dataset.
diff --git a/tutorials/roofit/roofit/rf501_simultaneouspdf.py b/tutorials/roofit/roofit/rf501_simultaneouspdf.py
index 4c645b91a922c..6f4f81da61807 100644
--- a/tutorials/roofit/roofit/rf501_simultaneouspdf.py
+++ b/tutorials/roofit/roofit/rf501_simultaneouspdf.py
@@ -13,7 +13,6 @@
 
 import ROOT
 
-
 # Create model for physics sample
 # -------------------------------------------------------------
 
@@ -96,28 +95,28 @@
 frame1 = x.frame(Title="Physics sample")
 
 # Plot all data tagged as physics sample
-combData.plotOn(frame1, Cut="sample==sample::physics")
+slicedData1 = combData.reduce(Cut="sample==sample::physics")
+slicedData1.plotOn(frame1)
 
 # Plot "physics" slice of simultaneous pdf.
-# NB: You *must* project the sample index category with data using ProjWData as
-# a RooSimultaneous makes no prediction on the shape in the index category and
-# can thus not be integrated. In other words: Since the PDF doesn't know the
-# number of events in the different category states, it doesn't know how much
-# of each component it has to project out. This info is read from the data.
-simPdf.plotOn(frame1, Slice=(sample, "physics"), ProjWData=(sample, combData))
-simPdf.plotOn(frame1, Slice=(sample, "physics"), Components="px", ProjWData=(sample, combData), LineStyle="--")
+simPdf.getPdf("physics").plotOn(frame1)
+simPdf.getPdf("physics").plotOn(frame1, Components="px", LineStyle="--")
 
 # The same plot for the control sample slice. We do this with a different
-# approach this time, for illustration purposes. Here, we are slicing the
-# dataset and then use the data slice for the projection, because then the
-# RooFit::Slice() becomes unnecessary. This approach is more general,
-# because you can plot sums of slices by using logical or in the Cut()
-# command.
+# approach, using the data slice for the projection. This approach is more
+# general, because you can plot sums of slices by using logical or in the
+# Cut() command.
+# NBL You _must_ project the sample index category with data using ProjWData
+# as a RooSimultaneous makes no prediction on the shape in the index category
+# and can thus not be integrated.
+# In other words: Since the PDF doesn't know the number of events in the different
+# category states, it doesn't know how much of each component it has to project out.
+# This information is read from the data.
 frame2 = x.frame(Title="Control sample")
-slicedData = combData.reduce(Cut="sample==sample::control")
-slicedData.plotOn(frame2)
-simPdf.plotOn(frame2, ProjWData=(sample, slicedData))
-simPdf.plotOn(frame2, Components="px_ctl", ProjWData=(sample, slicedData), LineStyle="--")
+slicedData2 = combData.reduce(Cut="sample==sample::control")
+slicedData2.plotOn(frame2)
+simPdf.plotOn(frame2, ProjWData=(sample, slicedData2))
+simPdf.plotOn(frame2, Components="px_ctl", ProjWData=(sample, slicedData2), LineStyle="--")
 
 # The same plot for all the phase space. Here, we can just use the original
 # combined dataset.

From 99a80117920d79ca457bbabb4139d5574639050f Mon Sep 17 00:00:00 2001
From: Jonas Rembser <jonas.rembser@cern.ch>
Date: Tue, 2 Dec 2025 14:37:07 +0100
Subject: [PATCH 3/6] [RF] Don't double-scale pdfs plotted by
 RooSimultaneous::plotOn()

After a8ef8b0, when plotting an extended pdf it automatically scales
itself to the number of expected events.

However, when plotting slices of a RooSimultaneous, the normalization is
already precomputed in `RooSimultanous::plotOn()` and should not be
overridden again. To prevent this, flag the calculated normatization as
the final number of events, and not some relative factor, which it is by
default.

A unit test that covers this case was also implemented.

Closes #20383.
---
 roofit/roofitcore/src/RooSimultaneous.cxx     |   3 +-
 roofit/roofitcore/test/gtest_wrapper.h        |  24 ++-
 .../roofitcore/test/testRooSimultaneous.cxx   | 149 +++++++++++++++++-
 3 files changed, 171 insertions(+), 5 deletions(-)

diff --git a/roofit/roofitcore/src/RooSimultaneous.cxx b/roofit/roofitcore/src/RooSimultaneous.cxx
index 0b241f961c493..7169024564fde 100644
--- a/roofit/roofitcore/src/RooSimultaneous.cxx
+++ b/roofit/roofitcore/src/RooSimultaneous.cxx
@@ -777,7 +777,8 @@ RooPlot* RooSimultaneous::plotOn(RooPlot *frame, RooLinkedList& cmdList) const
        projDataTmp(projData->reduce(RooFit::SelectVars(projDataVars), RooFit::Cut(cutString.c_str())));
 
     // Override normalization and projection dataset
-    RooCmdArg tmp1 = RooFit::Normalization(scaleFactor*wTable->getFrac(idxCatClone->getCurrentLabel()),stype) ;
+    RooCmdArg tmp1 =
+       RooFit::Normalization(scaleFactor * wTable->get(idxCatClone->getCurrentLabel()), RooAbsReal::NumEvent);
     RooCmdArg tmp2 = RooFit::ProjWData(*projDataSet,*projDataTmp) ;
 
     // WVE -- do not adjust normalization for asymmetry plots
diff --git a/roofit/roofitcore/test/gtest_wrapper.h b/roofit/roofitcore/test/gtest_wrapper.h
index de31e562ad9f5..f342c55575ab8 100644
--- a/roofit/roofitcore/test/gtest_wrapper.h
+++ b/roofit/roofitcore/test/gtest_wrapper.h
@@ -32,6 +32,28 @@
 
 #define ROOFIT_EVAL_BACKENDS ROOFIT_EVAL_BACKEND_LEGACY ROOFIT_EVAL_BACKEND_CUDA RooFit::EvalBackend::Cpu()
 
-#define ROOFIT_EVAL_BACKENDS_WITH_CODEGEN ROOFIT_EVAL_BACKENDS, ROOFIT_EVAL_BACKEND_CODEGEN RooFit::EvalBackend::CodegenNoGrad()
+#define ROOFIT_EVAL_BACKENDS_WITH_CODEGEN \
+   ROOFIT_EVAL_BACKENDS, ROOFIT_EVAL_BACKEND_CODEGEN RooFit::EvalBackend::CodegenNoGrad()
+
+#include <gtest/gtest.h>
+#include <gmock/gmock.h>
+
+#include <cmath>
+
+MATCHER_P2(RelativeNear, expected, rel_tol,
+           "is within relative tolerance around ref=" + ::testing::PrintToString(expected) +
+              " (tol=" + ::testing::PrintToString(rel_tol) + ")")
+{
+   const double diff = std::fabs(arg - expected);
+   const double scale = std::fabs(expected);
+   // One could also also consider the target value as scale, or an adaptive scale:
+   // const double scale = std::max(std::fabs(arg), std::fabs(expected));
+
+   if (diff <= rel_tol * scale)
+      return true;
+
+   *result_listener << "error relative to ref = " << diff / scale << ", absolute diff = " << diff;
+   return false;
+}
 
 #endif // RooFit_gtest_wrapper_h
diff --git a/roofit/roofitcore/test/testRooSimultaneous.cxx b/roofit/roofitcore/test/testRooSimultaneous.cxx
index 14ac5e25ac074..548d695dc3ff3 100644
--- a/roofit/roofitcore/test/testRooSimultaneous.cxx
+++ b/roofit/roofitcore/test/testRooSimultaneous.cxx
@@ -5,10 +5,12 @@
 #include <RooAddPdf.h>
 #include <RooAddition.h>
 #include <RooCategory.h>
+#include <RooChebychev.h>
 #include <RooConstVar.h>
 #include <RooDataSet.h>
 #include <RooExponential.h>
 #include <RooFitResult.h>
+#include <RooGaussian.h>
 #include <RooGenericPdf.h>
 #include <RooHelpers.h>
 #include <RooMinimizer.h>
@@ -417,6 +419,8 @@ TEST(RooSimultaneous, ConditionalProdPdf)
 // channels can be extended. Also check if the likelihood can be created.
 TEST(RooSimultaneous, PartiallyExtendedPdfs)
 {
+   RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING);
+
    RooWorkspace ws;
    ws.factory("Gaussian::pdfA(x_a[-10, 10], mu_a[0, -10, 10], sigma_a[2.0, 0.1, 10.0])");
    ws.factory("Gaussian::pdfB(x_b[-10, 10], mu_b[0, -10, 10], sigma_b[2.0, 0.1, 10.0])");
@@ -427,7 +431,6 @@ TEST(RooSimultaneous, PartiallyExtendedPdfs)
    RooArgSet observables{*ws.var("x_a"), *ws.var("x_b"), *ws.cat("cat")};
 
    auto &simPdf = *ws.pdf("simPdf");
-   std::cout << simPdf.getVal() << std::endl;
 
    // A completely extended pdf, just to easily create a toy dataset
    ws.factory("ExtendPdf::pdfAext(pdfA, n_b[1000., 100., 10000.])");
@@ -518,6 +521,7 @@ TEST_P(TestStatisticTest, RooSimultaneousSingleChannelCrossCheckWithCondVar)
 /// RooSimultaneous with the new CPU backend.
 TEST(RooSimultaneous, RangedExtendedRooAddPdf)
 {
+   RooHelpers::LocalChangeMsgLevel changeMsgLevel{RooFit::WARNING};
 
    const double nBkgA_nom = 9000;
    const double nBkgB_nom = 10000;
@@ -549,11 +553,12 @@ TEST(RooSimultaneous, RangedExtendedRooAddPdf)
    simPdf.getParameters(combData->get(), params);
    params.snapshot(paramsSnap);
 
-   Res fitResult{simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(RooFit::EvalBackend::Cpu()))};
+   Res fitResult{simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(EvalBackend::Cpu()), PrintLevel(-1))};
 
    params.assign(paramsSnap);
 
-   Res fitResultRef{simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(RooFit::EvalBackend::Legacy()))};
+   Res fitResultRef{
+      simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(EvalBackend::Legacy()), PrintLevel(-1))};
 
    EXPECT_TRUE(fitResult->isIdentical(*fitResultRef));
 }
@@ -563,6 +568,8 @@ TEST(RooSimultaneous, RangedExtendedRooAddPdf)
 /// with a projection dataset.
 TEST(RooSimultaneous, PlotProjWData)
 {
+   RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING);
+
    RooRealVar x("x", "x", -8, 8);
    x.setBins(1);
 
@@ -590,11 +597,147 @@ TEST(RooSimultaneous, PlotProjWData)
    EXPECT_DOUBLE_EQ(frame->getCurve()->interpolate(0.), combData.sumEntries());
 }
 
+/// Second part of GitHub issue #20383.
+/// Check that the the simultaneous pdf is normalized correctly to the data
+/// when plotting with a projection dataset, in the extended and non-extended
+/// case, based on the reproducer provided by the user who opened the issue.
+TEST(RooSimultaneous, PlotProjWDataExtended)
+{
+   using namespace RooFit;
+
+   RooHelpers::LocalChangeMsgLevel changeMsgLevel{RooFit::WARNING};
+
+   RooRealVar xvar1("x1", "", 0.0, 10.0);
+   RooRealVar xvar2("x2", "", 0.0, 10.0);
+
+   RooRealVar mean1("mean1", "", 3.0, 1.0, 4.0);
+   RooRealVar mean2("mean2", "", 5.0, 4.0, 6.0);
+   RooRealVar sigma("sigma", "", 1.0, 0.01, 2.0);
+
+   RooGaussian gauss1("gauss1", "", xvar1, mean1, sigma);
+   RooGaussian gauss2("gauss2", "", xvar2, mean2, sigma);
+
+   RooRealVar a0("a0", "a0", -0.1, -1.0, 1.0);
+   RooChebychev bkg1("bkg1", "b1", xvar1, {a0});
+
+   RooRealVar c0("c0", "c0", -0.1, -1.0, 1.0);
+   RooChebychev bkg2("bkg2", "b2", xvar2, {c0});
+
+   RooRealVar s1("s1", "s1", 75.0, 0.0, 100000.0);
+   RooRealVar b1("b1", "b1", 25.0, 0.0, 100000.0);
+
+   // Extended models
+   RooAddPdf model1e("model1e", "", {gauss1, bkg1}, {s1, b1});
+
+   RooRealVar s2("s2", "s2", 50.0, 0.0, 100000.0);
+   RooRealVar b2("b2", "b2", 50.0, 0.0, 100000.0);
+
+   RooAddPdf model2e("model2e", "model2e", {gauss2, bkg2}, {s2, b2});
+
+   // Non-extended models
+   RooRealVar f1("f1", "f1", 0.75, 0.0, 1.0);
+   RooAddPdf model1n("model1n", "model1n", {gauss1, bkg1}, {f1});
+
+   RooRealVar f2("f2", "f2", 0.50, 0.0, 1.0);
+   RooAddPdf model2n("model2n", "model2n", RooArgList(gauss2, bkg2), RooArgList(f2));
+
+   // Case handling
+   enum class Case {
+      Gaussian,
+      NonExtended,
+      Extended
+   };
+
+   auto caseName = [](Case c) {
+      switch (c) {
+      case Case::Gaussian: return "Gaussian";
+      case Case::NonExtended: return "NonExtended";
+      case Case::Extended: return "Extended";
+      }
+      return "Unknown";
+   };
+
+   const std::vector<Case> cases = {Case::Gaussian, Case::NonExtended, Case::Extended};
+
+   // Helpers
+   auto integrateLastCurve = [](RooPlot *plot) {
+      const double xmin = plot->getPlotVar()->getMin();
+      const double xmax = plot->getPlotVar()->getMax();
+      return plot->getCurve()->average(xmin, xmax) * plot->getPlotVar()->numBins();
+   };
+
+   constexpr double tol = 0.01; // tolerate 1 % sampling error
+
+   // Test body
+   auto runCase = [&](Case c) {
+      SCOPED_TRACE(std::string("Case = ") + caseName(c));
+
+      RooAbsPdf *model1 = nullptr;
+      RooAbsPdf *model2 = nullptr;
+
+      switch (c) {
+      case Case::Gaussian:
+         model1 = &gauss1;
+         model2 = &gauss2;
+         break;
+      case Case::NonExtended:
+         model1 = &model1n;
+         model2 = &model2n;
+         break;
+      case Case::Extended:
+         model1 = &model1e;
+         model2 = &model2e;
+         break;
+      }
+
+      ASSERT_NE(model1, nullptr);
+      ASSERT_NE(model2, nullptr);
+
+      // Generate data
+      std::unique_ptr<RooDataSet> data1{model1->generate(xvar1, 10000)};
+      std::unique_ptr<RooDataSet> data2{model2->generate(xvar2, 1000)};
+
+      // Category
+      RooCategory sample("sample", "");
+      sample.defineType("Fit1");
+      sample.defineType("Fit2");
+
+      RooDataSet data("combinedData", "", {xvar1, xvar2}, Index(sample),
+                      Import({{"Fit1", data1.get()}, {"Fit2", data2.get()}}));
+
+      // Simultaneous PDF
+      RooSimultaneous sim_pdf("sim_pdf", "", sample);
+      sim_pdf.addPdf(*model1, "Fit1");
+      sim_pdf.addPdf(*model2, "Fit2");
+
+      std::unique_ptr<RooFitResult> result{sim_pdf.fitTo(data, Save(true), PrintLevel(-1))};
+
+      // Plot + checks
+      RooPlot *frame1 = xvar1.frame();
+      data.plotOn(frame1, Cut("sample==sample::Fit1"));
+      sim_pdf.plotOn(frame1, Slice(sample, "Fit1"), ProjWData(sample, data));
+
+      EXPECT_THAT(integrateLastCurve(frame1), RelativeNear(data1->sumEntries(), tol));
+
+      RooPlot *frame2 = xvar2.frame();
+      data.plotOn(frame2, Cut("sample==sample::Fit2"));
+      sim_pdf.plotOn(frame2, Slice(sample, "Fit2"), ProjWData(sample, data));
+
+      EXPECT_THAT(integrateLastCurve(frame2), RelativeNear(data2->sumEntries(), tol));
+   };
+
+   // Execute
+   for (Case c : cases) {
+      runCase(c);
+   }
+}
+
 /// JIRA ticket https://its.cern.ch/jira/browse/ROOT-7499
 /// Check that we can also generate Asimov datasets with non-integer weights
 /// via RooSimultaneous.
 TEST(RooSimultaneous, ExpectedDataWithNonIntegerWeights)
 {
+   RooHelpers::LocalChangeMsgLevel changeMsgLevel{RooFit::WARNING};
 
    RooWorkspace ws{"ws"};
    ws.factory("dummy_obs_a[0,1]");

From ca07a642791400f77e78f9fb7318392320854cdb Mon Sep 17 00:00:00 2001
From: Jonas Rembser <jonas.rembser@cern.ch>
Date: Sun, 21 Dec 2025 14:05:27 +0100
Subject: [PATCH 4/6] [RF] Some code modernization of RooFit result

  * smart pointers for transient members

  * more `std::string` instead of TString
---
 roofit/roofitcore/inc/RooFitResult.h   |  8 ++--
 roofit/roofitcore/src/RooFitResult.cxx | 52 ++++++++++----------------
 2 files changed, 23 insertions(+), 37 deletions(-)

diff --git a/roofit/roofitcore/inc/RooFitResult.h b/roofit/roofitcore/inc/RooFitResult.h
index 2d2181ee4adec..a13c9af5ad6da 100644
--- a/roofit/roofitcore/inc/RooFitResult.h
+++ b/roofit/roofitcore/inc/RooFitResult.h
@@ -33,9 +33,7 @@
 class RooArgSet ;
 class RooAbsPdf ;
 class RooPlot;
-class TObject ;
 class TH2 ;
-typedef RooArgSet* pRooArgSet ;
 
 class RooFitResult : public TNamed, public RooPrintable, public RooDirItem {
 public:
@@ -185,11 +183,11 @@ class RooFitResult : public TNamed, public RooPrintable, public RooDirItem {
   RooArgList* _initPars = nullptr;     ///< List of floating parameters with initial values
   RooArgList* _finalPars = nullptr;    ///< List of floating parameters with final values
 
-  mutable RooArgList* _globalCorr = nullptr;   ///<! List of global correlation coefficients
+  mutable std::unique_ptr<RooArgList> _globalCorr;   ///<! List of global correlation coefficients
   mutable TList       _corrMatrix ;   ///<! Correlation matrix (list of RooArgLists)
 
-  mutable RooArgList *_randomPars = nullptr; ///<! List of floating parameters with most recent random perturbation applied
-  mutable TMatrixF* _Lt = nullptr;           ///<! triangular matrix used for generate random perturbations
+  mutable std::unique_ptr<RooArgList> _randomPars; ///<! List of floating parameters with most recent random perturbation applied
+  mutable std::unique_ptr<TMatrixF> _Lt;           ///<! triangular matrix used for generate random perturbations
 
   TMatrixDSym* _CM = nullptr;  ///< Correlation matrix
   TMatrixDSym* _VM = nullptr;  ///< Covariance matrix
diff --git a/roofit/roofitcore/src/RooFitResult.cxx b/roofit/roofitcore/src/RooFitResult.cxx
index b86417a372c6c..8a13639a1500a 100644
--- a/roofit/roofitcore/src/RooFitResult.cxx
+++ b/roofit/roofitcore/src/RooFitResult.cxx
@@ -92,10 +92,10 @@ RooFitResult::RooFitResult(const RooFitResult &other)
 
   other._finalPars->snapshot(*_finalPars);
   if (other._randomPars) {
-    _randomPars = new RooArgList;
+    _randomPars = std::make_unique<RooArgList>();
     other._randomPars->snapshot(*_randomPars);
   }
-  if (other._Lt) _Lt = new TMatrix(*other._Lt);
+  if (other._Lt) _Lt = std::make_unique<TMatrix>(*other._Lt);
   if (other._VM) _VM = new TMatrixDSym(*other._VM) ;
   if (other._CM) _CM = new TMatrixDSym(*other._CM) ;
   if (other._GC) _GC = new TVectorD(*other._GC) ;
@@ -114,9 +114,6 @@ RooFitResult::~RooFitResult()
   if (_constPars) delete _constPars ;
   if (_initPars)  delete _initPars ;
   if (_finalPars) delete _finalPars ;
-  if (_globalCorr) delete _globalCorr;
-  if (_randomPars) delete _randomPars;
-  if (_Lt) delete _Lt;
   if (_CM) delete _CM ;
   if (_VM) delete _VM ;
   if (_GC) delete _GC ;
@@ -334,10 +331,10 @@ RooPlot *RooFitResult::plotOn(RooPlot *frame, const char *parName1, const char *
 const RooArgList& RooFitResult::randomizePars() const
 {
   Int_t nPar= _finalPars->size();
-  if(nullptr == _randomPars) { // first-time initialization
+  if(!_randomPars) { // first-time initialization
     assert(nullptr != _finalPars);
     // create the list of random values to fill
-    _randomPars = new RooArgList;
+    _randomPars = std::make_unique<RooArgList>();
     _finalPars->snapshot(*_randomPars);
     // calculate the elements of the upper-triangular matrix L that gives Lt*L = C
     // where Lt is the transpose of L (the "square-root method")
@@ -360,7 +357,7 @@ const RooArgList& RooFitResult::randomizePars() const
       }
     }
     // remember Lt
-    _Lt= new TMatrix(TMatrix::kTransposed,L);
+    _Lt= std::make_unique<TMatrix>(TMatrix::kTransposed,L);
   }
   else {
     // reset to the final fit values
@@ -455,7 +452,7 @@ const RooArgList* RooFitResult::globalCorr()
     fillLegacyCorrMatrix() ;
   }
 
-  return _globalCorr ;
+  return _globalCorr.get();
 }
 
 
@@ -622,39 +619,28 @@ void RooFitResult::fillLegacyCorrMatrix() const
   if (!_CM) return ;
 
   // Delete eventual previous correlation data holders
-  if (_globalCorr) delete _globalCorr ;
   _corrMatrix.Delete();
 
   // Build holding arrays for correlation coefficients
-  _globalCorr = new RooArgList("globalCorrelations") ;
+  _globalCorr = std::make_unique<RooArgList>("globalCorrelations") ;
 
   for(RooAbsArg* arg : *_initPars) {
     // Create global correlation value holder
-    TString gcName("GC[") ;
-    gcName.Append(arg->GetName()) ;
-    gcName.Append("]") ;
-    TString gcTitle(arg->GetTitle()) ;
-    gcTitle.Append(" Global Correlation") ;
-    _globalCorr->addOwned(std::make_unique<RooRealVar>(gcName.Data(),gcTitle.Data(),0.));
+    std::string argName = arg->GetName();
+    std::string argTitle = arg->GetTitle();
+    std::string gcName = "GC[" + argName + "]";
+    std::string gcTitle = argTitle + " Global Correlation";
+    _globalCorr->addOwned(std::make_unique<RooRealVar>(gcName.c_str(),gcTitle.c_str(),0.));
 
     // Create array with correlation holders for this parameter
-    TString name("C[") ;
-    name.Append(arg->GetName()) ;
-    name.Append(",*]") ;
-    RooArgList* corrMatrixRow = new RooArgList(name.Data()) ;
+    RooArgList* corrMatrixRow = new RooArgList(("C[" + argName + ",*]").c_str()) ;
     _corrMatrix.Add(corrMatrixRow) ;
     for(RooAbsArg* arg2 : *_initPars) {
 
-      TString cName("C[") ;
-      cName.Append(arg->GetName()) ;
-      cName.Append(",") ;
-      cName.Append(arg2->GetName()) ;
-      cName.Append("]") ;
-      TString cTitle("Correlation between ") ;
-      cTitle.Append(arg->GetName()) ;
-      cTitle.Append(" and ") ;
-      cTitle.Append(arg2->GetName()) ;
-      corrMatrixRow->addOwned(std::make_unique<RooRealVar>(cName.Data(),cTitle.Data(),0.));
+      std::string arg2Name = arg2->GetName();
+      std::string cName = "C[" + argName + "," + arg2Name + "]";
+      std::string cTitle = "Correlation between " + argName + " and " + arg2Name;
+      corrMatrixRow->addOwned(std::make_unique<RooRealVar>(cName.c_str(),cTitle.c_str(),0.));
     }
   }
 
@@ -1315,7 +1301,9 @@ void RooFitResult::Streamer(TBuffer &R__b)
       R__b >> _constPars;
       R__b >> _initPars;
       R__b >> _finalPars;
-      R__b >> _globalCorr;
+      RooArgList* globalCorr;
+      R__b >> globalCorr;
+      _globalCorr.reset(globalCorr);
       _corrMatrix.Streamer(R__b);
       R__b.CheckByteCount(R__s, R__c, RooFitResult::IsA());
 

From cef89dda50f0e736586690617e3146b3dc80efdf Mon Sep 17 00:00:00 2001
From: Jonas Rembser <jonas.rembser@cern.ch>
Date: Sun, 21 Dec 2025 14:07:27 +0100
Subject: [PATCH 5/6] [RF] Apply `git clang-format` to RooFitResult

---
 roofit/roofitcore/inc/RooFitResult.h   |  9 +--
 roofit/roofitcore/src/RooFitResult.cxx | 87 +++++++++++++-------------
 2 files changed, 49 insertions(+), 47 deletions(-)

diff --git a/roofit/roofitcore/inc/RooFitResult.h b/roofit/roofitcore/inc/RooFitResult.h
index a13c9af5ad6da..0eabae7df175e 100644
--- a/roofit/roofitcore/inc/RooFitResult.h
+++ b/roofit/roofitcore/inc/RooFitResult.h
@@ -33,7 +33,7 @@
 class RooArgSet ;
 class RooAbsPdf ;
 class RooPlot;
-class TH2 ;
+class TH2;
 
 class RooFitResult : public TNamed, public RooPrintable, public RooDirItem {
 public:
@@ -183,11 +183,12 @@ class RooFitResult : public TNamed, public RooPrintable, public RooDirItem {
   RooArgList* _initPars = nullptr;     ///< List of floating parameters with initial values
   RooArgList* _finalPars = nullptr;    ///< List of floating parameters with final values
 
-  mutable std::unique_ptr<RooArgList> _globalCorr;   ///<! List of global correlation coefficients
+  mutable std::unique_ptr<RooArgList> _globalCorr; ///<! List of global correlation coefficients
   mutable TList       _corrMatrix ;   ///<! Correlation matrix (list of RooArgLists)
 
-  mutable std::unique_ptr<RooArgList> _randomPars; ///<! List of floating parameters with most recent random perturbation applied
-  mutable std::unique_ptr<TMatrixF> _Lt;           ///<! triangular matrix used for generate random perturbations
+  mutable std::unique_ptr<RooArgList>
+     _randomPars;                        ///<! List of floating parameters with most recent random perturbation applied
+  mutable std::unique_ptr<TMatrixF> _Lt; ///<! triangular matrix used for generate random perturbations
 
   TMatrixDSym* _CM = nullptr;  ///< Correlation matrix
   TMatrixDSym* _VM = nullptr;  ///< Covariance matrix
diff --git a/roofit/roofitcore/src/RooFitResult.cxx b/roofit/roofitcore/src/RooFitResult.cxx
index 8a13639a1500a..a1900eb46f6aa 100644
--- a/roofit/roofitcore/src/RooFitResult.cxx
+++ b/roofit/roofitcore/src/RooFitResult.cxx
@@ -92,10 +92,11 @@ RooFitResult::RooFitResult(const RooFitResult &other)
 
   other._finalPars->snapshot(*_finalPars);
   if (other._randomPars) {
-    _randomPars = std::make_unique<RooArgList>();
-    other._randomPars->snapshot(*_randomPars);
+     _randomPars = std::make_unique<RooArgList>();
+     other._randomPars->snapshot(*_randomPars);
   }
-  if (other._Lt) _Lt = std::make_unique<TMatrix>(*other._Lt);
+  if (other._Lt)
+     _Lt = std::make_unique<TMatrix>(*other._Lt);
   if (other._VM) _VM = new TMatrixDSym(*other._VM) ;
   if (other._CM) _CM = new TMatrixDSym(*other._CM) ;
   if (other._GC) _GC = new TVectorD(*other._GC) ;
@@ -113,7 +114,8 @@ RooFitResult::~RooFitResult()
 {
   if (_constPars) delete _constPars ;
   if (_initPars)  delete _initPars ;
-  if (_finalPars) delete _finalPars ;
+  if (_finalPars)
+     delete _finalPars;
   if (_CM) delete _CM ;
   if (_VM) delete _VM ;
   if (_GC) delete _GC ;
@@ -331,37 +333,36 @@ RooPlot *RooFitResult::plotOn(RooPlot *frame, const char *parName1, const char *
 const RooArgList& RooFitResult::randomizePars() const
 {
   Int_t nPar= _finalPars->size();
-  if(!_randomPars) { // first-time initialization
-    assert(nullptr != _finalPars);
-    // create the list of random values to fill
-    _randomPars = std::make_unique<RooArgList>();
-    _finalPars->snapshot(*_randomPars);
-    // calculate the elements of the upper-triangular matrix L that gives Lt*L = C
-    // where Lt is the transpose of L (the "square-root method")
-    TMatrix L(nPar,nPar);
-    for(Int_t iPar= 0; iPar < nPar; iPar++) {
-      // calculate the diagonal term first
-      L(iPar,iPar)= covariance(iPar,iPar);
-      for(Int_t k= 0; k < iPar; k++) {
-   double tmp= L(k,iPar);
-   L(iPar,iPar)-= tmp*tmp;
-      }
-      L(iPar,iPar)= sqrt(L(iPar,iPar));
-      // then the off-diagonal terms
-      for(Int_t jPar= iPar+1; jPar < nPar; jPar++) {
-   L(iPar,jPar)= covariance(iPar,jPar);
-   for(Int_t k= 0; k < iPar; k++) {
-     L(iPar,jPar)-= L(k,iPar)*L(k,jPar);
-   }
-   L(iPar,jPar)/= L(iPar,iPar);
-      }
-    }
-    // remember Lt
-    _Lt= std::make_unique<TMatrix>(TMatrix::kTransposed,L);
-  }
-  else {
-    // reset to the final fit values
-    _randomPars->assign(*_finalPars);
+  if (!_randomPars) { // first-time initialization
+     assert(nullptr != _finalPars);
+     // create the list of random values to fill
+     _randomPars = std::make_unique<RooArgList>();
+     _finalPars->snapshot(*_randomPars);
+     // calculate the elements of the upper-triangular matrix L that gives Lt*L = C
+     // where Lt is the transpose of L (the "square-root method")
+     TMatrix L(nPar, nPar);
+     for (Int_t iPar = 0; iPar < nPar; iPar++) {
+        // calculate the diagonal term first
+        L(iPar, iPar) = covariance(iPar, iPar);
+        for (Int_t k = 0; k < iPar; k++) {
+           double tmp = L(k, iPar);
+           L(iPar, iPar) -= tmp * tmp;
+        }
+        L(iPar, iPar) = sqrt(L(iPar, iPar));
+        // then the off-diagonal terms
+        for (Int_t jPar = iPar + 1; jPar < nPar; jPar++) {
+           L(iPar, jPar) = covariance(iPar, jPar);
+           for (Int_t k = 0; k < iPar; k++) {
+              L(iPar, jPar) -= L(k, iPar) * L(k, jPar);
+           }
+           L(iPar, jPar) /= L(iPar, iPar);
+        }
+     }
+     // remember Lt
+     _Lt = std::make_unique<TMatrix>(TMatrix::kTransposed, L);
+  } else {
+     // reset to the final fit values
+     _randomPars->assign(*_finalPars);
   }
 
   // create a vector of unit Gaussian variables
@@ -622,7 +623,7 @@ void RooFitResult::fillLegacyCorrMatrix() const
   _corrMatrix.Delete();
 
   // Build holding arrays for correlation coefficients
-  _globalCorr = std::make_unique<RooArgList>("globalCorrelations") ;
+  _globalCorr = std::make_unique<RooArgList>("globalCorrelations");
 
   for(RooAbsArg* arg : *_initPars) {
     // Create global correlation value holder
@@ -630,17 +631,17 @@ void RooFitResult::fillLegacyCorrMatrix() const
     std::string argTitle = arg->GetTitle();
     std::string gcName = "GC[" + argName + "]";
     std::string gcTitle = argTitle + " Global Correlation";
-    _globalCorr->addOwned(std::make_unique<RooRealVar>(gcName.c_str(),gcTitle.c_str(),0.));
+    _globalCorr->addOwned(std::make_unique<RooRealVar>(gcName.c_str(), gcTitle.c_str(), 0.));
 
     // Create array with correlation holders for this parameter
-    RooArgList* corrMatrixRow = new RooArgList(("C[" + argName + ",*]").c_str()) ;
+    RooArgList *corrMatrixRow = new RooArgList(("C[" + argName + ",*]").c_str());
     _corrMatrix.Add(corrMatrixRow) ;
     for(RooAbsArg* arg2 : *_initPars) {
 
-      std::string arg2Name = arg2->GetName();
-      std::string cName = "C[" + argName + "," + arg2Name + "]";
-      std::string cTitle = "Correlation between " + argName + " and " + arg2Name;
-      corrMatrixRow->addOwned(std::make_unique<RooRealVar>(cName.c_str(),cTitle.c_str(),0.));
+       std::string arg2Name = arg2->GetName();
+       std::string cName = "C[" + argName + "," + arg2Name + "]";
+       std::string cTitle = "Correlation between " + argName + " and " + arg2Name;
+       corrMatrixRow->addOwned(std::make_unique<RooRealVar>(cName.c_str(), cTitle.c_str(), 0.));
     }
   }
 
@@ -1301,7 +1302,7 @@ void RooFitResult::Streamer(TBuffer &R__b)
       R__b >> _constPars;
       R__b >> _initPars;
       R__b >> _finalPars;
-      RooArgList* globalCorr;
+      RooArgList *globalCorr;
       R__b >> globalCorr;
       _globalCorr.reset(globalCorr);
       _corrMatrix.Streamer(R__b);

From 1d1b9d86963470637d8fc39226a0f14c9402c2b8 Mon Sep 17 00:00:00 2001
From: Jonas Rembser <jonas.rembser@cern.ch>
Date: Sun, 21 Dec 2025 14:46:51 +0100
Subject: [PATCH 6/6] [RF][HF] Some code cleanup for measurement classes

  * delete commented-out obsolete code

  * use more range-based loops

  * inline definitions where appropriate

  * initialize `Measurement` data member in class declaration
---
 .../inc/RooStats/HistFactory/Measurement.h    |  17 +-
 roofit/histfactory/src/Measurement.cxx        | 223 ++----------------
 2 files changed, 28 insertions(+), 212 deletions(-)

diff --git a/roofit/histfactory/inc/RooStats/HistFactory/Measurement.h b/roofit/histfactory/inc/RooStats/HistFactory/Measurement.h
index 8c8cf30dc0bd8..8907b200b8f3c 100644
--- a/roofit/histfactory/inc/RooStats/HistFactory/Measurement.h
+++ b/roofit/histfactory/inc/RooStats/HistFactory/Measurement.h
@@ -650,9 +650,9 @@ extern Channel BadChannel;
 class Measurement : public TNamed {
 
 public:
-   Measurement();
-   ///  Measurement( const Measurement& other ); // Copy
-   Measurement(const char *Name, const char *Title = "");
+   Measurement() = default;
+   /// Constructor specifying name and title of measurement
+   Measurement(const char *Name, const char *Title = "") : TNamed(Name, Title) {}
 
    ///  set output prefix
    void SetOutputFilePrefix(const std::string &prefix) { fOutputFilePrefix = prefix; }
@@ -724,7 +724,8 @@ class Measurement : public TNamed {
    std::vector<RooStats::HistFactory::Channel> &GetChannels() { return fChannels; }
    const std::vector<RooStats::HistFactory::Channel> &GetChannels() const { return fChannels; }
    RooStats::HistFactory::Channel &GetChannel(std::string);
-   void AddChannel(RooStats::HistFactory::Channel chan);
+   /// Add a completely configured channel.
+   void AddChannel(RooStats::HistFactory::Channel chan) { fChannels.push_back(chan); }
 
    bool HasChannel(std::string);
    void writeToFile(TFile *file);
@@ -752,10 +753,10 @@ class Measurement : public TNamed {
    /// Configurables of this measurement
    std::string fOutputFilePrefix;
    std::vector<std::string> fPOI;
-   double fLumi;
-   double fLumiRelErr;
-   int fBinLow;
-   int fBinHigh;
+   double fLumi = 1.0;
+   double fLumiRelErr = 0.1;
+   int fBinLow = 0;
+   int fBinHigh = 1;
    bool fExportOnly = true;
    std::string fInterpolationScheme;
 
diff --git a/roofit/histfactory/src/Measurement.cxx b/roofit/histfactory/src/Measurement.cxx
index bc7b84e43e673..9955131003768 100644
--- a/roofit/histfactory/src/Measurement.cxx
+++ b/roofit/histfactory/src/Measurement.cxx
@@ -41,29 +41,6 @@ uncertainty or the functional form of constraints on nuisance parameters.
 
 namespace RooStats::HistFactory {
 
-/// Standard constructor
-Measurement::Measurement() : fLumi(1.0), fLumiRelErr(.10), fBinLow(0), fBinHigh(1) {}
-
-/*
-Measurement::Measurement(const Measurement& other) :
-  POI( other.POI ), Lumi( other.Lumi ), LumiRelErr( other.LumiRelErr ),
-  BinLow( other.BinLow ), BinHigh( other.BinHigh ), ExportOnly( other.ExportOnly ),
-  channels( other.channels ), OutputFilePrefix( other.outputFilePrefix ),
-  constantParams( other.constantParams ), { ; }
-*/
-
-/// Standard constructor specifying name and title of measurement
-Measurement::Measurement(const char *Name, const char *Title)
-   : TNamed(Name, Title), fLumi(1.0), fLumiRelErr(.10), fBinLow(0), fBinHigh(1)
-{
-}
-
-/// add a completely configured channel
-void Measurement::AddChannel(Channel chan)
-{
-   fChannels.push_back(chan);
-}
-
 /// Set a parameter in the model to be constant.
 /// the parameter does not have to exist yet, the information will be used when
 /// the model is actually created.
@@ -183,13 +160,6 @@ Channel &Measurement::GetChannel(std::string ChanName)
    // return BadChannel;
 }
 
-/*
-  void Measurement::Print( Option_t* option ) const {
-  Measurement::Print( std::cout );
-  return;
-  }
-*/
-
 /// Print information about measurement object in tree-like structure to given stream
 void Measurement::PrintTree(std::ostream &stream)
 {
@@ -270,14 +240,7 @@ void Measurement::PrintXML(std::string directory, std::string newOutputPrefix)
    // Add the time
    xml << "<!--" << std::endl;
    xml << "This xml file created automatically on: " << std::endl;
-   /*
-     time_t t = time(0);   // get time now
-     struct tm * now = localtime( &t );
-     xml << (now->tm_year + 1900) << '-'
-         << (now->tm_mon + 1) << '-'
-         << now->tm_mday
-         << std::endl;
-   */
+
    // LM: use TTimeStamp
    TTimeStamp t;
    UInt_t year = 0;
@@ -301,12 +264,6 @@ void Measurement::PrintXML(std::string directory, std::string newOutputPrefix)
    for (unsigned int i = 0; i < fFunctionObjects.size(); ++i) {
       PreprocessFunction func = fFunctionObjects.at(i);
       func.PrintXML(xml);
-      /*
-      xml << "<Function Name=\"" << func.GetName() << "\" "
-     << "Expression=\""     << func.GetExpression() << "\" "
-     << "Dependents=\""     << func.GetDependents() << "\" "
-     << "/>" << std::endl;
-      */
    }
 
    xml << std::endl;
@@ -418,15 +375,13 @@ void Measurement::writeToFile(TFile *file)
    // Collect all histograms from file:
    // HistCollector collector;
 
-   for (unsigned int chanItr = 0; chanItr < outMeas.fChannels.size(); ++chanItr) {
-
+   for (Channel &channel : fChannels) {
       // Go to the main directory
       // in the file
       file->cd();
       file->Flush();
 
       // Get the name of the channel:
-      Channel &channel = outMeas.fChannels.at(chanItr);
       std::string chanName = channel.GetName();
 
       if (!channel.CheckHistograms()) {
@@ -461,39 +416,21 @@ void Measurement::writeToFile(TFile *file)
 
       channel.fData.writeToFile(OutputFileName, GetDirPath(dataDir));
 
-      /*
-      // Write the data file to this directory
-      TH1* hData = channel.data.GetHisto();
-      hData->Write();
-
-      // Set the location of the data
-      // in the output measurement
-
-      channel.data.InputFile = OutputFileName;
-      channel.data.HistoName = hData->GetName();
-      channel.data.HistoPath = GetDirPath( dataDir );
-      */
-
       // Loop over samples:
-
-      for (unsigned int sampItr = 0; sampItr < channel.GetSamples().size(); ++sampItr) {
-
-         Sample &sample = channel.GetSamples().at(sampItr);
-         std::string sampName = sample.GetName();
-
-         cxcoutPHF << "Writing sample: " << sampName << std::endl;
+      for (Sample &sample : channel.GetSamples()) {
+         cxcoutPHF << "Writing sample: " << sample.GetName() << std::endl;
 
          file->cd();
          chanDir->cd();
-         TDirectory *sampleDir = chanDir->mkdir(sampName.c_str());
+         TDirectory *sampleDir = chanDir->mkdir(sample.GetName().c_str());
          if (sampleDir == nullptr) {
-            cxcoutEHF << "Error: Directory " << sampName << " not created properly" << std::endl;
+            cxcoutEHF << "Error: Directory " << sample.GetName() << " not created properly" << std::endl;
             throw hf_exc();
          }
          std::string sampleDirPath = GetDirPath(sampleDir);
 
          if (!sampleDir) {
-            cxcoutEHF << "Error making directory: " << sampName << " in directory: " << chanName << std::endl;
+            cxcoutEHF << "Error making directory: " << sample.GetName() << " in directory: " << chanName << std::endl;
             throw hf_exc();
          }
 
@@ -501,59 +438,6 @@ void Measurement::writeToFile(TFile *file)
          sampleDir->cd();
 
          sample.writeToFile(OutputFileName, sampleDirPath);
-         /*
-         TH1* hSample = sample.GetHisto();
-         if( ! hSample ) {
-      std::cout << "Error getting histogram for sample: "
-           << sampName << std::endl;
-      throw -1;
-         }
-         sampleDir->cd();
-         hSample->Write();
-
-         sample.InputFile = OutputFileName;
-         sample.HistoName = hSample->GetName();
-         sample.HistoPath = sampleDirPath;
-         */
-
-         // Write the histograms associated with
-         // systematics
-
-         /*  THIS IS WHAT I"M COMMENTING
-         sample.GetStatError().writeToFile( OutputFileName, sampleDirPath );
-
-         // Must write all systematics that contain internal histograms
-         // (This is not all systematics)
-
-         for( unsigned int i = 0; i < sample.GetHistoSysList().size(); ++i ) {
-      sample.GetHistoSysList().at(i).writeToFile( OutputFileName, sampleDirPath );
-         }
-         for( unsigned int i = 0; i < sample.GetHistoFactorList().size(); ++i ) {
-      sample.GetHistoFactorList().at(i).writeToFile( OutputFileName, sampleDirPath );
-         }
-         for( unsigned int i = 0; i < sample.GetShapeSysList().size(); ++i ) {
-      sample.GetShapeSysList().at(i).writeToFile( OutputFileName, sampleDirPath );
-         }
-         END COMMENT  */
-         /*
-         sample.statError.writeToFile( OutputFileName, sampleDirPath );
-
-         // Now, get the Stat config histograms
-         if( sample.statError.HistoName != "" ) {
-      TH1* hStatError = sample.statError.GetErrorHist();
-      if( ! hStatError ) {
-        std::cout << "Error getting stat error histogram for sample: "
-             << sampName << std::endl;
-        throw -1;
-      }
-      hStatError->Write();
-
-      sample.statError.InputFile = OutputFileName;
-      sample.statError.HistoName = hStatError->GetName();
-      sample.statError.HistoPath = sampleDirPath;
-
-         }
-         */
       }
    }
 
@@ -571,7 +455,6 @@ void Measurement::writeToFile(TFile *file)
 /// stripped of unnecessary prefixes
 std::string Measurement::GetDirPath(TDirectory *dir)
 {
-
    std::string path = dir->GetPath();
 
    if (path.find(':') != std::string::npos) {
@@ -579,18 +462,7 @@ std::string Measurement::GetDirPath(TDirectory *dir)
       path.replace(0, index + 1, "");
    }
 
-   path = path + "/";
-
-   return path;
-
-   /*
-       if( path.find(":") != std::string::npos ) {
-    size_t index = path.find(":");
-    SampleName.replace( 0, index, "" );
-       }
-
-       // Remove the file:
-       */
+   return path + "/";
 }
 
 /// The most common way to add histograms to channels is to have them
@@ -603,11 +475,7 @@ std::string Measurement::GetDirPath(TDirectory *dir)
 /// save all necessary histograms.
 void Measurement::CollectHistograms()
 {
-
-   for (unsigned int chanItr = 0; chanItr < fChannels.size(); ++chanItr) {
-
-      Channel &chan = fChannels.at(chanItr);
-
+   for (Channel &chan : fChannels) {
       chan.CollectHistograms();
    }
 }
@@ -727,8 +595,6 @@ void Sample::writeToFile(std::string OutputFileName, std::string DirName)
    for (unsigned int i = 0; i < GetShapeFactorList().size(); ++i) {
       GetShapeFactorList().at(i).writeToFile(OutputFileName, DirName);
    }
-
-   return;
 }
 
 void Sample::SetValue(double val)
@@ -1200,10 +1066,7 @@ void Channel::CollectHistograms()
    }
 
    // Get the histograms for the samples:
-   for (unsigned int sampItr = 0; sampItr < fSamples.size(); ++sampItr) {
-
-      Sample &sample = fSamples.at(sampItr);
-
+   for (Sample &sample : fSamples) {
       // Get the nominal histogram:
       cxcoutDHF << "Collecting Nominal Histogram" << std::endl;
       TH1 *Nominal = GetHistogram(sample.GetInputFile(), sample.GetHistoPath(), sample.GetHistoName(), fileHandles);
@@ -1218,10 +1081,7 @@ void Channel::CollectHistograms()
       }
 
       // Get the HistoSys Variations:
-      for (unsigned int histoSysItr = 0; histoSysItr < sample.GetHistoSysList().size(); ++histoSysItr) {
-
-         HistoSys &histoSys = sample.GetHistoSysList().at(histoSysItr);
-
+      for (HistoSys &histoSys : sample.GetHistoSysList()) {
          histoSys.SetHistoLow(GetHistogram(histoSys.GetInputFileLow(), histoSys.GetHistoPathLow(),
                                            histoSys.GetHistoNameLow(), fileHandles));
 
@@ -1230,10 +1090,7 @@ void Channel::CollectHistograms()
       } // End Loop over HistoSys
 
       // Get the HistoFactor Variations:
-      for (unsigned int histoFactorItr = 0; histoFactorItr < sample.GetHistoFactorList().size(); ++histoFactorItr) {
-
-         HistoFactor &histoFactor = sample.GetHistoFactorList().at(histoFactorItr);
-
+      for (HistoFactor &histoFactor : sample.GetHistoFactorList()) {
          histoFactor.SetHistoLow(GetHistogram(histoFactor.GetInputFileLow(), histoFactor.GetHistoPathLow(),
                                               histoFactor.GetHistoNameLow(), fileHandles));
 
@@ -1242,19 +1099,13 @@ void Channel::CollectHistograms()
       } // End Loop over HistoFactor
 
       // Get the ShapeSys Variations:
-      for (unsigned int shapeSysItr = 0; shapeSysItr < sample.GetShapeSysList().size(); ++shapeSysItr) {
-
-         ShapeSys &shapeSys = sample.GetShapeSysList().at(shapeSysItr);
-
+      for (ShapeSys &shapeSys : sample.GetShapeSysList()) {
          shapeSys.SetErrorHist(
             GetHistogram(shapeSys.GetInputFile(), shapeSys.GetHistoPath(), shapeSys.GetHistoName(), fileHandles));
       } // End Loop over ShapeSys
 
       // Get any initial shape for a ShapeFactor
-      for (unsigned int shapeFactorItr = 0; shapeFactorItr < sample.GetShapeFactorList().size(); ++shapeFactorItr) {
-
-         ShapeFactor &shapeFactor = sample.GetShapeFactorList().at(shapeFactorItr);
-
+      for (ShapeFactor &shapeFactor : sample.GetShapeFactorList()) {
          // Check if we need an InitialShape
          if (shapeFactor.HasInitialShape()) {
             TH1 *hist = GetHistogram(shapeFactor.GetInputFile(), shapeFactor.GetHistoPath(), shapeFactor.GetHistoName(),
@@ -1279,10 +1130,7 @@ bool Channel::CheckHistograms() const
    }
 
    // Get the histograms for the samples:
-   for (unsigned int sampItr = 0; sampItr < fSamples.size(); ++sampItr) {
-
-      const Sample &sample = fSamples.at(sampItr);
-
+   for (Sample const &sample : fSamples) {
       // Get the nominal histogram:
       if (sample.GetHisto() == nullptr) {
          cxcoutEHF << "Error: Nominal Histogram for sample " << sample.GetName() << " is nullptr." << std::endl;
@@ -1322,9 +1170,7 @@ bool Channel::CheckHistograms() const
       }
 
       // Get the HistoSys Variations:
-      for (unsigned int histoSysItr = 0; histoSysItr < sample.GetHistoSysList().size(); ++histoSysItr) {
-
-         const HistoSys &histoSys = sample.GetHistoSysList().at(histoSysItr);
+      for (const HistoSys &histoSys : sample.GetHistoSysList()) {
 
          if (histoSys.GetHistoLow() == nullptr) {
             cxcoutEHF << "Error: HistoSyst Low for Systematic " << histoSys.GetName() << " in sample "
@@ -1340,9 +1186,7 @@ bool Channel::CheckHistograms() const
       } // End Loop over HistoSys
 
       // Get the HistoFactor Variations:
-      for (unsigned int histoFactorItr = 0; histoFactorItr < sample.GetHistoFactorList().size(); ++histoFactorItr) {
-
-         const HistoFactor &histoFactor = sample.GetHistoFactorList().at(histoFactorItr);
+      for (const HistoFactor &histoFactor : sample.GetHistoFactorList()) {
 
          if (histoFactor.GetHistoLow() == nullptr) {
             cxcoutEHF << "Error: HistoSyst Low for Systematic " << histoFactor.GetName() << " in sample "
@@ -1358,16 +1202,12 @@ bool Channel::CheckHistograms() const
       } // End Loop over HistoFactor
 
       // Get the ShapeSys Variations:
-      for (unsigned int shapeSysItr = 0; shapeSysItr < sample.GetShapeSysList().size(); ++shapeSysItr) {
-
-         const ShapeSys &shapeSys = sample.GetShapeSysList().at(shapeSysItr);
-
+      for (const ShapeSys &shapeSys : sample.GetShapeSysList()) {
          if (shapeSys.GetErrorHist() == nullptr) {
             cxcoutEHF << "Error: HistoSyst High for Systematic " << shapeSys.GetName() << " in sample "
                       << sample.GetName() << " is nullptr." << std::endl;
             return false;
          }
-
       } // End Loop over ShapeSys
 
    } // End Loop over Samples
@@ -1614,10 +1454,7 @@ void Asimov::ConfigureWorkspace(RooWorkspace *wspace)
    // Then, we set any variables to constant
    //
 
-   for (std::map<std::string, bool>::iterator itr = fParamsToFix.begin(); itr != fParamsToFix.end(); ++itr) {
-
-      std::string param = itr->first;
-      bool isConstant = itr->second;
+   for (auto &[param, isConstant] : fParamsToFix) {
 
       // Try to get the variable in the workspace
       RooRealVar *var = wspace->var(param);
@@ -1630,8 +1467,6 @@ void Asimov::ConfigureWorkspace(RooWorkspace *wspace)
       std::cout << "Configuring Asimov Dataset: Setting " << param << " to constant " << std::endl;
       var->setConstant(isConstant);
    }
-
-   return;
 }
 
 /** \class RooStats::HistFactory::HistRef
@@ -1737,8 +1572,6 @@ Constraint::Type Constraint::GetType(const std::string &Name)
    }
 }
 
-// Norm Factor
-
 void NormFactor::Print(std::ostream &stream) const
 {
    stream << "\t \t Name: " << fName << "\t Val: " << fVal << "\t Low: " << fLow << "\t High: " << fHigh << std::endl;
@@ -1753,7 +1586,6 @@ void NormFactor::PrintXML(std::ostream &xml) const
        << "  /> " << std::endl;
 }
 
-// Overall Sys
 void OverallSys::Print(std::ostream &stream) const
 {
    stream << "\t \t Name: " << fName << "\t Low: " << fLow << "\t High: " << fHigh << std::endl;
@@ -1814,8 +1646,6 @@ void HistogramUncertaintyBase::writeToFile(const std::string &FileName, const st
    fInputFileHigh = FileName;
    fHistoPathHigh = DirName;
    fHistoNameHigh = histHigh->GetName();
-
-   return;
 }
 
 void HistoSys::PrintXML(std::ostream &xml) const
@@ -1831,8 +1661,6 @@ void HistoSys::PrintXML(std::ostream &xml) const
        << "  /> " << std::endl;
 }
 
-// Shape Sys
-
 void ShapeSys::SetErrorHist(TH1 *hError)
 {
    fhHigh.reset(hError);
@@ -1856,7 +1684,6 @@ void ShapeSys::PrintXML(std::ostream &xml) const
 
 void ShapeSys::writeToFile(const std::string &FileName, const std::string &DirName)
 {
-
    auto histError = GetErrorHist();
    if (histError == nullptr) {
       std::cout << "Error: Cannot write " << GetName() << " to file: " << FileName << " ErrorHist is nullptr"
@@ -1867,12 +1694,8 @@ void ShapeSys::writeToFile(const std::string &FileName, const std::string &DirNa
    fInputFileHigh = FileName;
    fHistoPathHigh = DirName;
    fHistoNameHigh = histError->GetName();
-
-   return;
 }
 
-// HistoFactor
-
 void HistoFactor::PrintXML(std::ostream &xml) const
 {
    xml << "      <HistoFactor Name=\"" << GetName() << "\" "
@@ -1887,8 +1710,6 @@ void HistoFactor::PrintXML(std::ostream &xml) const
        << "  /> " << std::endl;
 }
 
-// Shape Factor
-
 void ShapeFactor::SetInitialShape(TH1 *shape)
 {
    fhHigh.reset(shape);
@@ -1925,8 +1746,6 @@ void ShapeFactor::writeToFile(const std::string &FileName, const std::string &Di
       fHistoPathHigh = DirName;
       fHistoNameHigh = histInitialShape->GetName();
    }
-
-   return;
 }
 
 void ShapeFactor::PrintXML(std::ostream &xml) const
@@ -1945,7 +1764,6 @@ void StatError::SetErrorHist(TH1 *Error)
    fhHigh.reset(Error);
 }
 
-// Stat Error Config
 void StatErrorConfig::Print(std::ostream &stream) const
 {
    stream << "\t \t RelErrorThreshold: " << fRelErrorThreshold
@@ -1960,7 +1778,6 @@ void StatErrorConfig::PrintXML(std::ostream &xml) const
        << std::endl;
 }
 
-// Stat Error
 void StatError::Print(std::ostream &stream) const
 {
    stream << "\t \t Activate: " << fActivate << "\t InputFile: " << fInputFileHigh << "\t HistoName: " << fHistoNameHigh
@@ -1997,8 +1814,6 @@ void StatError::writeToFile(const std::string &OutputFileName, const std::string
       fHistoNameHigh = statErrorHistName;
       fHistoPathHigh = DirName;
    }
-
-   return;
 }
 
 HistogramUncertaintyBase::HistogramUncertaintyBase(const HistogramUncertaintyBase &oth)