[RF] Always save generated parameters in RooMCStudy

Fix three related problems around RooMCStudy::genParDataSet(), reported
in GitHub issue #9490 and its discussion thread:

1. The generated parameter values were only recorded when internal
   constraints were used via Constrain(), despite the documentation
   promising that they are saved. The _genParData dataset is now always
   created and filled with one entry per generated toy.

2. Parameters constrained with ExternalConstraints() were not sampled
   from the constraint p.d.f.s for each toy, because only
   RooAbsPdf::getAllConstraints() was consulted, which cannot see
   constraint terms that are not part of the model. The external
   constraint p.d.f.s are now included in the per-toy sampling, with the
   constrained parameters found by intersecting the constraint
   observables with the model parameters. Passing the constraints inside
   FitOptions() remains available to apply them in the fit only.

3. When some toy fits failed, merging the generated parameter values
   into the fit parameter dataset failed with "ERROR: datasets have
   different size", because _genParData has one entry per generated toy
   while _fitParData only gets entries for converged fits. The merge now
   goes through a copy that only contains the toys whose fit converged,
   so the "<name>_gen" columns stay row-aligned with fitParDataSet(),
   while genParDataSet() keeps the entries of all generated toys
   (including the failed ones, so no information is dropped).

Further consequences of the implementation:

- genParDataSet() columns are no longer renamed in place to
  "<name>_gen"; only the merged copies are. This also fixes repeated
  runs accumulating "_gen" suffixes.
- The merge only happens for runs that both generate and fit, and only
  when the parameters are sampled from constraint p.d.f.s. This avoids
  clobbering the "<name>_gen" columns published by
  RooRandomizeParamMCSModule and keeps the fitParDataSet() schema
  unchanged for unconstrained studies. Fit-only runs previously merged
  the columns in reverse toy order, silently producing wrong pulls; they
  now emit a warning instead.
- _genParData is filled after the processBeforeGen() module hook, so it
  records the values actually used for generation.
- generate() and generateAndFit() reset _genParData like the other
  result containers.

Closes #9490

🤖 Done with the help of AI

Author: guitargeek

roofit/roofitcore/inc/RooMCStudy.h

