Resolve "Create a separate tracking library"

JugReco/CMakeLists.txt

@@ -9,7 +9,7 @@ find_package(EICD REQUIRED)
 #find_package(podio 0.11.01 REQUIRED)
 #include_directories(${podio_INCLUDE_DIR})
 
-find_package(Acts REQUIRED COMPONENTS Core PluginIdentification PluginTGeo PluginDD4hep )
+#find_package(Acts REQUIRED COMPONENTS Core PluginIdentification PluginTGeo PluginDD4hep )
 find_package(ROOT COMPONENTS RIO Tree Core REQUIRED)
 find_package(DD4hep COMPONENTS DDG4 DDG4IO DDRec REQUIRED)
 
@@ -29,20 +29,7 @@ gaudi_depends_on_subdirs(GaudiAlg GaudiKernel)
 
 gaudi_add_module(JugRecoPlugins
   src/components/TrackerHitReconstruction.cpp 
-  src/components/TrackerSourceLinker.cpp 
-  src/components/Tracker2SourceLinker.cpp 
-  src/components/TrackerSourcesLinker.cpp 
-  #src/components/TrackingHitsSourceLinker.cpp 
-  src/components/TrackFindingAlgorithm.cpp 
-  src/components/TrackFindingAlgorithmFunction.cpp 
-  #src/components/TrackFittingAlgorithm.cpp 
-  #src/components/TrackFittingFunction.cpp 
-  src/components/TestACTSLogger.cpp 
-  src/components/TrackParamTruthInit.cpp 
-  src/components/TrackParamClusterInit.cpp 
-  src/components/TrackParamVertexClusterInit.cpp 
   src/components/SimpleClustering.cpp 
-  src/components/ParticlesFromTrackFit.cpp 
   src/components/CalorimeterIslandCluster.cpp 
   src/components/CrystalEndcapsReco.cpp 
   src/components/EcalTungstenSamplingReco.cpp
@@ -55,17 +42,17 @@ gaudi_add_module(JugRecoPlugins
   src/components/SamplingECalHitsMerger.cpp 
   src/components/TopologicalCellCluster.cpp 
   src/components/ImagingClusterReco.cpp 
-  LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO EICD::eicd DDRec Acts)
+  LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO EICD::eicd DDRec )
 
 target_compile_options(JugRecoPlugins PRIVATE -Wno-suggest-override)
 
 #gaudi_add_test(ProduceForReadTest
 #               WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
 #               FRAMEWORK tests/options/simple_producer.py)
-gaudi_add_test(loggerTest
-               WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
-               FRAMEWORK tests/options/testLogger.py
-               )
+#gaudi_add_test(loggerTest
+#               WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
+#               FRAMEWORK tests/options/testLogger.py
+#               )
 #gaudi_add_test(CheckReadCollectionSize
 #               ENVIRONMENT PYTHONPATH+=${PODIO_PYTHON_DIR}
 #               WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}

JugTrack/CMakeLists.txt

