diff --git a/roofit/roofitcore/inc/RooFit/Detail/RooNormalizedPdf.h b/roofit/roofitcore/inc/RooFit/Detail/RooNormalizedPdf.h index f2424b1cb3051..6a18d65faac6a 100644 --- a/roofit/roofitcore/inc/RooFit/Detail/RooNormalizedPdf.h +++ b/roofit/roofitcore/inc/RooFit/Detail/RooNormalizedPdf.h @@ -55,12 +55,20 @@ class RooNormalizedPdf : public RooAbsPdf { Int_t getAnalyticalIntegralWN(RooArgSet &allVars, RooArgSet &analVars, const RooArgSet *normSet, const char *rangeName = nullptr) const override { - return _pdf->getAnalyticalIntegralWN(allVars, analVars, normSet ? normSet : &_normSet, rangeName); + return _pdf->getAnalyticalIntegralWN(allVars, analVars, hasNormSet(normSet) ? normSet : nullptr, rangeName); } /// Forward calculation of analytical integrals to input p.d.f double analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName = nullptr) const override { - return _pdf->analyticalIntegralWN(code, normSet ? normSet : &_normSet, rangeName); + if (hasNormSet(normSet)) { + return _pdf->analyticalIntegralWN(code, normSet, rangeName); + } + // Without normalization set (e.g. in a RooGenProdProj), integrate the + // normalized value as it is. This assumes that the normalization integral + // doesn't depend on the integrated observables. It is taken from this + // wrapper and not from the input pdf, whose cached normalization integral + // goes stale in the RooFit::Evaluator. + return _pdf->analyticalIntegralWN(code, nullptr, rangeName) / _normIntegral; } ExtendMode extendMode() const override { return static_cast(*_pdf).extendMode(); } @@ -87,6 +95,7 @@ class RooNormalizedPdf : public RooAbsPdf { double getValV(const RooArgSet * normSet) const override; private: + static bool hasNormSet(const RooArgSet *normSet) { return normSet && !normSet->empty(); } void logEvalErrorCounts() const; RooTemplateProxy _pdf; diff --git a/roofit/roofitcore/test/testRooProdPdf.cxx b/roofit/roofitcore/test/testRooProdPdf.cxx index 88a33ba9da456..3a25828850467 100644 --- a/roofit/roofitcore/test/testRooProdPdf.cxx +++ b/roofit/roofitcore/test/testRooProdPdf.cxx @@ -19,6 +19,7 @@ #include "gtest_wrapper.h" +#include #include #include #include @@ -177,6 +178,74 @@ TEST(RooProdPdf, DISABLED_ChangeServerNormSetForProdPdfInAddPdf) EXPECT_FLOAT_EQ(val3, val1); } +// A RooGenProdProj normalizes factors with overlapping normalization +// observables together, integrating them with an empty normalization set. The +// cpu backend must agree with the scalar evaluation also if these integrals +// reach normalized pdfs inside composite factors, like a nested RooProdPdf or +// a RooExtendPdf. A jointly normalized product nested in another one is +// compared to the equivalent flat product instead. +TEST(RooProdPdf, JointlyNormalizedCompositeFactors) +{ + 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("Gaussian::gaussz(z[-5, 5], 0.0, 1.0)"); + ws.factory("PROD::nested(gaussy, PROD::inner(gaussx, gaussz))"); + ws.factory("PROD::extended(gaussy, ExtendPdf::gaussxExt(gaussx, n[1000]))"); + ws.factory("Gaussian::gaussw(w[-5, 5], PolyVar::fw(y, {0.0, b[0.3, -1, 1]}), 1.0)"); + ws.factory("PROD::flat(gaussy, gaussx, gaussw)"); + ws.factory("PROD::nestedFirst(PROD::innerFirst(gaussy, gaussx), gaussw)"); + + RooRealVar &x = *ws.var("x"); + RooRealVar &y = *ws.var("y"); + RooRealVar &z = *ws.var("z"); + RooRealVar &w = *ws.var("w"); + RooRealVar &a0 = *ws.var("a0"); + RooRealVar &sigma = *ws.var("sigma"); + + auto check = [&](RooAbsPdf &model, RooArgSet const &obs) { + std::unique_ptr data{model.generate(obs, 1000)}; + std::unique_ptr nll{model.createNLL(*data, RooFit::EvalBackend::Cpu(), RooFit::Extended(false))}; + + auto scalarNll = [&]() { + double out = 0.0; + for (int i = 0; i < data->numEntries(); ++i) { + RooArgSet const *row = data->get(i); + for (RooAbsArg *arg : obs) { + static_cast(arg)->setVal(row->getRealValue(arg->GetName())); + } + 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-6 * std::abs(ref)) << model.GetName() << ", a0 = " << a0Val; + } + }; + + check(*ws.pdf("nested"), {x, y, z}); + check(*ws.pdf("extended"), {x, y}); + + // The nested product only agrees with the flat product to about 1e-3. This + // tolerance still catches integrals of the inner product that are + // inconsistent with its value, which double the NLL. + std::unique_ptr data{ws.pdf("flat")->generate({x, y, w}, 1000)}; + std::unique_ptr nllFlat{ws.pdf("flat")->createNLL(*data, RooFit::EvalBackend::Cpu())}; + std::unique_ptr nllNested{ws.pdf("nestedFirst")->createNLL(*data, RooFit::EvalBackend::Cpu())}; + for (double sigmaVal : {0.5, 1.0}) { + sigma.setVal(sigmaVal); + const double ref = nllFlat->getVal(); + EXPECT_NEAR(nllNested->getVal(), ref, 1e-2 * std::abs(ref)) << "sigma = " << sigmaVal; + } +} + // Make sure that the pdf that provides the expected number of events in a // RooProdPdf is not considered as a constraint. Covers JIRA issue ROOT-7604. TEST(RooProdPdf, RooProdPdfWithExtendedTerm)