From 550a9bd0aa370bc64dde4835bb5c74cf15f51143 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Fri, 6 Nov 2026 17:21:30 +0000 Subject: [PATCH] [RF] Don't wrap jointly normalized RooProdPdf factors for normalization MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Factors of a RooProdPdf with overlapping normalization observables are normalized together by a RooGenProdProj. RooProdPdf::compileForNormSet() now creates the computation graph of the RooFixedProdPdf from a shallow clone of the product where these factors are unwrapped from their RooNormalizedPdf, so the RooGenProdProj integrates the unnormalized pdfs like the scalar evaluation. The RooFixedProdPdf still forwards analytical integrals to the product with the wrapped pdfs: with the unwrapped ones, these integrals would not match its normalized value, which breaks nested products. This fixes two problems with the wrapped pdfs: * The RooGenProdProj clones its pdfs, and the clone of a RooNormalizedPdf owns a copy of its normalization integral, which the RooFit::Evaluator doesn't keep up to date. If the normalization depends on an observable, for example for a pdf conditional on it, the likelihood in the new test was wrong by 36 %. * Integrals of normalized pdfs are smaller, so the absolute tolerance of the numerical integration costs more precision. For the rf305_condcorrprod model of the HS3 test suite, the relative error of 1e-6 exceeded the tolerance of the suite. The new test requires agreement with the scalar evaluation to 1e-10. 🤖 Done with the help of AI --- roofit/roofitcore/inc/RooProdPdf.h | 3 +- roofit/roofitcore/src/RooProdPdf.cxx | 56 +++++++++++++++++++++-- roofit/roofitcore/test/testRooProdPdf.cxx | 48 +++++++++++++++++++ 3 files changed, 102 insertions(+), 5 deletions(-) diff --git a/roofit/roofitcore/inc/RooProdPdf.h b/roofit/roofitcore/inc/RooProdPdf.h index b789dac176e0e..1fb2bb9597ffc 100644 --- a/roofit/roofitcore/inc/RooProdPdf.h +++ b/roofit/roofitcore/inc/RooProdPdf.h @@ -208,7 +208,8 @@ namespace RooFit::Detail { /// normalization set. class RooFixedProdPdf : public RooAbsPdf { public: - RooFixedProdPdf(std::unique_ptr &&prodPdf, RooArgSet const &normSet); + RooFixedProdPdf(std::unique_ptr &&prodPdf, RooArgSet const &normSet, + RooProdPdf const *cacheProdPdf = nullptr); RooFixedProdPdf(const RooFixedProdPdf &other, const char *name = nullptr); inline TObject *clone(const char *newname) const override { return new RooFixedProdPdf(*this, newname); } diff --git a/roofit/roofitcore/src/RooProdPdf.cxx b/roofit/roofitcore/src/RooProdPdf.cxx index a844bb0d70047..1ecc11461e0c1 100644 --- a/roofit/roofitcore/src/RooProdPdf.cxx +++ b/roofit/roofitcore/src/RooProdPdf.cxx @@ -59,6 +59,7 @@ have to appear in any specific place in the list. #include "RooCustomizer.h" #include "RooRealIntegral.h" #include "RooFitImplHelpers.h" +#include "RooFit/Detail/RooNormalizedPdf.h" #include "strtok.h" #include @@ -67,6 +68,9 @@ have to appear in any specific place in the list. #include #include #include +#include +#include +#include #ifndef _WIN32 #include @@ -2073,6 +2077,26 @@ RooProdPdf::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileC } } + // Pdfs with overlapping normalization observables end up in the same term, + // which a RooGenProdProj normalizes as a whole. Like in the scalar + // evaluation, it should integrate the unnormalized pdfs. The RooGenProdProj + // clones its pdfs, and the clone of a RooNormalizedPdf owns a copy of the + // normalization integral that the RooFit::Evaluator doesn't keep up to + // date. Also, the integrals of normalized pdfs are smaller, so the absolute + // tolerance of the numerical integration costs more precision. + std::unordered_set jointlyNormalized; + if (!normSet.empty()) { + Factorized factorized; + prodPdfClone->factorizeProduct(normSet, RooArgSet{}, factorized); + for (auto *term : static_range_cast(factorized.terms)) { + if (term->size() > 1) { + for (RooAbsArg *pdf : *term) { + jointlyNormalized.insert(pdf->GetName()); + } + } + } + } + for (const auto server : prodPdfClone->servers()) { auto nsetForServer = fillNormSetForServer(normSet, *server); RooArgSet const &nset = nsetForServer ? *nsetForServer : normSet; @@ -2087,14 +2111,31 @@ RooProdPdf::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileC pdf->setNormRange(oldRange.empty() ? nullptr : oldRange.c_str()); } + // The computation graph is created from a shallow clone with the jointly + // normalized pdfs unwrapped. The RooFixedProdPdf forwards analytical + // integrals to prodPdfClone, which keeps them wrapped so that these + // integrals match the normalized value, e.g. in a nested RooProdPdf. + std::unordered_map unwrapped; + for (RooAbsArg *pdf : prodPdfClone->_pdfList) { + auto *normalized = dynamic_cast(pdf); + if (normalized && jointlyNormalized.count(normalized->pdf().GetName())) { + unwrapped[normalized] = const_cast(&normalized->pdf()); + } + } + std::unique_ptr cacheProdPdf; + if (!unwrapped.empty()) { + cacheProdPdf.reset(static_cast(prodPdfClone->Clone())); + cacheProdPdf->redirectServers(unwrapped); + } + // The RooFixedProdPdf keeps the normalization set around to create its // cache elements and to evaluate the expected number of events. It must // refer to the observables in the compiled computation graph, because the // original ones can get out of sync with it: RooSimultaneous prefixes the // observable names of each channel after compilation, and the stale // observables would then not be found anymore by name. - auto fixedProdPdf = - std::make_unique(std::move(prodPdfClone), ctx.mapToCompiled(normSet)); + auto fixedProdPdf = std::make_unique( + std::move(prodPdfClone), ctx.mapToCompiled(normSet), cacheProdPdf.get()); ctx.markAsCompiled(*fixedProdPdf); return fixedProdPdf; @@ -2102,13 +2143,20 @@ RooProdPdf::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileC namespace RooFit::Detail { -RooFixedProdPdf::RooFixedProdPdf(std::unique_ptr &&prodPdf, RooArgSet const &normSet) +/// \param[in] prodPdf The product with the fixed normalization set. +/// \param[in] normSet The normalization set. +/// \param[in] cacheProdPdf If not null, the cache elements that make up the +/// computation graph are created from this product instead of +/// `prodPdf`. It must have the same observables and the same pdfs, +/// but they may be normalized differently. +RooFixedProdPdf::RooFixedProdPdf(std::unique_ptr &&prodPdf, RooArgSet const &normSet, + RooProdPdf const *cacheProdPdf) : RooAbsPdf(prodPdf->GetName(), prodPdf->GetTitle()), _normSet{normSet}, _servers("!servers", "List of servers", this), _prodPdf{std::move(prodPdf)} { - auto cache = _prodPdf->createCacheElem(&_normSet, nullptr); + auto cache = (cacheProdPdf ? cacheProdPdf : _prodPdf.get())->createCacheElem(&_normSet, nullptr); _isRearranged = cache->_isRearranged; // We don't want to carry the full cache object around, so we let it go out diff --git a/roofit/roofitcore/test/testRooProdPdf.cxx b/roofit/roofitcore/test/testRooProdPdf.cxx index 88a33ba9da456..7cf39dab81525 100644 --- a/roofit/roofitcore/test/testRooProdPdf.cxx +++ b/roofit/roofitcore/test/testRooProdPdf.cxx @@ -119,6 +119,54 @@ TEST(RooProdPdf, TestDepsAreCond) #endif } +// The RooGenProdProj that normalizes factors with overlapping normalization +// observables must integrate the unnormalized pdfs. Then, the cpu backend +// agrees with the scalar evaluation to high precision, also if the +// normalization of a factor depends on an observable. Covers the HS3 test +// suite model rf305_condcorrprod. +TEST(RooProdPdf, NonFactorizingProductLikelihood) +{ + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + RooRandom::randomGenerator()->SetSeed(1337); + + RooWorkspace ws; + ws.factory("Gaussian::gaussy(y[-5, 5], 0.0, 3.0)"); + ws.factory("Gaussian::gaussx(x[-5, 5], PolyVar::fy(y, {a0[-1.5, -5, 5], a1[-0.5, -1, 1]}), sigma[0.5, 0.1, 3])"); + ws.factory("PROD::joint(gaussy, gaussx)"); + // Here, gaussx is normalized over y only, so its normalization depends on x. + ws.factory("PROD::condOnX(gaussy, gaussx|x)"); + + RooRealVar &x = *ws.var("x"); + RooRealVar &y = *ws.var("y"); + RooRealVar &a0 = *ws.var("a0"); + + for (const char *name : {"joint", "condOnX"}) { + RooAbsPdf &model = *ws.pdf(name); + RooArgSet obs{x, y}; + + std::unique_ptr data{model.generate(obs, 1000)}; + std::unique_ptr nll{model.createNLL(*data, RooFit::EvalBackend::Cpu())}; + + auto scalarNll = [&]() { + double out = 0.0; + for (int i = 0; i < data->numEntries(); ++i) { + RooArgSet const *row = data->get(i); + x.setVal(row->getRealValue("x")); + y.setVal(row->getRealValue("y")); + out -= std::log(model.getVal(obs)); + } + return out; + }; + + for (double a0Val : {-1.5, 0.0, -1.5}) { + a0.setVal(a0Val); + const double ref = scalarNll(); + EXPECT_NEAR(nll->getVal(), ref, 1e-10 * std::abs(ref)) << name << ", a0 = " << a0Val; + } + } +} + /// This test covers a potential problem with the custom normalization ranges /// of PDFs: the value of a RooProdPdf should not depend on the normalization /// ranges of its clients. TODO: this is still a problem! It needs to be fixed,