diff --git a/roofit/roofitcore/src/RooProdPdf.cxx b/roofit/roofitcore/src/RooProdPdf.cxx index 2538a4c7bf88a..ee9095ec71449 100644 --- a/roofit/roofitcore/src/RooProdPdf.cxx +++ b/roofit/roofitcore/src/RooProdPdf.cxx @@ -2061,6 +2061,41 @@ RooProdPdf::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileC std::unique_ptr prodPdfClone{static_cast(this->Clone())}; ctx.markAsCompiled(*prodPdfClone); + // If this RooProdPdf has a normalization range (e.g. the fit range in a + // ranged fit), propagate it to the component pdfs while they are compiled. + // This is required so that a RooAddPdf nested in a RooProdPdf reinterprets + // its coefficients with respect to the full range, exactly like a top-level + // RooAddPdf would. Without this, e.g. a RooProdPdf wrapping an extended + // RooAddPdf (a very common way to attach constraint terms) gives a different + // yield than the bare RooAddPdf in a ranged fit (GitHub issue #16673). + // + // The range must only be set on components for which it is actually defined: + // constraint pdfs of nuisance parameters are normalized over observables that + // don't know about the fit range, and forcing it on them would be wrong (and + // for a multi-range even throws because the undefined sub-ranges collapse to + // the full range and overlap). + std::vector> restoreNormRanges; + if (const char *prodNormRange = normRange()) { + std::vector rangeTokens = ROOT::Split(prodNormRange, ",", /*skipEmpty=*/true); + for (auto *pdf : static_range_cast(prodPdfClone->_pdfList)) { + RooArgSet pdfObs; + pdf->getObservables(&normSet, pdfObs); + bool rangeDefined = !pdfObs.empty(); + for (RooAbsArg *obs : pdfObs) { + auto *lval = dynamic_cast(obs); + for (auto const &token : rangeTokens) { + if (!lval || !lval->hasRange(token.c_str())) { + rangeDefined = false; + } + } + } + if (rangeDefined && std::string(pdf->normRange() ? pdf->normRange() : "") != prodNormRange) { + restoreNormRanges.emplace_back(pdf, pdf->normRange() ? pdf->normRange() : ""); + pdf->setNormRange(prodNormRange); + } + } + } + for (const auto server : prodPdfClone->servers()) { auto nsetForServer = fillNormSetForServer(normSet, *server); RooArgSet const &nset = nsetForServer ? *nsetForServer : normSet; @@ -2071,6 +2106,10 @@ RooProdPdf::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileC ctx.compileServer(*server, *prodPdfClone, depList); } + for (auto const &[pdf, oldRange] : restoreNormRanges) { + pdf->setNormRange(oldRange.empty() ? nullptr : oldRange.c_str()); + } + auto fixedProdPdf = std::make_unique(std::move(prodPdfClone), normSet); ctx.markAsCompiled(*fixedProdPdf); diff --git a/roofit/roofitcore/test/testRooProdPdf.cxx b/roofit/roofitcore/test/testRooProdPdf.cxx index 2b7d0cbd9dfed..c32d97fecf680 100644 --- a/roofit/roofitcore/test/testRooProdPdf.cxx +++ b/roofit/roofitcore/test/testRooProdPdf.cxx @@ -13,6 +13,7 @@ #include #include #include +#include #include @@ -248,3 +249,65 @@ TEST(RooProdPdf, RooProdPdfWithExtendedTerm) double nll1Val = nll1->getVal(); EXPECT_FLOAT_EQ(nll1Val, refVal); } + +// A RooProdPdf that just wraps a (extended) RooAddPdf must give the same result +// in a ranged fit as the bare RooAddPdf. In particular, the RooAddPdf yield must +// be reinterpreted with respect to the full range in both cases. This is the +// common workspace pattern of attaching constraint terms to a model via a +// RooProdPdf. Regression test for GitHub issue #16673. +TEST(RooProdPdf, RangedFitWrappingRooAddPdf) +{ + using namespace RooFit; + RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING); + + RooRandom::randomGenerator()->SetSeed(42); + + RooRealVar x("x", "x", 0, 100); + x.setRange("LEFT", 0, 20); + x.setRange("RIGHT", 60, 100); + + RooRealVar alphaGen("alpha", "alpha", -0.04, -0.1, -0.0); + RooGenericPdf modelGen("model", "exp(alpha*x)", {x, alphaGen}); + std::unique_ptr data{modelGen.generate(x, 10000)}; + + // A nuisance parameter with a Gaussian constraint. Its observable does not + // know about the "LEFT"/"RIGHT" ranges, so the fit range must not be forced + // onto the constraint pdf (otherwise the ranged normalization throws). + RooRealVar theta("theta", "theta", 0, -5, 5); + RooRealVar thetaGlob("thetaGlob", "thetaGlob", 0); + RooRealVar thetaErr("thetaErr", "thetaErr", 1.0); + RooGenericPdf constraint("constraint", "exp(-0.5*(theta-thetaGlob)^2/thetaErr^2)", {theta, thetaGlob, thetaErr}); + + auto fitYield = [&](bool wrapInProd, bool withConstraint) { + RooRealVar alpha("alpha", "alpha", -0.04, -0.1, -0.0); + RooGenericPdf model("model", "exp(alpha*x)", {x, alpha}); + RooRealVar nBkg("Nbkg", "Nbkg", 10000, 0, 20000); + RooAddPdf add("pdfadd", "", RooArgList{model}, RooArgList{nBkg}); + std::unique_ptr owned; + RooAbsPdf *pdf = &add; + if (wrapInProd) { + RooArgList factors{add}; + if (withConstraint) { + factors.add(constraint); + } + owned = std::make_unique("pdfprod", "", factors); + pdf = owned.get(); + } + pdf->fitTo(*data, Range("LEFT,RIGHT"), PrintLevel(-1)); + return nBkg.getVal(); + }; + + const double yieldAdd = fitYield(false, false); + const double yieldProd = fitYield(true, false); + const double yieldProdConstr = fitYield(true, true); + + // The yield must refer to the full range, i.e. be clearly larger than the + // number of events actually inside the fit range. + const double nInFitRange = std::unique_ptr{data->reduce(CutRange("LEFT,RIGHT"))}->sumEntries(); + EXPECT_GT(yieldAdd, 1.2 * nInFitRange); + + // Wrapping in a RooProdPdf (with or without a constraint term) must not + // change the fitted yield. + EXPECT_NEAR(yieldProd, yieldAdd, 1e-3 * yieldAdd); + EXPECT_NEAR(yieldProdConstr, yieldAdd, 1e-3 * yieldAdd); +}