diff --git a/roofit/roofitcore/inc/RooHistPdf.h b/roofit/roofitcore/inc/RooHistPdf.h index 544bc46704396..7c1d0f560a54c 100644 --- a/roofit/roofitcore/inc/RooHistPdf.h +++ b/roofit/roofitcore/inc/RooHistPdf.h @@ -145,9 +145,16 @@ class RooHistPdf : public RooAbsPdf { inline void initializeOwnedDataHist(std::unique_ptr &&dataHist) { - _ownedDataHist = std::move(dataHist); + // The constructor may already have taken ownership of a sanitized clone + // of the input histogram (see clampNegativeBins()). In that case, the + // clone stays and the original histogram can be disposed of. + if (!_ownedDataHist) { + _ownedDataHist = std::move(dataHist); + } } + void clampNegativeBins(); + ClassDefOverride(RooHistPdf,4) // Histogram based PDF }; diff --git a/roofit/roofitcore/src/RooHistPdf.cxx b/roofit/roofitcore/src/RooHistPdf.cxx index 7ac8eff0ea9ae..c2aee28ba16e2 100644 --- a/roofit/roofitcore/src/RooHistPdf.cxx +++ b/roofit/roofitcore/src/RooHistPdf.cxx @@ -23,6 +23,10 @@ A probability density function sampled from a multidimensional histogram. The histogram distribution is explicitly normalized by RooHistPdf and can have an arbitrary number of real or discrete dimensions. + +A p.d.f. cannot be negative. If the input histogram contains bins with +negative content, the bin contents are clipped to zero and the bin errors are kept the same. +The input histogram is not modified. **/ #include "Riostream.h" @@ -41,14 +45,16 @@ discrete dimensions. #include "TError.h" #include "TBuffer.h" - - +#include +#include //////////////////////////////////////////////////////////////////////////////// /// Constructor from a RooDataHist. RooDataHist dimensions /// can be either real or discrete. See RooDataHist::RooDataHist for details on the binning. /// RooHistPdf neither owns or clone 'dhist' and the user must ensure the input histogram exists /// for the entire life span of this PDF. +/// The only exception is a 'dhist' that contains bins with negative content: those are set to +/// zero in an internally-owned clone that is used instead (see clampNegativeBins()). RooHistPdf::RooHistPdf(const char *name, const char *title, const RooArgSet& vars, const RooDataHist& dhist, Int_t intOrder) : @@ -86,6 +92,7 @@ RooHistPdf::RooHistPdf(const char *name, const char *title, const RooArgSet& var } } + clampNegativeBins(); } @@ -139,6 +146,8 @@ RooHistPdf::RooHistPdf(const char *name, const char *title, const RooArgList& pd (static_cast(hobs))->setRange(dhreal->getMin(),dhreal->getMax()) ; } } + + clampNegativeBins(); } RooHistPdf::RooHistPdf(const char *name, const char *title, const RooArgSet &vars, std::unique_ptr dhist, @@ -177,6 +186,51 @@ RooDataHist* RooHistPdf::cloneAndOwnDataHist(const char* newname) { return _dataHist; } +void RooHistPdf::clampNegativeBins() +{ + const std::size_t nBins = _dataHist->numEntries(); + + std::size_t nNegative = 0; + double sumNegative = 0.; + for (std::size_t i = 0; i < nBins; ++i) { + if (_dataHist->weight(i) < 0.) { + ++nNegative; + sumNegative += _dataHist->weight(i); + } + } + if (nNegative == 0) { + return; + } + + coutW(InputArguments) << "RooHistPdf::ctor(" << GetName() << ") WARNING: input histogram \"" << _dataHist->GetName() + << "\" contains " << nNegative + << " bins with negative content (sum of negative contents: " << sumNegative + << "). A p.d.f. cannot be negative, so these bins contents are clipped to zero while " + "preserving the error. The input " + "histogram is not modified. To avoid this message, remove the negative bin contents " + "before constructing the RooHistPdf." + << std::endl; + + RooDataHist *dh = cloneAndOwnDataHist(); + const bool hasSumW2 = dh->sumW2Array() != nullptr; + for (std::size_t i = 0; i < nBins; ++i) { + if (dh->weight(i) < 0.) { + // Keep the original bin error: clamping the content is a + // normalization-consistency measure, not a statement that the bin is + // now known exactly. The error still quantifies the statistical + // uncertainty of the original bin content estimate (e.g. whether the + // negative content is compatible with a fluctuation around zero), + // and setting it to zero would introduce undercoverage, which is + // always undesired. It would also irreversibly discard information + // for anyone retrieving the histogram via dataHist(), including any + // future per-bin MC-stat treatment, where a zero error would wrongly + // fix the bin at exactly zero. + const double wgtErr = hasSumW2 ? std::sqrt(std::max(dh->weightSquared(i), 0.)) : 0.; + dh->set(i, 0., wgtErr); + } + } +} + void RooHistPdf::doEval(RooFit::EvalContext &ctx) const { std::span output = ctx.output(); diff --git a/roofit/roofitcore/test/testRooHistPdf.cxx b/roofit/roofitcore/test/testRooHistPdf.cxx index ae4c32edea56a..3dc5c870e4db9 100644 --- a/roofit/roofitcore/test/testRooHistPdf.cxx +++ b/roofit/roofitcore/test/testRooHistPdf.cxx @@ -12,6 +12,7 @@ #include +#include #include #include @@ -39,6 +40,54 @@ TEST(RooHistPdf, AnalyticIntWithRooLinearVar) EXPECT_EQ(integ.anaIntVars().size(), 1); } +// A RooHistPdf clamps negative bin contents to zero when evaluating, but the +// normalization used to be based on the raw sum of weights including the +// negative ones, so the pdf did not integrate to unity and extended fits were +// silently biased (ROOT-10825). Now, negative bins are zeroed in an +// internally-owned clone of the histogram at construction time, making the +// shape and the normalization consistent. +TEST(RooHistPdf, NegativeBinsAreClampedConsistently) +{ + RooRealVar x{"x", "x", 0, 4}; + x.setBins(4); + + RooDataHist dataHist{"dataHist", "dataHist", x}; + double const contents[4] = {4., 3., -2., 5.}; + double const sumClamped = 12.; + for (int i = 0; i < x.numBins(); ++i) { + dataHist.set(i, contents[i], 1.0); + } + + RooHelpers::HijackMessageStream hijack{RooFit::WARNING, RooFit::InputArguments, "pdf"}; + RooHistPdf pdf{"pdf", "pdf", x, dataHist}; + EXPECT_FALSE(hijack.str().empty()) << "constructing from a histogram with negative bins should warn"; + + // The negative bin is dropped from both the pdf value and the + // normalization, so the pdf values are max(w, 0) / sum(max(w, 0)) and the + // pdf integrates to unity. + RooArgSet normSet{x}; + double integral = 0.; + for (int i = 0; i < x.numBins(); ++i) { + x.setBin(i); + EXPECT_DOUBLE_EQ(pdf.getVal(normSet), std::max(contents[i], 0.) / sumClamped) << "wrong value in bin " << i; + integral += pdf.getVal(normSet) * x.getBinWidth(i); + } + EXPECT_DOUBLE_EQ(integral, 1.); + + // The input histogram must not be modified. + for (int i = 0; i < x.numBins(); ++i) { + EXPECT_DOUBLE_EQ(dataHist.weight(i), contents[i]) << "input histogram was modified in bin " << i; + } + + // Same behavior for the constructor where the RooHistPdf takes ownership + // of the input histogram. + RooHistPdf pdf2{"pdf2", "pdf2", x, std::make_unique(dataHist, "dataHist2")}; + for (int i = 0; i < x.numBins(); ++i) { + x.setBin(i); + EXPECT_DOUBLE_EQ(pdf2.getVal(normSet), std::max(contents[i], 0.) / sumClamped) << "wrong value in bin " << i; + } +} + namespace { // Narrow Gaussian peak (sigma smaller than the bin width) on a steeply falling