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 <limits>
+#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<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));
+}