+################################################################################
+# Package: JugTrack
+################################################################################
+gaudi_subdir(JugTrack v1r0)
+
+find_package(EICD REQUIRED)
+#set(PODIO $ENV{PODIO})
+#set(CMAKE_MODULE_PATH CMAKE_MODULE_PATH PODIO)
+#find_package(podio 0.11.01 REQUIRED)
+#include_directories(${podio_INCLUDE_DIR})
+
+find_package(Acts REQUIRED COMPONENTS Core PluginIdentification PluginTGeo PluginDD4hep )
+find_package(ROOT COMPONENTS RIO Tree Core REQUIRED)
+find_package(DD4hep COMPONENTS DDG4 DDG4IO DDRec REQUIRED)
+
+# this declaration will not be needed in the future
+gaudi_depends_on_subdirs(GaudiAlg GaudiKernel)
+#gaudi_install_scripts()
+#gaudi_install_python_modules()
+
+#gaudi_add_library(JugBase
+#  src/*.cpp 
+#                  INCLUDE_DIRS EICD PODIO ROOT $ENV{HOME}/stow/podio/include
+#                  LINK_LIBRARIES GaudiAlgLib GaudiKernel ROOT DD4hep::DDG4IO
+#                  PUBLIC_HEADERS JugBase)
+#target_link_libraries(JugBase 
+#  podio::podioRootIO
+#  )
+
+gaudi_add_module(JugTrackPlugins
+  src/components/TrackerSourceLinker.cpp 
+  src/components/Tracker2SourceLinker.cpp 
+  src/components/TrackerSourcesLinker.cpp 
+  src/components/TrackFindingAlgorithm.cpp 
+  src/components/TrackFindingAlgorithmFunction.cpp 
+  #src/components/TrackFittingAlgorithm.cpp 
+  #src/components/TrackFittingFunction.cpp 
+  src/components/TestACTSLogger.cpp 
+  src/components/TrackParamTruthInit.cpp 
+  src/components/TrackParamClusterInit.cpp 
+  src/components/TrackParamVertexClusterInit.cpp 
+  src/components/ParticlesFromTrackFit.cpp 
+  LINK_LIBRARIES GaudiAlgLib GaudiKernel JugBase ROOT NPDet::DD4podIO EICD::eicd DDRec Acts)
+
+target_compile_options(JugTrackPlugins PRIVATE -Wno-suggest-override)
+
+#gaudi_add_test(loggerTest
+#               WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
+#               FRAMEWORK tests/options/testLogger.py
+#               )
+#gaudi_add_test(CheckReadCollectionSize
+#               ENVIRONMENT PYTHONPATH+=${PODIO_PYTHON_DIR}
+#               WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
+#               COMMAND python FWCore/tests/scripts/check_coll_after_read.py
+#               DEPENDS ReadTest)

JugTrack/JugTrack/ACTSLogger.h

+
+
+#ifndef RECTRACKER_ACTSLOGGER_H
+#define RECTRACKER_ACTSLOGGER_H
+
+#include "Acts/Utilities/Logger.hpp"
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/IMessageSvc.h"
+#include "GaudiKernel/MsgStream.h"
+
+class GaudiFilterPolicy : public Acts::Logging::OutputFilterPolicy {
+public:
+  GaudiFilterPolicy(IMessageSvc* owner) : m_messenger(owner), m_currentLevel(m_messenger.currentLevel()) {}
+
+  bool doPrint(const Acts::Logging::Level& lvl) const {
+    MSG::Level l = MSG::VERBOSE;
+    switch (lvl) {
+    case Acts::Logging::VERBOSE:
+      l = MSG::VERBOSE;
+      break;
+    case Acts::Logging::DEBUG:
+      l = MSG::DEBUG;
+      break;
+    case Acts::Logging::INFO:
+      l = MSG::INFO;
+      break;
+    case Acts::Logging::WARNING:
+      l = MSG::WARNING;
+      break;
+    case Acts::Logging::ERROR:
+      l = MSG::ERROR;
+      break;
+    case Acts::Logging::FATAL:
+      l = MSG::FATAL;
+      break;
+    }
+    MSG::Level cl = m_currentLevel;
+    return l < cl;
+  }
+
+private:
+  MsgStream m_messenger;
+  MSG::Level m_currentLevel;
+};
+
+class GaudiPrintPolicy : public Acts::Logging::OutputPrintPolicy {
+public:
+  GaudiPrintPolicy(IMessageSvc* owner) : m_messenger(owner) {}
+
+  void flush(const Acts::Logging::Level& lvl, const std::ostringstream& input) {
+    MSG::Level l = MSG::VERBOSE;
+    switch (lvl) {
+    case Acts::Logging::VERBOSE:
+      l = MSG::VERBOSE;
+      break;
+    case Acts::Logging::DEBUG:
+      l = MSG::DEBUG;
+      break;
+    case Acts::Logging::INFO:
+      l = MSG::INFO;
+      break;
+    case Acts::Logging::WARNING:
+      l = MSG::WARNING;
+      break;
+    case Acts::Logging::ERROR:
+      l = MSG::ERROR;
+      break;
+    case Acts::Logging::FATAL:
+      l = MSG::FATAL;
+      break;
+    }
+
+    m_messenger << l << input.str() << endmsg;
+  }
+
+private:
+  MsgStream m_messenger;
+};
+
+#endif  // RECTRACKER_ACTSLOGGER_H

