From c01e078bac830d98ca82cd409c96a71e5da89498 Mon Sep 17 00:00:00 2001 From: Timothy Gray Date: Mon, 20 Jul 2026 19:51:52 +0000 Subject: [PATCH] [hist][math] Allow multiple ranges in HFitInterface::FillData ROOT::Fit::FillData only honoured the first range interval of each coordinate and warned that "support only one range interval" for any further one. Fitting a single function over several disjoint regions (e.g. a background with the peak region excluded) was therefore not possible. Collect the bins by looping over every combination of the per-coordinate ranges. A coordinate without any range contributes a single iteration spanning the full axis, so the behaviour without ranges is unchanged. DataRange::CleanRangeSet() now also merges partially overlapping ranges instead of only dropping fully contained ones, so that AddRange() keeps the range set a collection of disjoint intervals. It takes xmin and xmax by reference to widen them to the merged interval. Adds unit tests for the overlapping-range handling in DataRange and for the range handling of FillData itself. --- hist/hist/src/HFitInterface.cxx | 75 +++++++------ hist/hist/test/CMakeLists.txt | 1 + hist/hist/test/test_FillDataRange.cxx | 150 ++++++++++++++++++++++++++ math/mathcore/inc/Fit/DataRange.h | 7 +- math/mathcore/src/DataRange.cxx | 24 ++++- math/mathcore/test/CMakeLists.txt | 2 + math/mathcore/test/test_Range.cxx | 57 ++++++++++ 7 files changed, 279 insertions(+), 37 deletions(-) create mode 100644 hist/hist/test/test_FillDataRange.cxx create mode 100644 math/mathcore/test/test_Range.cxx diff --git a/hist/hist/src/HFitInterface.cxx b/hist/hist/src/HFitInterface.cxx index aa45fc95e71bf..6b2bc269f5c94 100644 --- a/hist/hist/src/HFitInterface.cxx +++ b/hist/hist/src/HFitInterface.cxx @@ -136,36 +136,6 @@ void FillData(BinData & dv, const TH1 * hfit, TF1 * func) // get the range (add the function range ??) // to check if inclusion/exclusion at end/point const DataRange & range = dv.Range(); - if (range.Size(0) != 0) { - HFitInterface::ExamineRange( hfit->GetXaxis(), range(0), hxfirst, hxlast); - if (range.Size(0) > 1 ) { - Warning("ROOT::Fit::FillData","support only one range interval for X coordinate"); - } - } - - if (hfit->GetDimension() > 1 && range.Size(1) != 0) { - HFitInterface::ExamineRange( hfit->GetYaxis(), range(1), hyfirst, hylast); - if (range.Size(1) > 1 ) - Warning("ROOT::Fit::FillData","support only one range interval for Y coordinate"); - } - - if (hfit->GetDimension() > 2 && range.Size(2) != 0) { - HFitInterface::ExamineRange( hfit->GetZaxis(), range(2), hzfirst, hzlast); - if (range.Size(2) > 1 ) - Warning("ROOT::Fit::FillData","support only one range interval for Z coordinate"); - } - - - int n = (hxlast-hxfirst+1)*(hylast-hyfirst+1)*(hzlast-hzfirst+1); - -#ifdef DEBUG - std::cout << "THFitInterface: ifirst = " << hxfirst << " ilast = " << hxlast - << " total bins " << n - << std::endl; -#endif - - // reserve n for more efficient usage - //dv.Data().reserve(n); int hdim = hfit->GetDimension(); int ndim = hdim; @@ -175,7 +145,17 @@ void FillData(BinData & dv, const TH1 * hfit, TF1 * func) assert( ndim > 0 ); //typedef BinPoint::CoordData CoordData; //CoordData x = CoordData( hfit->GetDimension() ); - dv.Initialize(n,ndim, (fitOpt.fErrors1) ? ROOT::Fit::BinData::kNoError : ROOT::Fit::BinData::kValueError); + + const ROOT::Fit::BinData::ErrorType errorType = + (fitOpt.fErrors1) ? ROOT::Fit::BinData::kNoError : ROOT::Fit::BinData::kValueError; + + // Several disjoint ranges can be defined for each coordinate. The bins are + // collected by looping over every combination of the per-coordinate ranges. + // A coordinate without any range contributes a single iteration spanning the + // full axis. + const unsigned int nRangesX = (range.Size(0) != 0) ? range.Size(0) : 1; + const unsigned int nRangesY = (hdim > 1 && range.Size(1) != 0) ? range.Size(1) : 1; + const unsigned int nRangesZ = (hdim > 2 && range.Size(2) != 0) ? range.Size(2) : 1; double x[3]; double s[3]; @@ -188,6 +168,35 @@ void FillData(BinData & dv, const TH1 * hfit, TF1 * func) const TAxis *yaxis = hfit->GetYaxis(); const TAxis *zaxis = hfit->GetZaxis(); + for (unsigned int ix = 0; ix < nRangesX; ++ix) { + hxfirst = xaxis->GetFirst(); + hxlast = xaxis->GetLast(); + if (range.Size(0) != 0) + HFitInterface::ExamineRange( xaxis, range(0,ix), hxfirst, hxlast); + + for (unsigned int iy = 0; iy < nRangesY; ++iy) { + hyfirst = yaxis->GetFirst(); + hylast = yaxis->GetLast(); + if (hdim > 1 && range.Size(1) != 0) + HFitInterface::ExamineRange( yaxis, range(1,iy), hyfirst, hylast); + + for (unsigned int iz = 0; iz < nRangesZ; ++iz) { + hzfirst = zaxis->GetFirst(); + hzlast = zaxis->GetLast(); + if (hdim > 2 && range.Size(2) != 0) + HFitInterface::ExamineRange( zaxis, range(2,iz), hzfirst, hzlast); + + int n = (hxlast-hxfirst+1)*(hylast-hyfirst+1)*(hzlast-hzfirst+1); + +#ifdef DEBUG + std::cout << "THFitInterface: ifirst = " << hxfirst << " ilast = " << hxlast + << " total bins " << n + << std::endl; +#endif + + // reserve n for more efficient usage + dv.Append(n,ndim,errorType); + for ( binx = hxfirst; binx <= hxlast; ++binx) { if (useBinEdges) { x[0] = xaxis->GetBinLowEdge(binx); @@ -251,6 +260,10 @@ void FillData(BinData & dv, const TH1 * hfit, TF1 * func) } // end loop on y bins } // end loop on x axis + } // end loop on z ranges + } // end loop on y ranges + } // end loop on x ranges + #ifdef DEBUG std::cout << "THFitInterface::FillData: Hist FitData size is " << dv.Size() << std::endl; diff --git a/hist/hist/test/CMakeLists.txt b/hist/hist/test/CMakeLists.txt index 56e0a4fd2519a..78cfbcce14ae7 100644 --- a/hist/hist/test/CMakeLists.txt +++ b/hist/hist/test/CMakeLists.txt @@ -6,6 +6,7 @@ ROOT_ADD_GTEST(testTProfile2Poly test_tprofile2poly.cxx LIBRARIES Hist Matrix MathCore RIO) ROOT_ADD_GTEST(testTFractionFitter test_TFractionFitter.cxx LIBRARIES Hist MathCore) +ROOT_ADD_GTEST(testFillDataRange test_FillDataRange.cxx LIBRARIES Hist MathCore) ROOT_ADD_GTEST(testTH2PolyBinError test_TH2Poly_BinError.cxx LIBRARIES Hist Matrix MathCore RIO) ROOT_ADD_GTEST(testTH2PolyAdd test_TH2Poly_Add.cxx LIBRARIES Hist Matrix MathCore RIO) ROOT_ADD_GTEST(testTH2PolyGetNumberOfBins test_TH2Poly_GetNumberOfBins.cxx LIBRARIES Hist Matrix MathCore RIO) diff --git a/hist/hist/test/test_FillDataRange.cxx b/hist/hist/test/test_FillDataRange.cxx new file mode 100644 index 0000000000000..5fd4f6fb2c3ca --- /dev/null +++ b/hist/hist/test/test_FillDataRange.cxx @@ -0,0 +1,150 @@ +#include "gtest/gtest.h" + +#include "Fit/BinData.h" +#include "Fit/DataRange.h" +#include "Fit/Fitter.h" +#include "HFitInterface.h" +#include "Math/WrappedMultiTF1.h" +#include "TF1.h" +#include "TH1.h" +#include "TH2.h" + +#include +#include + +namespace { + +// linearly increasing bin contents, so that every bin is distinguishable +std::unique_ptr MakeH1() +{ + auto h = std::make_unique("h1_filldata", "h1", 100, -10., 10.); + for (int i = 1; i <= 100; ++i) + h->SetBinContent(i, 1. + 0.01 * i); + return h; +} + +std::vector CollectedCoords(const ROOT::Fit::BinData &data) +{ + std::vector x; + x.reserve(data.Size()); + for (unsigned int i = 0; i < data.Size(); ++i) + x.push_back(*data.Coords(i)); + return x; +} + +} // namespace + +// Without any range the whole axis must be used. This is the default path of +// TH1::Fit and must not be affected by the multiple-range support. +TEST(FillDataRange, NoRange) +{ + auto h = MakeH1(); + + ROOT::Fit::DataOptions opt; + ROOT::Fit::DataRange range; + ROOT::Fit::BinData data(opt, range); + ROOT::Fit::FillData(data, h.get()); + + EXPECT_EQ(data.Size(), 100u); +} + +// A single range behaves as before. +TEST(FillDataRange, SingleRange) +{ + auto h = MakeH1(); + + ROOT::Fit::DataOptions opt; + ROOT::Fit::DataRange range; + range.AddRange(0, -10., -5.); // bins 1..25 + ROOT::Fit::BinData data(opt, range); + ROOT::Fit::FillData(data, h.get()); + + EXPECT_EQ(data.Size(), 25u); + for (double x : CollectedCoords(data)) { + EXPECT_GE(x, -10.); + EXPECT_LE(x, -5.); + } +} + +// Two disjoint ranges must both contribute, and nothing outside them. +TEST(FillDataRange, TwoDisjointRanges) +{ + auto h = MakeH1(); + + ROOT::Fit::DataOptions opt; + ROOT::Fit::DataRange range; + range.AddRange(0, -10., -5.); // bins 1..25 + range.AddRange(0, 5., 10.); // bins 76..100 + ROOT::Fit::BinData data(opt, range); + ROOT::Fit::FillData(data, h.get()); + + EXPECT_EQ(data.Size(), 50u); + for (double x : CollectedCoords(data)) + EXPECT_TRUE((x >= -10. && x <= -5.) || (x >= 5. && x <= 10.)); +} + +// Overlapping ranges are merged by DataRange, so no bin is counted twice. +TEST(FillDataRange, OverlappingRangesAreNotDoubleCounted) +{ + auto h = MakeH1(); + + ROOT::Fit::DataOptions opt; + ROOT::Fit::DataRange range; + range.AddRange(0, -10., 0.); + range.AddRange(0, -5., 5.); // overlaps the previous one + ROOT::Fit::BinData data(opt, range); + ROOT::Fit::FillData(data, h.get()); + + ASSERT_EQ(range.Size(0), 1u); // merged into [-10,5] + EXPECT_EQ(data.Size(), 75u); + + auto x = CollectedCoords(data); + EXPECT_EQ(std::adjacent_find(x.begin(), x.end()), x.end()); // no duplicates +} + +// Two ranges along both axes of a 2D histogram give the product of the bins. +TEST(FillDataRange, TwoDimensions) +{ + TH2F h("h2_filldata", "h2", 10, 0., 10., 10, 0., 10.); + for (int i = 1; i <= 10; ++i) + for (int j = 1; j <= 10; ++j) + h.SetBinContent(i, j, 1. + i + 10. * j); + + ROOT::Fit::DataOptions opt; + ROOT::Fit::DataRange range; + range.AddRange(0, 0., 2.); // x bins 1..2 + range.AddRange(0, 8., 10.);// x bins 9..10 + range.AddRange(1, 0., 3.); // y bins 1..3 + ROOT::Fit::BinData data(opt, range); + ROOT::Fit::FillData(data, &h); + + EXPECT_EQ(data.Size(), 4u * 3u); +} + +// End-to-end: fit a linear background over two ranges chosen to exclude a peak. +TEST(FillDataRange, BackgroundFitExcludingPeak) +{ + TF1 comb("comb_filldata", "[0] + x*[1] + [2]*exp(-((x-[3])**2)/[4])", -10., 10.); + comb.SetParameters(10., -0.5, 15., 1., 1.5); + + TH1F h("h1_peak", "h1", 100, -10., 10.); + for (int i = 1; i <= 100; ++i) + h.SetBinContent(i, comb.Eval(h.GetBinCenter(i))); + + ROOT::Fit::DataOptions opt; + ROOT::Fit::DataRange range; + range.AddRange(0, -10., -2.5); + range.AddRange(0, 4.5, 10.); + ROOT::Fit::BinData data(opt, range); + ROOT::Fit::FillData(data, &h); + + TF1 back("back_filldata", "[0] + x*[1]", -10., 10.); + back.SetParameters(1., 0.); + ROOT::Math::WrappedMultiTF1 fitFunc(back, back.GetNdim()); + ROOT::Fit::Fitter fitter; + fitter.SetFunction(fitFunc, false); + + ASSERT_TRUE(fitter.Fit(data)); + EXPECT_NEAR(fitter.Result().Parameter(0), 10., 1e-3); + EXPECT_NEAR(fitter.Result().Parameter(1), -0.5, 1e-3); +} diff --git a/math/mathcore/inc/Fit/DataRange.h b/math/mathcore/inc/Fit/DataRange.h index 5761388212224..adba241e6648d 100644 --- a/math/mathcore/inc/Fit/DataRange.h +++ b/math/mathcore/inc/Fit/DataRange.h @@ -219,10 +219,11 @@ class DataRange { protected: /** - internal function to remove all the existing ranges between xmin and xmax - called when a new range is inserted + internal function to remove all the existing ranges overlapping with + [xmin,xmax], called when a new range is inserted. xmin and xmax are + widened to cover the removed ranges. */ - void CleanRangeSet(unsigned int icoord, double xmin, double xmax); + void CleanRangeSet(unsigned int icoord, double & xmin, double & xmax); // get the full range (-inf, +inf) static void GetInfRange(double &x1, double &x2); diff --git a/math/mathcore/src/DataRange.cxx b/math/mathcore/src/DataRange.cxx index 89980b00babbd..f254972068e3d 100644 --- a/math/mathcore/src/DataRange.cxx +++ b/math/mathcore/src/DataRange.cxx @@ -162,9 +162,10 @@ void DataRange::Clear(unsigned int icoord ) { } -void DataRange::CleanRangeSet(unsigned int icoord, double xmin, double xmax) { - // remove all the existing ranges between xmin and xmax - // called when a new range is inserted +void DataRange::CleanRangeSet(unsigned int icoord, double & xmin, double & xmax) { + // remove all the existing ranges overlapping with [xmin,xmax] and widen + // [xmin,xmax] to cover them. Called when a new range is inserted, so that + // the resulting range set stays a set of disjoint intervals. // loop on existing ranges RangeSet & ranges = fRanges[icoord]; @@ -174,6 +175,23 @@ void DataRange::CleanRangeSet(unsigned int icoord, double xmin, double xmax) { itr = ranges.erase(itr); // itr goes to next element, so go back before adding --itr; + continue; + } + // the new range is contained in an existing one: keep the existing one + if ( xmin >= itr->first && xmax <= itr->second) { + xmin = itr->first; + xmax = itr->second; + itr = ranges.erase(itr); + --itr; + continue; + } + // partial overlap on either side: merge the two ranges + if ( (itr->first <= xmin && xmin <= itr->second) || + (itr->first <= xmax && xmax <= itr->second) ) { + xmin = std::min(itr->first, xmin); + xmax = std::max(itr->second, xmax); + itr = ranges.erase(itr); + --itr; } } diff --git a/math/mathcore/test/CMakeLists.txt b/math/mathcore/test/CMakeLists.txt index a0b9dc01f8b7c..4ba098f6afafc 100644 --- a/math/mathcore/test/CMakeLists.txt +++ b/math/mathcore/test/CMakeLists.txt @@ -76,6 +76,8 @@ ROOT_ADD_GTEST(MulmodUnitNoInt128 mulmod_noint128.cxx) ROOT_ADD_GTEST(MulmodUnitOpt mulmod_opt.cxx) ROOT_ADD_GTEST(RanluxLCGUnit ranlux_lcg.cxx) ROOT_ADD_GTEST(RanluxppEngineTests RanluxppEngine.cxx LIBRARIES Core MathCore) + +ROOT_ADD_GTEST(testRange test_Range.cxx LIBRARIES MathCore) ROOT_ADD_GTEST(testDelaunay2D testDelaunay2D.cxx LIBRARIES Core MathCore) ROOT_ADD_GTEST(testFitter testFitter.cxx LIBRARIES Core MathCore Hist) ROOT_ADD_GTEST(testKNNDensity testKNNDensity.cxx LIBRARIES Core MathCore Hist) diff --git a/math/mathcore/test/test_Range.cxx b/math/mathcore/test/test_Range.cxx new file mode 100644 index 0000000000000..292fcd02ed9cb --- /dev/null +++ b/math/mathcore/test/test_Range.cxx @@ -0,0 +1,57 @@ +#include "gtest/gtest.h" + +#include "Fit/DataRange.h" + +//Overlap Range +TEST(Range, Overlap) +{ + ROOT::Fit::DataRange range; + range.AddRange(0,5); + range.AddRange(2,3); + + EXPECT_EQ(range.Size(), 1); + + EXPECT_EQ(range(0,0).first, 0); + EXPECT_EQ(range(0,0).second, 5); + + range.AddRange(-1,6); + EXPECT_EQ(range.Size(), 1); + + EXPECT_EQ(range(0,0).first, -1); + EXPECT_EQ(range(0,0).second, 6); + + range.AddRange(-2,4); + EXPECT_EQ(range.Size(), 1); + + EXPECT_EQ(range(0,0).first, -2); + EXPECT_EQ(range(0,0).second, 6); + + range.AddRange(5,7); + EXPECT_EQ(range.Size(), 1); + + EXPECT_EQ(range(0,0).first, -2); + EXPECT_EQ(range(0,0).second, 7); + + range.AddRange(20,25); + EXPECT_EQ(range.Size(), 2); + + EXPECT_EQ(range(0,0).first, -2); + EXPECT_EQ(range(0,0).second, 7); + + EXPECT_EQ(range(0,1).first, 20); + EXPECT_EQ(range(0,1).second, 25); + + range.AddRange(24,26); + EXPECT_EQ(range.Size(), 2); + EXPECT_EQ(range(0,1).first, 20); + EXPECT_EQ(range(0,1).second, 26); + + range.AddRange(19,20); + EXPECT_EQ(range(0,1).first, 19); + EXPECT_EQ(range(0,1).second, 26); + + range.AddRange(6,20); + EXPECT_EQ(range.Size(), 1); + EXPECT_EQ(range(0,0).first, -2); + EXPECT_EQ(range(0,0).second, 26); +}