@@ -60,8 +60,11 @@ class RooMCStudy : public TNamed {
   const RooFitResult* fitResult(Int_t sampleNum) const ;
         RooAbsData* genData(Int_t sampleNum) const ;
   const RooDataSet& fitParDataSet() ;
-  /// Return dataset with generator parameters for each toy. When constraints are used these
-  /// may generally not be the same as the fitted parameters.
+  /// Return dataset with the generator parameter values used for each toy, including any
+  /// modification by constraint-p.d.f. sampling or by study modules. When constraints are
+  /// used, the values are sampled from the constraint p.d.f.s for each toy and thus generally
+  /// differ from the fitted parameter values. Unlike fitParDataSet(), this dataset has one
+  /// entry per generated toy, including toys for which the fit did not converge.
   const RooDataSet* genParDataSet() const {
     return _genParData.get();
   }

roofit/roofitcore/src/RooMCStudy.cxx

@@ -95,6 +95,8 @@ fitting the PDF to data and accumulating the fit statistics.
 <tr><td> Verbose(bool flag)              <td> Activate informational messages in event generation phase
 <tr><td> Extended(bool flag)             <td> Determine number of events for each sample anew from a Poisson distribution
 <tr><td> Constrain(const RooArgSet& pars)  <td> Apply internal constraints on given parameters in fit and sample constrained parameter values from constraint p.d.f for each toy.
+<tr><td> ExternalConstraints(const RooArgSet& cpdfs) <td> Apply given external constraint p.d.f.s in fit and sample values of the parameters they constrain from them for each toy.
+                                                To apply the constraints in the fit only, without the per-toy sampling, pass them inside FitOptions() instead.
 <tr><td> ProtoData(const RooDataSet&, bool randOrder)
          <td> Prototype data for the event generation. If the randOrder flag is set, the order of the dataset will be re-randomized for each generation
               cycle to protect against systematic biases if the number of generated events does not exactly match the number of events in the prototype dataset
@@ -169,28 +171,31 @@ RooMCStudy::RooMCStudy(const RooAbsPdf& model, const RooArgSet& observables,
     _fitOptList.Add(RooFit::ExternalConstraints(*extCons).Clone()) ;
   }
 
-  // Make list of all constraints
+  // Make list of all constraints and of the parameters they constrain
   RooArgSet allConstraints ;
   RooArgSet consPars ;
   if (cPars) {
     if (std::unique_ptr<RooArgSet> constraints{model.getAllConstraints(observables,*cPars,true)}) {
       allConstraints.add(*constraints) ;
     }
+    consPars.add(*cPars) ;
+  }
+  if (extCons) {
+    // External constraint p.d.f.s are not part of the model, so the parameters
+    // they constrain are found among their observables instead
+    allConstraints.add(*extCons) ;
+    RooArgSet params;
+    model.getParameters(&observables, params);
+    for (RooAbsArg const* con : *extCons) {
+      RooArgSet cparams;
+      con->getObservables(&params, cparams);
+      consPars.add(cparams, /*silent=*/true) ;
+    }
   }
 
   // Construct constraint p.d.f
   if (!allConstraints.empty()) {
     _constrPdf = std::make_unique<RooProdPdf>("mcs_constr_prod","RooMCStudy constraints product",allConstraints);
-
-    if (cPars) {
-      consPars.add(*cPars) ;
-    } else {
-      RooArgSet params;
-      model.getParameters(&observables, params);
-      RooArgSet cparams;
-      _constrPdf->getObservables(&params, cparams);
-      consPars.add(cparams) ;
-    }
     _constrGenContext.reset(_constrPdf->genContext(consPars,nullptr,nullptr,_verboseGen));
 
     _perExptGenParams = true ;
@@ -261,9 +266,7 @@ RooMCStudy::RooMCStudy(const RooAbsPdf& model, const RooArgSet& observables,
   tmp2.setAttribAll("StoreError",false) ;
   tmp2.setAttribAll("StoreAsymError",false) ;
 
-  if (_perExptGenParams) {
-    _genParData = std::make_unique<RooDataSet>("genParData","Generated Parameters dataset",_genParams);
-  }
+  _genParData = std::make_unique<RooDataSet>("genParData","Generated Parameters dataset",_genParams);
 
   // Append proto variables to allDependents
   if (_genProtoData) {
@@ -309,6 +312,11 @@ void RooMCStudy::addModule(RooAbsMCStudyModule& module)
 /// If keepGenData is set, all generated data sets will be kept in memory and can be accessed
 /// later via genData().
 ///
+/// When generating, the generator parameter values used for each sample are recorded in the
+/// dataset returned by genParDataSet(). When constraints are used and the run both generates
+/// and fits, the sampled parameter values are in addition merged into the fit parameter
+/// dataset as `<name>_gen` columns for each toy whose fit converged.
+///
 /// When generating, data sets will be written out in ascii form if the pattern string is supplied
 /// The pattern, which is a template for snprintf, should look something like "data/toymc_%04d.dat"
 /// and should contain one integer field that encodes the sample serial number.
@@ -329,8 +337,25 @@ bool RooMCStudy::run(bool doGenerate, bool DoFit, Int_t nSamples, Int_t nEvtPerS
     mod->initializeRun(nSamples) ;
   }
 
+  if (DoFit && !doGenerate && _perExptGenParams) {
+    coutW(Generation) << "RooMCStudy::run: WARNING fitting previously generated samples in a separate run:"
+        " the per-toy sampled generator parameters are not merged into the fit parameter dataset,"
+        " so pulls are computed with respect to the initial parameter values instead of the sampled ones" << std::endl ;
+  }
+
   int prescale = nSamples>100 ? int(nSamples/100) : 1 ;
 
+  // Generator parameter values of the toys whose fit converged, filled in the
+  // same order as _fitParData so that the two datasets can be merged after the
+  // loop. Only done when the parameters are sampled from constraint p.d.f.s:
+  // otherwise the values are the constant initial ones, and study modules like
+  // RooRandomizeParamMCSModule publish their own "<name>_gen" columns that
+  // must not be overwritten by the merge.
+  std::unique_ptr<RooDataSet> genParDataForMerge;
+  if (doGenerate && DoFit && _perExptGenParams && _genParData) {
+    genParDataForMerge = std::make_unique<RooDataSet>("genParDataForMerge","Generated Parameters dataset",*_genParData->get());
+  }
+
   while(nSamples--) {
 
     if (nSamples%prescale==0) {
@@ -356,16 +381,17 @@ bool RooMCStudy::run(bool doGenerate, bool DoFit, Int_t nSamples, Int_t nEvtPerS
         _genParams.assign(*std::unique_ptr<RooDataSet>{_constrGenContext->generate(1)}->get());
       }
 
-      // Save generated parameters if required
-      if (_genParData) {
-   _genParData->add(_genParams) ;
-      }
-
       // Call module before-generation hook
       for (RooAbsMCStudyModule *mod : _modList) {
          mod->processBeforeGen(nSamples) ;
       }
 
+      // Save the generator parameters used for this toy, including any
+      // modification applied by the study modules above
+      if (_genParData) {
+   _genParData->add(_genParams) ;
+      }
+
       if (_binGenData) {
 
    // Calculate the number of (extended) events for this run
@@ -436,6 +462,14 @@ bool RooMCStudy::run(bool doGenerate, bool DoFit, Int_t nSamples, Int_t nEvtPerS
     bool fitOk = true;
     if (DoFit) fitOk = !fitSample(_genSample) ;
 
+    // Keep the generator parameters of this toy for merging into the fit
+    // parameter dataset. Only converged fits get an entry in _fitParData, so
+    // the toys with failed fits have to be skipped here as well. The values
+    // are taken from _genParData because the fit changes the parameters.
+    if (genParDataForMerge && _genParData && fitOk) {
+      genParDataForMerge->add(*_genParData->get(_genParData->numEntries()-1)) ;
+    }
+
     // Call module between generation and fitting hook
     for (RooAbsMCStudyModule *mod : _modList) {
       mod->processAfterFit(fitOk) ;
@@ -468,12 +502,14 @@ bool RooMCStudy::run(bool doGenerate, bool DoFit, Int_t nSamples, Int_t nEvtPerS
 
   _canAddFitResults = false ;
 
-  if (_genParData) {
-    for(RooAbsArg * arg : *_genParData->get()) {
-      _genParData->changeObservableName(arg->GetName(),(std::string(arg->GetName()) + "_gen").c_str());
+  if (genParDataForMerge) {
+    // Append the generator parameter values as additional "<name>_gen"
+    // columns to the fit parameter dataset
+    for(RooAbsArg * arg : *genParDataForMerge->get()) {
+      genParDataForMerge->changeObservableName(arg->GetName(),(std::string(arg->GetName()) + "_gen").c_str());
     }
 
-    _fitParData->merge(_genParData.get());
+    _fitParData->merge(genParDataForMerge.get());
   }
 
   if (DoFit) calcPulls() ;
@@ -506,6 +542,7 @@ bool RooMCStudy::generateAndFit(Int_t nSamples, Int_t nEvtPerSample, bool keepGe
   _fitResList.Delete() ; // even though the fit results are owned by gROOT, we still want to scratch them here.
   _genDataList.Delete() ;
   _fitParData->reset() ;
+  if (_genParData) _genParData->reset() ;
 
   return run(true,true,nSamples,nEvtPerSample,keepGenData,asciiFilePat) ;
 }
@@ -526,6 +563,7 @@ bool RooMCStudy::generate(Int_t nSamples, Int_t nEvtPerSample, bool keepGenData,
 {
   // Clear any previous data in memory
   _genDataList.Delete() ;
+  if (_genParData) _genParData->reset() ;
 
   return run(true,false,nSamples,nEvtPerSample,keepGenData,asciiFilePat) ;
 }

roofit/roofitcore/test/CMakeLists.txt

@@ -82,6 +82,7 @@ if(clad)
 endif()
 ROOT_ADD_GTEST(testNaNPacker testNaNPacker.cxx LIBRARIES RooFitCore)
 ROOT_ADD_GTEST(testRooExtendedBinding testRooExtendedBinding.cxx LIBRARIES RooFitCore RooFit)
+ROOT_ADD_GTEST(testRooMCStudy testRooMCStudy.cxx LIBRARIES RooFitCore RooFit)
 ROOT_ADD_GTEST(testRooMinimizer testRooMinimizer.cxx LIBRARIES RooFitCore RooFit)
 ROOT_ADD_GTEST(testRooMulti testRooMulti.cxx LIBRARIES RooFitCore RooFit)
 ROOT_ADD_GTEST(testRooRombergIntegrator testRooRombergIntegrator.cxx LIBRARIES MathCore RooFitCore)

roofit/roofitcore/test/testRooMCStudy.cxx

@@ -0,0 +1,175 @@
+// Tests for RooMCStudy
+// Authors: Jonas Rembser, CERN 2026
+
+#include <RooAbsPdf.h>
+#include <RooArgSet.h>
+#include <RooDataSet.h>
+#include <RooHelpers.h>
+#include <RooMCStudy.h>
+#include <RooRandom.h>
+#include <RooRealVar.h>
+#include <RooWorkspace.h>
+
+#include "gtest_wrapper.h"
+
+#include <algorithm>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+namespace {
+
+void fillModel(RooWorkspace &ws)
+{
+   ws.factory("x[-10, 10]");
+   ws.factory("Gaussian::pdf1(x, m[-1, 1], s[5, 10])");
+   ws.factory("Gaussian::pdf2(x, m, s2[1, 3])");
+   ws.factory("SUM::pdf(N1[0, 100] * pdf1, N2[0, 100] * pdf2)");
+}
+
+std::vector<double> getColumn(RooDataSet const &data, const char *name)
+{
+   auto *var = static_cast<RooRealVar const *>(data.get()->find(name));
+   if (var == nullptr) {
+      throw std::runtime_error(std::string{"dataset has no column named \""} + name + "\"");
+   }
+   std::vector<double> out;
+   out.reserve(data.numEntries());
+   for (int i = 0; i < data.numEntries(); ++i) {
+      data.get(i);
+      out.push_back(var->getVal());
+   }
+   return out;
+}
+
+} // namespace
+
+/// Covers GitHub issue #9490: the generated parameter values must be saved
+/// also when no constraints are used.
+TEST(RooMCStudy, GenParDataSetNoConstraints)
+{
+   RooRandom::randomGenerator()->SetSeed(4357);
+   RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING};
+
+   RooWorkspace ws;
+   fillModel(ws);
+
+   RooMCStudy mcstudy{*ws.pdf("pdf"), *ws.var("x"), RooFit::Silence()};
+   mcstudy.generate(3, 100, true);
+
+   RooDataSet const *genParData = mcstudy.genParDataSet();
+   ASSERT_NE(genParData, nullptr);
+   EXPECT_EQ(genParData->numEntries(), 3);
+   // The columns keep the original parameter names
+   EXPECT_NE(genParData->get()->find("s"), nullptr);
+}
+
+/// Covers GitHub issue #9490: parameters constrained by external constraint
+/// p.d.f.s must be sampled from them for each toy, like for internal
+/// constraints.
+TEST(RooMCStudy, GenParDataSetExternalConstraints)
+{
+   RooRandom::randomGenerator()->SetSeed(4357);
+   RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING};
+
+   RooWorkspace ws;
+   fillModel(ws);
+   ws.factory("Gaussian::constraint(s, cm[7], cs[0.5])");
+   RooArgSet extCons{*ws.pdf("constraint")};
+
+   RooMCStudy mcstudy{*ws.pdf("pdf"), *ws.var("x"), RooFit::ExternalConstraints(extCons), RooFit::Silence()};
+   const int nToys = 20;
+   mcstudy.generate(nToys, 100, true);
+
+   RooDataSet const *genParData = mcstudy.genParDataSet();
+   ASSERT_NE(genParData, nullptr);
+   ASSERT_EQ(genParData->numEntries(), nToys);
+
+   std::vector<double> svals = getColumn(*genParData, "s");
+   double smin = svals[0];
+   double smax = svals[0];
+   double ssum = 0.0;
+   for (double v : svals) {
+      smin = std::min(smin, v);
+      smax = std::max(smax, v);
+      ssum += v;
+   }
+   // The values of "s" are sampled from Gaussian(s | 7, 0.5) for each toy
+   EXPECT_GT(smax, smin);
+   EXPECT_NEAR(ssum / nToys, 7.0, 1.0);
+}
+
+/// The internal-constraints behavior that worked before GitHub issue #9490
+/// must be unchanged: per-toy sampled parameters in genParDataSet, and
+/// "<name>_gen" columns merged into fitParDataSet.
+TEST(RooMCStudy, GenParDataSetInternalConstraints)
+{
+   RooRandom::randomGenerator()->SetSeed(4357);
+   RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING};
+
+   RooWorkspace ws;
+   fillModel(ws);
+   ws.factory("Gaussian::constraint(s, cm[7], cs[0.5])");
+   ws.factory("PROD::prodpdf({pdf,constraint})");
+
+   const int nToys = 3;
+   RooMCStudy mcstudy{*ws.pdf("prodpdf"), *ws.var("x"), RooFit::Constrain(*ws.var("s")), RooFit::Silence(),
+                      RooFit::FitOptions(RooFit::PrintLevel(-1))};
+   mcstudy.generateAndFit(nToys, 200);
+
+   RooDataSet const *genParData = mcstudy.genParDataSet();
+   ASSERT_NE(genParData, nullptr);
+   EXPECT_EQ(genParData->numEntries(), nToys);
+   EXPECT_NE(genParData->get()->find("s"), nullptr);
+   EXPECT_NE(mcstudy.fitParDataSet().get()->find("s_gen"), nullptr);
+}
+
+/// Covers the fitParData/genParData size mismatch from the discussion in
+/// GitHub issue #9490: when some fits fail, the "<name>_gen" columns merged
+/// into fitParDataSet must stay aligned with the successful fits, while
+/// genParDataSet keeps the entries of all generated toys.
+TEST(RooMCStudy, FailedFitsMergeConsistency)
+{
+   RooRandom::randomGenerator()->SetSeed(4357);
+   RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING};
+
+   RooWorkspace ws;
+   ws.factory("x[-10, 10]");
+   ws.factory("Gaussian::gauss(x, m[0, -1, 1], s[7, 5, 10])");
+   // With only 0.7 expected events, some toys have zero events and their fits fail
+   ws.factory("ExtendPdf::model(gauss, nev[0.7])");
+   ws.factory("Gaussian::constraint(s, cm[7], cs[0.5])");
+   RooArgSet extCons{*ws.pdf("constraint")};
+
+   const int nToys = 30;
+   RooMCStudy mcstudy{*ws.pdf("model"),
+                      *ws.var("x"),
+                      RooFit::ExternalConstraints(extCons),
+                      RooFit::Extended(),
+                      RooFit::Silence(),
+                      RooFit::FitOptions(RooFit::PrintLevel(-1))};
+   mcstudy.generateAndFit(nToys, 0);
+
+   RooDataSet const *genParData = mcstudy.genParDataSet();
+   ASSERT_NE(genParData, nullptr);
+   EXPECT_EQ(genParData->numEntries(), nToys);
+
+   RooDataSet const &fitParData = mcstudy.fitParDataSet();
+   const int nFit = fitParData.numEntries();
+   // Make sure the test setup is meaningful: some toys must have failed
+   ASSERT_LT(nFit, nToys);
+   ASSERT_GT(nFit, 0);
+   ASSERT_NE(fitParData.get()->find("s_gen"), nullptr);
+
+   // The merged values must be an ordered subsequence of the generated ones
+   std::vector<double> fitGenVals = getColumn(fitParData, "s_gen");
+   std::vector<double> allGenVals = getColumn(*genParData, "s");
+   std::size_t iAll = 0;
+   for (double val : fitGenVals) {
+      while (iAll < allGenVals.size() && allGenVals[iAll] != val) {
+         ++iAll;
+      }
+      EXPECT_LT(iAll, allGenVals.size()) << "merged s_gen value " << val << " misaligned with generated toys";
+      ++iAll;
+   }
+}