From 0876d396d051855fbe9a01a6a590492468e21edb Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Mon, 14 Sep 2026 08:00:43 +0000 Subject: [PATCH 1/5] [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 a233715e4835a..89d310a0b6f05 100644 --- a/README/ReleaseNotes/v642/index.md +++ b/README/ReleaseNotes/v642/index.md @@ -201,6 +201,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) From 4b87ceb333300e93a888a723216f12667e538e4a Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Mon, 12 Oct 2026 20:13:28 +0000 Subject: [PATCH 2/5] [RF] Address review comments on makeResidHist()/makePullHist() - Declare makeResidHist() and makePullHist() as variadic templates that forward the command arguments as a RooLinkedList to the implementation in RooFit::Detail, like RooAbsPdf::fitTo(), instead of eight RooCmdArg parameters. - Fix the memory leak of the raw pointers returned by createIntegral(). The integral object is intentionally not reused across bins: the integration strategy of a RooRealIntegral (e.g. the analytical integral code of RooHistPdf, which checks if the range is the full range) is fixed at creation time and can depend on the range bounds. - Don't leave the scratch range for the per-bin integration on the observable: a free name is found with a local counter instead of a static atomic, and the range is removed at the end of the call. - Use a plain static bool for the once-per-process warning in RooPlot::residHist(). - Document multi-range selection with the usual comma-separated Range("a,b") convention of fitTo() and plotOn(). Passing Range() twice triggers the "argument is duplicated" warning of RooCmdConfig. - Doc wording suggestions, \copydetails for RooHist::makePullHist(), and RooFit:: qualification in the rf109 tutorial. - Tests: compare against a RooHistPdf with analytic integrals at 1e-10 tolerance in addition to the numerically integrated RooGenericPdf, check that no scratch range is left behind, use a named bin count, and replace the loose |pull| < 5 bound with a check of the pull definition, |pull| < 3.5 and chi2/ndf in [0.5, 1.5]. --- roofit/roofitcore/inc/RooGlobalFunc.h | 43 +++--- roofit/roofitcore/inc/RooHist.h | 2 + roofit/roofitcore/src/RooGlobalFunc.cxx | 87 +++++------ roofit/roofitcore/src/RooPlot.cxx | 7 +- roofit/roofitcore/test/testRooHist.cxx | 143 +++++++++++++----- tutorials/roofit/roofit/rf109_chi2residpull.C | 4 +- 6 files changed, 173 insertions(+), 113 deletions(-) diff --git a/roofit/roofitcore/inc/RooGlobalFunc.h b/roofit/roofitcore/inc/RooGlobalFunc.h index bf131a096d55b..4f88926465d61 100644 --- a/roofit/roofitcore/inc/RooGlobalFunc.h +++ b/roofit/roofitcore/inc/RooGlobalFunc.h @@ -490,17 +490,17 @@ RooConstVar& RooConst(double val) ; /// 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 bin (with RooAbsReal::createIntegral()), normalized to /// 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. +/// - `Range(lo, hi)` / `Range("name")` : compute residuals only in this range. +/// The model expectation is normalized in this range. Multiple +/// comma-separated range names can be specified 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. @@ -519,9 +519,6 @@ RooConstVar& RooConst(double val) ; /// (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 @@ -539,17 +536,13 @@ RooConstVar& RooConst(double val) ; /// 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 = {}); +template +std::unique_ptr makeResidHist(RooAbsReal &fitModel, RooAbsData const &data, CmdArgs_t const &...cmdArgs); /// 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 = {}); +template +std::unique_ptr makePullHist(RooAbsReal &fitModel, RooAbsData const &data, CmdArgs_t const &...cmdArgs); /** * @} @@ -607,8 +600,24 @@ RooCmdArg ImportFlatMap(FlatMap const &args); RooCmdArg ImportFlatMap(FlatMap const &args); RooCmdArg LinkFlatMap(FlatMap const &args); +// Implementation of RooFit::makeResidHist() and RooFit::makePullHist(). +std::unique_ptr +makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData const &data, bool normalize, RooLinkedList const &cmdList); + } // namespace Detail +template +std::unique_ptr makeResidHist(RooAbsReal &fitModel, RooAbsData const &data, CmdArgs_t const &...cmdArgs) +{ + return Detail::makeResidOrPullHist(fitModel, data, false, *Detail::createCmdList(&cmdArgs...)); +} + +template +std::unique_ptr makePullHist(RooAbsReal &fitModel, RooAbsData const &data, CmdArgs_t const &...cmdArgs) +{ + return Detail::makeResidOrPullHist(fitModel, data, true, *Detail::createCmdList(&cmdArgs...)); +} + } // namespace RooFit namespace RooFitShortHand { diff --git a/roofit/roofitcore/inc/RooHist.h b/roofit/roofitcore/inc/RooHist.h index dbe58b601ef47..7e6bf15ffc33c 100644 --- a/roofit/roofitcore/inc/RooHist.h +++ b/roofit/roofitcore/inc/RooHist.h @@ -81,6 +81,8 @@ class RooHist : public TGraphAsymmErrors, public RooPlotable { /// 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; + /// Compute pulls with respect to a curve, i.e. residuals divided by the data uncertainty. + /// \copydetails makeResidHist() RooHist* makePullHist(const RooCurve& curve, bool useAverage=false) const {return makeResidHist(curve,true,useAverage); } diff --git a/roofit/roofitcore/src/RooGlobalFunc.cxx b/roofit/roofitcore/src/RooGlobalFunc.cxx index f970ddae959bf..4f96d1be2120d 100644 --- a/roofit/roofitcore/src/RooGlobalFunc.cxx +++ b/roofit/roofitcore/src/RooGlobalFunc.cxx @@ -42,11 +42,9 @@ #include #include -#include #include #include #include -#include namespace RooFit { @@ -1164,13 +1162,11 @@ RooConstVar &RooConst(double val) return RooRealConstant::value(val); } -namespace { +namespace Detail { /// 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) +std::unique_ptr +makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData const &data, bool normalize, RooLinkedList const &cmdList) { const char *funcName = normalize ? "RooFit::makePullHist" : "RooFit::makeResidHist"; @@ -1190,7 +1186,7 @@ std::unique_ptr makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData co 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); + pc.process(cmdList); if (!pc.ok(true)) { throw std::invalid_argument(std::string(funcName) + "(): unrecognized arguments"); } @@ -1323,10 +1319,8 @@ std::unique_ptr makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData co // 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. + // the data weight in these ranges. Multiple comma-separated range names + // select a union of ranges (sidebands), by the usual RooFit convention. std::vector> ranges; if (pc.hasProcessed("RangeWithName")) { std::string rangeNames = pc.getString("rangeName", nullptr, true); @@ -1351,25 +1345,43 @@ std::unique_ptr makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData co // 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"; + // A scratch range on the model observable that is redefined before each + // integration, used to integrate the model exactly over the normalization + // ranges and over each bin (same mechanism as RooBinSamplingPdf). To not + // clobber an existing range of the user's variable, the first name that + // is not in use yet is taken, and the range is removed again at the end + // of this function. + std::string scratchRange; + for (unsigned int i = 0;; ++i) { + scratchRange = "_binInt_" + std::to_string(i); + if (!intVar->hasRange(scratchRange.c_str())) + break; + } + struct ScratchRangeGuard { + RooRealVar &var; + std::string name; + ~ScratchRangeGuard() { var.removeBinning(name.c_str()); } + } scratchRangeGuard{*intVar, scratchRange}; + + // Integrates the model over the current bounds of the scratch range. A new + // integral object is created for each integration on purpose: the + // integration strategy of a RooRealIntegral (e.g. the analytical integral + // code of the integrand) is fixed at creation and may depend on the range + // bounds. RooHistPdf for instance only sums bin contents in a sub-range if + // the range was not the full range at creation time. + auto integrate = [&](double lo, double hi) { + intVar->setRange(scratchRange.c_str(), lo, hi); + std::unique_ptr integral{fitModel.createIntegral(obsSet, normSet, scratchRange.c_str())}; + return integral->getVal(); + }; 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(); + for (auto const &range : ranges) { + normIntegral += integrate(range.first, range.second); } if (normIntegral == 0.) { throw std::invalid_argument(std::string(funcName) + "(): model has no support in the normalization range(s)"); @@ -1410,8 +1422,7 @@ std::unique_ptr makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData co if (!inRange) continue; - intVar->setRange(rangeNameBin.c_str(), binLo, binHi); - const double fraction = fitModel.createIntegral(obsSet, normSet, rangeNameBin.c_str())->getVal() / normIntegral; + const double fraction = integrate(binLo, binHi) / normIntegral; const double expected = scaleFactor * nData * fraction; const double actual = binData->weight(iBin); @@ -1437,28 +1448,6 @@ std::unique_ptr makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData co 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) { return processFlatMap("SliceCatMany", processSliceItem, args); diff --git a/roofit/roofitcore/src/RooPlot.cxx b/roofit/roofitcore/src/RooPlot.cxx index 3986ed9da3d9e..a31671ff167b2 100644 --- a/roofit/roofitcore/src/RooPlot.cxx +++ b/roofit/roofitcore/src/RooPlot.cxx @@ -43,8 +43,6 @@ object onto a one-dimensional plot. #include "RooPlot.h" -#include - #include "RooAbsReal.h" #include "RooAbsRealLValue.h" #include "RooPlotable.h" @@ -1073,8 +1071,9 @@ double RooPlot::chiSquare(const char* curvename, const char* histname, int nFitP /// 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)) { + static bool warned = false; + if (!warned) { + warned = 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 " diff --git a/roofit/roofitcore/test/testRooHist.cxx b/roofit/roofitcore/test/testRooHist.cxx index adf1c7636eb3a..e647de05bb36f 100644 --- a/roofit/roofitcore/test/testRooHist.cxx +++ b/roofit/roofitcore/test/testRooHist.cxx @@ -6,8 +6,6 @@ #include #include #include -#include -#include #include #include #include @@ -21,8 +19,12 @@ #include +#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. @@ -146,16 +148,6 @@ double gaussIntegral(double lo, double hi) 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 @@ -178,10 +170,42 @@ void fillExactGauss(TH1D &hist, double nEvents, std::vector 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 residuals of the exact model must be zero up to the precision of + // the model integration: floating-point precision for the analytic + // integrals of the RooHistPdf, and the numeric integration precision for + // the RooGenericPdf. + for (auto const &[pdf, tol] : {std::pair{&model.gaussAnalytic, kTolAnalytic}, + std::pair{&model.gaussNumeric, kTolNumeric}}) { + std::unique_ptr resid{makeResidHist(*pdf, binData)}; + std::unique_ptr pull{makePullHist(*pdf, binData)}; + ASSERT_EQ(resid->GetN(), kNBins); + ASSERT_EQ(pull->GetN(), kNBins); + expectExactResiduals(*resid, hist, tol); + for (int i = 0; i < pull->GetN(); ++i) { + EXPECT_LT(std::abs(pull->GetPointY(i)), tol) << pdf->GetName() << " bin " << i; + } } // The Name() and Normalization() command args - std::unique_ptr scaled{makeResidHist(model.gauss, binData, Name("scaled"), Normalization(0.5))}; + std::unique_ptr scaled{makeResidHist(model.gaussAnalytic, 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; + EXPECT_NEAR(scaled->GetPointY(i), 0.5 * content, kTolAnalytic * std::max(content, 1e-6)) << "bin " << i; + } + + // The scratch range used for the per-bin integration must not be left on + // the observable + for (auto const &name : model.x.getBinningNames()) { + EXPECT_EQ(name.find("_binInt_"), std::string::npos) << "leftover range " << name; } } @@ -251,20 +287,42 @@ TEST(RooHist, MakeResidAndPullHistUnbinned) 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); + const int nBins = 40; + std::unique_ptr resid{makeResidHist(gauss, *data, Binning(nBins))}; + std::unique_ptr pull{makePullHist(gauss, *data, Binning(nBins))}; + EXPECT_EQ(resid->GetN(), nBins); + EXPECT_EQ(pull->GetN(), nBins); - // The correct model should give small pulls + // The pulls must be the residuals divided by the data error of the + // corresponding residual point + for (int i = 0; i < pull->GetN(); ++i) { + EXPECT_DOUBLE_EQ(pull->GetPointX(i), resid->GetPointX(i)); + const double r = resid->GetPointY(i); + const double err = r > 0 ? resid->GetErrorYlow(i) : resid->GetErrorYhigh(i); + if (err == 0.) { + EXPECT_DOUBLE_EQ(pull->GetPointY(i), 0.) << "point " << i; + continue; + } + EXPECT_NEAR(pull->GetPointY(i), r / err, 1e-9) << "point " << i; + } + + // The pulls of the true model must be standard-normal like: no outliers + // (the seed is fixed), and a chi2/ndf around one. The chi2 is a strong + // check of the per-bin expectations, since with a wrong normalization or + // bin-center evaluation it would be far off. int nChecked = 0; + double chi2 = 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; + const double p = pull->GetPointY(i); + EXPECT_LT(std::abs(p), 3.5) << "point " << i; + chi2 += p * p; ++nChecked; } EXPECT_GT(nChecked, 30); + EXPECT_GT(chi2 / nChecked, 0.5); + EXPECT_LT(chi2 / nChecked, 1.5); } /// The Range() argument of RooFit::makeResidHist()/makePullHist(): only bins @@ -289,16 +347,19 @@ TEST(RooHist, MakeResidAndPullHistRange) 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))}; + // With coordinate range + std::unique_ptr resid{makeResidHist(model.gaussAnalytic, binData, Range(rlo, rhi))}; EXPECT_EQ(resid->GetN(), 10); - expectExactResiduals(*resid, hist, 5e-3); + expectExactResiduals(*resid, hist, kTolAnalytic); + + // Same with the numerically integrated model + std::unique_ptr residNumeric{makeResidHist(model.gaussNumeric, binData, Range(rlo, rhi))}; + EXPECT_EQ(residNumeric->GetN(), 10); + expectExactResiduals(*residNumeric, hist, kTolNumeric); // With named range model.x.setRange("central", rlo, rhi); - std::unique_ptr resid2{makeResidHist(model.gauss, binData, Range("central"))}; + std::unique_ptr resid2{makeResidHist(model.gaussAnalytic, 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)); @@ -306,11 +367,11 @@ TEST(RooHist, MakeResidAndPullHistRange) } // An unknown named range must throw - EXPECT_THROW(makeResidHist(model.gauss, binData, Range("unknown")), std::invalid_argument); + EXPECT_THROW(makeResidHist(model.gaussAnalytic, 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). +/// The Range() argument of RooFit::makeResidHist() accepts multiple +/// comma-separated range names for a union of named ranges (sidebands). TEST(RooHist, MakeResidAndPullHistMultiRange) { using namespace RooFit; @@ -327,9 +388,9 @@ TEST(RooHist, MakeResidAndPullHistMultiRange) 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"))}; + std::unique_ptr resid{makeResidHist(model.gaussAnalytic, binData, Range("sb1,sb2"))}; ASSERT_EQ(resid->GetN(), 8); - expectExactResiduals(*resid, hist, 5e-3); + expectExactResiduals(*resid, hist, kTolAnalytic); } /// The documented failure modes of RooFit::makeResidHist() and diff --git a/tutorials/roofit/roofit/rf109_chi2residpull.C b/tutorials/roofit/roofit/rf109_chi2residpull.C index d548f660ea1c5..c132aefa47c06 100644 --- a/tutorials/roofit/roofit/rf109_chi2residpull.C +++ b/tutorials/roofit/roofit/rf109_chi2residpull.C @@ -65,10 +65,10 @@ void rf109_chi2residpull() // 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))}; + std::unique_ptr hresid{RooFit::makeResidHist(gauss, *data, Binning(40))}; // Construct a histogram with the pulls of the data w.r.t the model - std::unique_ptr hpull{makePullHist(gauss, *data, Binning(40))}; + std::unique_ptr hpull{RooFit::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. From 4116c368c7cb9181e787864208b9fa8f6a61736e Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Mon, 12 Oct 2026 21:50:01 +0000 Subject: [PATCH 3/5] [RF] Fix normalization and caching issues in makeResidHist()/makePullHist() Found in a review of the branch: - Composite models like RooProdPdf and RooProjectedPdf cache the integrals of their components by integration range name, and the integration strategy of a RooHistPdf (full sum vs. partial sum) is fixed when the integral is created. Redefining one scratch range per bin therefore gave wrong results for a RooHistPdf wrapped in a RooProdPdf or projected with createProjection(), which is exactly what the documentation recommends for multi-dimensional fits. Now one scratch range per used bin is defined, each only once, and the integration is done on a deep clone of the model so that no cache entries for the removed ranges are left in the user's model. - Binning(nbins, xlo, xhi) with a range narrower than the observable clamped the entries outside of the range into the edge bins (that is what RooDataHist::add() does). They are now dropped, like in RooAbsData::plotOn(). - Without Binning(), unbinned data was binned with the binning snapshot stored in the dataset instead of the current binning of the observable as documented. The default binning of the model's observable is used now. - Ranges that are not aligned with the bin edges biased the expectation, because the model integral over the partially covered bins entered the normalization while their data did not, and overlapping named ranges double counted the model integral. The expectation is now normalized to the data in the used bins (the bins fully inside the ranges), which is self-consistent for any range and also saves the normalization integrals. A range that doesn't fully contain any bin throws. The rf109 tutorial now uses the same DataError(SumW2) for the residuals and pulls as for the plotted data. --- roofit/roofitcore/inc/RooGlobalFunc.h | 26 +-- roofit/roofitcore/src/RooGlobalFunc.cxx | 192 ++++++++++-------- roofit/roofitcore/test/testRooHist.cxx | 103 ++++++++++ tutorials/roofit/roofit/rf109_chi2residpull.C | 7 +- .../roofit/roofit/rf109_chi2residpull.py | 7 +- 5 files changed, 233 insertions(+), 102 deletions(-) diff --git a/roofit/roofitcore/inc/RooGlobalFunc.h b/roofit/roofitcore/inc/RooGlobalFunc.h index 4f88926465d61..f47d2bc5b9671 100644 --- a/roofit/roofitcore/inc/RooGlobalFunc.h +++ b/roofit/roofitcore/inc/RooGlobalFunc.h @@ -490,17 +490,19 @@ RooConstVar& RooConst(double val) ; /// 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 data inside the normalization range. Unbinned input data is binned -/// internally. The supported RooCmdArgs are: +/// the bin (with RooAbsReal::createIntegral()), normalized to the data in +/// the used bins. 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")` : compute residuals only in this range. -/// The model expectation is normalized in this range. Multiple -/// comma-separated range names can be specified 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. +/// binning of the model's observable). Entries outside of the binning +/// range are dropped. Not supported for binned input data. +/// - `Range(lo, hi)` / `Range("name")` : compute residuals only for the +/// bins fully inside this range. The model expectation is normalized to +/// the data in these bins. Multiple comma-separated range names can be +/// specified for a union of named ranges (sidebands). Use it for sideband +/// fits or blinded data, so that the model is normalized to the data +/// outside of the blinded region. /// - `Normalization(double)` : additional scale factor for the expectation /// (default 1.0). The scale types of Normalization(scale, type) are not /// supported. @@ -521,9 +523,9 @@ RooConstVar& RooConst(double val) ; /// /// 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: +/// ranges or binnings, ranges that don't fully contain any bin, or when the +/// model has no support in the used bins. 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})}; diff --git a/roofit/roofitcore/src/RooGlobalFunc.cxx b/roofit/roofitcore/src/RooGlobalFunc.cxx index 4f96d1be2120d..ee44facb1baf1 100644 --- a/roofit/roofitcore/src/RooGlobalFunc.cxx +++ b/roofit/roofitcore/src/RooGlobalFunc.cxx @@ -1225,6 +1225,20 @@ makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData const &data, bool normalize throw std::invalid_argument(std::string(funcName) + "(): observable must be a RooRealVar"); } + // 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}; + + // Work on a deep clone of the model. Composite models like RooProdPdf or + // RooProjectedPdf cache the integrals of their components by integration + // range name, and the scratch ranges for the per-bin integration below + // are removed again at the end of this function. Integrating the user's + // model would leave stale cache entries behind that are reused in later + // calls with different bin boundaries, and would also grow its caches by + // one entry per bin and call. + std::unique_ptr model{static_cast(fitModel.cloneTree())}; + // 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 @@ -1235,7 +1249,7 @@ makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData const &data, bool normalize // data-side one and the model integral is analytic in that case. RooRealVar *intVar = obs; { - std::unique_ptr modelObs{fitModel.getObservables(*data.get())}; + std::unique_ptr modelObs{model->getObservables(*data.get())}; if (auto *modelVar = dynamic_cast(modelObs->find(obs->GetName()))) { intVar = modelVar; } @@ -1271,21 +1285,33 @@ makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData const &data, bool normalize 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. + // If the input data is unbinned, bin it. Without Binning() arguments, the + // current default binning of the model's observable is used (the data + // only stores a snapshot of the binning at the time of its creation). The + // binning is realized on a clone of the data's observable, so the binning + // of the user's variables is never modified. Entries outside of the + // binning range are skipped, as RooDataHist::add() would clamp them into + // the edge bins. std::unique_ptr ownedBinData; - std::unique_ptr ownedObs; + 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(); + if (!binning) { + binning = &intVar->getBinning(); } + ownedObs.reset(static_cast(obs->Clone())); + ownedObs->setBinning(*binning); std::string binnedName = std::string(data.GetName()) + "_binned"; - ownedBinData = std::make_unique(binnedName, data.GetTitle(), RooArgSet{*usedObs}, data); + ownedBinData = std::make_unique(binnedName, data.GetTitle(), RooArgSet{*ownedObs}); + const double xlo = binning->lowBound(); + const double xhi = binning->highBound(); + for (int i = 0; i < data.numEntries(); ++i) { + const RooArgSet *row = data.get(i); + const double xval = obs->getVal(); + if (xval < xlo || xval > xhi) + continue; + ownedBinData->add(*row, data.weight(), data.weightSquared()); + } binData = ownedBinData.get(); } else if (binning) { throw std::invalid_argument(std::string(funcName) + @@ -1316,11 +1342,11 @@ makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData const &data, bool normalize 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. Multiple comma-separated range names - // select a union of ranges (sidebands), by the usual RooFit convention. + // Resolve the optional Range() arguments. Only the bins fully inside one + // of the ranges are used (default: all bins), and the model expectation + // is normalized to the data in these bins. Multiple comma-separated range + // names select a union of ranges (sidebands), by the usual RooFit + // convention. std::vector> ranges; if (pc.hasProcessed("RangeWithName")) { std::string rangeNames = pc.getString("rangeName", nullptr, true); @@ -1336,55 +1362,6 @@ makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData const &data, bool normalize } } 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}; - - // A scratch range on the model observable that is redefined before each - // integration, used to integrate the model exactly over the normalization - // ranges and over each bin (same mechanism as RooBinSamplingPdf). To not - // clobber an existing range of the user's variable, the first name that - // is not in use yet is taken, and the range is removed again at the end - // of this function. - std::string scratchRange; - for (unsigned int i = 0;; ++i) { - scratchRange = "_binInt_" + std::to_string(i); - if (!intVar->hasRange(scratchRange.c_str())) - break; - } - struct ScratchRangeGuard { - RooRealVar &var; - std::string name; - ~ScratchRangeGuard() { var.removeBinning(name.c_str()); } - } scratchRangeGuard{*intVar, scratchRange}; - - // Integrates the model over the current bounds of the scratch range. A new - // integral object is created for each integration on purpose: the - // integration strategy of a RooRealIntegral (e.g. the analytical integral - // code of the integrand) is fixed at creation and may depend on the range - // bounds. RooHistPdf for instance only sums bin contents in a sub-range if - // the range was not the full range at creation time. - auto integrate = [&](double lo, double hi) { - intVar->setRange(scratchRange.c_str(), lo, hi); - std::unique_ptr integral{fitModel.createIntegral(obsSet, normSet, scratchRange.c_str())}; - return integral->getVal(); - }; - - 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 (auto const &range : ranges) { - normIntegral += integrate(range.first, range.second); - } - 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 @@ -1395,38 +1372,85 @@ makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData const &data, bool normalize } 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). + // The used bins and their data weight + std::vector usedBins; 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; - }); + const bool inRange = ranges.empty() || std::any_of(ranges.begin(), ranges.end(), [&](auto const &r) { + return binLo >= r.first - eps && binHi <= r.second + eps; + }); if (inRange) { + usedBins.push_back(iBin); nData += binData->weight(iBin); } } + if (usedBins.empty()) { + throw std::invalid_argument(std::string(funcName) + "(): no bin of the data is fully inside the range(s)"); + } - for (int iBin = 0; iBin < nBins; ++iBin) { - binData->get(iBin); + // Scratch ranges on the model observable, one per used bin, to integrate + // the model exactly over each bin (same mechanism as RooBinSamplingPdf). + // Each range is defined only once, because composite models cache the + // integrals of their components by range name, and the integration + // strategy of a RooRealIntegral (e.g. the analytical integral code of a + // RooHistPdf, which depends on whether the range is the full range) is + // fixed at creation time. To not clobber existing ranges of the user's + // variable, the first set of names not in use yet is taken, and the + // ranges are removed again at the end of this function (named ranges are + // shared between the model clone and the user's variable). + std::vector scratchRanges(usedBins.size()); + for (unsigned int k = 0;; ++k) { + bool allFree = true; + for (std::size_t j = 0; j < usedBins.size(); ++j) { + scratchRanges[j] = "_binInt_" + std::to_string(k) + "_" + std::to_string(j); + if (intVar->hasRange(scratchRanges[j].c_str())) { + allFree = false; + break; + } + } + if (allFree) + break; + } + struct ScratchRangeGuard { + RooRealVar &var; + std::vector const &names; + ~ScratchRangeGuard() + { + for (auto const &name : names) { + var.removeBinning(name.c_str()); + } + } + } scratchRangeGuard{*intVar, scratchRanges}; + + // The model integrals over the used bins + std::vector binIntegrals(usedBins.size()); + double sumBinIntegrals = 0.; + for (std::size_t j = 0; j < usedBins.size(); ++j) { + intVar->setRange(scratchRanges[j].c_str(), dataBinning.binLow(usedBins[j]), dataBinning.binHigh(usedBins[j])); + std::unique_ptr integral{model->createIntegral(obsSet, normSet, scratchRanges[j].c_str())}; + binIntegrals[j] = integral->getVal(); + sumBinIntegrals += binIntegrals[j]; + } + if (sumBinIntegrals == 0.) { + throw std::invalid_argument(std::string(funcName) + "(): model has no support in the used bins"); + } - 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; + // The expectation is normalized to the data in the used bins + const double expectedNorm = scaleFactor * nData / sumBinIntegrals; - const double fraction = integrate(binLo, binHi) / normIntegral; - const double expected = scaleFactor * nData * fraction; + RooAbsData::ErrorType errType = + etype == RooAbsData::Auto ? (binData->isNonPoissonWeighted() ? RooAbsData::SumW2 : RooAbsData::Poisson) : etype; + + for (std::size_t j = 0; j < usedBins.size(); ++j) { + const int iBin = usedBins[j]; + binData->get(iBin); + const double expected = expectedNorm * binIntegrals[j]; const double actual = binData->weight(iBin); - double eylo, eyhi; + double eylo; + double eyhi; binData->weightError(eylo, eyhi, errType); double resid = actual - expected; diff --git a/roofit/roofitcore/test/testRooHist.cxx b/roofit/roofitcore/test/testRooHist.cxx index e647de05bb36f..86770005c3038 100644 --- a/roofit/roofitcore/test/testRooHist.cxx +++ b/roofit/roofitcore/test/testRooHist.cxx @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -323,6 +324,87 @@ TEST(RooHist, MakeResidAndPullHistUnbinned) EXPECT_GT(nChecked, 30); EXPECT_GT(chi2 / nChecked, 0.5); EXPECT_LT(chi2 / nChecked, 1.5); + + // Without Binning(), the current binning of the observable is used, not + // the binning stored in the dataset at the time of its creation + x.setBins(25); + std::unique_ptr residDefault{makeResidHist(gauss, *data)}; + EXPECT_EQ(residDefault->GetN(), 25); +} + +/// With Binning(nbins, xlo, xhi) narrower than the observable range, the +/// entries outside of the binning range must be dropped (and not clamped +/// into the edge bins, as RooDataHist::add() does). +TEST(RooHist, MakeResidAndPullHistBinningRange) +{ + 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)}; + + const int nBins = 20; + const double xlo = -2.; + const double xhi = 2.; + + // With Normalization(0.), the model expectation is zero and the residuals + // are the binned data contents + std::unique_ptr resid{makeResidHist(gauss, *data, Binning(nBins, xlo, xhi), Normalization(0.))}; + ASSERT_EQ(resid->GetN(), nBins); + + TH1D ref("ref", "ref", nBins, xlo, xhi); + for (int i = 0; i < data->numEntries(); ++i) { + ref.Fill(static_cast(*data->get(i)->find("x")).getVal()); + } + EXPECT_LT(ref.Integral(), 0.6 * data->numEntries()); // a large fraction of the data is outside the range + for (int i = 0; i < nBins; ++i) { + EXPECT_DOUBLE_EQ(resid->GetPointX(i), ref.GetBinCenter(i + 1)) << "bin " << i; + EXPECT_DOUBLE_EQ(resid->GetPointY(i), ref.GetBinContent(i + 1)) << "bin " << i; + } +} + +/// Composite models cache the integrals of their components by range name, +/// and the integration strategy of a RooHistPdf (full sum vs. partial sum) +/// is fixed when the integral is created. The per-bin integration must not +/// be affected by this: a RooHistPdf wrapped in a RooProdPdf or projected +/// with createProjection() must give the same residuals as the bare +/// RooHistPdf. +TEST(RooHist, MakeResidAndPullHistWrappedModel) +{ + using namespace RooFit; + + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING}; + + GaussModel model; + RooRealVar y("y", "y", -1., 1.); + RooGenericPdf gy("gy", "gy", "1 + y * y", RooArgSet{y}); + RooProdPdf prod("prod", "prod", RooArgList{model.gaussAnalytic, gy}); + std::unique_ptr proj{prod.createProjection(RooArgSet{y})}; + + TH1D hist("hist", "hist", kNBins, kGaussX0, kGaussX1); + fillExactGauss(hist, 1e6, {{kGaussX0, kGaussX1}}); + RooDataHist binData("binData", "binData", RooArgSet(model.x), &hist); + + for (RooAbsPdf *pdf : {static_cast(&prod), proj.get()}) { + std::unique_ptr resid{makeResidHist(*pdf, binData)}; + ASSERT_EQ(resid->GetN(), kNBins) << pdf->GetName(); + expectExactResiduals(*resid, hist, kTolAnalytic); + // Also with a sub-range, and in a second call with the same model + std::unique_ptr residRange{makeResidHist(*pdf, binData, Range(119.6, 123.6))}; + ASSERT_EQ(residRange->GetN(), 10) << pdf->GetName(); + expectExactResiduals(*residRange, hist, kTolAnalytic); + } + + // The user's variable must not keep any scratch range + for (auto const &name : model.x.getBinningNames()) { + EXPECT_EQ(name.find("_binInt_"), std::string::npos) << "leftover range " << name; + } } /// The Range() argument of RooFit::makeResidHist()/makePullHist(): only bins @@ -357,6 +439,14 @@ TEST(RooHist, MakeResidAndPullHistRange) EXPECT_EQ(residNumeric->GetN(), 10); expectExactResiduals(*residNumeric, hist, kTolNumeric); + // With a range that is not aligned with the bin edges, only the bins + // fully inside the range are used, and the expectation is normalized to + // the data in these bins (so the residuals of the exact data are still + // zero) + std::unique_ptr residUnaligned{makeResidHist(model.gaussAnalytic, binData, Range(rlo - 0.1, rhi + 0.1))}; + EXPECT_EQ(residUnaligned->GetN(), 10); + expectExactResiduals(*residUnaligned, hist, kTolAnalytic); + // With named range model.x.setRange("central", rlo, rhi); std::unique_ptr resid2{makeResidHist(model.gaussAnalytic, binData, Range("central"))}; @@ -391,6 +481,16 @@ TEST(RooHist, MakeResidAndPullHistMultiRange) std::unique_ptr resid{makeResidHist(model.gaussAnalytic, binData, Range("sb1,sb2"))}; ASSERT_EQ(resid->GetN(), 8); expectExactResiduals(*resid, hist, kTolAnalytic); + + // Overlapping ranges must not double count the model expectation: bins + // 13..18 with the overlapping ranges sb1 and [116.0, 117.6] + model.x.setRange("sbo", 116.0, 117.6); + TH1D hist2("hist2", "hist2", kNBins, kGaussX0, kGaussX1); + fillExactGauss(hist2, 1e6, {{115.2, 117.6}}, 3.e6); + RooDataHist binData2("binData2", "binData2", RooArgSet(model.x), &hist2); + std::unique_ptr residOverlap{makeResidHist(model.gaussAnalytic, binData2, Range("sb1,sbo"))}; + ASSERT_EQ(residOverlap->GetN(), 6); + expectExactResiduals(*residOverlap, hist2, kTolAnalytic); } /// The documented failure modes of RooFit::makeResidHist() and @@ -428,6 +528,9 @@ TEST(RooHist, MakeResidAndPullHistErrors) // Unknown named ranges EXPECT_THROW(makeResidHist(gauss, *data, Range("unknown")), std::invalid_argument); + // A range that doesn't fully contain any bin (bins are 0.5 wide) + EXPECT_THROW(makeResidHist(gauss, *data, Range(0.1, 0.4)), std::invalid_argument); + // Unsupported Normalization scale types EXPECT_THROW(makeResidHist(gauss, *data, Normalization(0.5, RooAbsReal::NumEvent)), std::invalid_argument); diff --git a/tutorials/roofit/roofit/rf109_chi2residpull.C b/tutorials/roofit/roofit/rf109_chi2residpull.C index c132aefa47c06..66a68e6c4ed19 100644 --- a/tutorials/roofit/roofit/rf109_chi2residpull.C +++ b/tutorials/roofit/roofit/rf109_chi2residpull.C @@ -64,11 +64,12 @@ void rf109_chi2residpull() // that curve-based residuals show for sharply peaked pdfs. // 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{RooFit::makeResidHist(gauss, *data, Binning(40))}; + // The Binning() argument selects how the unbinned dataset is binned, and + // DataError() selects the same error model as in the data plot above. + std::unique_ptr hresid{RooFit::makeResidHist(gauss, *data, Binning(40), DataError(RooAbsData::SumW2))}; // Construct a histogram with the pulls of the data w.r.t the model - std::unique_ptr hpull{RooFit::makePullHist(gauss, *data, Binning(40))}; + std::unique_ptr hpull{RooFit::makePullHist(gauss, *data, Binning(40), DataError(RooAbsData::SumW2))}; // 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. diff --git a/tutorials/roofit/roofit/rf109_chi2residpull.py b/tutorials/roofit/roofit/rf109_chi2residpull.py index aa4b262f05125..f478241a2c525 100644 --- a/tutorials/roofit/roofit/rf109_chi2residpull.py +++ b/tutorials/roofit/roofit/rf109_chi2residpull.py @@ -53,13 +53,14 @@ # 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. +# argument selects how the unbinned dataset is binned, and DataError() +# selects the same error model as in the data plot above. # Construct a histogram with the residuals of the data w.r.t. the model -hresid = ROOT.RooFit.makeResidHist(gauss, data, ROOT.RooFit.Binning(40)) +hresid = ROOT.RooFit.makeResidHist(gauss, data, ROOT.RooFit.Binning(40), ROOT.RooFit.DataError(ROOT.RooAbsData.SumW2)) # Construct a histogram with the pulls of the data w.r.t the model -hpull = ROOT.RooFit.makePullHist(gauss, data, ROOT.RooFit.Binning(40)) +hpull = ROOT.RooFit.makePullHist(gauss, data, ROOT.RooFit.Binning(40), ROOT.RooFit.DataError(ROOT.RooAbsData.SumW2)) # Create a frame to draw the residual distribution and add the # distribution to the frame. Note that addPlotable() transfers ownership From d412e353ee13917d64947c220036ed6d51f1a58f Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Mon, 12 Oct 2026 22:01:31 +0000 Subject: [PATCH 4/5] [RF] Support the Normalization() scale types in makeResidHist()/makePullHist() The Normalization(scale, type) argument now follows the RooAbsReal::ScaleType semantics of RooAbsPdf::plotOn() instead of rejecting all scale types but Relative: - Relative (default): the expectation is normalized to the data in the used bins, as before. - RelativeExpected: the expectation is the expected number of events of the extended model. This allows to compare the data with the fitted yield of an extended fit, e.g. for blinded or sideband fits, and matches the default normalization of RooAbsPdf::plotOn() for extended pdfs. - NumEvent: the expectation is normalized to the given number of events in the full range of the observable. - Raw: the expectation is the integral of the unnormalized model over each bin. --- README/ReleaseNotes/v642/index.md | 2 +- roofit/roofitcore/inc/RooGlobalFunc.h | 27 +++++++--- roofit/roofitcore/src/RooGlobalFunc.cxx | 53 ++++++++++++++----- roofit/roofitcore/test/testRooHist.cxx | 68 ++++++++++++++++++++++++- 4 files changed, 126 insertions(+), 24 deletions(-) diff --git a/README/ReleaseNotes/v642/index.md b/README/ReleaseNotes/v642/index.md index 89d310a0b6f05..bcee7d61f5743 100644 --- a/README/ReleaseNotes/v642/index.md +++ b/README/ReleaseNotes/v642/index.md @@ -205,7 +205,7 @@ the cut instead of being selected based on `sqrt(abs(x))`. 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. +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()` with the scale types of `RooAbsPdf::plotOn()` (e.g. `RooAbsReal::RelativeExpected` to compare with the fitted yield of an extended pdf), and `RooFit::DataError()` to select the point error model. ```cpp RooDataHist binData("binData", "binData", x, *unbinnedData); diff --git a/roofit/roofitcore/inc/RooGlobalFunc.h b/roofit/roofitcore/inc/RooGlobalFunc.h index f47d2bc5b9671..1f74845890bda 100644 --- a/roofit/roofitcore/inc/RooGlobalFunc.h +++ b/roofit/roofitcore/inc/RooGlobalFunc.h @@ -490,9 +490,9 @@ RooConstVar& RooConst(double val) ; /// 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 data in -/// the used bins. Unbinned input data is binned internally. The supported -/// RooCmdArgs are: +/// the bin (with RooAbsReal::createIntegral()), by default normalized to the +/// data in the used bins (see `Normalization()`). 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 model's observable). Entries outside of the binning @@ -503,9 +503,19 @@ RooConstVar& RooConst(double val) ; /// specified for a union of named ranges (sidebands). Use it for sideband /// fits or blinded data, so that the model is normalized to the data /// outside of the blinded region. -/// - `Normalization(double)` : additional scale factor for the expectation -/// (default 1.0). The scale types of Normalization(scale, type) are not -/// supported. +/// - `Normalization(scale, type)` : scale factor and scale type of the +/// expectation, with the RooAbsReal::ScaleType semantics of +/// RooAbsPdf::plotOn() (default: `Normalization(1.0, RooAbsReal::Relative)`): +/// - `Relative` : the expectation is normalized to the data in the used +/// bins, and multiplied by `scale` +/// - `RelativeExpected` : the expectation is the expected number of +/// events of the extended model, multiplied by `scale`. Use it to +/// compare the data with the fitted yield, e.g. for blinded or +/// sideband fits. Throws if the model is not an extendable pdf. +/// - `NumEvent` : the expectation is normalized to `scale` events in +/// the full range of the observable +/// - `Raw` : the expectation is the integral of the unnormalized model +/// over the bin, multiplied by `scale` /// - `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) @@ -515,7 +525,10 @@ RooConstVar& RooConst(double val) ; /// 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. +/// match the displayed data and curve. In particular, RooAbsPdf::plotOn() +/// normalizes an extended pdf to its expected number of events unless a +/// range is given, which corresponds to `Normalization(1.0, +/// RooAbsReal::RelativeExpected)` here. /// /// \note If the model depends on observables beyond the data's observable /// (e.g. a conditional pdf), they are evaluated at their current values. diff --git a/roofit/roofitcore/src/RooGlobalFunc.cxx b/roofit/roofitcore/src/RooGlobalFunc.cxx index ee44facb1baf1..88ed7da4a8bf7 100644 --- a/roofit/roofitcore/src/RooGlobalFunc.cxx +++ b/roofit/roofitcore/src/RooGlobalFunc.cxx @@ -1172,7 +1172,7 @@ makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData const &data, bool normalize RooCmdConfig pc(funcName); pc.defineDouble("normScale", "Normalization", 0, 1.0); - pc.defineInt("normType", "Normalization", 0, 0); + pc.defineInt("normType", "Normalization", 0, static_cast(RooAbsReal::Relative)); pc.defineInt("etype", "DataError", 0, static_cast(RooAbsData::Auto)); pc.defineString("name", "Name", 0, ""); pc.defineString("title", "Title", 0, ""); @@ -1195,12 +1195,17 @@ makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData const &data, bool normalize 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 stype = static_cast(pc.getInt("normType")); + if (stype != RooAbsReal::Relative && stype != RooAbsReal::RelativeExpected && stype != RooAbsReal::NumEvent && + stype != RooAbsReal::Raw) { + throw std::invalid_argument(std::string(funcName) + "(): unsupported scale type in Normalization()"); + } + if (stype == RooAbsReal::RelativeExpected) { + auto *pdf = dynamic_cast(&fitModel); + if (!pdf || !pdf->canBeExtended()) { + throw std::invalid_argument(std::string(funcName) + + "(): the RelativeExpected scale type requires an extendable pdf"); + } } const auto etype = static_cast(pc.getInt("etype")); @@ -1424,21 +1429,41 @@ makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData const &data, bool normalize } } scratchRangeGuard{*intVar, scratchRanges}; - // The model integrals over the used bins + // The model integrals over the used bins. For the Raw scale type, the + // integrals of the unnormalized model are used. std::vector binIntegrals(usedBins.size()); double sumBinIntegrals = 0.; for (std::size_t j = 0; j < usedBins.size(); ++j) { intVar->setRange(scratchRanges[j].c_str(), dataBinning.binLow(usedBins[j]), dataBinning.binHigh(usedBins[j])); - std::unique_ptr integral{model->createIntegral(obsSet, normSet, scratchRanges[j].c_str())}; + std::unique_ptr integral{stype == RooAbsReal::Raw + ? model->createIntegral(obsSet, scratchRanges[j].c_str()) + : model->createIntegral(obsSet, normSet, scratchRanges[j].c_str())}; binIntegrals[j] = integral->getVal(); sumBinIntegrals += binIntegrals[j]; } - if (sumBinIntegrals == 0.) { - throw std::invalid_argument(std::string(funcName) + "(): model has no support in the used bins"); - } - // The expectation is normalized to the data in the used bins - const double expectedNorm = scaleFactor * nData / sumBinIntegrals; + // The normalization of the expectation, following the scale types of + // RooAbsPdf::plotOn() + double expectedNorm = scaleFactor; + if (stype == RooAbsReal::Relative) { + // Normalized to the data in the used bins + if (sumBinIntegrals == 0.) { + throw std::invalid_argument(std::string(funcName) + "(): model has no support in the used bins"); + } + expectedNorm *= nData / sumBinIntegrals; + } else if (stype == RooAbsReal::NumEvent || stype == RooAbsReal::RelativeExpected) { + // Normalized to the given or the expected number of events in the + // full range of the observable + std::unique_ptr fullIntegral{model->createIntegral(obsSet, normSet)}; + const double full = fullIntegral->getVal(); + if (full == 0.) { + throw std::invalid_argument(std::string(funcName) + "(): model has no support in the observable range"); + } + expectedNorm /= full; + if (stype == RooAbsReal::RelativeExpected) { + expectedNorm *= static_cast(*model).expectedEvents(normSet); + } + } RooAbsData::ErrorType errType = etype == RooAbsData::Auto ? (binData->isNonPoissonWeighted() ? RooAbsData::SumW2 : RooAbsData::Poisson) : etype; diff --git a/roofit/roofitcore/test/testRooHist.cxx b/roofit/roofitcore/test/testRooHist.cxx index 86770005c3038..3a2645c1a3e99 100644 --- a/roofit/roofitcore/test/testRooHist.cxx +++ b/roofit/roofitcore/test/testRooHist.cxx @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -493,6 +494,69 @@ TEST(RooHist, MakeResidAndPullHistMultiRange) expectExactResiduals(*residOverlap, hist2, kTolAnalytic); } +/// The scale types of the Normalization() argument of +/// RooFit::makeResidHist()/makePullHist(), following RooAbsPdf::plotOn(). +TEST(RooHist, MakeResidAndPullHistNormalization) +{ + using namespace RooFit; + + RooHelpers::LocalChangeMsgLevel chmsglvl{RooFit::WARNING}; + + GaussModel model; + const double nEvents = 1e6; + TH1D hist("hist", "hist", kNBins, kGaussX0, kGaussX1); + fillExactGauss(hist, nEvents, {{kGaussX0, kGaussX1}}); + RooDataHist binData("binData", "binData", RooArgSet(model.x), &hist); + + // Checks that the residuals are (-factor) times the exact bin contents, + // i.e. the expectation is (1 + factor) times the data + auto expectScaledResiduals = [&](RooHist const &resid, double factor, double tol) { + 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_NEAR(resid.GetPointY(i), -factor * content, tol * std::max(content, 1e-6)) << "bin " << iBin; + } + }; + + // An extended pdf with a yield of 1.5 times the data + RooRealVar nExpected("nExpected", "nExpected", 1.5 * nEvents); + RooExtendPdf ext("ext", "ext", model.gaussAnalytic, nExpected); + + // The default Relative scale type normalizes to the data also for an + // extended pdf + std::unique_ptr resid{makeResidHist(ext, binData)}; + ASSERT_EQ(resid->GetN(), kNBins); + expectScaledResiduals(*resid, 0., kTolAnalytic); + + // RelativeExpected: the expectation is the expected number of events + std::unique_ptr residExp{makeResidHist(ext, binData, Normalization(1.0, RooAbsReal::RelativeExpected))}; + ASSERT_EQ(residExp->GetN(), kNBins); + expectScaledResiduals(*residExp, 0.5, kTolAnalytic); + + // Also in a sub-range, the expectation is not renormalized to the data + std::unique_ptr residExpRange{ + makeResidHist(ext, binData, Normalization(2.0, RooAbsReal::RelativeExpected), Range(119.6, 123.6))}; + ASSERT_EQ(residExpRange->GetN(), 10); + expectScaledResiduals(*residExpRange, 2., kTolAnalytic); + + // NumEvent: the expectation is normalized to the given number of events + std::unique_ptr residNum{ + makeResidHist(model.gaussAnalytic, binData, Normalization(2. * nEvents, RooAbsReal::NumEvent))}; + ASSERT_EQ(residNum->GetN(), kNBins); + expectScaledResiduals(*residNum, 1., kTolAnalytic); + + // Raw: the integral of the unnormalized model. The RooGenericPdf + // exp(-0.5 u^2) integrates to sigma * sqrt(2 pi) times the Gaussian + // probability in each bin. + const double rawScale = nEvents / (kGaussSigma * std::sqrt(2. * M_PI)); + std::unique_ptr residRaw{ + makeResidHist(model.gaussNumeric, binData, Normalization(rawScale, RooAbsReal::Raw))}; + ASSERT_EQ(residRaw->GetN(), kNBins); + expectScaledResiduals(*residRaw, 0., kTolNumeric); +} + /// The documented failure modes of RooFit::makeResidHist() and /// RooFit::makePullHist() that throw exceptions. TEST(RooHist, MakeResidAndPullHistErrors) @@ -531,8 +595,8 @@ TEST(RooHist, MakeResidAndPullHistErrors) // A range that doesn't fully contain any bin (bins are 0.5 wide) EXPECT_THROW(makeResidHist(gauss, *data, Range(0.1, 0.4)), std::invalid_argument); - // Unsupported Normalization scale types - EXPECT_THROW(makeResidHist(gauss, *data, Normalization(0.5, RooAbsReal::NumEvent)), std::invalid_argument); + // The RelativeExpected scale type requires an extendable pdf + EXPECT_THROW(makeResidHist(gauss, *data, Normalization(1.0, RooAbsReal::RelativeExpected)), std::invalid_argument); // Unsupported DataError types EXPECT_THROW(makeResidHist(gauss, *data, DataError(99)), std::invalid_argument); From 2d8a43e271eba7fb443efe2f772d4180b9071e01 Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Tue, 13 Oct 2026 19:18:40 +0000 Subject: [PATCH 5/5] [RF] Fix testRooHist on Windows with TMath::TwoPi() instead of M_PI MSVC's `` only defines `M_PI` when `_USE_MATH_DEFINES` is set, so the Windows CI jobs failed to compile the RooHist tests. --- roofit/roofitcore/test/testRooHist.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/roofit/roofitcore/test/testRooHist.cxx b/roofit/roofitcore/test/testRooHist.cxx index 3a2645c1a3e99..8eaed5fa46c24 100644 --- a/roofit/roofitcore/test/testRooHist.cxx +++ b/roofit/roofitcore/test/testRooHist.cxx @@ -18,6 +18,7 @@ #include #include +#include #include @@ -550,7 +551,7 @@ TEST(RooHist, MakeResidAndPullHistNormalization) // Raw: the integral of the unnormalized model. The RooGenericPdf // exp(-0.5 u^2) integrates to sigma * sqrt(2 pi) times the Gaussian // probability in each bin. - const double rawScale = nEvents / (kGaussSigma * std::sqrt(2. * M_PI)); + const double rawScale = nEvents / (kGaussSigma * std::sqrt(TMath::TwoPi())); std::unique_ptr residRaw{ makeResidHist(model.gaussNumeric, binData, Normalization(rawScale, RooAbsReal::Raw))}; ASSERT_EQ(residRaw->GetN(), kNBins);