Repository navigation
[RF] Add RooFit::makeResidHist() and makePullHist() with bin integration - #23357
guitargeek wants to merge 5 commits into
Conversation
Test Results 23 files 23 suites 3d 19h 6m 9s ⏱️ For more details on these failures, see this check. Results for commit 2d8a43e. ♻️ This comment has been updated with latest results. |
hageboeck
left a comment
There was a problem hiding this comment.
Hello, below a few things I was wondering about on the first pass.
| /// bin, instead of interpolating the plotted curve. | ||
| RooHist* RooPlot::residHist(const char* histname, const char* curvename, bool normalize, bool useAverage) const | ||
| { | ||
| static std::atomic<bool> warned{false}; |
There was a problem hiding this comment.
Is there a particular reason for using a static atomic here?
Although I'm sure it works, I doubt the plotting and integration code is thread safe in the way that it should be used from multiple threads, is it? So this could possibly be done without the atomic and the (small) risk that this warning is seen more than once.
| 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<RooHist> hresid{makeResidHist(gauss, *data, Binning(40))}; |
There was a problem hiding this comment.
One could use RooFit::makeResidHist just to clarify where it came from. Similar below.
What do you think?
| /// the bin (with RooAbsReal::createIntegral()), normalized to the weight of | ||
| /// the data inside the normalization range. Unbinned input data is binned |
There was a problem hiding this comment.
That it's a weight is an implementation detail. Most users probably think of data as counts.
| /// the bin (with RooAbsReal::createIntegral()), normalized to the weight of | |
| /// the data inside the normalization range. Unbinned input data is binned | |
| /// the bin (with RooAbsReal::createIntegral()), normalized to | |
| /// the data inside the normalization range. Unbinned input data is binned |
| /// - `Range(lo, hi)` / `Range("name")` : the range that the model | ||
| /// expectation is normalized in. Only bins inside this range get points |
There was a problem hiding this comment.
| /// - `Range(lo, hi)` / `Range("name")` : the range that the model | |
| /// expectation is normalized in. Only bins inside this range get points | |
| /// - `Range(lo, hi)` / `Range("name")` : compute residuals only in this range. | |
| /// The model expectation is normalized in this range. |
| /// \note For internal use of RooAbsReal::createIntegral(), a scratch named | ||
| /// range is left on the model's observable after the call. |
There was a problem hiding this comment.
scratch named range?
A range with the name "xxx" is declared?
| std::unique_ptr<RooHist> pull{makePullHist(model.gauss, binData)}; | ||
| ASSERT_EQ(resid->GetN(), kNBins); | ||
| ASSERT_EQ(pull->GetN(), kNBins); | ||
| expectExactResiduals(*resid, hist, 1e-3); |
There was a problem hiding this comment.
Why only 1.E-3? Does the algo not manage better?
There was a problem hiding this comment.
Well, on the other hand, for a residual plot its accurate enough ...
| std::unique_ptr<RooHist> resid{makeResidHist(gauss, *data, Binning(40))}; | ||
| std::unique_ptr<RooHist> pull{makePullHist(gauss, *data, Binning(40))}; | ||
| EXPECT_EQ(resid->GetN(), 40); | ||
| EXPECT_EQ(pull->GetN(), 40); |
There was a problem hiding this comment.
On the future, consider declaring the number of bins somewhere, so later changes change all instances.
| for (int i = 0; i < pull->GetN(); ++i) { | ||
| if (pull->GetErrorYhigh(i) == 0.) | ||
| continue; | ||
| EXPECT_LT(std::abs(pull->GetPointY(i)), 5.) << "point " << i; |
| // 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<unsigned long> nScratchRanges{0}; |
There was a problem hiding this comment.
You expect to call this from multiple threads?
There was a problem hiding this comment.
As a slightly slower algo that has the advantage of using less global memory and working without static init & destroy (remember that these statics need to be allocated and finalised under lock by the C++ runtime):
- Put a local counter
- Craft a range name using that number
- Check if that range exists for the var. If not, just that, otherwise increment counter and repeat.
| continue; | ||
|
|
||
| intVar->setRange(rangeNameBin.c_str(), binLo, binHi); | ||
| const double fraction = fitModel.createIntegral(obsSet, normSet, rangeNameBin.c_str())->getVal() / normIntegral; |
There was a problem hiding this comment.
I don't remember this, but can the integral be reused?
There was a problem hiding this comment.
Thanks for catching this! It is at least a memory leak, because createIntegral returns a raw pointer that needs to be deleted...
0c0dfe6 to
cd6d6bc
Compare
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
- 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].
…Hist() 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.
…ullHist() 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.
MSVC's `<cmath>` only defines `M_PI` when `_USE_MATH_DEFINES` is set, so the Windows CI jobs failed to compile the RooHist tests.
cd6d6bc to
2d8a43e
Compare
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.
This direction was also discussed with @alread2223 over email.
🤖 Done with the help of AI