[hist] Conversion of RBinWithError

Author: hahnjo

hist/histv7util/inc/ROOT/Hist/ConvertToTH1.hxx

@@ -5,6 +5,7 @@
 #ifndef ROOT_Hist_ConvertToTH1
 #define ROOT_Hist_ConvertToTH1
 
+#include <ROOT/RBinWithError.hxx>
 #include <ROOT/RHist.hxx>
 #include <ROOT/RHistEngine.hxx>
 
@@ -61,6 +62,11 @@ std::unique_ptr<TH1F> ConvertToTH1F(const RHistEngine<float> &engine);
 /// \copydetails ConvertToTH1I(const RHistEngine<int> &engine)
 std::unique_ptr<TH1D> ConvertToTH1D(const RHistEngine<double> &engine);
 
+/// Convert a one-dimensional histogram to TH1D.
+///
+/// \copydetails ConvertToTH1I(const RHistEngine<int> &engine)
+std::unique_ptr<TH1D> ConvertToTH1D(const RHistEngine<RBinWithError> &engine);
+
 /// Convert a one-dimensional histogram to TH1C.
 ///
 /// \copydetails ConvertToTH1I(const RHist<int> &hist)
@@ -99,6 +105,11 @@ std::unique_ptr<TH1F> ConvertToTH1F(const RHist<float> &hist);
 /// \copydetails ConvertToTH1I(const RHist<int> &hist)
 std::unique_ptr<TH1D> ConvertToTH1D(const RHist<double> &hist);
 
+/// Convert a one-dimensional histogram to TH1D.
+///
+/// \copydetails ConvertToTH1I(const RHist<int> &hist)
+std::unique_ptr<TH1D> ConvertToTH1D(const RHist<RBinWithError> &hist);
+
 } // namespace Hist
 } // namespace Experimental
 } // namespace ROOT

hist/histv7util/src/ConvertToTH1.cxx

@@ -13,6 +13,7 @@
 #include <cstddef>
 #include <memory>
 #include <stdexcept>
+#include <type_traits>
 #include <variant>
 #include <vector>
 
@@ -27,8 +28,20 @@ std::unique_ptr<Hist> ConvertToTH1Impl(const RHistEngine<T> &engine)
    }
 
    std::unique_ptr<Hist> ret;
-   auto copyBinContent = [&ret, &engine](Int_t i, RBinIndex index) {
-      ret->GetArray()[i] = engine.GetBinContent(index);
+   Double_t *sumw2 = nullptr;
+   auto copyBinContent = [&ret, &engine, &sumw2](Int_t i, RBinIndex index) {
+      if constexpr (std::is_same_v<T, RBinWithError>) {
+         if (sumw2 == nullptr) {
+            ret->Sumw2();
+            sumw2 = ret->GetSumw2()->GetArray();
+         }
+         const RBinWithError &c = engine.GetBinContent(index);
+         ret->GetArray()[i] = c.fSum;
+         sumw2[i] = c.fSum2;
+      } else {
+         (void)sumw2;
+         ret->GetArray()[i] = engine.GetBinContent(index);
+      }
    };
 
    const auto &axes = engine.GetAxes();
@@ -114,6 +127,11 @@ std::unique_ptr<TH1D> ConvertToTH1D(const RHistEngine<double> &engine)
    return ConvertToTH1Impl<TH1D>(engine);
 }
 
+std::unique_ptr<TH1D> ConvertToTH1D(const RHistEngine<RBinWithError> &engine)
+{
+   return ConvertToTH1Impl<TH1D>(engine);
+}
+
 std::unique_ptr<TH1C> ConvertToTH1C(const RHist<char> &hist)
 {
    auto ret = ConvertToTH1C(hist.GetEngine());
@@ -163,6 +181,13 @@ std::unique_ptr<TH1D> ConvertToTH1D(const RHist<double> &hist)
    return ret;
 }
 
+std::unique_ptr<TH1D> ConvertToTH1D(const RHist<RBinWithError> &hist)
+{
+   auto ret = ConvertToTH1D(hist.GetEngine());
+   ConvertGlobalStatistics(*ret, hist.GetStats());
+   return ret;
+}
+
 } // namespace Hist
 } // namespace Experimental
 } // namespace ROOT

hist/histv7util/test/hist_convert_TH1.cxx

@@ -208,17 +208,50 @@ TEST(ConvertToTH1F, RHist)
 TEST(ConvertToTH1D, RHistEngine)
 {
    static constexpr std::size_t Bins = 20;
-   const RHistEngine<double> engine(Bins, {0, Bins});
+   const RHistEngine<double> engineD(Bins, {0, Bins});
 
-   auto th1d = ConvertToTH1D(engine);
+   auto th1d = ConvertToTH1D(engineD);
    ASSERT_TRUE(th1d);
+
+   RHistEngine<RBinWithError> engineE(Bins, {0, Bins});
+   for (std::size_t i = 0; i < Bins; i++) {
+      engineE.Fill(i, RWeight(0.1 + i * 0.03));
+   }
+
+   th1d = ConvertToTH1D(engineE);
+   ASSERT_TRUE(th1d);
+   const Double_t *sumw2 = th1d->GetSumw2()->GetArray();
+   ASSERT_TRUE(sumw2 != nullptr);
+
+   for (std::size_t i = 1; i <= Bins; i++) {
+      const double weight = 0.1 + (i - 1) * 0.03;
+      EXPECT_EQ(th1d->GetBinContent(i), weight);
+      EXPECT_EQ(sumw2[i], weight * weight);
+   }
 }
 
 TEST(ConvertToTH1D, RHist)
 {
    static constexpr std::size_t Bins = 20;
-   const RHist<double> hist(Bins, {0, Bins});
+   const RHist<double> histD(Bins, {0, Bins});
 
-   auto th1d = ConvertToTH1D(hist);
+   auto th1d = ConvertToTH1D(histD);
    ASSERT_TRUE(th1d);
+
+   RHist<RBinWithError> histE(Bins, {0, Bins});
+   for (std::size_t i = 0; i < Bins; i++) {
+      histE.Fill(i, RWeight(0.1 + i * 0.03));
+   }
+
+   th1d = ConvertToTH1D(histE);
+   ASSERT_TRUE(th1d);
+
+   ASSERT_EQ(histE.GetNEntries(), Bins);
+   EXPECT_EQ(th1d->GetEntries(), Bins);
+   Double_t stats[4];
+   th1d->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);
 }

hist/histv7util/test/histutil_test.hxx

@@ -2,11 +2,13 @@
 #define histutil_test
 
 #include <ROOT/Hist/ConvertToTH1.hxx>
+#include <ROOT/RBinWithError.hxx>
 #include <ROOT/RHist.hxx>
 #include <ROOT/RHistEngine.hxx>
 #include <ROOT/RWeight.hxx>
 
 using namespace ROOT::Experimental::Hist;
+using ROOT::Experimental::RBinWithError;
 using ROOT::Experimental::RHist;
 using ROOT::Experimental::RHistEngine;
 using ROOT::Experimental::RRegularAxis;