Skip to content
Open
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
39 changes: 39 additions & 0 deletions roofit/roofitcore/src/RooProdPdf.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -2061,6 +2061,41 @@ 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.
// 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<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 @@ -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<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