From a243be6ce216d509e8751631561798e91c92c6e4 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Sun, 26 Jul 2026 20:36:49 +0000 Subject: [PATCH] [RF] Use analytic integral in RooBinSamplingPdf when available MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit RooBinSamplingPdf integrates a continuous pdf over each bin to avoid the bias of evaluating it at the bin centre. So far this was always done with the numeric adaptive integrator, even when the wrapped pdf provides an analytic integral over the observable. Now integrate() detects, lazily and once, whether the wrapped pdf can integrate the observable analytically (getAnalyticalIntegralWN). If so, it uses the exact integral per bin via a private named range, which is both faster and more accurate; otherwise it falls back to the numeric integrator as before. The analytic branch is only taken for RooRealVar observables, since selecting the bin boundaries requires setRange(). The analytic result is normalized over the same normSet via analyticalIntegralWN() to match the numeric path, which integrates the value returned by operator() (i.e. the pdf value normalized over normSet). Addresses one item of #6557. 🤖 Done with the help of AI. --- roofit/roofitcore/inc/RooBinSamplingPdf.h | 4 ++ roofit/roofitcore/src/RooBinSamplingPdf.cxx | 63 ++++++++++++++++++- roofit/roofitcore/test/CMakeLists.txt | 2 +- .../roofitcore/test/testRooBinSamplingPdf.cxx | 52 +++++++++++++++ 4 files changed, 118 insertions(+), 3 deletions(-) diff --git a/roofit/roofitcore/inc/RooBinSamplingPdf.h b/roofit/roofitcore/inc/RooBinSamplingPdf.h index 4a10badb53afa..ff6ef582e72fd 100644 --- a/roofit/roofitcore/inc/RooBinSamplingPdf.h +++ b/roofit/roofitcore/inc/RooBinSamplingPdf.h @@ -24,6 +24,7 @@ #include "Math/Integrator.h" #include +#include class RooBinSamplingPdf : public RooAbsPdf { public: @@ -122,6 +123,7 @@ class RooBinSamplingPdf : public RooAbsPdf { // Call operator for our internal integrator. double operator()(double x) const; double integrate(const RooArgSet* normSet, double low, double high) const; + void initializeAnalyticalIntegral(const RooArgSet* normSet) const; RooTemplateProxy _pdf; @@ -130,6 +132,8 @@ class RooBinSamplingPdf : public RooAbsPdf { mutable std::unique_ptr _integrator{nullptr}; /// _binBoundaries; ///(&*_observable); + if (!observable) { + _analyticalIntegralCode = 0; + return; + } + + _analyticalIntegralRangeName = std::string("_binSampling_") + GetName(); + + // Define the range so that PDFs that inspect it in getAnalyticalIntegral() + // find a valid one. The actual bin boundaries are filled in for each bin in + // integrate(). + observable->setRange(_analyticalIntegralRangeName.c_str(), observable->getMin(), observable->getMax()); + + RooArgSet allVars{*_observable}; + RooArgSet analVars; + _analyticalIntegralCode = _pdf->getAnalyticalIntegralWN(allVars, analVars, normSet, _analyticalIntegralRangeName.c_str()); + + // Only use analytical integration if the observable is really integrated + // analytically. Otherwise, fall back to the numeric integrator. + if (_analyticalIntegralCode != 0 && !analVars.contains(*_observable)) { + _analyticalIntegralCode = 0; + } + + if (_analyticalIntegralCode != 0) { + coutI(NumIntegration) << "RooBinSamplingPdf::integrate(" << GetName() + << "): using the analytical integral of " << _pdf->GetName() + << " to sample the bins instead of the numeric integrator." << std::endl; + } +} + + +//////////////////////////////////////////////////////////////////////////////// +/// Integrate the wrapped PDF over a single bin, with the given norm set and limits. +/// If the wrapped PDF supports analytical integration over the observable, the +/// exact integral is used. Otherwise, the numeric integrator is employed. +/// +/// The result must match the numeric path, which integrates the value returned by +/// operator(), i.e. the PDF value normalized over `normSet`. Therefore, the +/// analytical integral is normalized over `normSet` as well via analyticalIntegralWN(). +double RooBinSamplingPdf::integrate(const RooArgSet* normSet, double low, double high) const { + if (_analyticalIntegralCode == -1) { + initializeAnalyticalIntegral(normSet); + } + + if (_analyticalIntegralCode != 0) { + // The analytical path is only enabled for RooRealVar observables (see + // initializeAnalyticalIntegral()), so this static_cast is safe. + static_cast(*_observable).setRange(_analyticalIntegralRangeName.c_str(), low, high); + return _pdf->analyticalIntegralWN(_analyticalIntegralCode, normSet, _analyticalIntegralRangeName.c_str()); + } + return integrator()->Integral(low, high); } diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt index def1fa324c837..e57a75936781f 100644 --- a/roofit/roofitcore/test/CMakeLists.txt +++ b/roofit/roofitcore/test/CMakeLists.txt @@ -25,7 +25,6 @@ ROOT_ADD_GTEST(testRooFitCore testGenericPdf.cxx testGlobalObservables.cxx testRooAbsCollection.cxx - testRooBinSamplingPdf.cxx testRooCacheManager.cxx testRooCategory.cxx testRooCollectionProxy.cxx @@ -58,6 +57,7 @@ ROOT_ADD_GTEST(testRooFitCore ${CMAKE_CURRENT_SOURCE_DIR}/data/rf502_workspace_v6.14.root ) +ROOT_ADD_GTEST(testRooBinSamplingPdf testRooBinSamplingPdf.cxx LIBRARIES RooFitCore RooFit) ROOT_ADD_GTEST(testRooAddPdf testRooAddPdf.cxx LIBRARIES RooFitCore RooFit RooStats) ROOT_ADD_GTEST(testRooAbsPdf testRooAbsPdf.cxx LIBRARIES RooFitCore RooFit) diff --git a/roofit/roofitcore/test/testRooBinSamplingPdf.cxx b/roofit/roofitcore/test/testRooBinSamplingPdf.cxx index 479ee82b8e4a8..9a297225972ac 100644 --- a/roofit/roofitcore/test/testRooBinSamplingPdf.cxx +++ b/roofit/roofitcore/test/testRooBinSamplingPdf.cxx @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include @@ -113,6 +114,57 @@ TEST(RooBinSamplingPdf, CheckConsistentNormalization) EXPECT_FLOAT_EQ(int2->getVal(), int3->getVal()); } +// The RooBinSamplingPdf uses the analytical integral of the wrapped pdf to +// sample the bins when it is available, and falls back to the numeric +// integrator otherwise. This test checks that the analytical path (here +// exercised with a RooGaussian) gives the same result as the numeric path +// (exercised with an identical RooGenericPdf that has no analytical integral). +TEST(RooBinSamplingPdf, AnalyticalMatchesNumeric) +{ + using namespace RooFit; + + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + RooRealVar x("x", "x", -5, 5); + x.setBins(20); + RooRealVar mean("mean", "mean", 0.7, -5, 5); + RooRealVar sigma("sigma", "sigma", 1.3, 0.1, 5); + + // Has an analytical integral over x -> analytical bin sampling. + RooGaussian gaus("gaus", "gaus", x, mean, sigma); + // Identical shape, but no analytical integral -> numeric bin sampling. + RooGenericPdf gen("gen", "gen", "std::exp(-0.5*(x-mean)*(x-mean)/(sigma*sigma))", {x, mean, sigma}); + + RooBinSamplingPdf bsAna("bsAna", "bsAna", x, gaus); + RooBinSamplingPdf bsNum("bsNum", "bsNum", x, gen); + + RooArgSet normSet{x}; + + // Compare the sampled values bin by bin. + for (int i = 0; i < x.numBins(); ++i) { + x.setBin(i); + EXPECT_NEAR(bsAna.getVal(normSet), bsNum.getVal(normSet), 1e-5 * bsNum.getVal(normSet)) + << "mismatch in bin " << i; + } + + // The results must also agree when used in a fit, for both the legacy and + // the vectorizing "cpu" evaluation backend. + std::unique_ptr dataH(gaus.generateBinned(x, 20000)); + RooDataSet data("data", "data", x, RooFit::Import(*dataH)); + + for (auto backend : {EvalBackend::Legacy(), EvalBackend::Cpu()}) { + mean.setVal(0.7); + sigma.setVal(1.3); + std::unique_ptr nllAna{gaus.createNLL(data, IntegrateBins(1.E-3), backend)}; + mean.setVal(0.7); + sigma.setVal(1.3); + std::unique_ptr nllNum{gen.createNLL(data, IntegrateBins(1.E-3), backend)}; + + EXPECT_NEAR(nllAna->getVal(), nllNum->getVal(), 1e-5 * std::abs(nllNum->getVal())) + << "NLL mismatch for backend " << backend.name(); + } +} + INSTANTIATE_TEST_SUITE_P(RooBinSamplingPdf, ParamTest, testing::Values("Off", "Cpu"), [](testing::TestParamInfo const ¶mInfo) { std::stringstream ss;