From 0c0dfe67244f088a6f72a8687573cacf3e291921 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Mon, 14 Sep 2026 08:00:43 +0000 Subject: [PATCH] [RF] Add RooFit::makeResidHist() and makePullHist() with bin integration MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Residual and pull histograms of binned data are currently created from a plotted RooCurve (RooPlot::residHist()/pullHist()), which interpolates or averages a sampled polyline of the pdf. For strongly peaked pdfs, this produces a systematic 'wiggle' pattern in residual panels and unrealistically large chi-square values when binned data is presented against the result of an unbinned fit, even when the fit is perfectly fine. See this ROOT forum thread for a study of the problem: https://root-forum.cern.ch/t/comparison-of-binned-data-to-results-of-unbinned-fit/64567 The new free functions take the fit model (any RooAbsReal) and the data (RooAbsData) directly and integrate the model exactly over each bin with RooAbsReal::createIntegral(), avoiding the curve-based bias. The expectation is normalized to the data weight, also correctly restricted to the range if Range() is passed. Unbinned input data is binned internally, and the usual Binning(), Normalization(), and DataError() command arguments are supported. For projections or components of composite models, project the model and data first or pass the component pdf with Normalization(). A warning is printed once per process when the curve-based RooPlot::residHist()/pullHist() are used, pointing to the new functions. The RooPlot functions are scheduled for deprecation in ROOT v6.44. The rf109 tutorial and the RooHist unit tests are updated/extended to the new interface. 🤖 Done with the help of AI --- README/ReleaseNotes/v642/index.md | 13 + roofit/roofitcore/inc/RooGlobalFunc.h | 83 +++++ roofit/roofitcore/inc/RooHist.h | 5 + roofit/roofitcore/inc/RooPlot.h | 8 +- roofit/roofitcore/src/RooGlobalFunc.cxx | 302 ++++++++++++++++++ roofit/roofitcore/src/RooPlot.cxx | 16 + roofit/roofitcore/test/testRooHist.cxx | 257 +++++++++++++++ tutorials/roofit/roofit/rf109_chi2residpull.C | 20 +- .../roofit/roofit/rf109_chi2residpull.py | 18 +- 9 files changed, 709 insertions(+), 13 deletions(-) diff --git a/README/ReleaseNotes/v642/index.md b/README/ReleaseNotes/v642/index.md index ccea5f564c7b0..81d6130889070 100644 --- a/README/ReleaseNotes/v642/index.md +++ b/README/ReleaseNotes/v642/index.md @@ -162,6 +162,19 @@ the cut instead of being selected based on `sqrt(abs(x))`. ## RooFit +### Exact bin-integrated residual and pull histograms + +The new `RooFit::makeResidHist()` and `RooFit::makePullHist()` functions create residual and pull histograms of binned data with respect to a fitted model, where the model expectation per bin is computed by integrating the model exactly over the bin. +The existing `RooPlot::residHist()` and `RooPlot::pullHist()` instead interpolate or average the curve of the plotted pdf, which biases the comparison for strongly peaked pdfs: it produces a characteristic "wiggle" pattern in the residuals and unrealistically large chi-square values when the binned data is compared with the result of an unbinned fit. +The new functions take the model (any `RooAbsReal`) and the data (`RooAbsData`, which is binned internally if unbinned) and support the usual RooFit command arguments, such as `RooFit::Binning()` to select the binning (also for unbinned data), `RooFit::Range()` to restrict the comparison to (a union of) fit ranges and normalize the model accordingly, `RooFit::Normalization()` for an extra scale factor, and `RooFit::DataError()` to select the point error model. + +```cpp +RooDataHist binData("binData", "binData", x, *unbinnedData); +auto pull = RooFit::makePullHist(pdf, binData); +``` + +For comparisons against projected pdfs or components, project the model with `RooAbsPdf::createProjection()` and the data with `RooAbsData::reduce()`, or pass the component pdf with a matching `RooFit::Normalization()` argument. + ### RooFit::MultiProcess without ZeroMQ, now enabled by default The `RooFit::MultiProcess` package that implements the parallel gradient minimization with `fitTo(..., RooFit::Parallelize(n))` previously communicated between the forked processes with ZeroMQ sockets, which required building ROOT with `roofit_multiprocess=ON` and the ZeroMQ (with draft API) and cppzmq dependencies. diff --git a/roofit/roofitcore/inc/RooGlobalFunc.h b/roofit/roofitcore/inc/RooGlobalFunc.h index eb1f27b7915fc..bf131a096d55b 100644 --- a/roofit/roofitcore/inc/RooGlobalFunc.h +++ b/roofit/roofitcore/inc/RooGlobalFunc.h @@ -24,6 +24,7 @@ #include #include +#include #include class RooDataHist ; @@ -44,6 +45,7 @@ class RooConstVar ; class RooRealVar ; class RooAbsCategory ; class RooNumIntConfig ; +class RooHist; class TH1 ; class TTree ; @@ -472,6 +474,87 @@ RooConstVar& RooConst(double val) ; * @} */ +/** + * \defgroup Residuals Residuals and pulls of binned data vs. a fitted model + * + * RooFit::makeResidHist() and RooFit::makePullHist() construct per-bin + * residual or pull distributions of binned data with respect to a fitted + * model. In contrast to RooPlot::residHist() and RooPlot::pullHist(), + * which interpolate or average the curve of a plotted pdf, these functions + * integrate the model itself exactly over each bin. This avoids the + * biased residuals (the "residual wiggle" and inflated \f$\chi^2\f$ values) + * that appear when a sharply peaked pdf is compared with binned data. + * @{ + */ + +/// Creates a RooHist with the per-bin residuals `data` - `fitModel`. +/// +/// For each bin, the model expectation is the model integrated exactly over +/// the bin (with RooAbsReal::createIntegral()), normalized to the weight of +/// the data inside the normalization range. Unbinned input data is binned +/// internally. The supported RooCmdArgs are: +/// - `Binning(binning)` / `Binning(nbins)` / `Binning(nbins, xlo, xhi)` : +/// the binning used to bin unbinned input data (default: the current +/// binning of the observable) +/// - `Range(lo, hi)` / `Range("name")` : the range that the model +/// expectation is normalized in. Only bins inside this range get points +/// in the result. Range("name") can be passed multiple times for a union +/// of named ranges (sidebands). Use it for sideband fits or blinded data, +/// so that the model is normalized in the same range as the fitted data. +/// - `Normalization(double)` : additional scale factor for the expectation +/// (default 1.0). The scale types of Normalization(scale, type) are not +/// supported. +/// - `DataError(RooAbsData::ErrorType)` : error model used for the point errors, +/// one of Poisson, SumW2 or Auto +/// (default: Auto, meaning SumW2 for weighted and Poisson for unweighted data) +/// - `Name(const char*)`, `Title(const char*)` : name and title of the created RooHist +/// +/// \note The residuals are computed with respect to `fitModel` itself, not +/// with respect to the curve of a plotted pdf. Make sure that the arguments +/// are consistent with what is plotted (same normalization and range, same +/// error model via DataError()), or the residual/pull histogram will not +/// match the displayed data and curve. +/// +/// \note If the model depends on observables beyond the data's observable +/// (e.g. a conditional pdf), they are evaluated at their current values. +/// There is no exact equivalent for a projection with ProjWData(). +/// +/// \note For internal use of RooAbsReal::createIntegral(), a scratch named +/// range is left on the model's observable after the call. +/// +/// Only 1-dimensional data is supported, otherwise std::invalid_argument is +/// thrown; the same exception is thrown for empty input data, unknown named +/// ranges or binnings, or when the model has no support in the normalization +/// range. For higher-dimensional fits, project both the data and the +/// model onto the observable of interest: +/// ~~~{.cpp} +/// std::unique_ptr projData{data.reduce(RooArgSet{x})}; +/// std::unique_ptr projModel{model.createProjection(RooArgSet{y})}; +/// auto pulls = RooFit::makePullHist(*projModel, *projData); +/// ~~~ +/// For a residual/pull histogram against one component of a composite pdf, +/// pass the component pdf itself, scaled by its coefficient: +/// ~~~{.cpp} +/// auto& background = static_cast(*addPdf.pdfList().find("bkg")); +/// auto pulls = RooFit::makePullHist(background, data, +/// RooFit::Normalization(bkgFraction.getVal())); +/// ~~~ +std::unique_ptr +makeResidHist(RooAbsReal &fitModel, RooAbsData const &data, RooCmdArg const &arg1 = {}, RooCmdArg const &arg2 = {}, + RooCmdArg const &arg3 = {}, RooCmdArg const &arg4 = {}, RooCmdArg const &arg5 = {}, + RooCmdArg const &arg6 = {}, RooCmdArg const &arg7 = {}, RooCmdArg const &arg8 = {}); + +/// Like makeResidHist(), but the residuals are divided by the corresponding +/// data uncertainty, creating a pull distribution. +std::unique_ptr +makePullHist(RooAbsReal &fitModel, RooAbsData const &data, RooCmdArg const &arg1 = {}, RooCmdArg const &arg2 = {}, + RooCmdArg const &arg3 = {}, RooCmdArg const &arg4 = {}, RooCmdArg const &arg5 = {}, + RooCmdArg const &arg6 = {}, RooCmdArg const &arg7 = {}, RooCmdArg const &arg8 = {}); + +/** + * @} + */ + namespace Detail { // Function to pack an arbitrary number of RooCmdArgs into a RooLinkedList. Implementation detail of many high-level RooFit functions. diff --git a/roofit/roofitcore/inc/RooHist.h b/roofit/roofitcore/inc/RooHist.h index 05392be18a73d..dbe58b601ef47 100644 --- a/roofit/roofitcore/inc/RooHist.h +++ b/roofit/roofitcore/inc/RooHist.h @@ -75,6 +75,11 @@ class RooHist : public TGraphAsymmErrors, public RooPlotable { bool hasIdenticalBinning(const RooHist& other) const ; + /// Compute residuals with respect to a curve. + /// \note For residuals/pulls of binned data against a fitted model, the + /// more accurate RooFit::makeResidHist() and RooFit::makePullHist() are + /// the recommended interface: they integrate the model itself exactly over + /// each bin, instead of interpolating the curve. RooHist* makeResidHist(const RooCurve& curve,bool normalize=false, bool useAverage=false) const; RooHist* makePullHist(const RooCurve& curve, bool useAverage=false) const {return makeResidHist(curve,true,useAverage); } diff --git a/roofit/roofitcore/inc/RooPlot.h b/roofit/roofitcore/inc/RooPlot.h index a55a3cc86aa14..659298db284bf 100644 --- a/roofit/roofitcore/inc/RooPlot.h +++ b/roofit/roofitcore/inc/RooPlot.h @@ -174,8 +174,14 @@ class RooPlot : public TNamed, public RooPrintable { double chiSquare(int nFitParam=0) const { return chiSquare(nullptr,nullptr,nFitParam) ; } double chiSquare(const char* pdfname, const char* histname, int nFitParam=0) const ; + // NOTE for developers: residHist() and pullHist() are superseded by + // RooFit::makeResidHist() / RooFit::makePullHist(). They should be + // deprecated in ROOT v6.44 and removed in v6.46. RooHist* residHist(const char* histname=nullptr, const char* pdfname=nullptr,bool normalize=false, bool useAverage=true) const ; - ///Uses residHist() and sets normalize=true + ///Uses residHist() and sets normalize=true. + ///\note For pulls of binned data against a fitted model, the more accurate + ///RooFit::makePullHist() is the recommended interface: it integrates the + ///model itself exactly over each bin, instead of interpolating the curve. RooHist* pullHist(const char* histname=nullptr, const char* pdfname=nullptr, bool useAverage=true) const { return residHist(histname,pdfname,true,useAverage); } diff --git a/roofit/roofitcore/src/RooGlobalFunc.cxx b/roofit/roofitcore/src/RooGlobalFunc.cxx index f6c183b208005..f970ddae959bf 100644 --- a/roofit/roofitcore/src/RooGlobalFunc.cxx +++ b/roofit/roofitcore/src/RooGlobalFunc.cxx @@ -18,19 +18,23 @@ #include +#include #include #include #include #include +#include #include #include #include #include #include +#include #include #include #include #include +#include #include "RooEvaluatorWrapper.h" @@ -38,6 +42,11 @@ #include #include +#include +#include +#include +#include +#include namespace RooFit { @@ -1155,6 +1164,299 @@ RooConstVar &RooConst(double val) return RooRealConstant::value(val); } +namespace { + +/// Implementation of RooFit::makeResidHist() and RooFit::makePullHist(). +std::unique_ptr makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData const &data, bool normalize, + RooCmdArg const &arg1, RooCmdArg const &arg2, RooCmdArg const &arg3, + RooCmdArg const &arg4, RooCmdArg const &arg5, RooCmdArg const &arg6, + RooCmdArg const &arg7, RooCmdArg const &arg8) +{ + const char *funcName = normalize ? "RooFit::makePullHist" : "RooFit::makeResidHist"; + + RooCmdConfig pc(funcName); + pc.defineDouble("normScale", "Normalization", 0, 1.0); + pc.defineInt("normType", "Normalization", 0, 0); + pc.defineInt("etype", "DataError", 0, static_cast(RooAbsData::Auto)); + pc.defineString("name", "Name", 0, ""); + pc.defineString("title", "Title", 0, ""); + pc.defineObject("binning", "Binning", 0); + pc.defineString("binningName", "BinningName", 0, ""); + pc.defineInt("nbins", "BinningSpec", 0, 0); + pc.defineDouble("xlo", "BinningSpec", 0, 0.0); + pc.defineDouble("xhi", "BinningSpec", 1, 0.0); + pc.defineMutex("Binning", "BinningName", "BinningSpec"); + pc.defineDouble("rangeLo", "Range", 0, -999.); + pc.defineDouble("rangeHi", "Range", 1, -999.); + pc.defineString("rangeName", "RangeWithName", 0, "", true); + pc.defineMutex("Range", "RangeWithName"); + pc.process(arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8); + if (!pc.ok(true)) { + throw std::invalid_argument(std::string(funcName) + "(): unrecognized arguments"); + } + + if (data.numEntries() == 0) { + throw std::invalid_argument(std::string(funcName) + "(): input data is empty"); + } + + // Only the pure scale-factor form of Normalization() is supported; the + // scale types Normalization(scale, type) like Relative or NumEvent are + // RooPlot-specific. + if (pc.hasProcessed("Normalization") && pc.getInt("normType") != static_cast(RooAbsReal::Relative)) { + throw std::invalid_argument(std::string(funcName) + + "(): only Normalization(scaleFactor) is supported, without a scale type"); + } + + const auto etype = static_cast(pc.getInt("etype")); + if (etype != RooAbsData::Auto && etype != RooAbsData::Poisson && etype != RooAbsData::SumW2) { + throw std::invalid_argument(std::string(funcName) + "(): unsupported DataError(), only Auto, Poisson and " + "SumW2 are supported"); + } + + // The observable must be a single RooRealVar (we need named ranges on it + // for the per-bin integration). This check will also catch multi-column + // datasets and RooDataHists with category dimensions (like reduced + // simultaneous-fit data). + RooRealVar *obs = nullptr; + { + const RooArgSet &coords = *data.get(); + if (coords.size() != 1) { + throw std::invalid_argument(std::string(funcName) + "(): only 1-dimensional data is supported"); + } + obs = dynamic_cast(coords[0]); + } + if (!obs) { + throw std::invalid_argument(std::string(funcName) + "(): observable must be a RooRealVar"); + } + + // The observable instance found in the data is not necessarily the same + // object as the model's observable (a RooDataSet owns its own copies). + // For the model integration, we must use the model's own variable with + // the same name, otherwise the interpretation of the integration ranges + // is model-dependent (e.g. works for plain pdfs, but not for + // RooProjectedPdf). If the model doesn't have an observable with this + // name (e.g. it's a constant in the observable), we still use the + // data-side one and the model integral is analytic in that case. + RooRealVar *intVar = obs; + { + std::unique_ptr modelObs{fitModel.getObservables(*data.get())}; + if (auto *modelVar = dynamic_cast(modelObs->find(obs->GetName()))) { + intVar = modelVar; + } + } + + // Resolve the optional Binning()/BinningName() arguments to the binning + // that should be used for the binning of unbinned data. + std::unique_ptr ownedBinning; + const RooAbsBinning *binning = nullptr; + if (pc.hasProcessed("Binning")) { + binning = static_cast(pc.getObject("binning")); + } else if (pc.hasProcessed("BinningName")) { + const char *binningName = pc.getString("binningName", nullptr, true); + // The named binning may live on the model-side or the data-side + // variable + if (intVar->hasBinning(binningName)) { + binning = &intVar->getBinning(binningName); + } else if (obs->hasBinning(binningName)) { + binning = &obs->getBinning(binningName); + } else { + throw std::invalid_argument(std::string(funcName) + "(): unknown binning '" + binningName + "'"); + } + } else if (pc.hasProcessed("BinningSpec")) { + double xlo = pc.getDouble("xlo"); + double xhi = pc.getDouble("xhi"); + // The RooFit convention (see RooAbsData::plotOn) is that Binning(nbins) + // encodes equal xlo and xhi, meaning "take the range from the variable" + if (xlo == xhi) { + xlo = obs->getMin(); + xhi = obs->getMax(); + } + ownedBinning = std::make_unique(xlo, xhi, pc.getInt("nbins")); + binning = ownedBinning.get(); + } + + // If the input data is unbinned, bin it. A binning passed with the + // Binning() arguments is realized on a clone of the observable, so the + // binning of the user's variables are never modified. + std::unique_ptr ownedBinData; + std::unique_ptr ownedObs; + RooDataHist const *binData = dynamic_cast(&data); + if (!binData) { + const RooAbsArg *usedObs = obs; + if (binning) { + ownedObs.reset(static_cast(obs->Clone())); + static_cast(*ownedObs).setBinning(*binning); + usedObs = ownedObs.get(); + } + std::string binnedName = std::string(data.GetName()) + "_binned"; + ownedBinData = std::make_unique(binnedName, data.GetTitle(), RooArgSet{*usedObs}, data); + binData = ownedBinData.get(); + } else if (binning) { + throw std::invalid_argument(std::string(funcName) + + "(): Binning() arguments are not supported for already binned data"); + } + + const double scaleFactor = pc.getDouble("normScale"); + + std::string name = pc.getString("name"); + if (name.empty()) { + name = (normalize ? "pull_" : "resid_") + std::string(data.GetName()) + "_" + fitModel.GetName(); + } + std::string title = pc.getString("title"); + if (title.empty()) { + title = std::string(normalize ? "Pull of " : "Residual of ") + data.GetTitle() + " and " + fitModel.GetTitle() + + " (integrated per bin)"; + } + + // Use the binning that the bin contents were computed with, which can be + // different from the default binning of the observable. + const RooAbsBinning &dataBinning = *binData->getBinnings().front(); + const int nBins = dataBinning.numBins(); + + auto out = std::make_unique(); + out->SetName(name.c_str()); + out->SetTitle(title.c_str()); + + RooArgSet obsSet{*intVar}; + RooArgSet normSet{*intVar}; + + // Resolve the normalization ranges (default: the full range of the + // observable). Only bins fully inside one of the ranges get points in the + // residual/pull histogram, and the model expectation is normalized to + // the data weight in these ranges. Range("name") can be given multiple + // times for a union of ranges (sidebands); the names are concatenated + // into the comma-separated RangeWithName argument by the usual RooFit + // convention. + std::vector> ranges; + if (pc.hasProcessed("RangeWithName")) { + std::string rangeNames = pc.getString("rangeName", nullptr, true); + for (std::istringstream tn(rangeNames); tn.good();) { + std::string rname; + std::getline(tn, rname, ','); + if (rname.empty()) + continue; + if (!intVar->hasRange(rname.c_str())) { + throw std::invalid_argument(std::string(funcName) + "(): unknown range '" + rname + "'"); + } + ranges.push_back(intVar->getRange(rname.c_str())); + } + } else if (pc.hasProcessed("Range")) { + ranges.emplace_back(pc.getDouble("rangeLo"), pc.getDouble("rangeHi")); + } else { + ranges.emplace_back(obs->getMin(), obs->getMax()); + } + + // The numeric integration with createIntegral() is noisy on INFO level + // (it initializes an integrator for every bin); suppress it, but keep + // warnings visible. + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING}; + + // Scratch ranges on the model observable that are redefined for every + // bin, used to integrate the model exactly over each bin (same mechanism + // as RooBinSamplingPdf). The ranges are intentionally named uniquely per + // call to not clobber existing ranges of the same name on the user's + // variable; they are shared ranges and stay on the variable after the + // call. + static std::atomic nScratchRanges{0}; + const std::string scratchRange = "_binInt_" + std::to_string(nScratchRanges++); + const std::string rangeNameBin = scratchRange + "_bin"; + + RooAbsData::ErrorType errType = + etype == RooAbsData::Auto ? (binData->isNonPoissonWeighted() ? RooAbsData::SumW2 : RooAbsData::Poisson) : etype; + + // The model integral over the (union of) normalization ranges + double normIntegral = 0.; + for (std::size_t iRange = 0; iRange < ranges.size(); ++iRange) { + std::string rangeNameRange = scratchRange + "_range" + std::to_string(iRange); + intVar->setRange(rangeNameRange.c_str(), ranges[iRange].first, ranges[iRange].second); + normIntegral += fitModel.createIntegral(obsSet, normSet, rangeNameRange.c_str())->getVal(); + } + if (normIntegral == 0.) { + throw std::invalid_argument(std::string(funcName) + "(): model has no support in the normalization range(s)"); + } + + // Tolerate floating-point imprecision of bin boundaries in the range + // comparison + double epsMax = 0.; + for (auto const &range : ranges) { + epsMax = std::max(epsMax, std::max(std::abs(range.first), std::abs(range.second))); + } + const double eps = 10. * std::numeric_limits::epsilon() * epsMax; + + // The data weight inside the normalization ranges, summing only the kept + // bins (the convention for plotted curves is analogous: a curve is + // rescaled to the event count of the plotted, possibly cut, histogram). + double nData = 0.; + for (int iBin = 0; iBin < nBins; ++iBin) { + const double binLo = dataBinning.binLow(iBin); + const double binHi = dataBinning.binHigh(iBin); + const bool inRange = std::any_of(ranges.begin(), ranges.end(), [&](auto const &r) { + return binLo >= r.first - eps && binHi <= r.second + eps; + }); + if (inRange) { + nData += binData->weight(iBin); + } + } + + for (int iBin = 0; iBin < nBins; ++iBin) { + binData->get(iBin); + + const double binLo = dataBinning.binLow(iBin); + const double binHi = dataBinning.binHigh(iBin); + // Only show points for bins inside the (union of) normalization ranges + const bool inRange = std::any_of(ranges.begin(), ranges.end(), [&](auto const &r) { + return binLo >= r.first - eps && binHi <= r.second + eps; + }); + if (!inRange) + continue; + + intVar->setRange(rangeNameBin.c_str(), binLo, binHi); + const double fraction = fitModel.createIntegral(obsSet, normSet, rangeNameBin.c_str())->getVal() / normIntegral; + const double expected = scaleFactor * nData * fraction; + + const double actual = binData->weight(iBin); + double eylo, eyhi; + binData->weightError(eylo, eyhi, errType); + + double resid = actual - expected; + if (normalize) { + const double norm = resid > 0 ? eylo : eyhi; + if (norm == 0.) { + oocoutW(&fitModel, Plotting) << funcName << "(" << fitModel.GetName() << ") WARNING: point " << iBin + << " has zero error, setting residual to zero" << std::endl; + resid = eyhi = eylo = 0.; + } else { + resid /= norm; + eyhi /= norm; + eylo /= norm; + } + } + out->addBinWithError(dataBinning.binCenter(iBin), resid, eylo, eyhi); + } + + return out; +} + +} // namespace + +/// \see RooFit::makeResidHist() +std::unique_ptr makeResidHist(RooAbsReal &fitModel, RooAbsData const &data, RooCmdArg const &arg1, + RooCmdArg const &arg2, RooCmdArg const &arg3, RooCmdArg const &arg4, + RooCmdArg const &arg5, RooCmdArg const &arg6, RooCmdArg const &arg7, + RooCmdArg const &arg8) +{ + return makeResidOrPullHist(fitModel, data, false, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8); +} + +/// \see RooFit::makePullHist() +std::unique_ptr makePullHist(RooAbsReal &fitModel, RooAbsData const &data, RooCmdArg const &arg1, + RooCmdArg const &arg2, RooCmdArg const &arg3, RooCmdArg const &arg4, + RooCmdArg const &arg5, RooCmdArg const &arg6, RooCmdArg const &arg7, + RooCmdArg const &arg8) +{ + return makeResidOrPullHist(fitModel, data, true, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8); +} + namespace Detail { RooCmdArg SliceFlatMap(FlatMap const &args) diff --git a/roofit/roofitcore/src/RooPlot.cxx b/roofit/roofitcore/src/RooPlot.cxx index 9d0cf19e84f1a..3986ed9da3d9e 100644 --- a/roofit/roofitcore/src/RooPlot.cxx +++ b/roofit/roofitcore/src/RooPlot.cxx @@ -43,6 +43,8 @@ object onto a one-dimensional plot. #include "RooPlot.h" +#include + #include "RooAbsReal.h" #include "RooAbsRealLValue.h" #include "RooPlotable.h" @@ -1064,8 +1066,22 @@ double RooPlot::chiSquare(const char* curvename, const char* histname, int nFitP /// of the histogram, effectively returning a pull histogram. /// \param useAverage If true, the histogram is compared with the curve averaged in each bin. /// Otherwise, the curve is evaluated at the bin centres, which is not accurate for strongly curved distributions. +/// +/// \note For the comparison of binned data with a fitted model, the more +/// accurate RooFit::makeResidHist() and RooFit::makePullHist() are the +/// recommended interface: they integrate the model itself exactly over each +/// bin, instead of interpolating the plotted curve. RooHist* RooPlot::residHist(const char* histname, const char* curvename, bool normalize, bool useAverage) const { + static std::atomic warned{false}; + if (!warned.exchange(true)) { + coutW(Plotting) << "RooPlot::residHist/pullHist compare data with the interpolated curve, which is biased for " + "sharply peaked models (systematic residual 'wiggle'). For binned data compared to a fitted " + "model, use RooFit::makeResidHist() / RooFit::makePullHist(), which integrate the model exactly " + "over each bin. These methods will be deprecated in ROOT v6.44." + << std::endl; + } + // Find all curve objects with the name "curvename" or the name of the last // plotted curve (there might be multiple in the case of multi-range fits). std::vector curves; diff --git a/roofit/roofitcore/test/testRooHist.cxx b/roofit/roofitcore/test/testRooHist.cxx index 36b2d7d8e01e4..adf1c7636eb3a 100644 --- a/roofit/roofitcore/test/testRooHist.cxx +++ b/roofit/roofitcore/test/testRooHist.cxx @@ -4,6 +4,10 @@ #include #include #include +#include +#include +#include +#include #include #include #include @@ -13,10 +17,12 @@ #include #include +#include #include #include +#include /// Check that the values returned by `RooHist::getFitRangeNEvt(double xmin, /// double xmax)` are correct also for non-uniform binning. Covers ROOT-9649. @@ -120,3 +126,254 @@ TEST(RooHist, ResidualsAndPulls) // The bulk of the Gaussian sample must have populated bins EXPECT_GT(nChecked, 30); } + +namespace { + +// Shared setup for the exact-bin-integration tests of +// RooFit::makeResidHist()/makePullHist(): a narrowly-peaked Gaussian and a +// TH1D out of it with 50 bins + +constexpr double kGaussMean = 125.; +constexpr double kGaussSigma = 1.4; +constexpr double kGaussX0 = 110.; +constexpr double kGaussX1 = 130.; +constexpr int kNBins = 50; + +/// Integral of a Gaussian with (kGaussMean, kGaussSigma) between lo and hi. +double gaussIntegral(double lo, double hi) +{ + auto phi = [](double u) { return 0.5 * (1 + std::erf(u / std::sqrt(2.))); }; + return phi((hi - kGaussMean) / kGaussSigma) - phi((lo - kGaussMean) / kGaussSigma); +} + +/// The Gaussian as a model plus observable; members are in construction order. +struct GaussModel { + RooRealVar x{"x", "x", kGaussX0, kGaussX1}; + RooRealVar mean{"mean", "mean", kGaussMean}; + RooRealVar sigma{"sigma", "sigma", kGaussSigma}; + RooGenericPdf gauss{"gauss", "gauss", "exp(-0.5 * pow((x - mean) / sigma, 2))", RooArgSet{x, mean, sigma}}; + + GaussModel() { x.setBins(kNBins); } +}; + +/// Fill the histogram with the exact Gaussian bin contents on the given +/// ranges, normalized to nEvents on those ranges. Bins outside get the +/// (constant) otherContent fill, which must not influence the in-range +/// expectations once a Range() selection is made. +void fillExactGauss(TH1D &hist, double nEvents, std::vector> const &ranges, + double otherContent = 0.) +{ + double norm = 0.; + for (auto const &r : ranges) { + norm += gaussIntegral(r.first, r.second); + } + for (int iBin = 1; iBin <= hist.GetNbinsX(); ++iBin) { + const double lo = hist.GetBinLowEdge(iBin); + const double hi = lo + hist.GetBinWidth(iBin); + // epsilon in the boundary check to guard against floating point + // arithmetic on the bin boundaries + const bool inRange = std::any_of(ranges.begin(), ranges.end(), + [&](auto const &r) { return lo >= r.first - 1e-9 && hi <= r.second + 1e-9; }); + hist.SetBinContent(iBin, inRange ? nEvents / norm * gaussIntegral(lo, hi) : otherContent); + } +} + +/// Check that all residuals in the histogram are consistent with zero, given +/// the exact per-bin contents of `hist`. The tolerance is relative to the bin +/// content, with an absolute floor to tolerate the ~1e-3 precision of the +/// numeric integration in the far tails. +void expectExactResiduals(RooHist const &resid, TH1D const &hist, double tolerance) +{ + for (int i = 0; i < resid.GetN(); ++i) { + const int iBin = int((resid.GetPointX(i) - kGaussX0) * kNBins / (kGaussX1 - kGaussX0)); + const double content = hist.GetBinContent(iBin + 1); + if (content == 0.) + continue; + EXPECT_LT(std::abs(resid.GetPointY(i)), tolerance * std::max(content, 1e-6)) << "bin " << iBin; + } +} + +} // namespace + +/// RooFit::makeResidHist() and RooFit::makePullHist() must compare the data +/// with the model integrated exactly over each bin. If the model instead was +/// evaluated at the bin center, a sharply-peaked Gaussian would show large +/// residuals in the peak bins (O(10^3) events here) even for data that is +/// filled with the exact per-bin expectation values. +TEST(RooHist, MakeResidAndPullHist) +{ + using namespace RooFit; + + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING}; + + GaussModel model; + + TH1D hist("hist", "hist", kNBins, kGaussX0, kGaussX1); + fillExactGauss(hist, 1e6, {{kGaussX0, kGaussX1}}); + RooDataHist binData("binData", "binData", RooArgSet(model.x), &hist); + + // The residuals of the exact model must be (numerically) zero + std::unique_ptr resid{makeResidHist(model.gauss, binData)}; + std::unique_ptr pull{makePullHist(model.gauss, binData)}; + ASSERT_EQ(resid->GetN(), kNBins); + ASSERT_EQ(pull->GetN(), kNBins); + expectExactResiduals(*resid, hist, 1e-3); + for (int i = 0; i < pull->GetN(); ++i) { + EXPECT_LT(std::abs(pull->GetPointY(i)), 1e-3) << "bin " << i; + } + + // The Name() and Normalization() command args + std::unique_ptr scaled{makeResidHist(model.gauss, binData, Name("scaled"), Normalization(0.5))}; + ASSERT_EQ(scaled->GetN(), kNBins); + EXPECT_STREQ(scaled->GetName(), "scaled"); + for (int i = 0; i < scaled->GetN(); ++i) { + const double content = hist.GetBinContent(i + 1); + if (content == 0.) + continue; + EXPECT_NEAR(scaled->GetPointY(i), 0.5 * content, 1e-3 * std::max(content, 1e-6)) << "bin " << i; + } +} + +/// RooFit::makePullHist() and RooFit::makeResidHist() also accept unbinned +/// input data, which is binned internally with the binning of the observable. +TEST(RooHist, MakeResidAndPullHistUnbinned) +{ + using namespace RooFit; + + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING}; + RooRandom::randomGenerator()->SetSeed(1337); + + RooRealVar x("x", "x", -10, 10); + RooRealVar mean("mean", "mean", 0.); + RooRealVar sigma("sigma", "sigma", 3.); + RooGenericPdf gauss("gauss", "gauss", "exp(-0.5 * pow((x - mean) / sigma, 2))", RooArgSet{x, mean, sigma}); + + std::unique_ptr data{gauss.generate(x, 10000)}; + + // The Binning() argument overrides the default binning of the observable + std::unique_ptr resid{makeResidHist(gauss, *data, Binning(40))}; + std::unique_ptr pull{makePullHist(gauss, *data, Binning(40))}; + EXPECT_EQ(resid->GetN(), 40); + EXPECT_EQ(pull->GetN(), 40); + + // The correct model should give small pulls + int nChecked = 0; + for (int i = 0; i < pull->GetN(); ++i) { + if (pull->GetErrorYhigh(i) == 0.) + continue; + EXPECT_LT(std::abs(pull->GetPointY(i)), 5.) << "point " << i; + ++nChecked; + } + EXPECT_GT(nChecked, 30); +} + +/// The Range() argument of RooFit::makeResidHist()/makePullHist(): only bins +/// inside the range get points, and the model is normalized to the data +/// weight in the range (as needed for sideband fits). +TEST(RooHist, MakeResidAndPullHistRange) +{ + using namespace RooFit; + + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING}; + + GaussModel model; + + // [rlo, rhi] chosen on bin boundaries: bins 24..33 of 50. The asymmetric + // range contains only ~16% of the model, so a wrong range normalization + // would produce huge residuals. + const double rlo = 119.6, rhi = 123.6; + + TH1D hist("hist", "hist", kNBins, kGaussX0, kGaussX1); + // The constant content outside of the range must not influence the + // expectation normalization, which only counts the data weight in the range + fillExactGauss(hist, 1e6, {{rlo, rhi}}, 2.e6); + RooDataHist binData("binData", "binData", RooArgSet(model.x), &hist); + + // With coordinate range. RooGenericPdf is integrated numerically, whose + // precision is limited to ~1e-3 relative in the far tails, hence the 5e-3 + // tolerance. + std::unique_ptr resid{makeResidHist(model.gauss, binData, Range(rlo, rhi))}; + EXPECT_EQ(resid->GetN(), 10); + expectExactResiduals(*resid, hist, 5e-3); + + // With named range + model.x.setRange("central", rlo, rhi); + std::unique_ptr resid2{makeResidHist(model.gauss, binData, Range("central"))}; + ASSERT_EQ(resid2->GetN(), resid->GetN()); + for (int i = 0; i < resid2->GetN(); ++i) { + EXPECT_DOUBLE_EQ(resid2->GetPointX(i), resid->GetPointX(i)); + EXPECT_DOUBLE_EQ(resid2->GetPointY(i), resid->GetPointY(i)); + } + + // An unknown named range must throw + EXPECT_THROW(makeResidHist(model.gauss, binData, Range("unknown")), std::invalid_argument); +} + +/// The Range("name") argument of RooFit::makeResidHist() can be passed +/// multiple times for a union of named ranges (sidebands). +TEST(RooHist, MakeResidAndPullHistMultiRange) +{ + using namespace RooFit; + + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING}; + + GaussModel model; + + // Bins 13..16 and 33..36 (on bin boundaries) + model.x.setRange("sb1", 115.2, 116.8); + model.x.setRange("sb2", 123.2, 124.8); + + TH1D hist("hist", "hist", kNBins, kGaussX0, kGaussX1); + fillExactGauss(hist, 1e6, {{115.2, 116.8}, {123.2, 124.8}}, 3.e6); + RooDataHist binData("binData", "binData", RooArgSet(model.x), &hist); + + std::unique_ptr resid{makeResidHist(model.gauss, binData, Range("sb1"), Range("sb2"))}; + ASSERT_EQ(resid->GetN(), 8); + expectExactResiduals(*resid, hist, 5e-3); +} + +/// The documented failure modes of RooFit::makeResidHist() and +/// RooFit::makePullHist() that throw exceptions. +TEST(RooHist, MakeResidAndPullHistErrors) +{ + using namespace RooFit; + + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING}; + + RooRealVar x("x", "x", -10, 10); + RooRealVar y("y", "y", -10, 10); + x.setBins(40); + RooRealVar mean("mean", "mean", 0.); + RooRealVar sigma("sigma", "sigma", 3.); + RooGenericPdf gauss("gauss", "gauss", "exp(-0.5 * pow((x - mean) / sigma, 2))", RooArgSet{x, mean, sigma}); + + RooRandom::randomGenerator()->SetSeed(1337); + std::unique_ptr data{gauss.generate(x, 1000)}; + + // Multi-dimensional data + TH2D hist2("hist2", "hist2", 4, -10, 10, 4, -10, 10); + hist2.Fill(0., 0.); + RooDataHist hist2D("hist2D", "hist2D", RooArgSet(x, y), &hist2); + EXPECT_THROW(makeResidHist(gauss, hist2D), std::invalid_argument); + + // Empty datasets + RooDataSet emptyData("emptyData", "emptyData", RooArgSet(x)); + EXPECT_THROW(makeResidHist(gauss, emptyData), std::invalid_argument); + + // Binning() arguments with already binned data + RooDataHist binned("binned", "binned", RooArgSet(x), *data); + EXPECT_THROW(makeResidHist(gauss, binned, Binning(20)), std::invalid_argument); + + // Unknown named ranges + EXPECT_THROW(makeResidHist(gauss, *data, Range("unknown")), std::invalid_argument); + + // Unsupported Normalization scale types + EXPECT_THROW(makeResidHist(gauss, *data, Normalization(0.5, RooAbsReal::NumEvent)), std::invalid_argument); + + // Unsupported DataError types + EXPECT_THROW(makeResidHist(gauss, *data, DataError(99)), std::invalid_argument); + + // A model with no support in the range must throw + RooGenericPdf zeroModel("zeroModel", "zeroModel", "x * 0.", RooArgSet{x}); + EXPECT_THROW(makeResidHist(zeroModel, *data), std::invalid_argument); +} diff --git a/tutorials/roofit/roofit/rf109_chi2residpull.C b/tutorials/roofit/roofit/rf109_chi2residpull.C index cd3834a8e6a0e..d548f660ea1c5 100644 --- a/tutorials/roofit/roofit/rf109_chi2residpull.C +++ b/tutorials/roofit/roofit/rf109_chi2residpull.C @@ -59,19 +59,25 @@ void rf109_chi2residpull() // S h o w r e s i d u a l a n d p u l l d i s t s // ------------------------------------------------------- - // Construct a histogram with the residuals of the data w.r.t. the curve - RooHist *hresid = frame1->residHist(); + // For the residuals and pulls, we compare the binned data with the model + // integrated exactly over each bin. This avoids the systematic "wiggle" + // that curve-based residuals show for sharply peaked pdfs. - // Construct a histogram with the pulls of the data w.r.t the curve - RooHist *hpull = frame1->pullHist(); + // Construct a histogram with the residuals of the data w.r.t. the model. + // The Binning() argument selects how the unbinned dataset is binned. + std::unique_ptr hresid{makeResidHist(gauss, *data, Binning(40))}; - // Create a new frame to draw the residual distribution and add the distribution to the frame + // Construct a histogram with the pulls of the data w.r.t the model + std::unique_ptr hpull{makePullHist(gauss, *data, Binning(40))}; + + // Create a new frame to draw the residual distribution and add the distribution to the frame. + // Note that addPlotable() transfers ownership of the histogram to the frame. RooPlot *frame2 = x.frame(Title("Residual Distribution")); - frame2->addPlotable(hresid, "P"); + frame2->addPlotable(hresid.release(), "P"); // Create a new frame to draw the pull distribution and add the distribution to the frame RooPlot *frame3 = x.frame(Title("Pull Distribution")); - frame3->addPlotable(hpull, "P"); + frame3->addPlotable(hpull.release(), "P"); TCanvas *c = new TCanvas("rf109_chi2residpull", "rf109_chi2residpull", 900, 300); c->Divide(3); diff --git a/tutorials/roofit/roofit/rf109_chi2residpull.py b/tutorials/roofit/roofit/rf109_chi2residpull.py index 18c699b54afe0..aa4b262f05125 100644 --- a/tutorials/roofit/roofit/rf109_chi2residpull.py +++ b/tutorials/roofit/roofit/rf109_chi2residpull.py @@ -50,21 +50,29 @@ # Show residual and pull dists # ------------------------------------------------------- -# Construct a histogram with the residuals of the data w.r.t. the curve -hresid = frame1.residHist() +# For the residuals and pulls, we compare the binned data with the model +# integrated exactly over each bin. This avoids the systematic "wiggle" +# that curve-based residuals show for sharply peaked pdfs. The Binning() +# argument selects how the unbinned dataset is binned. -# Construct a histogram with the pulls of the data w.r.t the curve -hpull = frame1.pullHist() +# Construct a histogram with the residuals of the data w.r.t. the model +hresid = ROOT.RooFit.makeResidHist(gauss, data, ROOT.RooFit.Binning(40)) + +# Construct a histogram with the pulls of the data w.r.t the model +hpull = ROOT.RooFit.makePullHist(gauss, data, ROOT.RooFit.Binning(40)) # Create a frame to draw the residual distribution and add the -# distribution to the frame +# distribution to the frame. Note that addPlotable() transfers ownership +# of the histogram to the frame. frame2 = x.frame(Title="Residual Distribution") frame2.addPlotable(hresid, "P") +ROOT.SetOwnership(hresid, False) # Create a frame to draw the pull distribution and add the distribution to # the frame frame3 = x.frame(Title="Pull Distribution") frame3.addPlotable(hpull, "P") +ROOT.SetOwnership(hpull, False) c = ROOT.TCanvas("rf109_chi2residpull", "rf109_chi2residpull", 900, 300) c.Divide(3)