From a17cc1572481f30a9249b53c54c3cb41bfc6bc2c Mon Sep 17 00:00:00 2001 From: Aapo Date: Tue, 28 Oct 2025 15:22:19 +0200 Subject: [PATCH] [hist] Don't use FindBin/GetBinCenter in TH1 derived projections Since outgoing bin indices can be computed depending only on the option to keep original binning and set axis ranges, there is no need to use FindBin and GetBinCenter to get the output bin where integrated contents are to be placed. This PR removes the usage of GetBinCenter and FindBin from the output bin index calculation, instead using simple index arithmetic. This fixes the corner-case where non-finite bin edges cause projected bins to be empty due to GetBinCenter returning a bin index that is not actually inside the bin. A slight performance improvement could also be expected, but this is unimportant since I don't expect a histogram projection to be a bottleneck in anyone's analysis. I tried adding more tests to verify that the new logic is sound and equivalent to the existing implementation in cases other than the mentioned infinite bin edges. I think someone more familiar with the internals should make sure I have covered all use cases so that this change doesn't introduce any new bugs. Fixes https://github.com/root-project/root/issues/20174 --- hist/hist/src/TH2.cxx | 24 ++- hist/hist/src/TH3.cxx | 63 ++++++-- hist/hist/test/test_projections.cxx | 242 ++++++++++++++++++++++++++++ 3 files changed, 311 insertions(+), 18 deletions(-) diff --git a/hist/hist/src/TH2.cxx b/hist/hist/src/TH2.cxx index 43d773a7c8311..35405dd2dd3c9 100644 --- a/hist/hist/src/TH2.cxx +++ b/hist/hist/src/TH2.cxx @@ -2272,10 +2272,25 @@ TH1D *TH2::DoProjection(bool onX, const char *name, Int_t firstbin, Int_t lastbi // if the out axis has labels and is extendable, temporary make it non-extendable to avoid adding extra bins Bool_t extendable = outAxis->CanExtend(); if ( labels && extendable ) h1->GetXaxis()->SetCanExtend(kFALSE); - for ( Int_t outbin = 0; outbin <= outAxis->GetNbins() + 1; ++outbin) { + + // compute loop indices + Int_t loop_start = 0; + Int_t loop_end = outAxis->GetNbins() + 1; + Int_t binOut = 0; + if (outAxis->TestBit(TAxis::kAxisRange)) { + loop_start = firstOutBin; + loop_end = lastOutBin; + if (originalRange) + // place results at original bin indices + binOut = firstOutBin; + else + // start from the first bin + binOut = 1; + } + + for (Int_t outbin = loop_start; outbin <= loop_end; ++outbin, ++binOut) { err2 = 0; cont = 0; - if (outAxis->TestBit(TAxis::kAxisRange) && ( outbin < firstOutBin || outbin > lastOutBin )) continue; for (Int_t inbin = firstbin ; inbin <= lastbin ; ++inbin) { Int_t binx, biny; @@ -2292,9 +2307,8 @@ TH1D *TH2::DoProjection(bool onX, const char *name, Int_t firstbin, Int_t lastbi err2 += exy*exy; } } - // find corresponding bin number in h1 for outbin - Int_t binOut = h1->GetXaxis()->FindBin( outAxis->GetBinCenter(outbin) ); - h1->SetBinContent(binOut ,cont); + + h1->SetBinContent(binOut, cont); if (computeErrors) h1->SetBinError(binOut,TMath::Sqrt(err2)); // sum all content totcont += cont; diff --git a/hist/hist/src/TH3.cxx b/hist/hist/src/TH3.cxx index 7cda8abb37676..7fff0dc496394 100644 --- a/hist/hist/src/TH3.cxx +++ b/hist/hist/src/TH3.cxx @@ -2038,8 +2038,23 @@ TH1D *TH3::DoProject1D(const char* name, const char * title, const TAxis* projX, // if the out axis has labels and is extendable, temporary make it non-extendable to avoid adding extra bins Bool_t extendable = projX->CanExtend(); if ( labels && extendable ) h1->GetXaxis()->SetCanExtend(kFALSE); - for (ixbin=0;ixbin<=1+projX->GetNbins();ixbin++) { - if ( projX->TestBit(TAxis::kAxisRange) && ( ixbin < ixmin || ixbin > ixmax )) continue; + + // compute loop indices + Int_t loop_start = 0; + Int_t loop_end = projX->GetNbins() + 1; + Int_t ix = 0; + if (projX->TestBit(TAxis::kAxisRange)) { + loop_start = ixmin; + loop_end = ixmax; + if (originalRange) + // place results at original bin indices + ix = ixmin; + else + // start from the first bin + ix = 1; + } + + for (ixbin = loop_start; ixbin <= loop_end; ++ixbin, ++ix) { Double_t cont = 0; Double_t err2 = 0; @@ -2058,12 +2073,10 @@ TH1D *TH3::DoProject1D(const char* name, const char * title, const TAxis* projX, } } } - Int_t ix = h1->FindBin( projX->GetBinCenter(ixbin) ); h1->SetBinContent(ix ,cont); if (computeErrors) h1->SetBinError(ix, TMath::Sqrt(err2) ); // sum all content totcont += cont; - } if ( labels ) h1->GetXaxis()->SetCanExtend(extendable); @@ -2244,7 +2257,7 @@ TH2D *TH3::DoProject2D(const char* name, const char * title, const TAxis* projX, } Int_t *refX = nullptr, *refY = nullptr, *refZ = nullptr; - Int_t ixbin, iybin, outbin; + Int_t ixbin, iybin, outbin, iy; if ( projX == GetXaxis() && projY == GetYaxis() ) { refX = &ixbin; refY = &iybin; refZ = &outbin; } if ( projX == GetYaxis() && projY == GetXaxis() ) { refX = &iybin; refY = &ixbin; refZ = &outbin; } if ( projX == GetXaxis() && projY == GetZaxis() ) { refX = &ixbin; refY = &outbin; refZ = &iybin; } @@ -2265,13 +2278,38 @@ TH2D *TH3::DoProject2D(const char* name, const char * title, const TAxis* projX, if (useUF && !out->TestBit(TAxis::kAxisRange) ) outmin -= 1; if (useOF && !out->TestBit(TAxis::kAxisRange) ) outmax += 1; - for (ixbin=0;ixbin<=1+projX->GetNbins();ixbin++) { - if ( projX->TestBit(TAxis::kAxisRange) && ( ixbin < ixmin || ixbin > ixmax )) continue; - Int_t ix = h2->GetYaxis()->FindBin( projX->GetBinCenter(ixbin) ); - - for (iybin=0;iybin<=1+projY->GetNbins();iybin++) { - if ( projY->TestBit(TAxis::kAxisRange) && ( iybin < iymin || iybin > iymax )) continue; - Int_t iy = h2->GetXaxis()->FindBin( projY->GetBinCenter(iybin) ); + // compute outer loop indices + Int_t ixbin_start = 0; + Int_t ixbin_end = projX->GetNbins() + 1; + Int_t ix = 0; + if (projX->TestBit(TAxis::kAxisRange)) { + ixbin_start = ixmin; + ixbin_end = ixmax; + if (originalRange) + // place results at original bin indices + ix = ixmin; + else + // start from the first bin + ix = 1; + } + // compute inner loop indices + Int_t iybin_start = 0; + Int_t iybin_end = projY->GetNbins() + 1; + Int_t iy_start = 0; + if (projY->TestBit(TAxis::kAxisRange)) { + iybin_start = iymin; + iybin_end = iymax; + if (originalRange) + // place results at original bin indices + iy_start = iymin; + else + // start from the first bin + iy_start = 1; + } + + for (ixbin = ixbin_start; ixbin <= ixbin_end; ++ixbin, ++ix) { + + for (iybin = iybin_start, iy = iy_start; iybin <= iybin_end; ++iybin, ++iy) { Double_t cont = 0; Double_t err2 = 0; @@ -2295,7 +2333,6 @@ TH2D *TH3::DoProject2D(const char* name, const char * title, const TAxis* projX, if (computeErrors) h2->SetBinError(iy, ix, TMath::Sqrt(err2) ); // sum all content totcont += cont; - } } diff --git a/hist/hist/test/test_projections.cxx b/hist/hist/test/test_projections.cxx index d57d4f325491d..70c27294fd626 100644 --- a/hist/hist/test/test_projections.cxx +++ b/hist/hist/test/test_projections.cxx @@ -3,6 +3,8 @@ #include "TProfile.h" // ProjectionX #include "TProfile2D.h" // ProjectionX #include "THashList.h" // GetLabels +#include +#include "TMath.h" #include "gtest/gtest.h" @@ -85,3 +87,243 @@ TEST(Projections, Issue_6658_Profile2D) EXPECT_EQ(xaxis_2d_nbins, xaxis_pxy_nbins); expect_list_eq_names(*labels_2d, *labels_pxy); } + +// Test projection from TH3D for correct output on user range on projected axis +TEST(Projections, RangesAndOptionO) +{ + TH3D h("h", "h", 3, 0., 3., 3, 0., 3., 3, 0., 3.); + + for (int ix = 1; ix <= 3; ++ix) { + for (int iy = 1; iy <= 3; ++iy) { + for (int iz = 1; iz <= 3; ++iz) { + auto bin = h.GetBin(ix, iy, iz); + h.SetBinContent(bin, 100 * ix + 10 * iy + iz); + } + } + } + + h.GetXaxis()->SetRange(2, 3); + auto expectedForX = [](int ix) { + double s = 0.; + for (int iy = 1; iy <= 3; ++iy) + for (int iz = 1; iz <= 3; ++iz) + s += 100 * ix + 10 * iy + iz; + return s; + }; + auto x2 = expectedForX(2); + auto x3 = expectedForX(3); + + { + auto px = h.Project3D("x"); + EXPECT_EQ(px->GetNbinsX(), 2); // selected length + EXPECT_DOUBLE_EQ(px->GetBinContent(1), x2); + EXPECT_DOUBLE_EQ(px->GetBinContent(2), x3); + } + + { + auto pxo = h.Project3D("xo"); + ASSERT_NE(pxo, nullptr); + EXPECT_EQ(pxo->GetNbinsX(), 3); // original length + EXPECT_DOUBLE_EQ(pxo->GetBinContent(1), 0.0); // outside selection + EXPECT_DOUBLE_EQ(pxo->GetBinContent(2), x2); + EXPECT_DOUBLE_EQ(pxo->GetBinContent(3), x3); + } +} + +// Test projection from TH3D for correct output on user range on integrated axis +TEST(Projections, SelectionAcrossNonTargetAxis) +{ + TH3D h("h", "h", 3, 0., 3., 3, 0., 3., 3, 0., 3.); + + for (int ix = 1; ix <= 3; ++ix) { + for (int iy = 1; iy <= 3; ++iy) { + for (int iz = 1; iz <= 3; ++iz) { + auto bin = h.GetBin(ix, iy, iz); + h.SetBinContent(bin, 100 * ix + 10 * iy + iz); + } + } + } + + h.GetYaxis()->SetRange(2, 3); + + auto px = h.Project3D("x"); + + auto expectedForX = [](int ix) { + double s = 0.; + for (int iy = 2; iy <= 3; ++iy) + for (int iz = 1; iz <= 3; ++iz) + s += 100 * ix + 10 * iy + iz; + return s; + }; + + EXPECT_DOUBLE_EQ(px->GetBinContent(1), expectedForX(1)); + EXPECT_DOUBLE_EQ(px->GetBinContent(2), expectedForX(2)); + EXPECT_DOUBLE_EQ(px->GetBinContent(3), expectedForX(3)); +} + +// Test TH2D projection for correctness for user ranges on both axes +TEST(Projections, ProjectionYRange) +{ + Double_t xedges[] = {0, 1, 2}; + Double_t yedges[] = {-2, -1, 0, 1, 2}; + TH2D h("h", "h;X;Y", 2, xedges, 4, yedges); + + for (int ix = 1; ix <= 3; ++ix) { + for (int iy = 1; iy <= 4; ++iy) { + auto bin = h.GetBin(ix, iy); + h.SetBinContent(bin, 10 * ix + iy); + } + } + + h.GetXaxis()->SetRange(1, 2); + h.GetYaxis()->SetRange(2, 4); + + auto py = h.ProjectionY(); + + auto expectedForY = [](int iy) { + double s = 0.; + for (int ix = 1; ix <= 2; ++ix) + s += 10 * ix + iy; + return s; + }; + + EXPECT_DOUBLE_EQ(py->GetBinContent(0), 0.0); + EXPECT_DOUBLE_EQ(py->GetBinContent(1), expectedForY(2)); + EXPECT_DOUBLE_EQ(py->GetBinContent(2), expectedForY(3)); + EXPECT_DOUBLE_EQ(py->GetBinContent(3), expectedForY(4)); + EXPECT_DOUBLE_EQ(py->GetBinContent(4), 0.0); +} +// Test TH2D projection for correct flow inclusion for default options +TEST(Projections, UFOF) +{ + TH2D h2("h2", "", 3, 0, 3, 4, 0, 4); + h2.Sumw2(); + for (int bx = 0; bx <= 4; ++bx) { + for (int by = 0; by <= 5; ++by) { + h2.SetBinContent(bx, by, 10 * bx + by); + h2.SetBinError(bx, by, 1.0); + } + } + + auto hpx = h2.ProjectionX(); + for (int bx = 0; bx <= 4; ++bx) { + const double exp = 60 * bx + 15; + double got = hpx->GetBinContent(bx); + EXPECT_DOUBLE_EQ(got, exp); + double e = hpx->GetBinError(bx); + EXPECT_DOUBLE_EQ(e, TMath::Sqrt(6.0)); + } +} +// Test TH2D projection for correct flow exclusion for specified user range +TEST(Projections, UFOFWithRange) +{ + TH2D h2("h2", "", 5, 0, 5, 4, 0, 4); + h2.Sumw2(); + for (int bx = 0; bx <= 6; ++bx) + for (int by = 0; by <= 5; ++by) { + h2.SetBinContent(bx, by, 100 * bx + by); + h2.SetBinError(bx, by, 2.0); + } + + h2.GetXaxis()->SetRange(2, 4); + + auto hpx = h2.ProjectionX("hpx_ranged", 1, 4, ""); + EXPECT_EQ(hpx->GetXaxis()->GetNbins(), 3); + for (int i = 0; i < 3; ++i) { + int outbin = 2 + i; + double exp = 0; + for (int by = 1; by <= 4; ++by) + exp += 100 * outbin + by; + double got = hpx->GetBinContent(i + 1); + EXPECT_DOUBLE_EQ(got, exp); + EXPECT_DOUBLE_EQ(hpx->GetBinError(i + 1), 4.0); + } + EXPECT_DOUBLE_EQ(hpx->GetBinContent(0), 0); + EXPECT_DOUBLE_EQ(hpx->GetBinContent(4), 0); +} + +// Test TH2D projection correctness with option "o" +TEST(Projections, OriginalRange) +{ + TH2D h2("h2", "", 6, 0, 6, 3, 0, 3); + for (int bx = 0; bx <= 7; ++bx) + for (int by = 0; by <= 4; ++by) + h2.SetBinContent(bx, by, bx + 10 * by); + + h2.GetXaxis()->SetRange(2, 5); + + auto hpx = h2.ProjectionX("h_o", 1, 3, "o"); + EXPECT_EQ(hpx->GetXaxis()->GetNbins(), 6); + for (int bx = 1; bx <= 6; ++bx) { + double got = hpx->GetBinContent(bx); + if (bx < 2 || bx > 5) { + EXPECT_EQ(got, 0.0); + continue; + } + double exp = 0; + for (int by = 1; by <= 3; ++by) + exp += bx + 10 * by; + EXPECT_DOUBLE_EQ(got, exp); + } +} + +// Test TH2D projection with variable bins and user range along projected axis +TEST(Projections, VarBinsRange) +{ + double edgesX[] = {0, 0.5, 1.5, 3.0, 5.0}; + TH2D h("h", "", 4, edgesX, 2, 0, 2); + + for (int bx = 0; bx <= 5; ++bx) + for (int by = 0; by <= 3; ++by) + h.SetBinContent(bx, by, 100 * bx + by); + + h.GetXaxis()->SetRange(2, 3); + + auto hpx = h.ProjectionX("hpx_var", 1, 2, ""); + EXPECT_EQ(hpx->GetXaxis()->GetNbins(), 2); + EXPECT_DOUBLE_EQ(hpx->GetXaxis()->GetBinLowEdge(1), edgesX[1]); + EXPECT_DOUBLE_EQ(hpx->GetXaxis()->GetBinUpEdge(2), edgesX[3]); + + for (int i = 0; i < 2; ++i) { + int outbin = 2 + i; + double exp = 0; + for (int by = 1; by <= 2; ++by) + exp += 100 * outbin + by; + double got = hpx->GetBinContent(i + 1); + EXPECT_DOUBLE_EQ(got, exp); + } +} + +// https://github.com/root-project/issues/20174 +TEST(Projections, ProjectionYInfiniteUpperEdge) +{ + Double_t xedges[] = {0., 1.}; + Double_t yedges[] = {1., std::numeric_limits::infinity()}; + TH2D h("h_inf", "h_inf;X;Y", 1, xedges, 1, yedges); + h.SetBinContent(1, 1, 11.); + auto projY = h.ProjectionY(); + EXPECT_EQ(projY->GetBinContent(1), h.Integral(1, 1, 1, 1)); +} +TEST(Projections, ProjectionYInfiniteLowerEdge) +{ + Double_t xedges[] = {0, 1.}; + Double_t yedges[] = {-std::numeric_limits::infinity(), 1}; + TH2D h("h_inf", "h_inf;X;Y", 1, xedges, 1, yedges); + h.SetBinContent(1, 1, 11.); + auto projY = h.ProjectionY(); + EXPECT_EQ(projY->GetBinContent(1), h.Integral(1, 1, 1, 1)); +} + +TEST(Projections, Projection3DInfiniteEdges) +{ + Double_t xedges[] = {0, 1.}; + Double_t yedges[] = {-std::numeric_limits::infinity(), 1}; + Double_t zedges[] = {0, std::numeric_limits::infinity()}; + TH3D h("h_inf", "h_inf;X;Y;Z", 1, xedges, 1, yedges, 1, zedges); + h.SetBinContent(1, 1, 1, 11.); + auto projY = h.Project3D("zy"); + EXPECT_EQ(projY->GetBinContent(projY->GetBin(1, 1)), h.Integral(1, 1, 1, 1, 1, 1)); + + auto projx = h.Project3D("x"); + EXPECT_EQ(projx->GetBinContent(projx->GetBin(1)), h.Integral(1, 1, 1, 1, 1, 1)); +}