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