From bfad5a6a27acd16ba37ac49d5c0ea7df91de235d Mon Sep 17 00:00:00 2001 From: Jonas Hahnfeld Date: Tue, 4 Aug 2026 09:28:41 +0200 Subject: [PATCH] [hist] Conversion functions to TH3* --- hist/histv7util/CMakeLists.txt | 2 + .../histv7util/inc/ROOT/Hist/ConvertToTH3.hxx | 120 ++++++ hist/histv7util/src/ConvertToTH3.cxx | 246 +++++++++++ hist/histv7util/test/CMakeLists.txt | 1 + hist/histv7util/test/hist_convert_TH3.cxx | 389 ++++++++++++++++++ hist/histv7util/test/histutil_test.hxx | 1 + 6 files changed, 759 insertions(+) create mode 100644 hist/histv7util/inc/ROOT/Hist/ConvertToTH3.hxx create mode 100644 hist/histv7util/src/ConvertToTH3.cxx create mode 100644 hist/histv7util/test/hist_convert_TH3.cxx diff --git a/hist/histv7util/CMakeLists.txt b/hist/histv7util/CMakeLists.txt index 66ea3f74ca1a5..c5c8977450ca8 100644 --- a/hist/histv7util/CMakeLists.txt +++ b/hist/histv7util/CMakeLists.txt @@ -3,10 +3,12 @@ ROOT_STANDARD_LIBRARY_PACKAGE(ROOTHistUtil ROOT/Hist/ConversionUtils.hxx ROOT/Hist/ConvertToTH1.hxx ROOT/Hist/ConvertToTH2.hxx + ROOT/Hist/ConvertToTH3.hxx SOURCES ConversionUtils.cxx ConvertToTH1.cxx ConvertToTH2.cxx + ConvertToTH3.cxx DEPENDENCIES Core Hist diff --git a/hist/histv7util/inc/ROOT/Hist/ConvertToTH3.hxx b/hist/histv7util/inc/ROOT/Hist/ConvertToTH3.hxx new file mode 100644 index 0000000000000..9ec73a739e520 --- /dev/null +++ b/hist/histv7util/inc/ROOT/Hist/ConvertToTH3.hxx @@ -0,0 +1,120 @@ +/// \file +/// \warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes. +/// Feedback is welcome! + +#ifndef ROOT_Hist_ConvertToTH3 +#define ROOT_Hist_ConvertToTH3 + +#include +#include +#include + +class TH3C; +class TH3S; +class TH3I; +class TH3L; +class TH3F; +class TH3D; + +#include + +namespace ROOT { +namespace Experimental { +namespace Hist { + +/// Convert a three-dimensional histogram to TH3C. +/// +/// \copydetails ConvertToTH3I(const RHistEngine &engine) +std::unique_ptr ConvertToTH3C(const RHistEngine &engine); + +/// Convert a three-dimensional histogram to TH3S. +/// +/// \copydetails ConvertToTH3I(const RHistEngine &engine) +std::unique_ptr ConvertToTH3S(const RHistEngine &engine); + +/// Convert a three-dimensional histogram to TH3I. +/// +/// As RHistEngine does not have global statistics, the number of entries and the total sum of weights will be unset. +/// +/// Throws an exception if the histogram has more than one dimension. +/// +/// \param[in] engine the RHistEngine to convert +/// \return the converted TH3 +std::unique_ptr ConvertToTH3I(const RHistEngine &engine); + +/// Convert a three-dimensional histogram to TH3L. +/// +/// \copydetails ConvertToTH3I(const RHistEngine &engine) +std::unique_ptr ConvertToTH3L(const RHistEngine &engine); + +/// Convert a three-dimensional histogram to TH3L. +/// +/// \copydetails ConvertToTH3I(const RHistEngine &engine) +std::unique_ptr ConvertToTH3L(const RHistEngine &engine); + +/// Convert a three-dimensional histogram to TH3F. +/// +/// \copydetails ConvertToTH3I(const RHistEngine &engine) +std::unique_ptr ConvertToTH3F(const RHistEngine &engine); + +/// Convert a three-dimensional histogram to TH3D. +/// +/// \copydetails ConvertToTH3I(const RHistEngine &engine) +std::unique_ptr ConvertToTH3D(const RHistEngine &engine); + +/// Convert a three-dimensional histogram to TH3D. +/// +/// \copydetails ConvertToTH3I(const RHistEngine &engine) +std::unique_ptr ConvertToTH3D(const RHistEngine &engine); + +/// Convert a three-dimensional histogram to TH3C. +/// +/// \copydetails ConvertToTH3I(const RHist &hist) +std::unique_ptr ConvertToTH3C(const RHist &hist); + +/// Convert a three-dimensional histogram to TH3S. +/// +/// \copydetails ConvertToTH3I(const RHist &hist) +std::unique_ptr ConvertToTH3S(const RHist &hist); + +/// Convert a three-dimensional histogram to TH3I. +/// +/// If the RHistStats are tainted, for example after setting bin contents, the number of entries and the total sum of +/// weights will be unset. +/// +/// Throws an exception if the histogram has more than one dimension. +/// +/// \param[in] hist the RHist to convert +/// \return the converted TH3 +std::unique_ptr ConvertToTH3I(const RHist &hist); + +/// Convert a three-dimensional histogram to TH3L. +/// +/// \copydetails ConvertToTH3I(const RHist &hist) +std::unique_ptr ConvertToTH3L(const RHist &hist); + +/// Convert a three-dimensional histogram to TH3L. +/// +/// \copydetails ConvertToTH3I(const RHist &hist) +std::unique_ptr ConvertToTH3L(const RHist &hist); + +/// Convert a three-dimensional histogram to TH3F. +/// +/// \copydetails ConvertToTH3I(const RHist &hist) +std::unique_ptr ConvertToTH3F(const RHist &hist); + +/// Convert a three-dimensional histogram to TH3D. +/// +/// \copydetails ConvertToTH3I(const RHist &hist) +std::unique_ptr ConvertToTH3D(const RHist &hist); + +/// Convert a three-dimensional histogram to TH3D. +/// +/// \copydetails ConvertToTH3I(const RHist &hist) +std::unique_ptr ConvertToTH3D(const RHist &hist); + +} // namespace Hist +} // namespace Experimental +} // namespace ROOT + +#endif diff --git a/hist/histv7util/src/ConvertToTH3.cxx b/hist/histv7util/src/ConvertToTH3.cxx new file mode 100644 index 0000000000000..1bbaa2c3ffc87 --- /dev/null +++ b/hist/histv7util/src/ConvertToTH3.cxx @@ -0,0 +1,246 @@ +/// \file +/// \warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes. +/// Feedback is welcome! + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +using namespace ROOT::Experimental; + +namespace { +template +std::unique_ptr ConvertToTH3Impl(const RHistEngine &engine) +{ + if (engine.GetNDimensions() != 3) { + throw std::invalid_argument("TH3 requires three dimensions"); + } + + auto ret = std::make_unique(); + ret->SetDirectory(nullptr); + + const auto &axis0 = engine.GetAxes()[0]; + ROOT::Experimental::Hist::Internal::ConvertAxis(*ret->GetXaxis(), axis0); + const auto &axis1 = engine.GetAxes()[1]; + ROOT::Experimental::Hist::Internal::ConvertAxis(*ret->GetYaxis(), axis1); + const auto &axis2 = engine.GetAxes()[2]; + ROOT::Experimental::Hist::Internal::ConvertAxis(*ret->GetZaxis(), axis2); + ret->SetBinsLength(); + + Double_t *sumw2 = nullptr; + auto copyBinContent = [&ret, &engine, &sumw2](Int_t i, RBinIndex index0, RBinIndex index1, RBinIndex index2) { + if constexpr (std::is_same_v) { + if (sumw2 == nullptr) { + ret->Sumw2(); + sumw2 = ret->GetSumw2()->GetArray(); + } + const RBinWithError &c = engine.GetBinContent(index0, index1, index2); + ret->GetArray()[i] = c.fSum; + sumw2[i] = c.fSum2; + } else { + (void)sumw2; + ret->GetArray()[i] = engine.GetBinContent(index0, index1, index2); + } + }; + + // Copy the bin contents, accounting for TH3 numbering conventions. + for (auto index0 : axis0.GetFullRange()) { + Int_t i0 = 0; + if (index0.IsUnderflow()) { + i0 = 0; + } else if (index0.IsOverflow()) { + i0 = axis0.GetNNormalBins() + 1; + } else { + assert(index0.IsNormal()); + i0 = index0.GetIndex() + 1; + } + Int_t n0 = ret->GetXaxis()->GetNbins() + 2; + + for (auto index1 : axis1.GetFullRange()) { + Int_t i1 = 0; + if (index1.IsUnderflow()) { + i1 = 0; + } else if (index1.IsOverflow()) { + i1 = axis1.GetNNormalBins() + 1; + } else { + assert(index1.IsNormal()); + i1 = index1.GetIndex() + 1; + } + Int_t n1 = ret->GetYaxis()->GetNbins() + 2; + + for (auto index2 : axis2.GetFullRange()) { + Int_t i2 = 0; + if (index2.IsUnderflow()) { + i2 = 0; + } else if (index2.IsOverflow()) { + i2 = axis2.GetNNormalBins() + 1; + } else { + assert(index2.IsNormal()); + i2 = index2.GetIndex() + 1; + } + + Int_t i = i0 + n0 * (i1 + n1 * i2); + copyBinContent(i, index0, index1, index2); + } + } + } + + return ret; +} + +template +void ConvertGlobalStatistics(Hist &h, const RHistStats &stats) +{ + if (stats.IsTainted()) { + return; + } + + h.SetEntries(stats.GetNEntries()); + + Double_t hStats[11] = { + stats.GetSumW(), + stats.GetSumW2(), + 0, + 0, + 0, + 0, + // We do not have sumwxy + 0, + 0, + 0, + // We do not have sumwxz + 0, + // We do not have sumwyz + 0, + }; + if (stats.IsEnabled(0)) { + hStats[2] = stats.GetDimensionStats(0).fSumWX; + hStats[3] = stats.GetDimensionStats(0).fSumWX2; + } + if (stats.IsEnabled(1)) { + hStats[4] = stats.GetDimensionStats(1).fSumWX; + hStats[5] = stats.GetDimensionStats(1).fSumWX2; + } + if (stats.IsEnabled(2)) { + hStats[7] = stats.GetDimensionStats(2).fSumWX; + hStats[8] = stats.GetDimensionStats(2).fSumWX2; + } + h.PutStats(hStats); +} +} // namespace + +namespace ROOT { +namespace Experimental { +namespace Hist { + +std::unique_ptr ConvertToTH3C(const RHistEngine &engine) +{ + return ConvertToTH3Impl(engine); +} + +std::unique_ptr ConvertToTH3S(const RHistEngine &engine) +{ + return ConvertToTH3Impl(engine); +} + +std::unique_ptr ConvertToTH3I(const RHistEngine &engine) +{ + return ConvertToTH3Impl(engine); +} + +std::unique_ptr ConvertToTH3L(const RHistEngine &engine) +{ + return ConvertToTH3Impl(engine); +} + +std::unique_ptr ConvertToTH3L(const RHistEngine &engine) +{ + return ConvertToTH3Impl(engine); +} + +std::unique_ptr ConvertToTH3F(const RHistEngine &engine) +{ + return ConvertToTH3Impl(engine); +} + +std::unique_ptr ConvertToTH3D(const RHistEngine &engine) +{ + return ConvertToTH3Impl(engine); +} + +std::unique_ptr ConvertToTH3D(const RHistEngine &engine) +{ + return ConvertToTH3Impl(engine); +} + +std::unique_ptr ConvertToTH3C(const RHist &hist) +{ + auto ret = ConvertToTH3C(hist.GetEngine()); + ConvertGlobalStatistics(*ret, hist.GetStats()); + return ret; +} + +std::unique_ptr ConvertToTH3S(const RHist &hist) +{ + auto ret = ConvertToTH3S(hist.GetEngine()); + ConvertGlobalStatistics(*ret, hist.GetStats()); + return ret; +} + +std::unique_ptr ConvertToTH3I(const RHist &hist) +{ + auto ret = ConvertToTH3I(hist.GetEngine()); + ConvertGlobalStatistics(*ret, hist.GetStats()); + return ret; +} + +std::unique_ptr ConvertToTH3L(const RHist &hist) +{ + auto ret = ConvertToTH3L(hist.GetEngine()); + ConvertGlobalStatistics(*ret, hist.GetStats()); + return ret; +} + +std::unique_ptr ConvertToTH3L(const RHist &hist) +{ + auto ret = ConvertToTH3L(hist.GetEngine()); + ConvertGlobalStatistics(*ret, hist.GetStats()); + return ret; +} + +std::unique_ptr ConvertToTH3F(const RHist &hist) +{ + auto ret = ConvertToTH3F(hist.GetEngine()); + ConvertGlobalStatistics(*ret, hist.GetStats()); + return ret; +} + +std::unique_ptr ConvertToTH3D(const RHist &hist) +{ + auto ret = ConvertToTH3D(hist.GetEngine()); + ConvertGlobalStatistics(*ret, hist.GetStats()); + return ret; +} + +std::unique_ptr ConvertToTH3D(const RHist &hist) +{ + auto ret = ConvertToTH3D(hist.GetEngine()); + ConvertGlobalStatistics(*ret, hist.GetStats()); + return ret; +} + +} // namespace Hist +} // namespace Experimental +} // namespace ROOT diff --git a/hist/histv7util/test/CMakeLists.txt b/hist/histv7util/test/CMakeLists.txt index 089c556239280..7e3ac21b1dbd3 100644 --- a/hist/histv7util/test/CMakeLists.txt +++ b/hist/histv7util/test/CMakeLists.txt @@ -1,3 +1,4 @@ ROOT_ADD_GTEST(hist_conversion_utils hist_conversion_utils.cxx LIBRARIES ROOTHist ROOTHistUtil) ROOT_ADD_GTEST(hist_convert_TH1 hist_convert_TH1.cxx LIBRARIES ROOTHist ROOTHistUtil) ROOT_ADD_GTEST(hist_convert_TH2 hist_convert_TH2.cxx LIBRARIES ROOTHist ROOTHistUtil) +ROOT_ADD_GTEST(hist_convert_TH3 hist_convert_TH3.cxx LIBRARIES ROOTHist ROOTHistUtil) diff --git a/hist/histv7util/test/hist_convert_TH3.cxx b/hist/histv7util/test/hist_convert_TH3.cxx new file mode 100644 index 0000000000000..886c330b8ae42 --- /dev/null +++ b/hist/histv7util/test/hist_convert_TH3.cxx @@ -0,0 +1,389 @@ +#include "histutil_test.hxx" + +#include + +#include +#include +#include +#include + +TEST(ConvertToTH3I, RHistEngine) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + RHistEngine engine(axis, axis, axis); + + engine.SetBinContent(RBinIndex::Underflow(), 0, 1, 1000); + engine.SetBinContent(RBinIndex::Underflow(), 2, 3, 2000); + for (std::size_t i = 0; i < Bins; i++) { + engine.SetBinContent(i, RBinIndex::Underflow(), 1, 100 * i); + for (std::size_t j = 0; j < Bins; j++) { + engine.SetBinContent(i, j, RBinIndex::Underflow(), 400 * j); + for (std::size_t k = 0; k < Bins; k++) { + engine.SetBinContent(i, j, k, (i * Bins + j) * Bins + k); + } + engine.SetBinContent(i, j, RBinIndex::Overflow(), 800 * j); + } + engine.SetBinContent(i, RBinIndex::Overflow(), 4, 200 * i); + } + engine.SetBinContent(RBinIndex::Overflow(), 4, 5, 3000); + engine.SetBinContent(RBinIndex::Overflow(), 6, 7, 4000); + + auto th3i = ConvertToTH3I(engine); + ASSERT_TRUE(th3i); + EXPECT_TRUE(th3i->GetDirectory() == nullptr); + ASSERT_EQ(th3i->GetDimension(), 3); + ASSERT_EQ(th3i->GetNbinsX(), Bins); + ASSERT_EQ(th3i->GetNbinsY(), Bins); + ASSERT_EQ(th3i->GetNbinsZ(), Bins); + + EXPECT_EQ(th3i->GetBinContent(0, 1, 2), 1000); + EXPECT_EQ(th3i->GetBinContent(0, 3, 4), 2000); + for (std::size_t i = 0; i < Bins; i++) { + EXPECT_EQ(th3i->GetBinContent(i + 1, 0, 2), 100 * i); + for (std::size_t j = 0; j < Bins; j++) { + EXPECT_EQ(th3i->GetBinContent(i + 1, j + 1, 0), 400 * j); + for (std::size_t k = 0; k < Bins; k++) { + EXPECT_EQ(th3i->GetBinContent(i + 1, j + 1, k + 1), (i * Bins + j) * Bins + k); + } + EXPECT_EQ(th3i->GetBinContent(i + 1, j + 1, Bins + 1), 800 * j); + } + EXPECT_EQ(th3i->GetBinContent(i + 1, Bins + 1, 5), 200 * i); + } + EXPECT_EQ(th3i->GetBinContent(Bins + 1, 5, 6), 3000); + EXPECT_EQ(th3i->GetBinContent(Bins + 1, 7, 8), 4000); + + EXPECT_EQ(th3i->GetEntries(), 0); + Double_t stats[11]; + th3i->GetStats(stats); + for (std::size_t i = 0; i < 11; i++) { + EXPECT_EQ(stats[i], 0); + } +} + +TEST(ConvertToTH3I, RHistEngineNoFlowBins) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}, /*enableFlowBins=*/false); + RHistEngine engine(axis, axis, axis); + + // Flow bins are disabled, so these fills will be discarded. + engine.Fill(-100, 0.5, 1.5); + engine.Fill(-100, 2.5, 3.5); + for (std::size_t i = 0; i < Bins; i++) { + // Flow bins are disabled, so this fill will be discarded. + engine.Fill(i + 0.5, -100, 1.5); + for (std::size_t j = 0; j < Bins; j++) { + // Flow bins are disabled, so this fill will be discarded. + engine.Fill(i + 0.5, j + 0.5, -100); + for (std::size_t k = 0; k < Bins; k++) { + engine.SetBinContent(i, j, k, (i * Bins + j) * Bins + k); + } + // Flow bins are disabled, so this fill will be discarded. + engine.Fill(i + 0.5, j + 0.5, 100); + } + // Flow bins are disabled, so this fill will be discarded. + engine.Fill(i + 0.5, 100, 4.5); + } + // Flow bins are disabled, so these fills will be discarded. + engine.Fill(100, 4.5, 5.5); + engine.Fill(100, 6.5, 7.5); + + auto th3i = ConvertToTH3I(engine); + ASSERT_TRUE(th3i); + + EXPECT_EQ(th3i->GetBinContent(0, 1, 2), 0); + EXPECT_EQ(th3i->GetBinContent(0, 3, 4), 0); + for (std::size_t i = 0; i < Bins; i++) { + EXPECT_EQ(th3i->GetBinContent(i + 1, 0, 1), 0); + for (std::size_t j = 0; j < Bins; j++) { + EXPECT_EQ(th3i->GetBinContent(i + 1, j + 1, 0), 0); + for (std::size_t k = 0; k < Bins; k++) { + EXPECT_EQ(th3i->GetBinContent(i + 1, j + 1, k + 1), (i * Bins + j) * Bins + k); + } + EXPECT_EQ(th3i->GetBinContent(i + 1, j + 1, Bins + 1), 0); + } + EXPECT_EQ(th3i->GetBinContent(i + 1, Bins + 1, 4), 0); + } + EXPECT_EQ(th3i->GetBinContent(Bins + 1, 5, 6), 0); + EXPECT_EQ(th3i->GetBinContent(Bins + 1, 7, 8), 0); +} + +TEST(ConvertToTH3I, RHistEngineInvalid) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + const RHistEngine engine(axis); + + EXPECT_THROW(ConvertToTH3I(engine), std::invalid_argument); +} + +TEST(ConvertToTH3I, RHist) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + RHist hist(axis, axis, axis); + + for (std::size_t i = 0; i < Bins; i++) { + hist.Fill(i, 2 * i, 3 * i); + } + + auto th3i = ConvertToTH3I(hist); + ASSERT_TRUE(th3i); + + ASSERT_EQ(hist.GetNEntries(), Bins); + EXPECT_EQ(th3i->GetEntries(), Bins); + Double_t stats[11]; + th3i->GetStats(stats); + EXPECT_EQ(stats[0], Bins); + EXPECT_EQ(stats[1], Bins); + EXPECT_EQ(stats[2], 190); + EXPECT_EQ(stats[3], 2470); + EXPECT_EQ(stats[4], 2 * 190); + EXPECT_EQ(stats[5], 4 * 2470); + EXPECT_EQ(stats[6], 0); + EXPECT_EQ(stats[7], 3 * 190); + EXPECT_EQ(stats[8], 9 * 2470); + EXPECT_EQ(stats[9], 0); + EXPECT_EQ(stats[10], 0); +} + +TEST(ConvertToTH3I, RHistSetBinContentTainted) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + RHist hist(axis, axis, axis); + const std::array indices = {1, 2, 3}; + hist.SetBinContent(indices, 42); + ASSERT_TRUE(hist.GetStats().IsTainted()); + + auto th3i = ConvertToTH3I(hist); + ASSERT_TRUE(th3i); + + EXPECT_EQ(th3i->GetBinContent(2 + (Bins + 2) * (3 + (Bins + 2) * 4)), 42); + + EXPECT_EQ(th3i->GetEntries(), 0); + Double_t stats[11]; + th3i->GetStats(stats); + for (std::size_t i = 0; i < 11; i++) { + EXPECT_EQ(stats[i], 0); + } +} + +TEST(ConvertToTH3I, RHistCategoricalAxis) +{ + const std::vector categories = {"a", "b", "c"}; + const RCategoricalAxis axis(categories); + RHist hist(axis, axis, axis); + ASSERT_FALSE(hist.GetStats().IsEnabled(0)); + ASSERT_FALSE(hist.GetStats().IsEnabled(1)); + ASSERT_FALSE(hist.GetStats().IsEnabled(2)); + + hist.Fill("a", "b", "c"); + + auto th3i = ConvertToTH3I(hist); + ASSERT_TRUE(th3i); + + EXPECT_EQ(th3i->GetBinContent(1 + 5 * (2 + 5 * 3)), 1); + + EXPECT_EQ(th3i->GetEntries(), 1); + Double_t stats[11]; + th3i->GetStats(stats); + EXPECT_EQ(stats[0], 1); + EXPECT_EQ(stats[1], 1); + for (std::size_t i = 2; i < 11; i++) { + EXPECT_EQ(stats[i], 0); + } +} + +TEST(ConvertToTH3C, RHistEngine) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + const RHistEngine engine(axis, axis, axis); + + auto th3c = ConvertToTH3C(engine); + ASSERT_TRUE(th3c); +} + +TEST(ConvertToTH3C, RHist) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + const RHist hist(axis, axis, axis); + + auto th3c = ConvertToTH3C(hist); + ASSERT_TRUE(th3c); +} + +TEST(ConvertToTH3S, RHistEngine) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + const RHistEngine engine(axis, axis, axis); + + auto th3s = ConvertToTH3S(engine); + ASSERT_TRUE(th3s); +} + +TEST(ConvertToTH3S, RHist) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + const RHist hist(axis, axis, axis); + + auto th3s = ConvertToTH3S(hist); + ASSERT_TRUE(th3s); +} + +TEST(ConvertToTH3L, RHistEngine) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + const RHistEngine engineL(axis, axis, axis); + + auto th3l = ConvertToTH3L(engineL); + ASSERT_TRUE(th3l); + + RHistEngine engineLL(axis, axis, axis); + + // Set one 64-bit long long value larger than what double can exactly represent. + static constexpr long long Large = (1LL << 60) - 1; + const std::array indices = {1, 2, 3}; + engineLL.SetBinContent(indices, Large); + + th3l = ConvertToTH3L(engineLL); + ASSERT_TRUE(th3l); + + // Get the value via TArrayL::At and store into a variable to be sure about the type. During direct comparison, a + // double return value may automatically promote Large to a double as well, introducing the truncation we want to + // test against. + const long long value = th3l->At(2 + (Bins + 2) * (3 + (Bins + 2) * 4)); + EXPECT_EQ(value, Large); +} + +TEST(ConvertToTH3L, RHist) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + const RHist histL(axis, axis, axis); + + auto th3l = ConvertToTH3L(histL); + ASSERT_TRUE(th3l); + + const RHist histLL(axis, axis, axis); + th3l = ConvertToTH3L(histLL); + ASSERT_TRUE(th3l); +} + +TEST(ConvertToTH3F, RHistEngine) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + RHistEngine engine(axis, axis, axis); + + engine.Fill(-100, 0.5, 1.5, RWeight(0.25)); + for (std::size_t i = 0; i < Bins / 3; i++) { + engine.Fill(i, 2 * i, 3 * i, RWeight(0.1 + i * 0.03)); + } + engine.Fill(100, Bins - 0.5, Bins - 1.5, RWeight(0.75)); + + auto th3f = ConvertToTH3F(engine); + ASSERT_TRUE(th3f); + + EXPECT_FLOAT_EQ(th3f->GetBinContent(0, 1, 2), 0.25); + for (std::size_t i = 0; i < Bins / 3; i++) { + EXPECT_FLOAT_EQ(th3f->GetBinContent(i + 1, 2 * i + 1, 3 * i + 1), 0.1 + i * 0.03); + } + EXPECT_EQ(th3f->GetBinContent(Bins + 1, Bins, Bins - 1), 0.75); +} + +TEST(ConvertToTH3F, RHist) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + RHist hist(axis, axis, axis); + + for (std::size_t i = 0; i < Bins; i++) { + hist.Fill(i, 2 * i, 3 * i, RWeight(0.1 + i * 0.03)); + } + + auto th3f = ConvertToTH3F(hist); + ASSERT_TRUE(th3f); + + ASSERT_EQ(hist.GetNEntries(), Bins); + EXPECT_EQ(th3f->GetEntries(), Bins); + Double_t stats[11]; + th3f->GetStats(stats); + EXPECT_DOUBLE_EQ(stats[0], 7.7); + EXPECT_DOUBLE_EQ(stats[1], 3.563); + EXPECT_DOUBLE_EQ(stats[2], 93.1); + EXPECT_DOUBLE_EQ(stats[3], 1330.0); + EXPECT_DOUBLE_EQ(stats[4], 2 * 93.1); + EXPECT_DOUBLE_EQ(stats[5], 4 * 1330.0); + EXPECT_DOUBLE_EQ(stats[6], 0); + EXPECT_DOUBLE_EQ(stats[7], 3 * 93.1); + EXPECT_DOUBLE_EQ(stats[8], 9 * 1330.0); + EXPECT_DOUBLE_EQ(stats[9], 0); + EXPECT_DOUBLE_EQ(stats[10], 0); +} + +TEST(ConvertToTH3D, RHistEngine) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + const RHistEngine engineD(axis, axis, axis); + + auto th3d = ConvertToTH3D(engineD); + ASSERT_TRUE(th3d); + + RHistEngine engineE(axis, axis, axis); + for (std::size_t i = 0; i < Bins / 3; i++) { + engineE.Fill(i, 2 * i, 3 * i, RWeight(0.1 + i * 0.03)); + } + + th3d = ConvertToTH3D(engineE); + ASSERT_TRUE(th3d); + const Double_t *sumw2 = th3d->GetSumw2()->GetArray(); + ASSERT_TRUE(sumw2 != nullptr); + + for (std::size_t i = 0; i < Bins / 3; i++) { + const double weight = 0.1 + i * 0.03; + EXPECT_EQ(th3d->GetBinContent(i + 1, 2 * i + 1, 3 * i + 1), weight); + EXPECT_EQ(sumw2[i + 1 + (Bins + 2) * (2 * i + 1 + (Bins + 2) * (3 * i + 1))], weight * weight); + } +} + +TEST(ConvertToTH3D, RHist) +{ + static constexpr std::size_t Bins = 20; + const RRegularAxis axis(Bins, {0, Bins}); + const RHist histD(axis, axis, axis); + + auto th3d = ConvertToTH3D(histD); + ASSERT_TRUE(th3d); + + RHist histE(axis, axis, axis); + for (std::size_t i = 0; i < Bins; i++) { + histE.Fill(i, 2 * i, 3 * i, RWeight(0.1 + i * 0.03)); + } + + th3d = ConvertToTH3D(histE); + ASSERT_TRUE(th3d); + + ASSERT_EQ(histE.GetNEntries(), Bins); + EXPECT_EQ(th3d->GetEntries(), Bins); + Double_t stats[11]; + th3d->GetStats(stats); + EXPECT_DOUBLE_EQ(stats[0], 7.7); + EXPECT_DOUBLE_EQ(stats[1], 3.563); + EXPECT_DOUBLE_EQ(stats[2], 93.1); + EXPECT_DOUBLE_EQ(stats[3], 1330.0); + EXPECT_DOUBLE_EQ(stats[4], 2 * 93.1); + EXPECT_DOUBLE_EQ(stats[5], 4 * 1330.0); + EXPECT_DOUBLE_EQ(stats[6], 0); + EXPECT_DOUBLE_EQ(stats[7], 3 * 93.1); + EXPECT_DOUBLE_EQ(stats[8], 9 * 1330.0); + EXPECT_DOUBLE_EQ(stats[9], 0); + EXPECT_DOUBLE_EQ(stats[10], 0); +} diff --git a/hist/histv7util/test/histutil_test.hxx b/hist/histv7util/test/histutil_test.hxx index 4015e60da3490..f0d0c11e6ee0f 100644 --- a/hist/histv7util/test/histutil_test.hxx +++ b/hist/histv7util/test/histutil_test.hxx @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include