JugTrack/JugTrack/BField.h

+#ifndef Jug_BFIELD_HH
+#define Jug_BFIELD_HH 1
+
+#include "Acts//Definitions/Units.hpp"
+#include "Acts/Utilities/detail/AxisFwd.hpp"
+#include "Acts/Utilities/detail/GridFwd.hpp"
+#include <memory>
+#include <tuple>
+#include <variant>
+
+#include "Acts/MagneticField/MagneticFieldContext.hpp"
+#include "Acts/Definitions/Common.hpp"
+
+// Forward declarations
+namespace Acts {
+  template <typename G>
+  struct InterpolatedBFieldMapper;
+  template <typename M>
+  class InterpolatedBFieldMap;
+  class ConstantBField;
+
+} // namespace Acts
+
+namespace Jug {
+  // namespace BField {
+  //  class ScalableBField;
+  //}
+  namespace BField {
+
+    /// The Context to be handed around
+    struct ScalableBFieldContext {
+      double scalor = 1.;
+    };
+
+    /// @ingroup MagneticField
+    ///
+    /// @brief returns a given constant field value at every point
+    ///
+    /// This class is based on the constant magnetic field class
+    /// but allows a event based context
+    class ScalableBField final {
+    public:
+      struct Cache {
+        double scalor = 1.;
+
+        /// @brief constructor with context
+        Cache(const Acts::MagneticFieldContext& mcfg)
+        {
+          scalor = std::any_cast<const ScalableBFieldContext>(mcfg).scalor;
+        }
+      };
+
+      /// @brief construct constant magnetic field from field vector
+      ///
+      /// @param [in] B magnetic field vector in global coordinate system
+      explicit ScalableBField(Acts::Vector3 B) : m_BField(std::move(B)) {}
+
+      /// @brief construct constant magnetic field from components
+      ///
+      /// @param [in] Bx magnetic field component in global x-direction
+      /// @param [in] By magnetic field component in global y-direction
+      /// @param [in] Bz magnetic field component in global z-direction
+      ScalableBField(double Bx = 0., double By = 0., double Bz = 0.) : m_BField(Bx, By, Bz) {}
+
+      /// @brief retrieve magnetic field value
+      ///
+      /// @param [in] position global position
+      /// @return magnetic field vector
+      ///
+      /// @note The @p position is ignored and only kept as argument to provide
+      ///       a consistent interface with other magnetic field services.
+      Acts::Vector3 getField(const Acts::Vector3& /*position*/) const { return m_BField; }
+
+      /// @brief retrieve magnetic field value
+      ///
+      /// @param [in] position global position
+      /// @param [in] cache Cache object (is ignored)
+      /// @return magnetic field vector
+      ///
+      /// @note The @p position is ignored and only kept as argument to provide
+      ///       a consistent interface with other magnetic field services.
+      Acts::Vector3 getField(const Acts::Vector3& /*position*/, Cache& cache) const
+      {
+        return m_BField * cache.scalor;
+      }
+
+      /// @brief retrieve magnetic field value & its gradient
+      ///
+      /// @param [in]  position   global position
+      /// @param [out] derivative gradient of magnetic field vector as (3x3)
+      /// matrix
+      /// @return magnetic field vector
+      ///
+      /// @note The @p position is ignored and only kept as argument to provide
+      ///       a consistent interface with other magnetic field services.
+      /// @note currently the derivative is not calculated
+      /// @todo return derivative
+      Acts::Vector3 getFieldGradient(const Acts::Vector3& /*position*/, Acts::ActsMatrix<3, 3>& /*derivative*/) const
+      {
+        return m_BField;
+      }
+
+      /// @brief retrieve magnetic field value & its gradient
+      ///
+      /// @param [in]  position   global position
+      /// @param [out] derivative gradient of magnetic field vector as (3x3)
+      /// matrix
+      /// @param [in] cache Cache object (is ignored)
+      /// @return magnetic field vector
+      ///
+      /// @note The @p position is ignored and only kept as argument to provide
+      ///       a consistent interface with other magnetic field services.
+      /// @note currently the derivative is not calculated
+      /// @todo return derivative
+      Acts::Vector3 getFieldGradient(const Acts::Vector3& /*position*/, Acts::ActsMatrix<3, 3>& /*derivative*/,
+                                      Cache& cache) const
+      {
+        return m_BField * cache.scalor;
+      }
+
+      /// @brief check whether given 3D position is inside look-up domain
+      ///
+      /// @param [in] position global 3D position
+      /// @return @c true if position is inside the defined look-up grid,
+      ///         otherwise @c false
+      /// @note The method will always return true for the constant B-Field
+      bool isInside(const Acts::Vector3& /*position*/) const { return true; }
+
+      /// @brief update magnetic field vector from components
+      ///
+      /// @param [in] Bx magnetic field component in global x-direction
+      /// @param [in] By magnetic field component in global y-direction
+      /// @param [in] Bz magnetic field component in global z-direction
+      void setField(double Bx, double By, double Bz) { m_BField << Bx, By, Bz; }
+
+      /// @brief update magnetic field vector
+      ///
+      /// @param [in] B magnetic field vector in global coordinate system
+      void setField(const Acts::Vector3& B) { m_BField = B; }
+
+    private:
+      /// magnetic field vector
+      Acts::Vector3 m_BField;
+    };
+
+}  // namespace BField
+} // namespace Jug
+
+using InterpolatedMapper2D = Acts::InterpolatedBFieldMapper<
+    Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis>>;
+
+using InterpolatedMapper3D =
+    Acts::InterpolatedBFieldMapper<Acts::detail::Grid<Acts::Vector3, Acts::detail::EquidistantAxis,
+                                                      Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis>>;
+using InterpolatedBFieldMap2D = Acts::InterpolatedBFieldMap<InterpolatedMapper2D>;
+using InterpolatedBFieldMap3D = Acts::InterpolatedBFieldMap<InterpolatedMapper3D>;
+
+using BFieldVariant = std::variant<std::shared_ptr<InterpolatedBFieldMap2D>, std::shared_ptr<InterpolatedBFieldMap3D>,
+                                   std::shared_ptr<Acts::ConstantBField>, std::shared_ptr<Jug::BField::ScalableBField>>;
+
+#endif

