diff --git a/hist/hist/src/TProfile3D.cxx b/hist/hist/src/TProfile3D.cxx
index 4756da07a4778..261d2f96a9b30 100644
--- a/hist/hist/src/TProfile3D.cxx
+++ b/hist/hist/src/TProfile3D.cxx
@@ -814,7 +814,8 @@ void TProfile3D::GetStats(Double_t *stats) const
    if (fBuffer) ((TProfile3D*)this)->BufferEmpty();
 
    // Loop on bins
-   if ( (fTsumw == 0 /* && fEntries > 0 */) || fXaxis.TestBit(TAxis::kAxisRange) || fYaxis.TestBit(TAxis::kAxisRange)) {
+   if ((fTsumw == 0 /* && fEntries > 0 */) || fXaxis.TestBit(TAxis::kAxisRange) || fYaxis.TestBit(TAxis::kAxisRange) ||
+       fZaxis.TestBit(TAxis::kAxisRange)) {
 
       // check for labels axis . In that case corresponding statistics do not make sense and it is set to zero
       Bool_t labelXaxis = ((const_cast<TAxis &>(fXaxis)).GetLabels() && fXaxis.CanExtend());
@@ -826,11 +827,38 @@ void TProfile3D::GetStats(Double_t *stats) const
       Double_t x,y,z;
       for (bin=0;bin<kNstat;bin++) stats[bin] = 0;
       if (!fBinEntries.fArray) return;
-      for (binz=fZaxis.GetFirst();binz<=fZaxis.GetLast();binz++) {
+      Int_t firstBinX = fXaxis.GetFirst();
+      Int_t lastBinX = fXaxis.GetLast();
+      Int_t firstBinY = fYaxis.GetFirst();
+      Int_t lastBinY = fYaxis.GetLast();
+      Int_t firstBinZ = fZaxis.GetFirst();
+      Int_t lastBinZ = fZaxis.GetLast();
+      // Include flow bins only on axes without an explicit range.
+      if (GetStatOverflowsBehaviour()) {
+         if (!fXaxis.TestBit(TAxis::kAxisRange)) {
+            if (firstBinX == 1)
+               firstBinX = 0;
+            if (lastBinX == fXaxis.GetNbins())
+               lastBinX += 1;
+         }
+         if (!fYaxis.TestBit(TAxis::kAxisRange)) {
+            if (firstBinY == 1)
+               firstBinY = 0;
+            if (lastBinY == fYaxis.GetNbins())
+               lastBinY += 1;
+         }
+         if (!fZaxis.TestBit(TAxis::kAxisRange)) {
+            if (firstBinZ == 1)
+               firstBinZ = 0;
+            if (lastBinZ == fZaxis.GetNbins())
+               lastBinZ += 1;
+         }
+      }
+      for (binz = firstBinZ; binz <= lastBinZ; binz++) {
          z = (!labelZaxis) ? fZaxis.GetBinCenter(binz) : 0;
-         for (biny=fYaxis.GetFirst();biny<=fYaxis.GetLast();biny++) {
+         for (biny = firstBinY; biny <= lastBinY; biny++) {
             y = (!labelYaxis) ? fYaxis.GetBinCenter(biny) : 0;
-            for (binx=fXaxis.GetFirst();binx<=fXaxis.GetLast();binx++) {
+            for (binx = firstBinX; binx <= lastBinX; binx++) {
                bin = GetBin(binx,biny,binz);
                w         = fBinEntries.fArray[bin];
                w2        = (fBinSumw2.fN ? fBinSumw2.fArray[bin] : w );
diff --git a/hist/hist/test/test_TProfile.cxx b/hist/hist/test/test_TProfile.cxx
index c08d0d2563608..b0ca19914a968 100644
--- a/hist/hist/test/test_TProfile.cxx
+++ b/hist/hist/test/test_TProfile.cxx
@@ -75,6 +75,32 @@ TEST(TProfile2D, Chi2Test)
    runTest(reference, sameDistr, differentDistr);
 }
 
+TEST(TProfile3D, GetStatsHonoursZRange)
+{
+   TProfile3D profile("profile", "", 4, 0, 4, 5, 0, 10, 3, 0, 30);
+   profile.Fill(0.5, 1., 5., 1.);
+   profile.Fill(3.5, 9., 25., 1.);
+   profile.GetZaxis()->SetRange(1, 1);
+
+   EXPECT_DOUBLE_EQ(0.5, profile.GetMean(1));
+   EXPECT_DOUBLE_EQ(1., profile.GetMean(2));
+   EXPECT_DOUBLE_EQ(5., profile.GetMean(3));
+}
+
+TEST(TProfile3D, GetStatsHonoursStatOverflowsWithZRange)
+{
+   TProfile3D profile("profile", "", 2, 0, 4, 3, 0, 12, 4, 0, 40);
+   profile.SetStatOverflows(TH1::kConsider);
+   profile.Fill(-1., 2., 5., 4.);
+   profile.Fill(1., 2., 5., 8.);
+   profile.GetZaxis()->SetRange(1, 1);
+
+   double stats[13] = {};
+   profile.GetStats(stats);
+   EXPECT_DOUBLE_EQ(2., stats[0]);
+   EXPECT_DOUBLE_EQ(0., profile.GetMean(1));
+}
+
 TEST(TProfile3D, Chi2Test)
 {
    TProfile3D reference("reference", "reference", 10, 0, 10, 11, 0, 11, 12, 0, 12);