From 2b5f5a96c7773defc8d27881d3aab9dced2b9e60 Mon Sep 17 00:00:00 2001 From: ferdymercury Date: Thu, 11 Apr 2024 23:11:49 +0200 Subject: [PATCH] [hist] THn: Add GetBinCenter and Integral helper methods (#15200) * [hist] THn: Add GetBinCenter and Integral helper methods Fixes https://github.com/root-project/root/issues/14173 * [hist] modernize C++ code thanks! Co-authored-by: Vincenzo Eduardo Padulano [hist,skip-ci] highlight difference wrt ComputeIntegral * [hist] add test functionality for THn Integral and GetBinCenter --- hist/hist/inc/THnBase.h | 3 +++ hist/hist/src/THnBase.cxx | 41 +++++++++++++++++++++++++++++++++++++++ hist/hist/test/THn.cxx | 33 +++++++++++++++++++++++++++++++ 3 files changed, 77 insertions(+) diff --git a/hist/hist/inc/THnBase.h b/hist/hist/inc/THnBase.h index dffcf4101cb71..47ac7ec4227c2 100644 --- a/hist/hist/inc/THnBase.h +++ b/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 GetBinCenter(const std::vector &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; diff --git a/hist/hist/src/THnBase.cxx b/hist/hist/src/THnBase.cxx index 92b4e5c6365b8..3e59a0d3fa10a 100644 --- a/hist/hist/src/THnBase.cxx +++ b/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 THnBase::GetBinCenter(const std::vector &idx) const +{ + if (idx.size() != static_cast(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 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 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. diff --git a/hist/hist/test/THn.cxx b/hist/hist/test/THn.cxx index 43d4ef121a3f6..f7c633aac9550 100644 --- a/hist/hist/test/THn.cxx +++ b/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); +}