JugTrack/JugTrack/GeometryContainers.hpp

+// This file is part of the Acts project.
+//
+// Copyright (C) 2017-2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "JugTrack/Utilities/GroupBy.hpp"
+#include "JugTrack/Utilities/Range.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+
+#include <algorithm>
+#include <cstddef>
+#include <utility>
+
+#include <boost/container/flat_map.hpp>
+#include <boost/container/flat_set.hpp>
+
+namespace Jug {
+namespace detail {
+// extract the geometry identifier from a variety of types
+struct GeometryIdGetter {
+  // explicit geometry identifier are just forwarded
+  constexpr Acts::GeometryIdentifier operator()(Acts::GeometryIdentifier geometryId) const {
+    return geometryId;
+  }
+  // encoded geometry ids are converted back to geometry identifiers.
+  constexpr Acts::GeometryIdentifier operator()(Acts::GeometryIdentifier::Value encoded) const {
+    return Acts::GeometryIdentifier(encoded);
+  }
+  // support elements in map-like structures.
+  template <typename T>
+  constexpr Acts::GeometryIdentifier operator()(
+      const std::pair<Acts::GeometryIdentifier, T>& mapItem) const {
+    return mapItem.first;
+  }
+  // support elements that implement `.geometryId()`.
+  template <typename T>
+  inline auto operator()(const T& thing) const
+      -> decltype(thing.geometryId(), Acts::GeometryIdentifier()) {
+    return thing.geometryId();
+  }
+};
+
+struct CompareGeometryId {
+  // indicate that comparisons between keys and full objects are allowed.
+  using is_transparent = void;
+  // compare two elements using the automatic key extraction.
+  template <typename Left, typename Right>
+  constexpr bool operator()(Left&& lhs, Right&& rhs) const {
+    return GeometryIdGetter()(lhs) < GeometryIdGetter()(rhs);
+  }
+};
+}  // namespace detail
+
+/// Store elements that know their detector geometry id, e.g. simulation hits.
+///
+/// @tparam T type to be stored, must be compatible with `CompareGeometryId`
+///
+/// The container stores an arbitrary number of elements for any geometry
+/// id. Elements can be retrieved via the geometry id; elements can be selected
+/// for a specific geometry id or for a larger range, e.g. a volume or a layer
+/// within the geometry hierachy using the helper functions below. Elements can
+/// also be accessed by index that uniquely identifies each element regardless
+/// of geometry id.
+template <typename T>
+using GeometryIdMultiset =
+    boost::container::flat_multiset<T, detail::CompareGeometryId>;
+
+/// Store elements indexed by an geometry id.
+///
+/// @tparam T type to be stored
+///
+/// The behaviour is the same as for the `GeometryIdMultiset` except that the
+/// stored elements do not know their geometry id themself. When iterating
+/// the iterator elements behave as for the `std::map`, i.e.
+///
+///     for (const auto& entry: elements) {
+///         auto id = entry.first; // geometry id
+///         const auto& el = entry.second; // stored element
+///     }
+///
+template <typename T>
+using GeometryIdMultimap = GeometryIdMultiset<std::pair<Acts::GeometryIdentifier, T>>;
+
+/// Select all elements within the given volume.
+template <typename T>
+inline Range<typename GeometryIdMultiset<T>::const_iterator> selectVolume(
+    const GeometryIdMultiset<T>& container, Acts::GeometryIdentifier::Value volume) {
+  auto cmp = Acts::GeometryIdentifier().setVolume(volume);
+  auto beg = std::lower_bound(container.begin(), container.end(), cmp,
+                              detail::CompareGeometryId{});
+  // WARNING overflows to volume==0 if the input volume is the last one
+  cmp = Acts::GeometryIdentifier().setVolume(volume + 1u);
+  // optimize search by using the lower bound as start point. also handles
+  // volume overflows since the geo id would be located before the start of
+  // the upper edge search window.
+  auto end =
+      std::lower_bound(beg, container.end(), cmp, detail::CompareGeometryId{});
+  return makeRange(beg, end);
+}
+template <typename T>
+inline auto selectVolume(const GeometryIdMultiset<T>& container,
+                         Acts::GeometryIdentifier id) {
+  return selectVolume(container, id.volume());
+}
+
+/// Select all elements within the given layer.
+template <typename T>
+inline Range<typename GeometryIdMultiset<T>::const_iterator> selectLayer(
+    const GeometryIdMultiset<T>& container, Acts::GeometryIdentifier::Value volume,
+    Acts::GeometryIdentifier::Value layer) {
+  auto cmp = Acts::GeometryIdentifier().setVolume(volume).setLayer(layer);
+  auto beg = std::lower_bound(container.begin(), container.end(), cmp,
+                              detail::CompareGeometryId{});
+  // WARNING resets to layer==0 if the input layer is the last one
+  cmp = Acts::GeometryIdentifier().setVolume(volume).setLayer(layer + 1u);
+  // optimize search by using the lower bound as start point. also handles
+  // volume overflows since the geo id would be located before the start of
+  // the upper edge search window.
+  auto end =
+      std::lower_bound(beg, container.end(), cmp, detail::CompareGeometryId{});
+  return makeRange(beg, end);
+}
+template <typename T>
+inline auto selectLayer(const GeometryIdMultiset<T>& container,
+                        Acts::GeometryIdentifier id) {
+  return selectLayer(container, id.volume(), id.layer());
+}
+
+/// Select all elements for the given module / sensitive surface.
+template <typename T>
+inline Range<typename GeometryIdMultiset<T>::const_iterator> selectModule(
+    const GeometryIdMultiset<T>& container, Acts::GeometryIdentifier geoId) {
+  // module is the lowest level and defines a single geometry id value
+  return makeRange(container.equal_range(geoId));
+}
+template <typename T>
+inline auto selectModule(const GeometryIdMultiset<T>& container,
+                         Acts::GeometryIdentifier::Value volume,
+                         Acts::GeometryIdentifier::Value layer,
+                         Acts::GeometryIdentifier::Value module) {
+  return selectModule(
+      container,
+      Acts::GeometryIdentifier().setVolume(volume).setLayer(layer).setSensitive(
+          module));
+}
+
+/// Iterate over groups of elements belonging to each module/ sensitive surface.
+template <typename T>
+inline GroupBy<typename GeometryIdMultiset<T>::const_iterator,
+               detail::GeometryIdGetter>
+groupByModule(const GeometryIdMultiset<T>& container) {
+  return makeGroupBy(container, detail::GeometryIdGetter());
+}
+
+}  // namespace FW

