diff --git a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp
index fd8937bd300..2144239a62c 100644
--- a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp
+++ b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp
@@ -8,18 +8,15 @@
 
 #pragma once
 
-#include "Acts/Definitions/Algebra.hpp"
-#include "Acts/Definitions/TrackParametrization.hpp"
-#include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Geometry/GeometryHierarchyMap.hpp"
 #include "Acts/Utilities/Logger.hpp"
 #include "ActsExamples/Digitization/DigitizationConfig.hpp"
 #include "ActsExamples/Digitization/MeasurementCreation.hpp"
 #include "ActsExamples/Digitization/SmearingConfig.hpp"
 #include "ActsExamples/EventData/Cluster.hpp"
-#include "ActsExamples/EventData/Index.hpp"
 #include "ActsExamples/EventData/Measurement.hpp"
 #include "ActsExamples/EventData/SimHit.hpp"
+#include "ActsExamples/EventData/SimParticle.hpp"
 #include "ActsExamples/Framework/DataHandle.hpp"
 #include "ActsExamples/Framework/IAlgorithm.hpp"
 #include "ActsExamples/Framework/ProcessCode.hpp"
@@ -29,17 +26,10 @@
 #include "ActsFatras/Digitization/UncorrelatedHitSmearer.hpp"
 
 #include <cstddef>
-#include <memory>
 #include <string>
-#include <tuple>
-#include <utility>
 #include <variant>
 #include <vector>
 
