diff --git a/tree/dataframe/inc/ROOT/RDF/ActionHelpers.hxx b/tree/dataframe/inc/ROOT/RDF/ActionHelpers.hxx
index c64870a0e0c48..cf717ed48a2bd 100644
--- a/tree/dataframe/inc/ROOT/RDF/ActionHelpers.hxx
+++ b/tree/dataframe/inc/ROOT/RDF/ActionHelpers.hxx
@@ -325,6 +325,44 @@ std::size_t GetSize(const T &val)
    }
 }
 
+// trait class to implement looping over data containers
+template <typename Helper>
+class R__CLING_PTRCHECK(off) ExecLoopTrait {
+private:
+   template <typename... Iterators>
+   void ExecLoop(unsigned int slot, std::size_t elements, Iterators... its)
+   {
+      for (std::size_t i = 0; i < elements; i++) {
+         Exec(slot, *its...);
+         (std::advance(its, 1), ...);
+      }
+   }
+
+public:
+   template <typename... ColumnTypes>
+   void Exec(unsigned int slot, const ColumnTypes &...columnValues)
+   {
+      if constexpr (std::disjunction_v<IsDataContainer<ColumnTypes>...>) {
+         constexpr std::array<bool, sizeof...(ColumnTypes)> isContainer{IsDataContainer<ColumnTypes>::value...};
+         constexpr std::size_t firstContainerIdx = FindIdxTrue(isContainer);
+         std::array<std::size_t, sizeof...(columnValues)> sizes = {{GetSize(columnValues)...}};
+         std::size_t elements = 0;
+         for (std::size_t i = 0; i < isContainer.size(); i++) {
+            if (isContainer[i]) {
+               if (i == firstContainerIdx) {
+                  elements = sizes[i];
+               } else if (elements != sizes[i]) {
+                  throw std::runtime_error("Cannot fill values in containers of different sizes.");
+               }
+            }
+         }
+         ExecLoop(slot, elements, MakeBegin(columnValues)...);
+      } else {
+         static_cast<Helper *>(this)->ExecSingle(slot, columnValues...);
+      }
+   }
+};
+
 // Helpers for dealing with histograms and similar:
 template <typename H, typename = decltype(std::declval<H>().Reset())>
 void ResetIfPossible(H *h)
@@ -485,8 +523,8 @@ public:
 
 #ifdef R__HAS_ROOT7
 template <typename BinContentType, bool WithWeight = false>
