Skip to content

Commit a17cc15

Browse files
committed
[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 #20174
1 parent 8152c4d commit a17cc15

File tree

3 files changed

+311
-18
lines changed

3 files changed

+311
-18
lines changed

hist/hist/src/TH2.cxx

Lines changed: 19 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -2272,10 +2272,25 @@ TH1D *TH2::DoProjection(bool onX, const char *name, Int_t firstbin, Int_t lastbi
22722272
// if the out axis has labels and is extendable, temporary make it non-extendable to avoid adding extra bins
22732273
Bool_t extendable = outAxis->CanExtend();
22742274
if ( labels && extendable ) h1->GetXaxis()->SetCanExtend(kFALSE);
2275-
for ( Int_t outbin = 0; outbin <= outAxis->GetNbins() + 1; ++outbin) {
2275+
2276+
// compute loop indices
2277+
Int_t loop_start = 0;
2278+
Int_t loop_end = outAxis->GetNbins() + 1;
2279+
Int_t binOut = 0;
2280+
if (outAxis->TestBit(TAxis::kAxisRange)) {
2281+
loop_start = firstOutBin;
2282+
loop_end = lastOutBin;
2283+
if (originalRange)
2284+
// place results at original bin indices
2285+
binOut = firstOutBin;
2286+
else
2287+
// start from the first bin
2288+
binOut = 1;
2289+
}
2290+
2291+
for (Int_t outbin = loop_start; outbin <= loop_end; ++outbin, ++binOut) {
22762292
err2 = 0;
22772293
cont = 0;
2278-
if (outAxis->TestBit(TAxis::kAxisRange) && ( outbin < firstOutBin || outbin > lastOutBin )) continue;
22792294

22802295
for (Int_t inbin = firstbin ; inbin <= lastbin ; ++inbin) {
22812296
Int_t binx, biny;
@@ -2292,9 +2307,8 @@ TH1D *TH2::DoProjection(bool onX, const char *name, Int_t firstbin, Int_t lastbi
22922307
err2 += exy*exy;
22932308
}
22942309
}
2295-
// find corresponding bin number in h1 for outbin
2296-
Int_t binOut = h1->GetXaxis()->FindBin( outAxis->GetBinCenter(outbin) );
2297-
h1->SetBinContent(binOut ,cont);
2310+
2311+
h1->SetBinContent(binOut, cont);
22982312
if (computeErrors) h1->SetBinError(binOut,TMath::Sqrt(err2));
22992313
// sum all content
23002314
totcont += cont;

hist/hist/src/TH3.cxx

Lines changed: 50 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -2038,8 +2038,23 @@ TH1D *TH3::DoProject1D(const char* name, const char * title, const TAxis* projX,
20382038
// if the out axis has labels and is extendable, temporary make it non-extendable to avoid adding extra bins
20392039
Bool_t extendable = projX->CanExtend();
20402040
if ( labels && extendable ) h1->GetXaxis()->SetCanExtend(kFALSE);
2041-
for (ixbin=0;ixbin<=1+projX->GetNbins();ixbin++) {
2042-
if ( projX->TestBit(TAxis::kAxisRange) && ( ixbin < ixmin || ixbin > ixmax )) continue;
2041+
2042+
// compute loop indices
2043+
Int_t loop_start = 0;
2044+
Int_t loop_end = projX->GetNbins() + 1;
2045+
Int_t ix = 0;
2046+
if (projX->TestBit(TAxis::kAxisRange)) {
2047+
loop_start = ixmin;
2048+
loop_end = ixmax;
2049+
if (originalRange)
2050+
// place results at original bin indices
2051+
ix = ixmin;
2052+
else
2053+
// start from the first bin
2054+
ix = 1;
2055+
}
2056+
2057+
for (ixbin = loop_start; ixbin <= loop_end; ++ixbin, ++ix) {
20432058

20442059
Double_t cont = 0;
20452060
Double_t err2 = 0;
@@ -2058,12 +2073,10 @@ TH1D *TH3::DoProject1D(const char* name, const char * title, const TAxis* projX,
20582073
}
20592074
}
20602075
}
2061-
Int_t ix = h1->FindBin( projX->GetBinCenter(ixbin) );
20622076
h1->SetBinContent(ix ,cont);
20632077
if (computeErrors) h1->SetBinError(ix, TMath::Sqrt(err2) );
20642078
// sum all content
20652079
totcont += cont;
2066-
20672080
}
20682081
if ( labels ) h1->GetXaxis()->SetCanExtend(extendable);
20692082

@@ -2244,7 +2257,7 @@ TH2D *TH3::DoProject2D(const char* name, const char * title, const TAxis* projX,
22442257
}
22452258

22462259
Int_t *refX = nullptr, *refY = nullptr, *refZ = nullptr;
2247-
Int_t ixbin, iybin, outbin;
2260+
Int_t ixbin, iybin, outbin, iy;
22482261
if ( projX == GetXaxis() && projY == GetYaxis() ) { refX = &ixbin; refY = &iybin; refZ = &outbin; }
22492262
if ( projX == GetYaxis() && projY == GetXaxis() ) { refX = &iybin; refY = &ixbin; refZ = &outbin; }
22502263
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,
22652278
if (useUF && !out->TestBit(TAxis::kAxisRange) ) outmin -= 1;
22662279
if (useOF && !out->TestBit(TAxis::kAxisRange) ) outmax += 1;
22672280

2268-
for (ixbin=0;ixbin<=1+projX->GetNbins();ixbin++) {
2269-
if ( projX->TestBit(TAxis::kAxisRange) && ( ixbin < ixmin || ixbin > ixmax )) continue;
2270-
Int_t ix = h2->GetYaxis()->FindBin( projX->GetBinCenter(ixbin) );
2271-
2272-
for (iybin=0;iybin<=1+projY->GetNbins();iybin++) {
2273-
if ( projY->TestBit(TAxis::kAxisRange) && ( iybin < iymin || iybin > iymax )) continue;
2274-
Int_t iy = h2->GetXaxis()->FindBin( projY->GetBinCenter(iybin) );
2281+
// compute outer loop indices
2282+
Int_t ixbin_start = 0;
2283+
Int_t ixbin_end = projX->GetNbins() + 1;
2284+
Int_t ix = 0;
2285+
if (projX->TestBit(TAxis::kAxisRange)) {
2286+
ixbin_start = ixmin;
2287+
ixbin_end = ixmax;
2288+
if (originalRange)
2289+
// place results at original bin indices
2290+
ix = ixmin;
2291+
else
2292+
// start from the first bin
2293+
ix = 1;
2294+
}
2295+
// compute inner loop indices
2296+
Int_t iybin_start = 0;
2297+
Int_t iybin_end = projY->GetNbins() + 1;
2298+
Int_t iy_start = 0;
2299+
if (projY->TestBit(TAxis::kAxisRange)) {
2300+
iybin_start = iymin;
2301+
iybin_end = iymax;
2302+
if (originalRange)
2303+
// place results at original bin indices
2304+
iy_start = iymin;
2305+
else
2306+
// start from the first bin
2307+
iy_start = 1;
2308+
}
2309+
2310+
for (ixbin = ixbin_start; ixbin <= ixbin_end; ++ixbin, ++ix) {
2311+
2312+
for (iybin = iybin_start, iy = iy_start; iybin <= iybin_end; ++iybin, ++iy) {
22752313

22762314
Double_t cont = 0;
22772315
Double_t err2 = 0;
@@ -2295,7 +2333,6 @@ TH2D *TH3::DoProject2D(const char* name, const char * title, const TAxis* projX,
22952333
if (computeErrors) h2->SetBinError(iy, ix, TMath::Sqrt(err2) );
22962334
// sum all content
22972335
totcont += cont;
2298-
22992336
}
23002337
}
23012338

hist/hist/test/test_projections.cxx

Lines changed: 242 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,8 @@
33
#include "TProfile.h" // ProjectionX
44
#include "TProfile2D.h" // ProjectionX
55
#include "THashList.h" // GetLabels
6+
#include <limits>
7+
#include "TMath.h"
68

79
#include "gtest/gtest.h"
810

@@ -85,3 +87,243 @@ TEST(Projections, Issue_6658_Profile2D)
8587
EXPECT_EQ(xaxis_2d_nbins, xaxis_pxy_nbins);
8688
expect_list_eq_names(*labels_2d, *labels_pxy);
8789
}
90+
91+
// Test projection from TH3D for correct output on user range on projected axis
92+
TEST(Projections, RangesAndOptionO)
93+
{
94+
TH3D h("h", "h", 3, 0., 3., 3, 0., 3., 3, 0., 3.);
95+
96+
for (int ix = 1; ix <= 3; ++ix) {
97+
for (int iy = 1; iy <= 3; ++iy) {
98+
for (int iz = 1; iz <= 3; ++iz) {
99+
auto bin = h.GetBin(ix, iy, iz);
100+
h.SetBinContent(bin, 100 * ix + 10 * iy + iz);
101+
}
102+
}
103+
}
104+
105+
h.GetXaxis()->SetRange(2, 3);
106+
auto expectedForX = [](int ix) {
107+
double s = 0.;
108+
for (int iy = 1; iy <= 3; ++iy)
109+
for (int iz = 1; iz <= 3; ++iz)
110+
s += 100 * ix + 10 * iy + iz;
111+
return s;
112+
};
113+
auto x2 = expectedForX(2);
114+
auto x3 = expectedForX(3);
115+
116+
{
117+
auto px = h.Project3D("x");
118+
EXPECT_EQ(px->GetNbinsX(), 2); // selected length
119+
EXPECT_DOUBLE_EQ(px->GetBinContent(1), x2);
120+
EXPECT_DOUBLE_EQ(px->GetBinContent(2), x3);
121+
}
122+
123+
{
124+
auto pxo = h.Project3D("xo");
125+
ASSERT_NE(pxo, nullptr);
126+
EXPECT_EQ(pxo->GetNbinsX(), 3); // original length
127+
EXPECT_DOUBLE_EQ(pxo->GetBinContent(1), 0.0); // outside selection
128+
EXPECT_DOUBLE_EQ(pxo->GetBinContent(2), x2);
129+
EXPECT_DOUBLE_EQ(pxo->GetBinContent(3), x3);
130+
}
131+
}
132+
133+
// Test projection from TH3D for correct output on user range on integrated axis
134+
TEST(Projections, SelectionAcrossNonTargetAxis)
135+
{
136+
TH3D h("h", "h", 3, 0., 3., 3, 0., 3., 3, 0., 3.);
137+
138+
for (int ix = 1; ix <= 3; ++ix) {
139+
for (int iy = 1; iy <= 3; ++iy) {
140+
for (int iz = 1; iz <= 3; ++iz) {
141+
auto bin = h.GetBin(ix, iy, iz);
142+
h.SetBinContent(bin, 100 * ix + 10 * iy + iz);
143+
}
144+
}
145+
}
146+
147+
h.GetYaxis()->SetRange(2, 3);
148+
149+
auto px = h.Project3D("x");
150+
151+
auto expectedForX = [](int ix) {
152+
double s = 0.;
153+
for (int iy = 2; iy <= 3; ++iy)
154+
for (int iz = 1; iz <= 3; ++iz)
155+
s += 100 * ix + 10 * iy + iz;
156+
return s;
157+
};
158+
159+
EXPECT_DOUBLE_EQ(px->GetBinContent(1), expectedForX(1));
160+
EXPECT_DOUBLE_EQ(px->GetBinContent(2), expectedForX(2));
161+
EXPECT_DOUBLE_EQ(px->GetBinContent(3), expectedForX(3));
162+
}
163+
164+
// Test TH2D projection for correctness for user ranges on both axes
165+
TEST(Projections, ProjectionYRange)
166+
{
167+
Double_t xedges[] = {0, 1, 2};
168+
Double_t yedges[] = {-2, -1, 0, 1, 2};
169+
TH2D h("h", "h;X;Y", 2, xedges, 4, yedges);
170+
171+
for (int ix = 1; ix <= 3; ++ix) {
172+
for (int iy = 1; iy <= 4; ++iy) {
173+
auto bin = h.GetBin(ix, iy);
174+
h.SetBinContent(bin, 10 * ix + iy);
175+
}
176+
}
177+
178+
h.GetXaxis()->SetRange(1, 2);
179+
h.GetYaxis()->SetRange(2, 4);
180+
181+
auto py = h.ProjectionY();
182+
183+
auto expectedForY = [](int iy) {
184+
double s = 0.;
185+
for (int ix = 1; ix <= 2; ++ix)
186+
s += 10 * ix + iy;
187+
return s;
188+
};
189+
190+
EXPECT_DOUBLE_EQ(py->GetBinContent(0), 0.0);
191+
EXPECT_DOUBLE_EQ(py->GetBinContent(1), expectedForY(2));
192+
EXPECT_DOUBLE_EQ(py->GetBinContent(2), expectedForY(3));
193+
EXPECT_DOUBLE_EQ(py->GetBinContent(3), expectedForY(4));
194+
EXPECT_DOUBLE_EQ(py->GetBinContent(4), 0.0);
195+
}
196+
// Test TH2D projection for correct flow inclusion for default options
197+
TEST(Projections, UFOF)
198+
{
199+
TH2D h2("h2", "", 3, 0, 3, 4, 0, 4);
200+
h2.Sumw2();
201+
for (int bx = 0; bx <= 4; ++bx) {
202+
for (int by = 0; by <= 5; ++by) {
203+
h2.SetBinContent(bx, by, 10 * bx + by);
204+
h2.SetBinError(bx, by, 1.0);
205+
}
206+
}
207+
208+
auto hpx = h2.ProjectionX();
209+
for (int bx = 0; bx <= 4; ++bx) {
210+
const double exp = 60 * bx + 15;
211+
double got = hpx->GetBinContent(bx);
212+
EXPECT_DOUBLE_EQ(got, exp);
213+
double e = hpx->GetBinError(bx);
214+
EXPECT_DOUBLE_EQ(e, TMath::Sqrt(6.0));
215+
}
216+
}
217+
// Test TH2D projection for correct flow exclusion for specified user range
218+
TEST(Projections, UFOFWithRange)
219+
{
220+
TH2D h2("h2", "", 5, 0, 5, 4, 0, 4);
221+
h2.Sumw2();
222+
for (int bx = 0; bx <= 6; ++bx)
223+
for (int by = 0; by <= 5; ++by) {
224+
h2.SetBinContent(bx, by, 100 * bx + by);
225+
h2.SetBinError(bx, by, 2.0);
226+
}
227+
228+
h2.GetXaxis()->SetRange(2, 4);
229+
230+
auto hpx = h2.ProjectionX("hpx_ranged", 1, 4, "");
231+
EXPECT_EQ(hpx->GetXaxis()->GetNbins(), 3);
232+
for (int i = 0; i < 3; ++i) {
233+
int outbin = 2 + i;
234+
double exp = 0;
235+
for (int by = 1; by <= 4; ++by)
236+
exp += 100 * outbin + by;
237+
double got = hpx->GetBinContent(i + 1);
238+
EXPECT_DOUBLE_EQ(got, exp);
239+
EXPECT_DOUBLE_EQ(hpx->GetBinError(i + 1), 4.0);
240+
}
241+
EXPECT_DOUBLE_EQ(hpx->GetBinContent(0), 0);
242+
EXPECT_DOUBLE_EQ(hpx->GetBinContent(4), 0);
243+
}
244+
245+
// Test TH2D projection correctness with option "o"
246+
TEST(Projections, OriginalRange)
247+
{
248+
TH2D h2("h2", "", 6, 0, 6, 3, 0, 3);
249+
for (int bx = 0; bx <= 7; ++bx)
250+
for (int by = 0; by <= 4; ++by)
251+
h2.SetBinContent(bx, by, bx + 10 * by);
252+
253+
h2.GetXaxis()->SetRange(2, 5);
254+
255+
auto hpx = h2.ProjectionX("h_o", 1, 3, "o");
256+
EXPECT_EQ(hpx->GetXaxis()->GetNbins(), 6);
257+
for (int bx = 1; bx <= 6; ++bx) {
258+
double got = hpx->GetBinContent(bx);
259+
if (bx < 2 || bx > 5) {
260+
EXPECT_EQ(got, 0.0);
261+
continue;
262+
}
263+
double exp = 0;
264+
for (int by = 1; by <= 3; ++by)
265+
exp += bx + 10 * by;
266+
EXPECT_DOUBLE_EQ(got, exp);
267+
}
268+
}
269+
270+
// Test TH2D projection with variable bins and user range along projected axis
271+
TEST(Projections, VarBinsRange)
272+
{
273+
double edgesX[] = {0, 0.5, 1.5, 3.0, 5.0};
274+
TH2D h("h", "", 4, edgesX, 2, 0, 2);
275+
276+
for (int bx = 0; bx <= 5; ++bx)
277+
for (int by = 0; by <= 3; ++by)
278+
h.SetBinContent(bx, by, 100 * bx + by);
279+
280+
h.GetXaxis()->SetRange(2, 3);
281+
282+
auto hpx = h.ProjectionX("hpx_var", 1, 2, "");
283+
EXPECT_EQ(hpx->GetXaxis()->GetNbins(), 2);
284+
EXPECT_DOUBLE_EQ(hpx->GetXaxis()->GetBinLowEdge(1), edgesX[1]);
285+
EXPECT_DOUBLE_EQ(hpx->GetXaxis()->GetBinUpEdge(2), edgesX[3]);
286+
287+
for (int i = 0; i < 2; ++i) {
288+
int outbin = 2 + i;
289+
double exp = 0;
290+
for (int by = 1; by <= 2; ++by)
291+
exp += 100 * outbin + by;
292+
double got = hpx->GetBinContent(i + 1);
293+
EXPECT_DOUBLE_EQ(got, exp);
294+
}
295+
}
296+
297+
// https://github.com/root-project/issues/20174
298+
TEST(Projections, ProjectionYInfiniteUpperEdge)
299+
{
300+
Double_t xedges[] = {0., 1.};
301+
Double_t yedges[] = {1., std::numeric_limits<double>::infinity()};
302+
TH2D h("h_inf", "h_inf;X;Y", 1, xedges, 1, yedges);
303+
h.SetBinContent(1, 1, 11.);
304+
auto projY = h.ProjectionY();
305+
EXPECT_EQ(projY->GetBinContent(1), h.Integral(1, 1, 1, 1));
306+
}
307+
TEST(Projections, ProjectionYInfiniteLowerEdge)
308+
{
309+
Double_t xedges[] = {0, 1.};
310+
Double_t yedges[] = {-std::numeric_limits<double>::infinity(), 1};
311+
TH2D h("h_inf", "h_inf;X;Y", 1, xedges, 1, yedges);
312+
h.SetBinContent(1, 1, 11.);
313+
auto projY = h.ProjectionY();
314+
EXPECT_EQ(projY->GetBinContent(1), h.Integral(1, 1, 1, 1));
315+
}
316+
317+
TEST(Projections, Projection3DInfiniteEdges)
318+
{
319+
Double_t xedges[] = {0, 1.};
320+
Double_t yedges[] = {-std::numeric_limits<double>::infinity(), 1};
321+
Double_t zedges[] = {0, std::numeric_limits<double>::infinity()};
322+
TH3D h("h_inf", "h_inf;X;Y;Z", 1, xedges, 1, yedges, 1, zedges);
323+
h.SetBinContent(1, 1, 1, 11.);
324+
auto projY = h.Project3D("zy");
325+
EXPECT_EQ(projY->GetBinContent(projY->GetBin(1, 1)), h.Integral(1, 1, 1, 1, 1, 1));
326+
327+
auto projx = h.Project3D("x");
328+
EXPECT_EQ(projx->GetBinContent(projx->GetBin(1)), h.Integral(1, 1, 1, 1, 1, 1));
329+
}

0 commit comments

Comments
 (0)