[RF] Don't randomize number of events when creating Asimov from RooSimultaneous

roofit/roofitcore/inc/RooSimSplitGenContext.h

@@ -54,7 +54,7 @@ class RooSimSplitGenContext : public RooAbsGenContext {
   std::vector<int>               _gcIndex ; ///< Index value corresponding to component
   TString _idxCatName ;                     ///< Name of index category
   Int_t _numPdf ;                           ///< Number of generated PDFs
-  double* _fracThresh ;                   ///< fraction thresholds
+  bool _expectedData = false;               ///< Asimov?
 
   RooArgSet _allVarsPdf ; ///< All pdf variables
 

roofit/roofitcore/src/RooSimSplitGenContext.cxx

@@ -78,12 +78,9 @@ RooSimSplitGenContext::RooSimSplitGenContext(const RooSimultaneous &model, const
 
   // Initialize fraction threshold array (used only in extended mode)
   _numPdf = model._pdfProxyList.GetSize() ;
-  _fracThresh = new double[_numPdf+1] ;
-  _fracThresh[0] = 0 ;
 
   // Generate index category and all registered PDFS
   _allVarsPdf.add(allPdfVars) ;
-  Int_t i(1) ;
   for(auto * proxy : static_range_cast<RooRealProxy*>(model._pdfProxyList)) {
     auto pdf = static_cast<RooAbsPdf*>(proxy->absArg());
 
@@ -96,14 +93,6 @@ RooSimSplitGenContext::RooSimSplitGenContext(const RooSimultaneous &model, const
     cx->SetName(proxy->name()) ;
     _gcList.push_back(cx) ;
     _gcIndex.push_back(state);
-
-    // Fill fraction threshold array
-    _fracThresh[i] = _fracThresh[i-1] + pdf->expectedEvents(&allPdfVars) ;
-    i++ ;
-  }
-
-  for(i=0 ; i<_numPdf ; i++) {
-    _fracThresh[i] /= _fracThresh[_numPdf] ;
   }
 
   // Clone the index category
@@ -121,7 +110,6 @@ RooSimSplitGenContext::RooSimSplitGenContext(const RooSimultaneous &model, const
 
 RooSimSplitGenContext::~RooSimSplitGenContext()
 {
-  delete[] _fracThresh ;
   for (RooAbsGenContext *item : _gcList) {
     delete item;
   }
@@ -193,39 +181,35 @@ RooDataSet* RooSimSplitGenContext::generate(double nEvents, bool skipInit, bool
   }
 
   // Generate lookup table from expected event counts
-  std::vector<double> nGen(_numPdf) ;
-  if (extendedMode ) {
-    Int_t i(0) ;
-    for(auto * proxy : static_range_cast<RooRealProxy*>(_pdf->_pdfProxyList)) {
-      RooAbsPdf* pdf=static_cast<RooAbsPdf*>(proxy->absArg()) ;
-      //nGen[i] = Int_t(pdf->expectedEvents(&_allVarsPdf)+0.5) ;
-      nGen[i] = pdf->expectedEvents(&_allVarsPdf) ;
-      i++ ;
-    }
+  std::vector<double> nGen(_numPdf);
+  std::vector<double> nExpected;
+  nExpected.reserve(_numPdf) ;
+  double nExpectedTotal = 0.;
+  for(auto * proxy : static_range_cast<RooRealProxy*>(_pdf->_pdfProxyList)) {
+    RooAbsPdf* pdf=static_cast<RooAbsPdf*>(proxy->absArg()) ;
+    nExpected.push_back(pdf->expectedEvents(&_allVarsPdf));
+    nExpectedTotal += nExpected.back();
+  }
 
+  // We don't randomize events in two cases:
+  //  1. When the generation is extended, each component pdf will already
+  //     randomize the expected number of events so we don't need to do it here.
+  //  2. If we want to create an expected Asimov dataset.
+  if (extendedMode || _expectedData) {
+    nGen = nExpected;
   } else {
-    Int_t i(1) ;
-    _fracThresh[0] = 0 ;
-    for(auto * proxy : static_range_cast<RooRealProxy*>(_pdf->_pdfProxyList)) {
-      RooAbsPdf* pdf=static_cast<RooAbsPdf*>(proxy->absArg()) ;
-      _fracThresh[i] = _fracThresh[i-1] + pdf->expectedEvents(&_allVarsPdf) ;
-      i++ ;
-    }
-    for(i=0 ; i<_numPdf ; i++) {
-      _fracThresh[i] /= _fracThresh[_numPdf] ;
-    }
-
     // Determine from that total number of events to be generated for each component
     double nGenSoFar(0) ;
     while (nGenSoFar<nEvents) {
-      double rand = RooRandom::uniform() ;
-      i=0 ;
-      for (i=0 ; i<_numPdf ; i++) {
-   if (rand>_fracThresh[i] && rand<_fracThresh[i+1]) {
-     nGen[i]++ ;
-     nGenSoFar++ ;
-     break ;
-   }
+      double rand = RooRandom::uniform() * nExpectedTotal;
+      double cumsum = 0;
+      for (int i=0 ; i<_numPdf ; i++) {
+        if (rand >= cumsum && rand < cumsum + nExpected[i]) {
+          nGen[i]++ ;
+          nGenSoFar++ ;
+          break ;
+        }
+        cumsum += nExpected[i];
       }
     }
   }
@@ -257,6 +241,7 @@ RooDataSet* RooSimSplitGenContext::generate(double nEvents, bool skipInit, bool
 
 void RooSimSplitGenContext::setExpectedData(bool flag)
 {
+  _expectedData = flag;
   for(RooAbsGenContext *elem : _gcList) {
     elem->setExpectedData(flag) ;
   }

roofit/roofitcore/test/testRooSimultaneous.cxx

@@ -1,6 +1,7 @@
 // Tests for the RooSimultaneous
 // Authors: Jonas Rembser, CERN  06/2021
 
+#include <Roo1DTable.h>
 #include <RooAddition.h>
 #include <RooCategory.h>
 #include <RooConstVar.h>
@@ -541,3 +542,35 @@ TEST(RooSimultaneous, PlotProjWData)
    // sum of data entries in the center.
    EXPECT_DOUBLE_EQ(frame->getCurve()->interpolate(0.), combData.sumEntries());
 }
+
+/// 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)
+{
+
+   RooWorkspace ws{"ws"};
+   ws.factory("dummy_obs_a[0,1]");
+   ws.factory("dummy_obs_b[0,1]");
+   ws.factory("Uniform::uniform_a(dummy_obs_a)");
+   ws.factory("Uniform::uniform_b(dummy_obs_b)");
+   ws.factory("SUM::model_a(coeff_a[3.5]*uniform_a)");
+   ws.factory("SUM::model_b(coeff_b[6.5]*uniform_b)");
+
+   RooRealVar &dummy_obs_a = *ws.var("dummy_obs_a");
+   RooRealVar &dummy_obs_b = *ws.var("dummy_obs_b");
+
+   ws.factory("dummy_cat[a]");
+   ws.factory("SIMUL::sim_model(dummy_cat, a = model_a, b = model_b)");
+   RooAbsCategory &dummy_cat = *ws.cat("dummy_cat");
+
+   // std::cout << "simultaneous expected = " << ws.pdf("sim_model")->expectedEvents(dummy_obs) << std::endl;
+   RooDataSet *data = ws.pdf("sim_model")->generate({dummy_obs_a, dummy_obs_b, dummy_cat}, RooFit::ExpectedData());
+
+   std::unique_ptr<Roo1DTable> tab{data->table(dummy_cat)};
+
+   // Check that the sum of entries for each category is as expected, matching
+   // the coefficients from the RooAddPdf.
+   EXPECT_FLOAT_EQ(tab->get("a"), ws.var("coeff_a")->getVal());
+   EXPECT_FLOAT_EQ(tab->get("b"), ws.var("coeff_b")->getVal());
+}