diff --git a/roofit/roofitcore/inc/RooSimSplitGenContext.h b/roofit/roofitcore/inc/RooSimSplitGenContext.h index 6df6bd9f00b51..5b48f399087ef 100644 --- a/roofit/roofitcore/inc/RooSimSplitGenContext.h +++ b/roofit/roofitcore/inc/RooSimSplitGenContext.h @@ -54,7 +54,7 @@ class RooSimSplitGenContext : public RooAbsGenContext { std::vector _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 diff --git a/roofit/roofitcore/src/RooSimSplitGenContext.cxx b/roofit/roofitcore/src/RooSimSplitGenContext.cxx index a416658a3d1e7..69d6c221bdadd 100644 --- a/roofit/roofitcore/src/RooSimSplitGenContext.cxx +++ b/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(model._pdfProxyList)) { auto pdf = static_cast(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 nGen(_numPdf) ; - if (extendedMode ) { - Int_t i(0) ; - for(auto * proxy : static_range_cast(_pdf->_pdfProxyList)) { - RooAbsPdf* pdf=static_cast(proxy->absArg()) ; - //nGen[i] = Int_t(pdf->expectedEvents(&_allVarsPdf)+0.5) ; - nGen[i] = pdf->expectedEvents(&_allVarsPdf) ; - i++ ; - } + std::vector nGen(_numPdf); + std::vector nExpected; + nExpected.reserve(_numPdf) ; + double nExpectedTotal = 0.; + for(auto * proxy : static_range_cast(_pdf->_pdfProxyList)) { + RooAbsPdf* pdf=static_cast(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(_pdf->_pdfProxyList)) { - RooAbsPdf* pdf=static_cast(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_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) ; } diff --git a/roofit/roofitcore/test/testRooSimultaneous.cxx b/roofit/roofitcore/test/testRooSimultaneous.cxx index ff58eb7275f90..944cd8e7b515b 100644 --- a/roofit/roofitcore/test/testRooSimultaneous.cxx +++ b/roofit/roofitcore/test/testRooSimultaneous.cxx @@ -1,6 +1,7 @@ // Tests for the RooSimultaneous // Authors: Jonas Rembser, CERN 06/2021 +#include #include #include #include @@ -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 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()); +}