JugTrack/JugTrack/Index.hpp

+// This file is part of the Acts project.
+//
+// Copyright (C) 2019-2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include <cstdint>
+
+#include <boost/container/flat_map.hpp>
+
+namespace Jug {
+
+/// Index type to reference elements in a container.
+///
+/// We do not expect to have more than 2^32 elements in any given container so a
+/// fixed sized integer type is sufficient.
+using Index = uint32_t;
+
+/// Store elements that are identified by an index, e.g. in another container.
+///
+/// Each index can have zero or more associated elements. A typical case could
+/// be to store all generating particles for a hit where the hit is identified
+/// by its index in the hit container.
+template <typename value_t>
+using IndexMultimap = boost::container::flat_multimap<Index, value_t>;
+
+/// Invert the multimap, i.e. from a -> {b...} to b -> {a...}.
+///
+/// @note This assumes that the value in the initial multimap is itself a
+///   sortable index-like object, as would be the case when mapping e.g.
+///   hit ids to particle ids/ barcodes.
+template <typename value_t>
+inline boost::container::flat_multimap<value_t, Index> invertIndexMultimap(
+    const IndexMultimap<value_t>& multimap) {
+  using InverseMultimap = boost::container::flat_multimap<value_t, Index>;
+
+  // switch key-value without enforcing the new ordering (linear copy)
+  typename InverseMultimap::sequence_type unordered;
+  unordered.reserve(multimap.size());
+  for (auto&& [index, value] : multimap) {
+    // value is now the key and the index is now the value
+    unordered.emplace_back(value, index);
+  }
+
+  // adopting the unordered sequence will reestablish the correct order
+  InverseMultimap inverse;
+  inverse.adopt_sequence(std::move(unordered));
+  return inverse;
+}
+
+}  // namespace ActsExamples

