[RF] Always save generated parameters in RooMCStudy

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;
+   }
+}