-namespace ActsFatras {
-class Barcode;
-}  // namespace ActsFatras
-
 namespace ActsExamples {
 struct AlgorithmContext;
 
@@ -98,17 +88,17 @@ class DigitizationAlgorithm final : public IAlgorithm {
   using CellsMap =
       std::map<Acts::GeometryIdentifier, std::vector<Cluster::Cell>>;
 
-  ReadDataHandle<SimHitContainer> m_simContainerReadHandle{this,
-                                                           "SimHitContainer"};
-
-  WriteDataHandle<MeasurementContainer> m_measurementWriteHandle{
-      this, "Measurements"};
-  WriteDataHandle<CellsMap> m_cellsWriteHandle{this, "Cells"};
-  WriteDataHandle<ClusterContainer> m_clusterWriteHandle{this, "Clusters"};
-  WriteDataHandle<IndexMultimap<ActsFatras::Barcode>>
-      m_measurementParticlesMapWriteHandle{this, "MeasurementParticlesMap"};
-  WriteDataHandle<IndexMultimap<Index>> m_measurementSimHitsMapWriteHandle{
-      this, "MeasurementSimHitsMap"};
+  ReadDataHandle<SimHitContainer> m_inputHits{this, "InputHits"};
+
+  WriteDataHandle<MeasurementContainer> m_outputMeasurements{
+      this, "OutputMeasurements"};
+  WriteDataHandle<CellsMap> m_outputCells{this, "OutputCells"};
+  WriteDataHandle<ClusterContainer> m_outputClusters{this, "OutputClusters"};
+
+  WriteDataHandle<IndexMultimap<SimBarcode>> m_outputMeasurementParticlesMap{
+      this, "OutputMeasurementParticlesMap"};
+  WriteDataHandle<IndexMultimap<Index>> m_outputMeasurementSimHitsMap{
+      this, "OutputMeasurementSimHitsMap"};
 
   /// Construct a fixed-size smearer from a configuration.
   ///
diff --git a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationConfig.hpp b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationConfig.hpp
index 6a933f50075..4772aa40405 100644
--- a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationConfig.hpp
+++ b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationConfig.hpp
@@ -10,13 +10,9 @@
 
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
-#include "Acts/Definitions/Units.hpp"
 #include "Acts/Geometry/GeometryHierarchyMap.hpp"
 #include "Acts/Utilities/BinUtility.hpp"
-#include "Acts/Utilities/BinningType.hpp"
-#include "Acts/Utilities/Logger.hpp"
 #include "Acts/Utilities/Result.hpp"
-#include "ActsExamples/Digitization/DigitizationConfig.hpp"
 #include "ActsExamples/Digitization/Smearers.hpp"
 #include "ActsExamples/Digitization/SmearingConfig.hpp"
 #include "ActsExamples/Framework/RandomNumbers.hpp"
@@ -24,7 +20,6 @@
 
 #include <algorithm>
 #include <cstddef>
-#include <functional>
 #include <map>
 #include <memory>
 #include <stdexcept>
@@ -145,6 +140,7 @@ class DigitizationConfig {
   std::string outputMeasurementParticlesMap = "measurement_particles_map";
   /// Output collection to map measured hits to simulated hits.
   std::string outputMeasurementSimHitsMap = "measurement_simhits_map";
+
   /// Map of surface by identifier to allow local - to global
   std::unordered_map<Acts::GeometryIdentifier, const Acts::Surface *>
       surfaceByIdentifier;
@@ -157,11 +153,11 @@ class DigitizationConfig {
   /// clusters, measurements, and sim-hit-maps are output.
   bool doClusterization = true;
   /// Do we merge hits or not
-  bool doMerge;
+  bool doMerge = false;
   /// How close do parameters have to be to consider merged
-  const double mergeNsigma;
+  double mergeNsigma = 3.0;
   /// Consider clusters that share a corner as merged (8-cell connectivity)
-  const bool mergeCommonCorner;
+  bool mergeCommonCorner = false;
   /// Energy deposit threshold for accepting a hit
   /// For a generic readout frontend we assume 1000 e/h pairs, in Si each
   /// e/h-pair requiers on average an energy of 3.65 eV (PDG  review 2023,
diff --git a/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp b/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp
index 72f3905ab57..b0c44a88e4a 100644
--- a/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp
+++ b/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp
@@ -11,27 +11,15 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/Geometry/GeometryIdentifier.hpp"
-#include "Acts/Geometry/TrackingGeometry.hpp"
-#include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/BinUtility.hpp"
-#include "Acts/Utilities/Result.hpp"
 #include "ActsExamples/Digitization/ModuleClusters.hpp"
-#include "ActsExamples/EventData/GeometryContainers.hpp"
-#include "ActsExamples/EventData/Index.hpp"
-#include "ActsExamples/EventData/SimHit.hpp"
 #include "ActsExamples/Framework/AlgorithmContext.hpp"
-#include "ActsExamples/Utilities/GroupBy.hpp"
-#include "ActsExamples/Utilities/Range.hpp"
 #include "ActsFatras/EventData/Barcode.hpp"
-#include "ActsFatras/EventData/Hit.hpp"
 
 #include <algorithm>
 #include <array>
-#include <cmath>
-#include <cstdint>
 #include <limits>
 #include <ostream>
-#include <set>
 #include <stdexcept>
 #include <string>
 #include <utility>
@@ -69,12 +57,11 @@ ActsExamples::DigitizationAlgorithm::DigitizationAlgorithm(
           "Missing hit-to-simulated-hits map output collection");
     }
 
-    m_measurementWriteHandle.initialize(m_cfg.outputMeasurements);
-    m_clusterWriteHandle.initialize(m_cfg.outputClusters);
-    m_measurementParticlesMapWriteHandle.initialize(
+    m_outputMeasurements.initialize(m_cfg.outputMeasurements);
+    m_outputClusters.initialize(m_cfg.outputClusters);
+    m_outputMeasurementParticlesMap.initialize(
         m_cfg.outputMeasurementParticlesMap);
-    m_measurementSimHitsMapWriteHandle.initialize(
-        m_cfg.outputMeasurementSimHitsMap);
+    m_outputMeasurementSimHitsMap.initialize(m_cfg.outputMeasurementSimHitsMap);
   }
 
   if (m_cfg.doOutputCells) {
@@ -82,10 +69,10 @@ ActsExamples::DigitizationAlgorithm::DigitizationAlgorithm(
       throw std::invalid_argument("Missing cell output collection");
     }
 
-    m_cellsWriteHandle.initialize(m_cfg.outputCells);
+    m_outputCells.initialize(m_cfg.outputCells);
   }
 
-  m_simContainerReadHandle.initialize(m_cfg.inputSimHits);
+  m_inputHits.initialize(m_cfg.inputSimHits);
 
   // Create the digitizers from the configuration
   std::vector<std::pair<Acts::GeometryIdentifier, Digitizer>> digitizerInput;
@@ -142,13 +129,14 @@ ActsExamples::DigitizationAlgorithm::DigitizationAlgorithm(
 ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute(
     const AlgorithmContext& ctx) const {
   // Retrieve input
-  const auto& simHits = m_simContainerReadHandle(ctx);
+  const auto& simHits = m_inputHits(ctx);
   ACTS_DEBUG("Loaded " << simHits.size() << " sim hits");
 
   // Prepare output containers
   // need list here for stable addresses
   MeasurementContainer measurements;
   ClusterContainer clusters;
+
   IndexMultimap<ActsFatras::Barcode> measurementParticlesMap;
   IndexMultimap<Index> measurementSimHitsMap;
   measurements.reserve(simHits.size());
@@ -271,7 +259,7 @@ ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute(
           // Store the cell data into a map.
           if (m_cfg.doOutputCells) {
             std::vector<Cluster::Cell> cells;
-            for (const auto& [dParameters, simhits] :
+            for (const auto& [dParameters, simHitsIdxs] :
                  digitizeParametersResult) {
               for (const auto& cell : dParameters.cluster.channels) {
                 cells.push_back(cell);
@@ -281,21 +269,24 @@ ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute(
           }
 
           if (m_cfg.doClusterization) {
-            for (auto& [dParameters, simhits] : digitizeParametersResult) {
-              // The measurement container is unordered and the index under
-              // which the measurement will be stored is known before adding it.
-              Index measurementIdx = measurements.size();
-
-              createMeasurement(measurements, moduleGeoId, dParameters);
+            for (auto& [dParameters, simHitsIdxs] : digitizeParametersResult) {
+              auto measurement =
+                  createMeasurement(measurements, moduleGeoId, dParameters);
               clusters.emplace_back(std::move(dParameters.cluster));
-              // this digitization does hit merging so there can be more than
-              // one mapping entry for each digitized hit.
-              for (auto simHitIdx : simhits) {
+
+              measurement.allocateHitIndices(simHitsIdxs.size());
+              measurement.allocateParticleIds(simHitsIdxs.size());
+              for (auto [i, simHitIdx] : Acts::enumerate(simHitsIdxs)) {
+                measurement.hitIndices()[i] = simHitIdx;
+                measurement.particleIds()[i] =
+                    simHits.nth(simHitIdx)->particleId();
+
                 measurementParticlesMap.emplace_hint(
-                    measurementParticlesMap.end(), measurementIdx,
+                    measurementParticlesMap.end(), measurement.index(),
                     simHits.nth(simHitIdx)->particleId());
                 measurementSimHitsMap.emplace_hint(measurementSimHitsMap.end(),
-                                                   measurementIdx, simHitIdx);
+                                                   measurement.index(),
+                                                   simHitIdx);
               }
             }
           }
@@ -309,16 +300,16 @@ ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute(
         << " skipped in Digitization. Enable DEBUG mode to see more details.");
   }
 
-  if (m_cfg.doOutputCells) {
-    m_cellsWriteHandle(ctx, std::move(cellsMap));
+  if (m_cfg.doClusterization) {
+    m_outputMeasurements(ctx, std::move(measurements));
+    m_outputClusters(ctx, std::move(clusters));
+
+    m_outputMeasurementParticlesMap(ctx, std::move(measurementParticlesMap));
+    m_outputMeasurementSimHitsMap(ctx, std::move(measurementSimHitsMap));
   }
 
-  if (m_cfg.doClusterization) {
-    m_measurementWriteHandle(ctx, std::move(measurements));
-    m_clusterWriteHandle(ctx, std::move(clusters));
-    m_measurementParticlesMapWriteHandle(ctx,
-                                         std::move(measurementParticlesMap));
-    m_measurementSimHitsMapWriteHandle(ctx, std::move(measurementSimHitsMap));
+  if (m_cfg.doOutputCells) {
+    m_outputCells(ctx, std::move(cellsMap));
   }
 
   return ProcessCode::SUCCESS;
diff --git a/Examples/Framework/include/ActsExamples/EventData/Cluster.hpp b/Examples/Framework/include/ActsExamples/EventData/Cluster.hpp
index 3eb850b1214..bfd81c9f130 100644
--- a/Examples/Framework/include/ActsExamples/EventData/Cluster.hpp
+++ b/Examples/Framework/include/ActsExamples/EventData/Cluster.hpp
@@ -8,10 +8,10 @@
 
 #pragma once
 
+#include "Acts/Definitions/Algebra.hpp"
 #include "ActsFatras/Digitization/Segmentizer.hpp"
 
 #include <numeric>
-#include <optional>
 #include <vector>
 
 namespace ActsExamples {
diff --git a/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp b/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp
index ab232628341..820507564c4 100644
--- a/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp
+++ b/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp
@@ -10,25 +10,17 @@
 
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
-#include "Acts/EventData/MeasurementHelpers.hpp"
 #include "Acts/EventData/SubspaceHelpers.hpp"
 #include "Acts/EventData/Types.hpp"
-#include "Acts/EventData/detail/ParameterTraits.hpp"
-#include "Acts/EventData/detail/PrintParameters.hpp"
 #include "Acts/Geometry/GeometryIdentifier.hpp"
 #include "Acts/Utilities/Iterator.hpp"
 #include "ActsExamples/EventData/GeometryContainers.hpp"
 #include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/MeasurementConcept.hpp"
+#include "ActsExamples/EventData/SimParticle.hpp"
 
-#include <array>
-#include <compare>
-#include <concepts>
 #include <cstddef>
-#include <iosfwd>
-#include <iterator>
 #include <type_traits>
-#include <variant>
 #include <vector>
 
 #include <boost/container/static_vector.hpp>
@@ -167,6 +159,8 @@ class MeasurementContainer {
     std::size_t subspaceIndexOffset{};
     std::size_t parameterOffset{};
     std::size_t covarianceOffset{};
+    std::pair<std::size_t, std::size_t> hitIndexRange{};
+    std::pair<std::size_t, std::size_t> particleIdRange{};
     std::uint8_t size{};
   };
 
@@ -176,6 +170,8 @@ class MeasurementContainer {
   std::vector<std::uint8_t> m_subspaceIndices;
   std::vector<double> m_parameters;
   std::vector<double> m_covariances;
+  std::vector<std::size_t> m_hitIndices;
+  std::vector<SimBarcode> m_particleIds;
 
   OrderedIndices m_orderedIndices;
 };
@@ -267,6 +263,90 @@ class MeasurementProxyBase {
     return self().subspaceHelper().expandMatrix(self().covariance());
   }
 
+  bool hasHitIndices() const {
+    const auto& entry = container().m_entries.at(m_index);
+    return entry.hitIndexRange.first != entry.hitIndexRange.second;
+  }
+
+  std::span<std::size_t> allocateHitIndices(std::size_t size)
+    requires(!ReadOnly)
+  {
+    if (hasHitIndices()) {
+      throw std::runtime_error("Hit indices already allocated");
+    }
+
+    auto& entry = container().m_entries.at(m_index);
+    entry.hitIndexRange.first = container().m_hitIndices.size();
+    entry.hitIndexRange.second = size;
+    container().m_hitIndices.resize(entry.hitIndexRange.first + size);
+
+    return hitIndices();
+  }
+
+  std::span<std::size_t> hitIndices()
+    requires(!ReadOnly)
+  {
+    if (!hasHitIndices()) {
+      throw std::runtime_error("Hit indices not allocated");
+    }
+
+    const auto& entry = container().m_entries.at(m_index);
+    return {container().m_hitIndices.data() + entry.hitIndexRange.first,
+            entry.hitIndexRange.second};
+  }
+
+  std::span<const std::size_t> hitIndices() const {
+    if (!hasHitIndices()) {
+      throw std::runtime_error("Hit indices not allocated");
+    }
+
+    const auto& entry = container().m_entries.at(m_index);
+    return {container().m_hitIndices.data() + entry.hitIndexRange.first,
+            entry.hitIndexRange.second};
+  }
+
+  bool hasParticleIds() const {
+    const auto& entry = container().m_entries.at(m_index);
+    return entry.particleIdRange.first != entry.particleIdRange.second;
+  }
+
+  std::span<SimBarcode> allocateParticleIds(std::size_t size)
+    requires(!ReadOnly)
+  {
+    if (hasParticleIds()) {
+      throw std::runtime_error("Particle IDs already allocated");
+    }
+
+    auto& entry = container().m_entries.at(m_index);
+    entry.particleIdRange.first = container().m_particleIds.size();
+    entry.particleIdRange.second = size;
+    container().m_particleIds.resize(entry.particleIdRange.first + size);
+
+    return particleIds();
+  }
+
+  std::span<SimBarcode> particleIds()
+    requires(!ReadOnly)
+  {
+    if (!hasParticleIds()) {
+      throw std::runtime_error("Particle IDs not allocated");
+    }
+
+    const auto& entry = container().m_entries.at(m_index);
+    return {container().m_particleIds.data() + entry.particleIdRange.first,
+            entry.particleIdRange.second};
+  }
+
+  std::span<const SimBarcode> particleIds() const {
+    if (!hasParticleIds()) {
+      throw std::runtime_error("Particle IDs not allocated");
+    }
+
+    const auto& entry = container().m_entries.at(m_index);
+    return {container().m_particleIds.data() + entry.particleIdRange.first,
+            entry.particleIdRange.second};
+  }
+
   /// @brief Construct the measurement from a subspace vector,
   /// parameters, and covariance.
   ///
@@ -291,6 +371,16 @@ class MeasurementProxyBase {
   {
     assert(size() == other.size() && "Size mismatch");
     fill(other.subspaceIndexVector(), other.parameters(), other.covariance());
+    if (other.hasHitIndices()) {
+      allocateHitIndices(other.hitIndices().size());
+      std::copy(other.hitIndices().begin(), other.hitIndices().end(),
+                hitIndices().begin());
+    }
+    if (other.hasParticleIds()) {
+      allocateParticleIds(other.particleIds().size());
+      std::copy(other.particleIds().begin(), other.particleIds().end(),
+                particleIds().begin());
+    }
   }
 
   /// @brief Copy the data from another measurement
diff --git a/Examples/Framework/include/ActsExamples/EventData/MeasurementConcept.hpp b/Examples/Framework/include/ActsExamples/EventData/MeasurementConcept.hpp
index 14627a78c53..e959b78bb90 100644
--- a/Examples/Framework/include/ActsExamples/EventData/MeasurementConcept.hpp
+++ b/Examples/Framework/include/ActsExamples/EventData/MeasurementConcept.hpp
@@ -8,7 +8,6 @@
 
 #pragma once
 
-#include "Acts/EventData/SourceLink.hpp"
 #include "Acts/Geometry/GeometryIdentifier.hpp"
 
 #include <concepts>
@@ -22,5 +21,10 @@ concept MeasurementConcept = requires(const T& m) {
   { m.subspaceIndexVector() };
   { m.parameters() };
   { m.covariance() };
+  { m.hasHitIndices() } -> std::same_as<bool>;
+  { m.hitIndices() };
+  { m.hasParticleIds() } -> std::same_as<bool>;
+  { m.particleIds() };
 };
+
 }  // namespace ActsExamples
diff --git a/Examples/Framework/src/EventData/Measurement.cpp b/Examples/Framework/src/EventData/Measurement.cpp
index ab2dad07489..e01b417e43d 100644
--- a/Examples/Framework/src/EventData/Measurement.cpp
+++ b/Examples/Framework/src/EventData/Measurement.cpp
@@ -28,8 +28,12 @@ void MeasurementContainer::reserve(std::size_t size) {
 
 std::size_t MeasurementContainer::addMeasurement(
     std::uint8_t size, Acts::GeometryIdentifier geometryId) {
-  m_entries.push_back({m_subspaceIndices.size(), m_parameters.size(),
-                       m_covariances.size(), size});
+  m_entries.push_back({m_subspaceIndices.size(),
+                       m_parameters.size(),
+                       m_covariances.size(),
+                       {},
+                       {},
+                       size});
   m_geometryIds.emplace_back(geometryId);
   m_subspaceIndices.resize(m_subspaceIndices.size() + size);
   m_parameters.resize(m_parameters.size() + size);
diff --git a/Examples/Python/src/Digitization.cpp b/Examples/Python/src/Digitization.cpp
index 948bcc9593b..38dd8573a8a 100644
--- a/Examples/Python/src/Digitization.cpp
+++ b/Examples/Python/src/Digitization.cpp
@@ -6,7 +6,6 @@
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
-#include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Geometry/GeometryHierarchyMap.hpp"
 #include "Acts/Plugins/Python/Utilities.hpp"
 #include "Acts/Utilities/Logger.hpp"
@@ -14,7 +13,6 @@
 #include "ActsExamples/Digitization/DigitizationConfig.hpp"
 #include "ActsExamples/Digitization/DigitizationConfigurator.hpp"
 #include "ActsExamples/Digitization/DigitizationCoordinatesConverter.hpp"
-#include "ActsExamples/Framework/AlgorithmContext.hpp"
 #include "ActsExamples/Io/Json/JsonDigitizationConfig.hpp"
 
 #include <array>
@@ -65,8 +63,6 @@ void addDigitization(Context& ctx) {
     ACTS_PYTHON_MEMBER(inputSimHits);
     ACTS_PYTHON_MEMBER(outputMeasurements);
     ACTS_PYTHON_MEMBER(outputClusters);
-    ACTS_PYTHON_MEMBER(outputMeasurementParticlesMap);
-    ACTS_PYTHON_MEMBER(outputMeasurementSimHitsMap);
     ACTS_PYTHON_MEMBER(surfaceByIdentifier);
     ACTS_PYTHON_MEMBER(randomNumbers);
     ACTS_PYTHON_MEMBER(doOutputCells);