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
9 changes: 8 additions & 1 deletion roofit/roofitcore/inc/RooHistPdf.h
Original file line number Diff line number Diff line change
Expand Up @@ -145,9 +145,16 @@ class RooHistPdf : public RooAbsPdf {

inline void initializeOwnedDataHist(std::unique_ptr<RooDataHist> &&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
};

Expand Down
58 changes: 56 additions & 2 deletions roofit/roofitcore/src/RooHistPdf.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -41,14 +45,16 @@ discrete dimensions.
#include "TError.h"
#include "TBuffer.h"



#include <algorithm>
#include <cmath>

////////////////////////////////////////////////////////////////////////////////
/// 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) :
Expand Down Expand Up @@ -86,6 +92,7 @@ RooHistPdf::RooHistPdf(const char *name, const char *title, const RooArgSet& var
}
}

clampNegativeBins();
}


Expand Down Expand Up @@ -139,6 +146,8 @@ RooHistPdf::RooHistPdf(const char *name, const char *title, const RooArgList& pd
(static_cast<RooRealVar*>(hobs))->setRange(dhreal->getMin(),dhreal->getMax()) ;
}
}

clampNegativeBins();
}

RooHistPdf::RooHistPdf(const char *name, const char *title, const RooArgSet &vars, std::unique_ptr<RooDataHist> dhist,
Expand Down Expand Up @@ -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<double> output = ctx.output();
Expand Down
49 changes: 49 additions & 0 deletions roofit/roofitcore/test/testRooHistPdf.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@

#include <gtest/gtest.h>

#include <algorithm>
#include <cmath>
#include <memory>

Expand Down Expand Up @@ -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<RooDataHist>(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
Expand Down
Loading