[hist] THn: Add GetBinCenter and Integral helper methods (#15200)

* [hist] THn: Add GetBinCenter and Integral helper methods

Fixes #14173

* [hist] modernize C++ code

thanks!

Co-authored-by: Vincenzo Eduardo Padulano <[email protected]>

[hist,skip-ci] highlight difference wrt ComputeIntegral

* [hist] add test functionality for THn Integral and GetBinCenter

hist/hist/inc/THnBase.h

@@ -192,6 +192,8 @@ class THnBase: public TNamed {
    void SetEntries(Double_t entries) { fEntries = entries; }
    void SetTitle(const char *title) override;
 
+   std::vector<Double_t> GetBinCenter(const std::vector<Int_t> &idx) const;
+
    Double_t GetBinContent(const Int_t *idx) const { return GetBinContent(GetBin(idx)); } // intentionally non-virtual
    virtual Double_t GetBinContent(Long64_t bin, Int_t* idx = nullptr) const = 0;
    virtual Double_t GetBinError2(Long64_t linidx) const = 0;
@@ -273,6 +275,7 @@ class THnBase: public TNamed {
 
    Double_t ComputeIntegral();
    void GetRandom(Double_t *rand, Bool_t subBinRandom = kTRUE);
+   Double_t Integral(Bool_t respectAxisRange) const;
 
    void Print(Option_t* option = "") const override;
    void PrintEntries(Long64_t from = 0, Long64_t howmany = -1, Option_t* options = nullptr) const;

hist/hist/src/THnBase.cxx

@@ -525,6 +525,29 @@ TFitResultPtr THnBase::Fit(TF1 *f ,Option_t *option ,Option_t *goption)
    return ROOT::Fit::FitObject(this, f , fitOption , minOption, goption, range);
 }
 
+////////////////////////////////////////////////////////////////////////////////
+/// \brief THnBase::GetBinCenter
+/// \param idx an array of bin index in each dimension.
+/// \return vector of bin centers in each dimension; empty in case of error.
+/// \note Throws error if size is different from nDimensions.
+/// \sa GetAxis(dim)::GetBinCenter(idx) as an alternative
+std::vector<Double_t> THnBase::GetBinCenter(const std::vector<Int_t> &idx) const
+{
+   if (idx.size() != static_cast<decltype(idx.size())>(fNdimensions)) {
+      Error("THnBase::GetBinCenter",
+            "Mismatched number of dimensions %d with bin index vector size %zu, returning an empty vector.",
+            fNdimensions, idx.size());
+      return {};
+   }
+   std::vector<Double_t> centers(fNdimensions);
+   std::generate(centers.begin(), centers.end(), [i = 0, &idx, this]() mutable {
+      auto bincenter = GetAxis(i)->GetBinCenter(idx[i]);
+      i++;
+      return bincenter;
+   });
+   return centers;
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 /// Generate an n-dimensional random tuple based on the histogrammed
 /// distribution. If subBinRandom, the returned tuple will be additionally
@@ -1312,6 +1335,24 @@ void THnBase::ResetBase(Option_t * /*option = ""*/)
    }
 }
 
+////////////////////////////////////////////////////////////////////////////////
+/// \brief Compute integral (sum of counts) of histogram in all dimensions
+/// \param respectAxisRange if false, count all bins including under/overflows,
+///                         if true, restrict sum to the user-set axis range
+/// \sa Projection(0)::Integral() as alternative
+/// \note this function is different from ComputeIntegral, that is a normalized
+/// cumulative sum
+Double_t THnBase::Integral(Bool_t respectAxisRange) const
+{
+   Long64_t myLinBin = 0;
+   std::unique_ptr<ROOT::Internal::THnBaseBinIter> iter{CreateIter(respectAxisRange)};
+   Double_t sum = 0.;
+   while ((myLinBin = iter->Next()) >= 0) {
+      sum += GetBinContent(myLinBin);
+   }
+   return sum;
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 ///  Compute integral (normalized cumulative sum of bins) w/o under/overflows
 ///  The result is stored in fIntegral and used by the GetRandom functions.

hist/hist/test/THn.cxx

@@ -111,3 +111,36 @@ TEST(THn, Projection) {
    }
 
 }
+
+TEST(THn, Integral)
+{
+   Int_t bins[2] = {2, 3};
+   Double_t xmin[2] = {0., -3.};
+   Double_t xmax[2] = {10., 3.};
+   THnD hn("hn", "hn", 2, bins, xmin, xmax);
+   Double_t x[2];
+   x[0] = 4;
+   x[1] = -0.1;
+   hn.Fill(x);
+   x[0] = 7;
+   x[1] = 1.0;
+   hn.Fill(x);
+   x[0] = 1;
+   x[1] = 2.6;
+   hn.Fill(x);
+   x[0] = 9;
+   x[1] = -0.9;
+   hn.Fill(x);
+   EXPECT_DOUBLE_EQ(hn.Integral(false), 4);
+}
+
+TEST(THn, GetBinCenter)
+{
+   Int_t bins[2] = {2, 2};
+   Double_t xmin[2] = {0., -3.};
+   Double_t xmax[2] = {10., 3.};
+   THnD hn("hn", "hn", 2, bins, xmin, xmax);
+   auto centers = hn.GetBinCenter({1, 1});
+   EXPECT_DOUBLE_EQ(centers.at(0), 2.5);
+   EXPECT_DOUBLE_EQ(centers.at(1), -1.5);
+}