From 1414f0427614e5b977c5a5fbcf946be56dc11bb2 Mon Sep 17 00:00:00 2001 From: Shuvam Pandey Date: Sat, 1 Aug 2026 21:51:03 +0545 Subject: [PATCH] [hist] Fix GetRandom with negative bins The negative-bin check in ComputeIntegral(true) used break, which only left the inner loop. In 2D and 3D, later rows could overwrite the NaN marker. GetRandom, GetRandom2 and GetRandom3 also compared the result directly to NaN, so their error paths were never taken. Return immediately without caching a partial integral, use std::isnan in the GetRandom methods, and make the TH2 and TH3 FillRandom overloads validate the source before its integral is cached. This also keeps TH2Poly::GetRandom2 from calling BinarySearch with a null integral. --- hist/hist/src/TH1.cxx | 11 ++- hist/hist/src/TH2.cxx | 16 +++- hist/hist/src/TH3.cxx | 17 ++++- hist/hist/test/test_TH1.cxx | 145 ++++++++++++++++++++++++++++++++++++ 4 files changed, 182 insertions(+), 7 deletions(-) diff --git a/hist/hist/src/TH1.cxx b/hist/hist/src/TH1.cxx index f4465a04f2862..371774ba27c25 100644 --- a/hist/hist/src/TH1.cxx +++ b/hist/hist/src/TH1.cxx @@ -2609,8 +2609,12 @@ Double_t TH1::ComputeIntegral(Bool_t onlyPositive, Option_t *option) if (onlyPositive && y < 0) { Error("ComputeIntegral","Bin content is negative - return a NaN value"); - fIntegral[nbins] = TMath::QuietNaN(); - break; + // Do not cache a partial integral: later rows could overwrite the + // failure marker, and other consumers assume fIntegral is a valid CDF. + // A subsequent request recomputes the integral and reports the error again. + delete[] fIntegral; + fIntegral = nullptr; + return TMath::QuietNaN(); } fIntegral[ibin] = fIntegral[ibin - 1] + y; } @@ -5204,7 +5208,8 @@ Double_t TH1::GetRandom(TRandom *rng, Option_t *option) const } if (integral == 0) return 0; // return a NaN in case some bins have negative content - if (integral == TMath::QuietNaN() ) return TMath::QuietNaN(); + if (std::isnan(integral)) + return TMath::QuietNaN(); Double_t r1 = (rng) ? rng->Rndm() : gRandom->Rndm(); Int_t ibin = TMath::BinarySearch(nbinsx,fIntegral,r1); diff --git a/hist/hist/src/TH2.cxx b/hist/hist/src/TH2.cxx index dfd5df93f51be..22eb6630033e8 100644 --- a/hist/hist/src/TH2.cxx +++ b/hist/hist/src/TH2.cxx @@ -28,6 +28,7 @@ #include "TVirtualHistPainter.h" #include "snprintf.h" +#include /** \addtogroup Histograms @{ @@ -765,7 +766,14 @@ void TH2::FillRandom(TH1 *h, Int_t ntimes, TRandom * rng) Error("FillRandom", "Histograms with different dimensions"); return; } - if (h->ComputeIntegral() == 0) return; + // Do not let GetRandom2 reuse an integral computed without the negative-bin check. + const Double_t integral = h->ComputeIntegral(true); + if (std::isnan(integral)) { + Error("FillRandom", "Histograms contains negative bins, does not represent probabilities"); + return; + } + if (integral == 0) + return; Int_t loop; Double_t x,y; @@ -1184,7 +1192,11 @@ void TH2::GetRandom2(Double_t &x, Double_t &y, TRandom *rng, Option_t *option) } if (integral == 0 ) { x = 0; y = 0; return;} // case histogram has negative bins - if (integral == TMath::QuietNaN() ) { x = TMath::QuietNaN(); y = TMath::QuietNaN(); return;} + if (std::isnan(integral)) { + x = TMath::QuietNaN(); + y = TMath::QuietNaN(); + return; + } if (!rng) rng = gRandom; Double_t r1 = rng->Rndm(); diff --git a/hist/hist/src/TH3.cxx b/hist/hist/src/TH3.cxx index de7a3de5e0aa7..8858093efa83d 100644 --- a/hist/hist/src/TH3.cxx +++ b/hist/hist/src/TH3.cxx @@ -27,6 +27,7 @@ #include #include +#include #include @@ -959,7 +960,14 @@ void TH3::FillRandom(TH1 *h, Int_t ntimes, TRandom * rng) Error("FillRandom", "Histograms with different dimensions"); return; } - if (h->ComputeIntegral() == 0) return; + // Do not let GetRandom3 reuse an integral computed without the negative-bin check. + const Double_t integral = h->ComputeIntegral(true); + if (std::isnan(integral)) { + Error("FillRandom", "Histograms contains negative bins, does not represent probabilities"); + return; + } + if (integral == 0) + return; TH3 *h3 = (TH3*)h; Int_t loop; @@ -1302,7 +1310,12 @@ void TH3::GetRandom3(Double_t &x, Double_t &y, Double_t &z, TRandom *rng, Option } if (integral == 0 ) { x = 0; y = 0; z = 0; return;} // case histogram has negative bins - if (integral == TMath::QuietNaN() ) { x = TMath::QuietNaN(); y = TMath::QuietNaN(); z = TMath::QuietNaN(); return;} + if (std::isnan(integral)) { + x = TMath::QuietNaN(); + y = TMath::QuietNaN(); + z = TMath::QuietNaN(); + return; + } if (!rng) rng = gRandom; Double_t r1 = rng->Rndm(); diff --git a/hist/hist/test/test_TH1.cxx b/hist/hist/test/test_TH1.cxx index dffcd66d679bb..8cceb9f7632d6 100644 --- a/hist/hist/test/test_TH1.cxx +++ b/hist/hist/test/test_TH1.cxx @@ -4,9 +4,11 @@ #include "TH2.h" #include "TH3.h" #include "TH1F.h" +#include "TH2Poly.h" #include "THLimitsFinder.h" #include "TDirectory.h" #include "TList.h" +#include "TRandom3.h" #include "TROOT.h" #include @@ -562,3 +564,146 @@ TEST(TH1, GetCumulativeErrors) } } } + +namespace { + +class TH1DIntegralAccess : public TH1D { +public: + using TH1D::TH1D; + bool HasCachedIntegral() const { return fIntegral != nullptr; } +}; + +} // namespace + +TEST(TH1, GetRandomNegativeBin) +{ + ROOT::TestSupport::CheckDiagsRAII diags; + diags.requiredDiag(kError, "TH1D::ComputeIntegral", "Bin content is negative - return a NaN value"); + + TH1DIntegralAccess h("hNegRandom", "", 2, 0., 2.); + h.SetBinContent(1, 1.); + h.SetBinContent(2, -0.25); + + TRandom3 rng(1); + TRandom3 referenceRng(1); + EXPECT_TRUE(std::isnan(h.GetRandom(&rng))); + EXPECT_FALSE(h.HasCachedIntegral()); + EXPECT_DOUBLE_EQ(rng.Rndm(), referenceRng.Rndm()); +} + +TEST(TH2, GetRandom2NegativeBin) +{ + ROOT::TestSupport::CheckDiagsRAII diags; + diags.requiredDiag(kError, "TH2D::ComputeIntegral", "Bin content is negative - return a NaN value"); + + TH2D h("hNegRandom2", "", 2, 0., 2., 2, 0., 2.); + h.SetBinContent(1, 1, 1.); + h.SetBinContent(2, 1, -0.25); + h.SetBinContent(1, 2, 1.); + h.SetBinContent(2, 2, 1.); + + TRandom3 rng(1); + TRandom3 referenceRng(1); + double x = 0., y = 0.; + h.GetRandom2(x, y, &rng); + EXPECT_TRUE(std::isnan(x)); + EXPECT_TRUE(std::isnan(y)); + EXPECT_DOUBLE_EQ(rng.Rndm(), referenceRng.Rndm()); +} + +TEST(TH3, GetRandom3NegativeBin) +{ + ROOT::TestSupport::CheckDiagsRAII diags; + diags.requiredDiag(kError, "TH3D::ComputeIntegral", "Bin content is negative - return a NaN value"); + + TH3D h("hNegRandom3", "", 2, 0., 2., 2, 0., 2., 2, 0., 2.); + h.SetBinContent(1, 1, 1, 1.); + h.SetBinContent(2, 2, 2, -0.25); + + TRandom3 rng(1); + TRandom3 referenceRng(1); + double x = 0., y = 0., z = 0.; + h.GetRandom3(x, y, z, &rng); + EXPECT_TRUE(std::isnan(x)); + EXPECT_TRUE(std::isnan(y)); + EXPECT_TRUE(std::isnan(z)); + EXPECT_DOUBLE_EQ(rng.Rndm(), referenceRng.Rndm()); +} + +TEST(TH2, FillRandomNegativeBin) +{ + ROOT::TestSupport::CheckDiagsRAII diags; + diags.requiredDiag(kError, "TH2D::ComputeIntegral", "Bin content is negative - return a NaN value"); + diags.requiredDiag(kError, "TH2D::FillRandom", + "Histograms contains negative bins, does not represent probabilities"); + + TH2D src("hFillRandom2Src", "", 2, 0., 2., 2, 0., 2.); + src.SetBinContent(1, 1, 1.); + src.SetBinContent(2, 1, -0.25); + src.SetBinContent(1, 2, 1.); + src.SetBinContent(2, 2, 1.); + + TH2D dest("hFillRandom2Dest", "", 2, 0., 2., 2, 0., 2.); + TRandom3 rng(1); + dest.FillRandom(&src, 1, &rng); + EXPECT_EQ(dest.GetEntries(), 0.); +} + +TEST(TH3, FillRandomNegativeBin) +{ + ROOT::TestSupport::CheckDiagsRAII diags; + diags.requiredDiag(kError, "TH3D::ComputeIntegral", "Bin content is negative - return a NaN value"); + diags.requiredDiag(kError, "TH3D::FillRandom", + "Histograms contains negative bins, does not represent probabilities"); + + TH3D src("hFillRandom3Src", "", 2, 0., 2., 2, 0., 2., 2, 0., 2.); + src.SetBinContent(1, 1, 1, 1.); + src.SetBinContent(2, 2, 2, -0.25); + + TH3D dest("hFillRandom3Dest", "", 2, 0., 2., 2, 0., 2., 2, 0., 2.); + TRandom3 rng(1); + dest.FillRandom(&src, 1, &rng); + EXPECT_EQ(dest.GetEntries(), 0.); +} + +TEST(TH2, FillRandomPositiveBins) +{ + ROOT::TestSupport::CheckDiagsRAII diags; + TRandom3 rng(1); + + TH2D src("hFillRandom2PositiveSrc", "", 1, 0., 1., 1, 0., 1.); + src.SetBinContent(1, 1, 1.); + TH2D dest("hFillRandom2PositiveDest", "", 1, 0., 1., 1, 0., 1.); + dest.FillRandom(&src, 1, &rng); + EXPECT_EQ(dest.GetEntries(), 1.); + EXPECT_EQ(dest.Integral(), 1.); +} + +TEST(TH3, FillRandomPositiveBins) +{ + ROOT::TestSupport::CheckDiagsRAII diags; + TRandom3 rng(1); + + TH3D src("hFillRandom3PositiveSrc", "", 1, 0., 1., 1, 0., 1., 1, 0., 1.); + src.SetBinContent(1, 1, 1, 1.); + TH3D dest("hFillRandom3PositiveDest", "", 1, 0., 1., 1, 0., 1., 1, 0., 1.); + dest.FillRandom(&src, 1, &rng); + EXPECT_EQ(dest.GetEntries(), 1.); + EXPECT_EQ(dest.Integral(), 1.); +} + +TEST(TH2Poly, GetRandom2NotImplemented) +{ + ROOT::TestSupport::CheckDiagsRAII diags; + diags.requiredDiag(kError, "TH2Poly::ComputeIntegral", "Not implemented for TH2Poly"); + + TH2Poly h("hPolyRandom", "", 0., 2., 0., 2.); + h.AddBin(0., 0., 1., 1.); + h.AddBin(1., 1., 2., 2.); + + TRandom3 rng(1); + double x = 0., y = 0.; + h.GetRandom2(x, y, &rng); + EXPECT_TRUE(std::isnan(x)); + EXPECT_TRUE(std::isnan(y)); +}