JugTrack/JugTrack/IndexContainers.hpp

+// This file is part of the Acts project.
+//
+// Copyright (C) 2019-2020 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include <cstddef>
+#include <utility>
+
+#include <boost/container/flat_map.hpp>
+
+namespace FW {
+
+/// Store elements that are identified by an index, e.g. in another container.
+///
+/// Each index can have zero or more associated elements. A typical case could
+/// be to store all generating particles for a hit where the hit is identified
+/// by its index in the hit container.
+template <typename Value, typename Key = std::size_t>
+using IndexMultimap = boost::container::flat_multimap<Key, Value>;
+
+/// Invert the multimap, i.e. from a -> {b...} to b -> {a...}.
+///
+/// @note This assumes that the value in the initial multimap is itself a
+///       sortable index-like object, as would be the case when mapping e.g.
+///       hit ids to particle ids/ barcodes.
+template <typename Value, typename Key>
+inline IndexMultimap<Key, Value> invertIndexMultimap(
+    const IndexMultimap<Value, Key>& multimap) {
+  // switch key-value without enforcing the new ordering (linear copy)
+  typename IndexMultimap<Key, Value>::sequence_type unordered;
+  unordered.reserve(multimap.size());
+  for (const auto& keyValue : multimap) {
+    // value is now the key and the key is now the value
+    unordered.emplace_back(keyValue.second, keyValue.first);
+  }
+  // adopting the unordered sequence will reestablish the correct order
+  IndexMultimap<Key, Value> inverse;
+  inverse.adopt_sequence(std::move(unordered));