Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 41 additions & 0 deletions roofit/roofitcore/src/RooProdPdf.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -2066,6 +2066,43 @@ RooProdPdf::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileC
std::unique_ptr<RooProdPdf> prodPdfClone{static_cast<RooProdPdf *>(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,
// so that a nested RooAddPdf reinterprets its coefficients with respect to
// the full range exactly like a top-level RooAddPdf would. Otherwise, a
// RooProdPdf wrapping an extended RooAddPdf (the common way to attach
// constraint terms) gives a different yield than the bare RooAddPdf in a
// ranged fit (GitHub issue #16673).
//
// This must not change how constraint terms are normalized, so the range is
// only set on components whose normalization observables define it.
// Constraint pdfs are normalized over the nuisance parameters, which don't
// define the fit range, so they are skipped. Forcing the range on them would
// be wrong, and for a multi-range fit it even throws, because each undefined
// sub-range falls back to the parameter's full range and the sub-ranges then
// overlap.
std::vector<std::pair<RooAbsPdf *, std::string>> restoreNormRanges;
if (const char *prodNormRange = normRange()) {
std::vector<std::string> rangeTokens = ROOT::Split(prodNormRange, ",", /*skipEmpty=*/true);
for (auto *pdf : static_range_cast<RooAbsPdf *>(prodPdfClone->_pdfList)) {
RooArgSet pdfObs;
pdf->getObservables(&normSet, pdfObs);
bool rangeDefined = !pdfObs.empty();
for (RooAbsArg *obs : pdfObs) {
auto *lval = dynamic_cast<RooAbsRealLValue *>(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;
Expand All @@ -2076,6 +2113,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<RooFit::Detail::RooFixedProdPdf>(std::move(prodPdfClone), normSet);
ctx.markAsCompiled(*fixedProdPdf);

Expand Down
63 changes: 63 additions & 0 deletions roofit/roofitcore/test/testRooProdPdf.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
#include <RooRealVar.h>
#include <RooWorkspace.h>
#include <RooHelpers.h>
#include <RooRandom.h>

#include <Math/PdfFuncMathCore.h>

Expand Down Expand Up @@ -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<RooDataSet> 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<RooAbsPdf> owned;
RooAbsPdf *pdf = &add;
if (wrapInProd) {
RooArgList factors{add};
if (withConstraint) {
factors.add(constraint);
}
owned = std::make_unique<RooProdPdf>("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<RooAbsData>{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);
}
Loading