-class R__CLING_PTRCHECK(off) RHistFillHelper
-   : public ROOT::Detail::RDF::RActionImpl<RHistFillHelper<BinContentType, WithWeight>> {
+class R__CLING_PTRCHECK(off) RHistFillHelper : public RActionImpl<RHistFillHelper<BinContentType, WithWeight>>,
+                                               public ExecLoopTrait<RHistFillHelper<BinContentType, WithWeight>> {
 public:
    using Result_t = ROOT::Experimental::RHist<BinContentType>;
 
@@ -524,7 +562,7 @@ public:
    }
 
    template <typename... ColumnTypes>
-   void Exec(unsigned int slot, const ColumnTypes &...columnValues)
+   void ExecSingle(unsigned int slot, const ColumnTypes &...columnValues)
    {
       if constexpr (WithWeight) {
          auto t = std::forward_as_tuple(columnValues...);
@@ -546,7 +584,8 @@ public:
 
 template <typename BinContentType, bool WithWeight = false>
 class R__CLING_PTRCHECK(off) RHistEngineFillHelper
-   : public ROOT::Detail::RDF::RActionImpl<RHistEngineFillHelper<BinContentType, WithWeight>> {
+   : public RActionImpl<RHistEngineFillHelper<BinContentType, WithWeight>>,
+     public ExecLoopTrait<RHistEngineFillHelper<BinContentType, WithWeight>> {
 public:
    using Result_t = ROOT::Experimental::RHistEngine<BinContentType>;
 
@@ -576,7 +615,7 @@ public:
    }
 
    template <typename... ColumnTypes>
-   void Exec(unsigned int, const ColumnTypes &...columnValues)
+   void ExecSingle(unsigned int, const ColumnTypes &...columnValues)
    {
       if constexpr (WithWeight) {
          auto t = std::forward_as_tuple(columnValues...);
diff --git a/tree/dataframe/inc/ROOT/RDF/RInterface.hxx b/tree/dataframe/inc/ROOT/RDF/RInterface.hxx
index aa3bb93285d00..8dbf31a98e446 100644
--- a/tree/dataframe/inc/ROOT/RDF/RInterface.hxx
+++ b/tree/dataframe/inc/ROOT/RDF/RInterface.hxx
@@ -2032,6 +2032,9 @@ public:
    /// \param[in] vName The name of the column that will fill the histogram.
    /// \return the histogram wrapped in a RResultPtr.
    ///
+   /// The column can be of a container type (e.g. `std::vector` or `RVec`), in which case the histogram is filled with
+   /// each one of the elements of the container, and the Container can also be nested.
+   ///
    /// This action is *lazy*: upon invocation of this method the calculation is
    /// booked but not executed. Also see RResultPtr.
    ///
@@ -2057,6 +2060,11 @@ public:
    /// \param[in] columnList A list containing the names of the columns that will be passed when calling `Fill`
    /// \return the histogram wrapped in a RResultPtr.
    ///
+   /// Columns can be of a container type (e.g. `std::vector` or `RVec`), in which case the histogram is filled with
+   /// each one of the elements of the container. In case multiple columns of container type are provided, they must
+   /// have the same length for each event (but possibly different lengths between events). Containers can be nested,
+   /// in which case their sizes must match recursively. Scalars are broadcasted to match container columns.
+   ///
    /// This action is *lazy*: upon invocation of this method the calculation is
    /// booked but not executed. Also see RResultPtr.
    ///
@@ -2086,6 +2094,11 @@ public:
    /// \param[in] columnList A list containing the names of the columns that will be passed when calling `Fill`
    /// \return the histogram wrapped in a RResultPtr.
    ///
+   /// Columns can be of a container type (e.g. `std::vector` or `RVec`), in which case the histogram is filled with
+   /// each one of the elements of the container. In case multiple columns of container type are provided, they must
+   /// have the same length for each event (but possibly different lengths between events). Containers can be nested,
+   /// in which case their sizes must match recursively. Scalars are broadcasted to match container columns.
+   ///
    /// This action is *lazy*: upon invocation of this method the calculation is
    /// booked but not executed. Also see RResultPtr.
    ///
@@ -2122,6 +2135,11 @@ public:
    /// \param[in] wName The name of the column that will provide the weights.
    /// \return the histogram wrapped in a RResultPtr.
    ///
+   /// Columns can be of a container type (e.g. `std::vector` or `RVec`), in which case the histogram is filled with
+   /// each one of the elements of the container. In case multiple columns of container type are provided, they must
+   /// have the same length for each event (but possibly different lengths between events). Containers can be nested,
+   /// in which case their sizes must match recursively. Scalars are broadcasted to match container columns.
+   ///
    /// This action is *lazy*: upon invocation of this method the calculation is
    /// booked but not executed. Also see RResultPtr.
    ///
@@ -2149,6 +2167,11 @@ public:
    /// \param[in] wName The name of the column that will provide the weights.
    /// \return the histogram wrapped in a RResultPtr.
    ///
+   /// Columns can be of a container type (e.g. `std::vector` or `RVec`), in which case the histogram is filled with
+   /// each one of the elements of the container. In case multiple columns of container type are provided, they must
+   /// have the same length for each event (but possibly different lengths between events). Containers can be nested,
+   /// in which case their sizes must match recursively. Scalars are broadcasted to match container columns.
+   ///
    /// This action is *lazy*: upon invocation of this method the calculation is
    /// booked but not executed. Also see RResultPtr.
    ///
@@ -2185,6 +2208,11 @@ public:
    /// \param[in] wName The name of the column that will provide the weights.
    /// \return the histogram wrapped in a RResultPtr.
    ///
+   /// Columns can be of a container type (e.g. `std::vector` or `RVec`), in which case the histogram is filled with
+   /// each one of the elements of the container. In case multiple columns of container type are provided, they must
+   /// have the same length for each event (but possibly different lengths between events). Containers can be nested,
+   /// in which case their sizes must match recursively. Scalars are broadcasted to match container columns.
+   ///
    /// This action is *lazy*: upon invocation of this method the calculation is
    /// booked but not executed. Also see RResultPtr.
    ///
@@ -2228,6 +2256,11 @@ public:
    /// \param[in] columnList A list containing the names of the columns that will be passed when calling `Fill`
    /// \return the histogram wrapped in a RResultPtr.
    ///
+   /// Columns can be of a container type (e.g. `std::vector` or `RVec`), in which case the histogram is filled with
+   /// each one of the elements of the container. In case multiple columns of container type are provided, they must
+   /// have the same length for each event (but possibly different lengths between events). Containers can be nested,
+   /// in which case their sizes must match recursively. Scalars are broadcasted to match container columns.
+   ///
    /// This action is *lazy*: upon invocation of this method the calculation is
    /// booked but not executed. Also see RResultPtr.
    ///
@@ -2262,6 +2295,11 @@ public:
    /// \param[in] wName The name of the column that will provide the weights.
    /// \return the histogram wrapped in a RResultPtr.
    ///
+   /// Columns can be of a container type (e.g. `std::vector` or `RVec`), in which case the histogram is filled with
+   /// each one of the elements of the container. In case multiple columns of container type are provided, they must
+   /// have the same length for each event (but possibly different lengths between events). Containers can be nested,
+   /// in which case their sizes must match recursively. Scalars are broadcasted to match container columns.
+   ///
    /// This action is *lazy*: upon invocation of this method the calculation is
    /// booked but not executed. Also see RResultPtr.
    ///
diff --git a/tree/dataframe/test/dataframe_helpers.cxx b/tree/dataframe/test/dataframe_helpers.cxx
index 97696841a27fc..89d7ec3ac61d4 100755
--- a/tree/dataframe/test/dataframe_helpers.cxx
+++ b/tree/dataframe/test/dataframe_helpers.cxx
@@ -10,8 +10,10 @@
 #include <algorithm>
 #include <deque>
 #include <iomanip>
-#include <vector>
+#include <stdexcept>
 #include <string>
+#include <utility>
+#include <vector>
 
 #include "gtest/gtest.h"
 
@@ -813,6 +815,147 @@ TEST(RDFHelpers, Cleanup_After_Exception)
       << "The Finalize method should have changed the value of testVal during the post-exception cleanup." << std::endl;
 }
 
+struct SimpleExecLoopAction : public ROOT::Internal::RDF::ExecLoopTrait<SimpleExecLoopAction> {
+   std::vector<double> fValues;
+
+   template <typename... ColumnTypes>
+   void ExecSingle(unsigned int, const ColumnTypes &...columnValues)
+   {
+      fValues.emplace_back(columnValues...);
+   }
+};
+
+TEST(RDFHelpers, ExecLoopTrait)
+{
+   SimpleExecLoopAction execLoop;
+
+   execLoop.Exec(/*slot=*/0, 42.0);
+   ASSERT_EQ(execLoop.fValues.size(), 1);
+   EXPECT_EQ(execLoop.fValues.front(), 42.0);
+   execLoop.fValues.clear();
+
+   // One-dimensional containers are processed in order.
+   std::vector<double> v = {43.0};
+   execLoop.Exec(/*slot=*/0, v);
+   EXPECT_EQ(execLoop.fValues.size(), 1);
+   EXPECT_EQ(execLoop.fValues, v);
+   execLoop.fValues.clear();
+
+   v = {44.0, 45.0};
+   execLoop.Exec(/*slot=*/0, v);
+   EXPECT_EQ(execLoop.fValues.size(), 2);
+   EXPECT_EQ(execLoop.fValues, v);
+   execLoop.fValues.clear();
+
+   v = {};
+   execLoop.Exec(/*slot=*/0, v);
+   EXPECT_TRUE(execLoop.fValues.empty());
+   execLoop.fValues.clear();
+
+   // Multi-dimensional containers are traversed recursively.
+   std::vector<std::vector<double>> v2 = {{1.0}, {}, {2.0, 3.0}, {4.0}};
+   execLoop.Exec(/*slot=*/0, v2);
+   EXPECT_EQ(execLoop.fValues.size(), 4);
+   v = {1.0, 2.0, 3.0, 4.0};
+   EXPECT_EQ(execLoop.fValues, v);
+   execLoop.fValues.clear();
+
+   std::vector<std::vector<std::vector<double>>> v3 = {{{1.0}, {2.0}}, {}, {{3.0, 4.0}}, {{5.0}}};
+   execLoop.Exec(/*slot=*/0, v3);
+   EXPECT_EQ(execLoop.fValues.size(), 5);
+   v = {1.0, 2.0, 3.0, 4.0, 5.0};
+   EXPECT_EQ(execLoop.fValues, v);
+   execLoop.fValues.clear();
+}
+
+struct ExecLoopActionTwo : public ROOT::Internal::RDF::ExecLoopTrait<ExecLoopActionTwo> {
+   std::vector<std::pair<double, double>> fValues;
+
+   template <typename... ColumnTypes>
+   void ExecSingle(unsigned int, const ColumnTypes &...columnValues)
+   {
+      fValues.emplace_back(columnValues...);
+   }
+};
+
+TEST(RDFHelpers, ExecLoopTraitTwo)
+{
+   ExecLoopActionTwo execLoop;
+
+   execLoop.Exec(/*slot=*/0, 42.0, 43.0);
+   ASSERT_EQ(execLoop.fValues.size(), 1);
+   EXPECT_EQ(execLoop.fValues.front(), std::make_pair(42.0, 43.0));
+   execLoop.fValues.clear();
+
+   // Two containers are traversed in sync.
+   std::vector<double> v1 = {1.0, 2.0};
+   std::vector<double> v2 = {3.0, 4.0};
+   execLoop.Exec(/*slot=*/0, v1, v2);
+   ASSERT_EQ(execLoop.fValues.size(), 2);
+   EXPECT_EQ(execLoop.fValues.at(0), std::make_pair(1.0, 3.0));
+   EXPECT_EQ(execLoop.fValues.at(1), std::make_pair(2.0, 4.0));
+   execLoop.fValues.clear();
+
+   // A scalar is broadcasted to match the containers.
+   execLoop.Exec(/*slot=*/0, v1, 42.0);
+   ASSERT_EQ(execLoop.fValues.size(), 2);
+   EXPECT_EQ(execLoop.fValues.at(0), std::make_pair(1.0, 42.0));
+   EXPECT_EQ(execLoop.fValues.at(1), std::make_pair(2.0, 42.0));
+   execLoop.fValues.clear();
+
+   execLoop.Exec(/*slot=*/0, 42.0, v2);
+   ASSERT_EQ(execLoop.fValues.size(), 2);
+   EXPECT_EQ(execLoop.fValues.at(0), std::make_pair(42.0, 3.0));
+   EXPECT_EQ(execLoop.fValues.at(1), std::make_pair(42.0, 4.0));
+   execLoop.fValues.clear();
+
+   // Two nested containers are traversed in sync.
+   std::vector<std::vector<double>> v3 = {{1.0, 2.0}, {3.0}};
+   std::vector<std::vector<double>> v4 = {{4.0, 5.0}, {6.0}};
+   execLoop.Exec(/*slot=*/0, v3, v4);
+   ASSERT_EQ(execLoop.fValues.size(), 3);
+   EXPECT_EQ(execLoop.fValues.at(0), std::make_pair(1.0, 4.0));
+   EXPECT_EQ(execLoop.fValues.at(1), std::make_pair(2.0, 5.0));
+   EXPECT_EQ(execLoop.fValues.at(2), std::make_pair(3.0, 6.0));
+   execLoop.fValues.clear();
+
+   // Broadcasting works recursively in multiple dimensions.
+   execLoop.Exec(/*slot=*/0, v3, 42.0);
+   ASSERT_EQ(execLoop.fValues.size(), 3);
+   EXPECT_EQ(execLoop.fValues.at(0), std::make_pair(1.0, 42.0));
+   EXPECT_EQ(execLoop.fValues.at(1), std::make_pair(2.0, 42.0));
+   EXPECT_EQ(execLoop.fValues.at(2), std::make_pair(3.0, 42.0));
+   execLoop.fValues.clear();
+
+   execLoop.Exec(/*slot=*/0, v1, v3);
+   ASSERT_EQ(execLoop.fValues.size(), 3);
+   EXPECT_EQ(execLoop.fValues.at(0), std::make_pair(1.0, 1.0));
+   EXPECT_EQ(execLoop.fValues.at(1), std::make_pair(1.0, 2.0));
+   EXPECT_EQ(execLoop.fValues.at(2), std::make_pair(2.0, 3.0));
+   execLoop.fValues.clear();
+}
+
+TEST(RDFHelpers, ExecLoopTraitInvalid)
+{
+   ExecLoopActionTwo execLoop;
+
+   // The number of elements has to match for all columns.
+   std::vector<double> v0;
+   std::vector<double> v1 = {1.0};
+   std::vector<double> v2 = {2.0, 3.0};
+   EXPECT_THROW(execLoop.Exec(/*slot=*/0, v0, v1), std::runtime_error);
+   EXPECT_THROW(execLoop.Exec(/*slot=*/0, v0, v2), std::runtime_error);
+   EXPECT_THROW(execLoop.Exec(/*slot=*/0, v2, v1), std::runtime_error);
+
+   // For multiple dimensions, the number of elements has to match recursively.
+   std::vector<std::vector<double>> v4 = {{1.0}, {2.0}};
+   std::vector<std::vector<double>> v5 = {{3.0, 4.0}};
+   std::vector<std::vector<double>> v6 = {{3.0, 4.0}, {}};
+   EXPECT_THROW(execLoop.Exec(/*slot=*/0, v4, v5), std::runtime_error);
+   EXPECT_THROW(execLoop.Exec(/*slot=*/0, v4, v6), std::runtime_error);
+   EXPECT_THROW(execLoop.Exec(/*slot=*/0, v6, v5), std::runtime_error);
+}
+
 TEST(RDFHelpers, ProgressBarRestorePrecision)
 {
    // Regression test for https://github.com/root-project/root/issues/21165
diff --git a/tree/dataframe/test/dataframe_hist.cxx b/tree/dataframe/test/dataframe_hist.cxx
index f5fe396b50de3..c4a003fb40fc3 100644
--- a/tree/dataframe/test/dataframe_hist.cxx
+++ b/tree/dataframe/test/dataframe_hist.cxx
@@ -8,6 +8,7 @@
 #include <ROOT/RHistEngine.hxx>
 #include <ROOT/RRegularAxis.hxx>
 #include <ROOT/RVariableBinAxis.hxx>
+#include <ROOT/RVec.hxx>
 
 #include <cstddef>
 #include <memory>
@@ -269,7 +270,7 @@ TEST_P(RDFHist, InvalidNumberOfArguments)
       // expected
    }
 
-   auto hist = std::make_shared<RHist<double>>(10, std::make_pair(5.0, 15.0));
+   auto hist = std::make_shared<RHist<double>>(axis);
    try {
       // Cannot use EXPECT_THROW because of template arguments...
       dfX.Hist<double, double>(hist, {"x", "x"});
@@ -278,7 +279,7 @@ TEST_P(RDFHist, InvalidNumberOfArguments)
       // expected
    }
 
-   auto engine = std::make_shared<RHistEngine<double>>(10, std::make_pair(5.0, 15.0));
+   auto engine = std::make_shared<RHistEngine<double>>(axis);
    try {
       // Cannot use EXPECT_THROW because of template arguments...
       dfX.Hist<double, double>(engine, {"x", "x"});
@@ -296,10 +297,10 @@ TEST_P(RDFHist, InvalidNumberOfArgumentsJit)
    const RRegularAxis axis(10, {5.0, 15.0});
    EXPECT_THROW(dfX.Hist({axis}, {"x", "x"}), std::invalid_argument);
 
-   auto hist = std::make_shared<RHist<double>>(10, std::make_pair(5.0, 15.0));
+   auto hist = std::make_shared<RHist<double>>(axis);
    EXPECT_THROW(dfX.Hist(hist, {"x", "x"}), std::invalid_argument);
 
-   auto engine = std::make_shared<RHistEngine<double>>(10, std::make_pair(5.0, 15.0));
+   auto engine = std::make_shared<RHistEngine<double>>(axis);
    EXPECT_THROW(dfX.Hist(engine, {"x", "x"}), std::invalid_argument);
 }
 
@@ -378,7 +379,7 @@ TEST_P(RDFHist, EngineWeight)
                   .Define("w", [](ULong64_t e) { return 0.1 + e * 0.03; }, {"rdfentry_"});
 
    const RRegularAxis axis(10, {5.0, 15.0});
-   auto hist = std::make_shared<RHistEngine<RBinWithError>>(10, std::make_pair(5.0, 15.0));
+   auto hist = std::make_shared<RHistEngine<RBinWithError>>(axis);
    auto resPtr = dfXW.Hist<double, double>(hist, {"x"}, "w");
    EXPECT_EQ(hist, resPtr.GetSharedPtr());
    for (auto index : axis.GetNormalRange()) {
@@ -395,7 +396,7 @@ TEST_P(RDFHist, EngineWeightJit)
    auto dfXW = df.Define("x", "rdfentry_ + 5.5").Define("w", "0.1 + rdfentry_ * 0.03");
 
    const RRegularAxis axis(10, {5.0, 15.0});
-   auto hist = std::make_shared<RHistEngine<RBinWithError>>(10, std::make_pair(5.0, 15.0));
+   auto hist = std::make_shared<RHistEngine<RBinWithError>>(axis);
    auto resPtr = dfXW.Hist(hist, {"x"}, "w");
    EXPECT_EQ(hist, resPtr.GetSharedPtr());
    for (auto index : axis.GetNormalRange()) {
@@ -421,7 +422,7 @@ TEST_P(RDFHist, WeightInvalidNumberOfArguments)
       // expected
    }
 
-   auto hist = std::make_shared<RHist<double>>(10, std::make_pair(5.0, 15.0));
+   auto hist = std::make_shared<RHist<double>>(axis);
    try {
       // Cannot use EXPECT_THROW because of template arguments...
       dfXW.Hist<double, double, double>(hist, {"x", "x"}, "w");
@@ -430,7 +431,7 @@ TEST_P(RDFHist, WeightInvalidNumberOfArguments)
       // expected
    }
 
-   auto engine = std::make_shared<RHistEngine<double>>(10, std::make_pair(5.0, 15.0));
+   auto engine = std::make_shared<RHistEngine<double>>(axis);
    try {
       // Cannot use EXPECT_THROW because of template arguments...
       dfXW.Hist<double, double, double>(engine, {"x", "x"}, "w");
@@ -448,13 +449,272 @@ TEST_P(RDFHist, WeightInvalidNumberOfArgumentsJit)
    const RRegularAxis axis(10, {5.0, 15.0});
    EXPECT_THROW(dfXW.Hist({axis}, {"x", "x"}, "w"), std::invalid_argument);
 
-   auto hist = std::make_shared<RHist<double>>(10, std::make_pair(5.0, 15.0));
+   auto hist = std::make_shared<RHist<double>>(axis);
    EXPECT_THROW(dfXW.Hist(hist, {"x", "x"}, "w"), std::invalid_argument);
 
-   auto engine = std::make_shared<RHistEngine<double>>(10, std::make_pair(5.0, 15.0));
+   auto engine = std::make_shared<RHistEngine<double>>(axis);
    EXPECT_THROW(dfXW.Hist(engine, {"x", "x"}, "w"), std::invalid_argument);
 }
 
+TEST_P(RDFHist, Container)
+{
+   RDataFrame df(10);
+   auto dfX = df.Define("x", [](ULong64_t e) { return ROOT::RVecD{e + 5.5, e + 15.5}; }, {"rdfentry_"});
+
+   const RRegularAxis axis(20, {5.0, 25.0});
+   auto hist = dfX.Hist</*BinContentType=*/double, ROOT::RVecD>({axis}, {"x"});
+   EXPECT_EQ(hist->GetNEntries(), 20);
+   for (auto index : axis.GetNormalRange()) {
+      EXPECT_EQ(hist->GetBinContent(index), 1.0);
+   }
+}
+
+TEST_P(RDFHist, ContainerNested)
+{
+   RDataFrame df(10);
+   auto dfX = df.Define("x", [](ULong64_t e) { return ROOT::RVec<ROOT::RVecD>{{e + 5.5}, {e + 15.5}}; }, {"rdfentry_"});
+
+   const RRegularAxis axis(20, {5.0, 25.0});
+   auto hist = dfX.Hist</*BinContentType=*/double, ROOT::RVec<ROOT::RVecD>>({axis}, {"x"});
+   EXPECT_EQ(hist->GetNEntries(), 20);
+   for (auto index : axis.GetNormalRange()) {
+      EXPECT_EQ(hist->GetBinContent(index), 1.0);
+   }
+}
+
+TEST_P(RDFHist, ContainerJit)
+{
+   RDataFrame df(10);
+   auto dfX = df.Define("x", "ROOT::RVecD{rdfentry_ + 5.5}");
+
+   const RRegularAxis axis(10, {5.0, 15.0});
+   auto hist = dfX.Hist({axis}, {"x"});
+   EXPECT_EQ(hist->GetNEntries(), 10);
+   for (auto index : axis.GetNormalRange()) {
+      EXPECT_EQ(hist->GetBinContent(index), 1.0);
+   }
+}
+
+TEST_P(RDFHist, ContainerMultiDim)
+{
+   RDataFrame df(10);
+   auto dfXY = df.Define("x", [](ULong64_t e) { return ROOT::RVecD{e + 5.5}; }, {"rdfentry_"})
+                  .Define("y", [](ULong64_t e) { return ROOT::RVecD{2 * e + 0.5}; }, {"rdfentry_"});
+
+   const RRegularAxis regularAxis(10, {5.0, 15.0});
+   static constexpr std::size_t BinsY = 20;
+   std::vector<double> bins;
+   for (std::size_t i = 0; i < BinsY; i++) {
+      bins.push_back(i);
+   }
+   bins.push_back(BinsY);
+   const RVariableBinAxis variableBinAxis(bins);
+
+   auto hist =
+      dfXY.Hist</*BinContentType=*/double, ROOT::RVecD, ROOT::RVecD>({regularAxis, variableBinAxis}, {"x", "y"});
+   EXPECT_EQ(hist->GetNEntries(), 10);
+   for (auto x : regularAxis.GetNormalRange()) {
+      for (auto y : variableBinAxis.GetNormalRange()) {
+         if (2 * x.GetIndex() == y.GetIndex()) {
+            EXPECT_EQ(hist->GetBinContent(x, y), 1.0);
+         } else {
+            EXPECT_EQ(hist->GetBinContent(x, y), 0.0);
+         }
+      }
+   }
+}
+
+TEST_P(RDFHist, ContainerMultiDimBroadcast)
+{
+   RDataFrame df(10);
+   auto dfXY = df.Define("x", [](ULong64_t e) { return e + 5.5; }, {"rdfentry_"})
+                  .Define("y", [](ULong64_t e) { return ROOT::RVecD{2 * e + 0.5}; }, {"rdfentry_"});
+
+   const RRegularAxis regularAxis(10, {5.0, 15.0});
+   static constexpr std::size_t BinsY = 20;
+   std::vector<double> bins;
+   for (std::size_t i = 0; i < BinsY; i++) {
+      bins.push_back(i);
+   }
+   bins.push_back(BinsY);
+   const RVariableBinAxis variableBinAxis(bins);
+
+   auto hist = dfXY.Hist</*BinContentType=*/double, double, ROOT::RVecD>({regularAxis, variableBinAxis}, {"x", "y"});
+   EXPECT_EQ(hist->GetNEntries(), 10);
+   for (auto x : regularAxis.GetNormalRange()) {
+      for (auto y : variableBinAxis.GetNormalRange()) {
+         if (2 * x.GetIndex() == y.GetIndex()) {
+            EXPECT_EQ(hist->GetBinContent(x, y), 1.0);
+         } else {
+            EXPECT_EQ(hist->GetBinContent(x, y), 0.0);
+         }
+      }
+   }
+}
+
+TEST_P(RDFHist, ContainerMultiDimNestedBroadcast)
+{
+   RDataFrame df(10);
+   auto dfXY = df.Define("x", [](ULong64_t e) { return e + 5.5; }, {"rdfentry_"})
+                  .Define("xV", [](ULong64_t e) { return ROOT::RVecD{e + 5.5}; }, {"rdfentry_"})
+                  .Define("y", [](ULong64_t e) { return ROOT::RVec<ROOT::RVecD>{{2 * e + 0.5}}; }, {"rdfentry_"});
+
+   const RRegularAxis regularAxis(10, {5.0, 15.0});
+   static constexpr std::size_t BinsY = 20;
+   std::vector<double> bins;
+   for (std::size_t i = 0; i < BinsY; i++) {
+      bins.push_back(i);
+   }
+   bins.push_back(BinsY);
+   const RVariableBinAxis variableBinAxis(bins);
+
+   auto hist =
+      dfXY.Hist</*BinContentType=*/double, double, ROOT::RVec<ROOT::RVecD>>({regularAxis, variableBinAxis}, {"x", "y"});
+   EXPECT_EQ(hist->GetNEntries(), 10);
+   for (auto x : regularAxis.GetNormalRange()) {
+      for (auto y : variableBinAxis.GetNormalRange()) {
+         if (2 * x.GetIndex() == y.GetIndex()) {
+            EXPECT_EQ(hist->GetBinContent(x, y), 1.0);
+         } else {
+            EXPECT_EQ(hist->GetBinContent(x, y), 0.0);
+         }
+      }
+   }
+
+   hist = dfXY.Hist</*BinContentType=*/double, ROOT::RVecD, ROOT::RVec<ROOT::RVecD>>({regularAxis, variableBinAxis},
+                                                                                     {"xV", "y"});
+   EXPECT_EQ(hist->GetNEntries(), 10);
+}
+
+TEST_P(RDFHist, ContainerInvalidElements)
+{
+   RDataFrame df(10);
+   auto dfXY = df.Define("x0", [] { return ROOT::RVecD{}; })
+                  .Define("x2", [](ULong64_t e) { return ROOT::RVecD{e + 5.5, e + 15.5}; }, {"rdfentry_"})
+                  .Define("y1", [](ULong64_t e) { return ROOT::RVecD{2 * e + 0.5}; }, {"rdfentry_"});
+
+   const RRegularAxis regularAxis(20, {5.0, 25.0});
+   static constexpr std::size_t BinsY = 20;
+   std::vector<double> bins;
+   for (std::size_t i = 0; i < BinsY; i++) {
+      bins.push_back(i);
+   }
+   bins.push_back(BinsY);
+   const RVariableBinAxis variableBinAxis(bins);
+
+   try {
+      // Cannot use EXPECT_THROW because of template arguments...
+      *dfXY.Hist</*BinContentType=*/double, ROOT::RVecD, ROOT::RVecD>({regularAxis, variableBinAxis}, {"x0", "y1"});
+      FAIL() << "expected std::runtime_error";
+   } catch (const std::runtime_error &) {
+      // expected
+   }
+
+   try {
+      // Cannot use EXPECT_THROW because of template arguments...
+      *dfXY.Hist</*BinContentType=*/double, ROOT::RVecD, ROOT::RVecD>({regularAxis, variableBinAxis}, {"x2", "y1"});
+      FAIL() << "expected std::runtime_error";
+   } catch (const std::runtime_error &e) {
+      // expected
+   }
+}
+
+TEST_P(RDFHist, EngineContainer)
+{
+   RDataFrame df(10);
+   auto dfX = df.Define("x", [](ULong64_t e) { return ROOT::RVecD{e + 5.5, e + 15.5}; }, {"rdfentry_"});
+
+   const RRegularAxis axis(20, {5.0, 25.0});
+   auto hist = std::make_shared<RHistEngine<double>>(axis);
+   auto resPtr = dfX.Hist<ROOT::RVecD>(hist, {"x"});
+   EXPECT_EQ(hist, resPtr.GetSharedPtr());
+   for (auto index : axis.GetNormalRange()) {
+      EXPECT_EQ(hist->GetBinContent(index), 1.0);
+   }
+}
+
+TEST_P(RDFHist, WeightContainer)
+{
+   RDataFrame df(10);
+   auto dfXW =
+      df.Define("xV", [](ULong64_t e) { return ROOT::RVecD{e + 5.5, e + 15.5}; }, {"rdfentry_"})
+         .Define("wV", [](ULong64_t e) { return ROOT::RVecD{0.1 + e * 0.03, 0.2 + e * e * 0.06}; }, {"rdfentry_"});
+
+   const RRegularAxis axis(20, {5.0, 25.0});
+   auto hist = dfXW.Hist</*BinContentType=*/RBinWithError, ROOT::RVecD, ROOT::RVecD>({axis}, {"xV"}, "wV");
+   EXPECT_EQ(hist->GetNEntries(), 20);
+   for (auto index : axis.GetNormalRange()) {
+      auto &bin = hist->GetBinContent(index);
+      auto i = index.GetIndex();
+      double weight = i >= 10 ? (0.2 + (i - 10) * (i - 10) * 0.06) : (0.1 + i * 0.03);
+      EXPECT_FLOAT_EQ(bin.fSum, weight);
+      EXPECT_FLOAT_EQ(bin.fSum2, weight * weight);
+   }
+}
+
+TEST_P(RDFHist, WeightContainerJit)
+{
+   RDataFrame df(10);
+   auto dfXW = df.Define("xV", "ROOT::RVecD{rdfentry_ + 5.5, rdfentry_ + 15.5}")
+                  .Define("wV", "ROOT::RVecD{0.1 + rdfentry_ * 0.03, 0.2 + rdfentry_ * rdfentry_ * 0.06}");
+
+   const RRegularAxis axis(20, {5.0, 25.0});
+   auto hist = dfXW.Hist({axis}, {"xV"}, "wV");
+   EXPECT_EQ(hist->GetNEntries(), 20);
+   for (auto index : axis.GetNormalRange()) {
+      auto &bin = hist->GetBinContent(index);
+      auto i = index.GetIndex();
+      double weight = i >= 10 ? (0.2 + (i - 10) * (i - 10) * 0.06) : (0.1 + i * 0.03);
+      EXPECT_FLOAT_EQ(bin.fSum, weight);
+      EXPECT_FLOAT_EQ(bin.fSum2, weight * weight);
+   }
+}
+
+TEST_P(RDFHist, WeightContainerBroadcast)
+{
+   RDataFrame df(10);
+   auto dfXW =
+      df.Define("x", [](ULong64_t e) { return e + 5.5; }, {"rdfentry_"})
+         .Define("xV", [](ULong64_t e) { return ROOT::RVecD{e + 5.5, e + 15.5}; }, {"rdfentry_"})
+         .Define("w", [](ULong64_t e) { return 0.1 + e * 0.03; }, {"rdfentry_"})
+         .Define("wV", [](ULong64_t e) { return ROOT::RVecD{0.1 + e * 0.03, 0.2 + e * e * 0.06}; }, {"rdfentry_"});
+
+   const RRegularAxis axis(20, {5.0, 25.0});
+   auto hist = dfXW.Hist</*BinContentType=*/RBinWithError, ROOT::RVecD, double>({axis}, {"xV"}, "w");
+   EXPECT_EQ(hist->GetNEntries(), 20);
+   for (auto index : axis.GetNormalRange()) {
+      auto &bin = hist->GetBinContent(index);
+      auto i = index.GetIndex();
+      double weight = 0.1 + (i >= 10 ? i - 10 : i) * 0.03;
+      EXPECT_FLOAT_EQ(bin.fSum, weight);
+      EXPECT_FLOAT_EQ(bin.fSum2, weight * weight);
+   }
+
+   hist = dfXW.Hist</*BinContentType=*/RBinWithError, double, ROOT::RVecD>({axis}, {"x"}, "wV");
+   EXPECT_EQ(hist->GetNEntries(), 20);
+   EXPECT_EQ(hist->GetBinContent(2).fSum, 0.1 + 2 * 0.03 + 0.2 + 2 * 2 * 0.06);
+}
+
+TEST_P(RDFHist, EngineWeightContainer)
+{
+   RDataFrame df(10);
+   auto dfXW =
+      df.Define("xV", [](ULong64_t e) { return ROOT::RVecD{e + 5.5, e + 15.5}; }, {"rdfentry_"})
+         .Define("wV", [](ULong64_t e) { return ROOT::RVecD{0.1 + e * 0.03, 0.2 + e * e * 0.06}; }, {"rdfentry_"});
+
+   const RRegularAxis axis(20, {5.0, 25.0});
+   auto hist = std::make_shared<RHistEngine<RBinWithError>>(axis);
+   auto resPtr = dfXW.Hist<ROOT::RVecD, ROOT::RVecD>(hist, {"xV"}, "wV");
+   EXPECT_EQ(hist, resPtr.GetSharedPtr());
+   for (auto index : axis.GetNormalRange()) {
+      auto &bin = hist->GetBinContent(index);
+      auto i = index.GetIndex();
+      double weight = i >= 10 ? (0.2 + (i - 10) * (i - 10) * 0.06) : (0.1 + i * 0.03);
+      EXPECT_FLOAT_EQ(bin.fSum, weight);
+      EXPECT_FLOAT_EQ(bin.fSum2, weight * weight);
+   }
+}
+
 INSTANTIATE_TEST_SUITE_P(Seq, RDFHist, ::testing::Values(false));
 
 #ifdef R__USE_IMT