[hist] Improve calculation of TAxis::FindBin(x) when x is at the bin edges

When x is at the bin edge values numerical error can cause the computation of the bin to be the one before (or after).
Correct for this case assuming the bin edge are computed as in `TH1::getBinLowEdge` and `TH1::GetBinUpEdge`. It is clear that the bin edge values are represented as floating point, so depending on how they are computed they can be sometime different.
However, it is better to have teh consistency to return the correct value when computed as internally in TAxis.

This should fix teh problem reported in https://root-forum.cern.ch/t/bug-in-taxis-findbin/57210 and #14091

Author: lmoneta

hist/hist/inc/TAxis.h

@@ -53,6 +53,7 @@ class TAxis : public TNamed, public TAttAxis {
    };
 
    Bool_t       HasBinWithoutLabel() const;
+   Int_t        DoFindFixBin(Double_t x) const;
 
 
    TAxisModLab *FindModLab(Int_t num, Double_t v = 0., Double_t eps = 0.) const;

hist/hist/src/TAxis.cxx

@@ -287,32 +287,25 @@ void TAxis::ExecuteEvent(Int_t event, Int_t px, Int_t py)
 
 Int_t TAxis::FindBin(Double_t x)
 {
-   Int_t bin;
    // NOTE: This should not be allowed for Alphanumeric histograms,
    // but it is heavily used (legacy) in the TTreePlayer to fill alphanumeric histograms.
    // but in case of alphanumeric do-not extend the axis. It makes no sense
    if (IsAlphanumeric() && gDebug) Info("FindBin","Numeric query on alphanumeric axis - Sorting the bins or extending the axes / rebinning can alter the correspondence between the label and the bin interval.");
    if (x < fXmin) {              //*-* underflow
-      bin = 0;
+      Int_t bin = 0;
       if (fParent == nullptr) return bin;
       if (!CanExtend() || IsAlphanumeric() ) return bin;
       ((TH1*)fParent)->ExtendAxis(x,this);
       return FindFixBin(x);
    } else  if ( !(x < fXmax)) {     //*-* overflow  (note the way to catch NaN)
-      bin = fNbins+1;
+      Int_t bin = fNbins+1;
       if (fParent == nullptr) return bin;
       if (!CanExtend() || IsAlphanumeric() ) return bin;
       ((TH1*)fParent)->ExtendAxis(x,this);
       return FindFixBin(x);
-   } else {
-      if (!fXbins.fN) {        //*-* fix bins
-         bin = 1 + int (fNbins*(x-fXmin)/(fXmax-fXmin) );
-      } else {                  //*-* variable bin sizes
-         //for (bin =1; x >= fXbins.fArray[bin]; bin++);
-         bin = 1 + TMath::BinarySearch(fXbins.fN,fXbins.fArray,x);
-      }
    }
-   return bin;
+   return DoFindFixBin(x);
+
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -404,6 +397,28 @@ Int_t TAxis::FindFixBin(const char *label) const
    return -1;
 }
 
+////////////////////////////////////////////////////////////////////////////////
+/// Internal function to find the fix bin
+////////////////////////////////////////////////////////////////////////////////
+Int_t TAxis::DoFindFixBin(Double_t x) const
+{
+   int bin = 0;
+   if (!fXbins.fN) {        //*-* fix bins
+      bin = 1 + int (fNbins * (x-fXmin)/(fXmax-fXmin) );
+      // apply correction when x  is at the bin edge value
+      // (computed as in GetBinLowEdge) a numerical error in the
+      // above division can cause a migration in a neighbour bin.
+      double binwidth = (fXmax - fXmin) / Double_t(fNbins);
+      double upEdge = fXmin + (bin) * binwidth;
+      double lowEdge = fXmin + (bin - 1) * binwidth;
+      if (upEdge <= x) bin += 1;
+      if (lowEdge > x) bin -= 1;
+   } else {                  //*-* variable bin sizes
+      //for (bin =1; x >= fXbins.fArray[bin]; bin++);
+      bin = 1 + TMath::BinarySearch(fXbins.fN,fXbins.fArray,x);
+   }
+   return bin;
+}
 
 ////////////////////////////////////////////////////////////////////////////////
 /// Find bin number corresponding to abscissa x
