Skip to content

[hist] Don't use FindBin/GetBinCenter in TH1 derived projections #20217

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
+311 −18
Open

[hist] Don't use FindBin/GetBinCenter in TH1 derived projections #20217

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 19 additions & 5 deletions hist/hist/src/TH2.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down
63 changes: 50 additions & 13 deletions hist/hist/src/TH3.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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);

Expand Down Expand Up @@ -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; }
Expand All @@ -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;
Expand All @@ -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;

}
}

Expand Down
242 changes: 242 additions & 0 deletions hist/hist/test/test_projections.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,8 @@
#include "TProfile.h" // ProjectionX
#include "TProfile2D.h" // ProjectionX
#include "THashList.h" // GetLabels
#include <limits>
#include "TMath.h"

#include "gtest/gtest.h"

Expand Down Expand Up @@ -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<double>::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<double>::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<double>::infinity(), 1};
Double_t zedges[] = {0, std::numeric_limits<double>::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));
}
Loading