diff --git a/README/ReleaseNotes/v642/index.md b/README/ReleaseNotes/v642/index.md index a233715e4835a..bcee7d61f5743 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()` 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); +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..1f74845890bda 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 ; @@ -468,6 +470,95 @@ RooCmdArg MultiArg(const RooCmdArg& arg1, const RooCmdArg& arg2, RooConstVar& RooConst(double val) ; // End group CmdArgs: +/** + * @} + */ + +/** + * \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()), 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 +/// 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(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) +/// - `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. 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. +/// There is no exact equivalent for a projection with ProjWData(). +/// +/// 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, 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})}; +/// 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())); +/// ~~~ +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. +template +std::unique_ptr makePullHist(RooAbsReal &fitModel, RooAbsData const &data, CmdArgs_t const &...cmdArgs); + /** * @} */ @@ -524,8 +615,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 05392be18a73d..7e6bf15ffc33c 100644 --- a/roofit/roofitcore/inc/RooHist.h +++ b/roofit/roofitcore/inc/RooHist.h @@ -75,7 +75,14 @@ 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; + /// 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/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..88ed7da4a8bf7 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,9 @@ #include #include +#include +#include +#include namespace RooFit { @@ -1157,6 +1164,339 @@ RooConstVar &RooConst(double val) namespace Detail { +/// Implementation of RooFit::makeResidHist() and RooFit::makePullHist(). +std::unique_ptr +makeResidOrPullHist(RooAbsReal &fitModel, RooAbsData const &data, bool normalize, RooLinkedList const &cmdList) +{ + const char *funcName = normalize ? "RooFit::makePullHist" : "RooFit::makeResidHist"; + + RooCmdConfig pc(funcName); + pc.defineDouble("normScale", "Normalization", 0, 1.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, ""); + 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(cmdList); + 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"); + } + + 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")); + 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 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 + // 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{model->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. 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; + RooDataHist const *binData = dynamic_cast(&data); + if (!binData) { + 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{*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) + + "(): 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 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); + 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")); + } + + // 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 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 = 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)"); + } + + // 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. 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{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]; + } + + // 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; + + 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; + double 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; +} + 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 9d0cf19e84f1a..a31671ff167b2 100644 --- a/roofit/roofitcore/src/RooPlot.cxx +++ b/roofit/roofitcore/src/RooPlot.cxx @@ -1064,8 +1064,23 @@ 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 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 " + "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..8eaed5fa46c24 100644 --- a/roofit/roofitcore/test/testRooHist.cxx +++ b/roofit/roofitcore/test/testRooHist.cxx @@ -4,19 +4,30 @@ #include #include #include +#include +#include +#include #include #include #include +#include #include #include #include #include #include +#include +#include #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. @@ -120,3 +131,478 @@ 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); +} + +/// 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); + } +} + +/// The exact Gaussian as a RooDataHist with kNBins bins, normalized to unity. +RooDataHist makeExactGaussDataHist(RooRealVar &x) +{ + TH1D hist("modelHist", "modelHist", kNBins, kGaussX0, kGaussX1); + fillExactGauss(hist, 1., {{kGaussX0, kGaussX1}}); + return RooDataHist{"modelDataHist", "modelDataHist", RooArgSet(x), &hist}; +} + +/// The Gaussian as a model plus observable; members are in construction order. +/// The model exists in two variants: a RooHistPdf with the exact per-bin +/// Gaussian integrals as bin contents, whose integrals over bins and ranges +/// are analytic (partial sums of the histogram), and a RooGenericPdf that is +/// integrated numerically. +struct GaussModel { + RooRealVar x{"x", "x", kGaussX0, kGaussX1}; + RooRealVar mean{"mean", "mean", kGaussMean}; + RooRealVar sigma{"sigma", "sigma", kGaussSigma}; + RooDataHist modelDataHist{makeExactGaussDataHist(x)}; + RooHistPdf gaussAnalytic{"gaussAnalytic", "gaussAnalytic", RooArgSet(x), modelDataHist}; + RooGenericPdf gaussNumeric{"gaussNumeric", "gaussNumeric", "exp(-0.5 * pow((x - mean) / sigma, 2))", + RooArgSet{x, mean, sigma}}; + + GaussModel() { x.setBins(kNBins); } +}; + +/// Tolerance (relative to the bin content) for the residuals of the exact +/// model. With the analytic integrals of the RooHistPdf, the comparison is +/// limited only by floating-point precision. +constexpr double kTolAnalytic = 1e-10; +/// For the numerically integrated RooGenericPdf, the relative precision of +/// the per-bin integrals is limited to ~1e-3 in the far tails. +constexpr double kTolNumeric = 1e-3; + +/// 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 for the nearly empty bins 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 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.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, 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; + } +} + +/// 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 + 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 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; + 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); + + // 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 +/// 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 + std::unique_ptr resid{makeResidHist(model.gaussAnalytic, binData, Range(rlo, rhi))}; + EXPECT_EQ(resid->GetN(), 10); + 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 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"))}; + 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.gaussAnalytic, binData, Range("unknown")), std::invalid_argument); +} + +/// 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; + + 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.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 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(TMath::TwoPi())); + 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) +{ + 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); + + // 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); + + // 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); + + // 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..66a68e6c4ed19 100644 --- a/tutorials/roofit/roofit/rf109_chi2residpull.C +++ b/tutorials/roofit/roofit/rf109_chi2residpull.C @@ -59,19 +59,26 @@ 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, 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))}; - // 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{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. 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..f478241a2c525 100644 --- a/tutorials/roofit/roofit/rf109_chi2residpull.py +++ b/tutorials/roofit/roofit/rf109_chi2residpull.py @@ -50,21 +50,30 @@ # 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, and DataError() +# selects the same error model as in the data plot above. -# 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), 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), ROOT.RooFit.DataError(ROOT.RooAbsData.SumW2)) # 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)