@@ -413,18 +428,13 @@ Int_t TAxis::FindFixBin(const char *label) const
 
 Int_t TAxis::FindFixBin(Double_t x) const
 {
-   Int_t bin;
+   Int_t bin = 0;
    if (x < fXmin) {              //*-* underflow
       bin = 0;
    } else  if ( !(x < fXmax)) {     //*-* overflow  (note the way to catch NaN
       bin = fNbins+1;
    } else {
-      if (!fXbins.fN) {        //*-* fix bins
-         bin = 1 + int (fNbins*(x-fXmin)/(fXmax-fXmin) );
-      } else {                  //*-* variable bin sizes
-//         for (bin =1; x >= fXbins.fArray[bin]; bin++);
-         bin = 1 + TMath::BinarySearch(fXbins.fN,fXbins.fArray,x);
-      }
+      bin = DoFindFixBin(x);
    }
    return bin;
 }

hist/hist/test/CMakeLists.txt

@@ -12,6 +12,7 @@ ROOT_ADD_GTEST(testTHn THn.cxx LIBRARIES Hist Matrix MathCore RIO)
 ROOT_ADD_GTEST(testTH1 test_TH1.cxx LIBRARIES Hist)
 ROOT_ADD_GTEST(testTHStack test_THStack.cxx LIBRARIES Hist)
 ROOT_ADD_GTEST(testProject3Dname test_Project3D_name.cxx LIBRARIES Hist)
+ROOT_ADD_GTEST(testTAxisFindBin test_TAxis_FindBin.cxx LIBRARIES Hist)
 ROOT_ADD_GTEST(testTFormula test_TFormula.cxx LIBRARIES Hist)
 ROOT_ADD_GTEST(testTKDE test_tkde.cxx LIBRARIES Hist)
 ROOT_ADD_GTEST(testTH1FindFirstBinAbove test_TH1_FindFirstBinAbove.cxx LIBRARIES Hist)

hist/hist/test/test_TAxis_FindBin.cxx

@@ -0,0 +1,34 @@
+#include "gtest/gtest.h"
+
+#include "TAxis.h"
+
+
+TEST(TAxis, FindBinExact)
+{
+   //Test the case where bin edges are exactly represented as floating points
+   TAxis ax(88, 1010, 1098);
+   for (int i = 1; i <= ax.GetNbins(); i++) {
+      double x = ax.GetBinLowEdge(i);
+      EXPECT_EQ(i, ax.FindBin(x));
+      EXPECT_EQ(i, ax.FindFixBin(x));
+      x = ax.GetBinUpEdge(i);
+      EXPECT_EQ(i+1, ax.FindBin(x));
+      EXPECT_EQ(i+1, ax.FindFixBin(x));
+      x -= x * std::numeric_limits<double>::epsilon();
+      EXPECT_EQ(i, ax.FindBin(x));
+   }
+}
+TEST(TAxis, FindBinApprox)
+{
+   TAxis ax(90, 0. , 10.);
+   for (int i = 1; i <= ax.GetNbins(); i++) {
+      double x = ax.GetBinLowEdge(i);
+      EXPECT_EQ(i, ax.FindBin(x));
+      EXPECT_EQ(i, ax.FindFixBin(x));
+      x = ax.GetBinUpEdge(i);
+      EXPECT_EQ(i+1, ax.FindBin(x));
+      EXPECT_EQ(i+1, ax.FindFixBin(x));
+      x -= x * std::numeric_limits<double>::epsilon();
+      EXPECT_EQ(i, ax.FindBin(x));
+   }
+}