[wip] add MeasurementND classes

edm4hep.yaml

@@ -418,6 +418,50 @@ datatypes :
       //TODO: edm4hep::Vertex  getEndVertex() { return  edm4hep::Vertex(  (getParticles(0).isAvailable() ? getParticles(0).getStartVertex() :  edm4hep::Vertex(0,0) ) ) ; }\n     
       "
 
+  edm4hep::Measurement2D:
+    Description: "Two dimensional measurement with covariance matrix"
+    Author: "T. Madlener, DESY"
+    Members:
+      - std::array<float, 2> vals    // The values of the measurement
+      - std::array<float, 3> covMatrix     // The (lower triangle) covariance matrix of the measurement
+      - std::array<uint16_t, 2> dims // The (encded) dimensions of the two values
+    ExtraCode:
+      includes: '#include "edm4hep/detail/enum_access.h"'
+      declaration: "
+      /// Get the value for the passed dimension\n
+      template<typename DimEnum>\n
+      float getValue(DimEnum e) const {\n
+        return detail::get_value(getVals(), e, getDims());\n
+      }\n
+      /// Get the value of the covariance matrix for the two passed dimensions\n
+      template<typename DimEnum>\n
+      float getCov(DimEnum ei, DimEnum ej) const {\n
+        return detail::get_cov(getCovMatrix(), ei, ej, getDims());\n
+      }\n
+      /// Get which dimensions are currently set for this measurement\n
+      template<typename DimEnum>\n
+      std::array<DimEnum, 2> getDimensions() const {\n
+        return {detail::to_enum<DimEnum>(getDims()[0]), detail::to_enum<DimEnum>(getDims()[1])};\n
+      }
+      "
+    MutableExtraCode:
+      includes: '#include "edm4hep/detail/enum_access.h"'
+      declaration: "
+      /// Set the value for the passed dimension\n
+      template<typename DimEnum>\n
+      void setValue(float val, DimEnum dim) {\n
+        detail::set_value(val, vals(), dim, getDims());\n
+      }\n
+      template<typename DimEnum>\n
+      void setCov(float val, DimEnum dimI, DimEnum dimJ) {\n
+        detail::set_cov(val, covMatrix(), dimI, dimJ, getDims());\n
+      }\n
+      template<typename DimEnum>\n
+      void setDimensions(DimEnum dim1, DimEnum dim2) {\n
+        setDims({detail::to_index(dim1), detail::to_index(dim2)});\n
+      }
+      "
+
   edm4hep::MCRecoParticleAssociation:
     Description: "Used to keep track of the correspondence between MC and reconstructed particles"
     Author: "C. Bernet, B. Hegner"

edm4hep/CMakeLists.txt

@@ -4,6 +4,9 @@
 PODIO_GENERATE_DATAMODEL(edm4hep ../edm4hep.yaml headers sources IO_BACKEND_HANDLERS ${PODIO_IO_HANDLERS})
 
 PODIO_ADD_DATAMODEL_CORE_LIB(edm4hep "${headers}" "${sources}")
+target_include_directories(edm4hep PUBLIC
+    $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/include>
+  )
 
 PODIO_ADD_ROOT_IO_DICT(edm4hepDict edm4hep "${headers}" src/selection.xml)
 add_library(edm4hep::edm4hepDict ALIAS edm4hepDict )

include/edm4hep/detail/enum_access.h

@@ -0,0 +1,112 @@
+#ifndef EDM4HEP_DETAIL_ENUM_ACCESS_H
+#define EDM4HEP_DETAIL_ENUM_ACCESS_H
+
+#include <algorithm>
+#include <array>
+#include <cstdint>
+#include <type_traits>
+
+namespace edm4hep {
+
+// Use 16 bits to encode the dimension
+// Could go to 8 bits, but would need a fix in podio first
+using DimType = std::uint16_t;
+
+namespace detail {
+// From c++23 this is functionality offerd by the STL
+#if __cpp_lib_to_underlying
+using to_index = std::to_underlying;
+#else
+// Otherwise it is simple enough to roll our own
+template <typename E> constexpr auto to_index(E e) {
+  return static_cast<std::underlying_type_t<E>>(e);
+}
+#endif
+
+/// Cast an index to an enum value
+template <typename DimEnum> constexpr DimEnum to_enum(DimType index) {
+  return static_cast<DimEnum>(index);
+}
+
+/// Decode which index in the array corresponds to the passed enum value.
+/// Assuming that they are unique!
+template <typename DimEnum, std::size_t N>
+constexpr int decode_index(const std::array<DimType, N> &dimIndices,
+                           DimEnum dim) {
+  for (auto i = 0u; i < N; ++i) {
+    if (dimIndices[i] == detail::to_index(dim)) {
+      return i;
+    }
+  }
+
+  return -1;
+}
+
+/// Get the value corresponding to the passed enum value from the array
+template <typename DimEnum, typename Scalar, std::size_t N>
+constexpr Scalar get_value(const std::array<Scalar, N> &values, DimEnum dim,
+                           const std::array<DimType, N> &dimIndices) {
+  const auto index = decode_index(dimIndices, dim);
+  if (index < 0) {
+    // TODO: error handling
+  }
+
+  return values[index];
+}
+
+template <typename DimEnum, typename Scalar, std::size_t N>
+constexpr void set_value(Scalar value, std::array<Scalar, N> &values,
+                         DimEnum dim,
+                         const std::array<DimType, N> &dimIndices) {
+  const auto index = decode_index(dimIndices, dim);
+  if (index < 0) {
+    // TODO: error handling
+  }
+
+  values[index] = value;
+}
+
+/// Convert from 2D matrix indices to 1d lower triangle indices
+constexpr int to_lower_tri(int i, int j, size_t N) {
+  if (j < i) {
+    std::swap(i, j);
+  }
+  return i * (2 * N - i - 1) / 2 + j;
+}
+
+/// Get the covariance matrix value corresponding to the passed enum values from
+/// the array
+template <typename DimEnum, typename Scalar, std::size_t N, std::size_t NDims>
+constexpr Scalar get_cov(const std::array<Scalar, N> &cov, DimEnum dimI,
+                         DimEnum dimJ,
+                         const std::array<DimType, NDims> &dimIndices) {
+  auto i = decode_index(dimIndices, dimI);
+  auto j = decode_index(dimIndices, dimJ);
+
+  if (i < 0 || j < 0) {
+    // TODO: error handling
+  }
+
+  return cov[to_lower_tri(i, j, NDims)];
+}
+
+template <typename DimEnum, typename Scalar, std::size_t N, std::size_t NDims>
+constexpr void set_cov(Scalar value, std::array<Scalar, N> &cov, DimEnum dimI,
+                       DimEnum dimJ,
+                       const std::array<DimType, NDims> &dimIndices) {
+  auto i = decode_index(dimIndices, dimI);
+  auto j = decode_index(dimIndices, dimJ);
+
+  if (i < 0 || j < 0) {
+    // TODO: error handling
+  }
+
+  // Covariance is in lower triangle this brings us from 2D indices to 1D
+  cov[to_lower_tri(i, j, NDims)] = value;
+}
+
+} // namespace detail
+
+} // namespace edm4hep
+
+#endif // EDM4HEP_DETAIL_ENUM_ACCESS_H

test/CMakeLists.txt

@@ -4,7 +4,7 @@ ENDIF()
 
 function(set_test_env _testname)
   set_property(TEST ${_testname} APPEND PROPERTY ENVIRONMENT
-      LD_LIBRARY_PATH=$<TARGET_FILE_DIR:edm4hep>:$<TARGET_FILE_DIR:podio::podio>:$ENV{LD_LIBRARY_PATH}
+      LD_LIBRARY_PATH=$<TARGET_FILE_DIR:edm4hep>:$<TARGET_FILE_DIR:podio::podio>:$<TARGET_FILE_DIR:ROOT::Tree>:$ENV{LD_LIBRARY_PATH}
       ROOT_INCLUDE_PATH=${PROJECT_SOURCE_DIR}/edm4hep:$ENV{ROOT_INCLUDE_PATH}
   )
   set_tests_properties(${_testname} PROPERTIES

test/dimensions.h

@@ -0,0 +1,10 @@
+#ifndef EDM4HEP_TEST_DIMENSIONS_H
+#define EDM4HEP_TEST_DIMENSIONS_H
+
+#include "edm4hep/detail/enum_access.h"
+
+enum Cartesian : edm4hep::DimType { X = 0, Y = 1, Z = 2 };
+
+enum class Polar : edm4hep::DimType { R = 0, PHI, THETA };
+
+#endif // EDM4HEP_TEST_DIMENSIONS_H

test/read_events.h

@@ -1,11 +1,14 @@
 #ifndef EDM4HEP_TEST_READ_EVENTS_H__
 #define EDM4HEP_TEST_READ_EVENTS_H__
 
+#include "dimensions.h"
+
 // test data model
 #include "edm4hep/MCParticleCollection.h"
 #include "edm4hep/SimTrackerHitCollection.h"
 #include "edm4hep/CaloHitContributionCollection.h"
 #include "edm4hep/SimCalorimeterHitCollection.h"
+#include "edm4hep/Measurement2DCollection.h"
 
 // podio specific includes
 #include "podio/EventStore.h"
@@ -16,6 +19,56 @@
 #include <exception>
 #include <cassert>
 
+
+bool checkMeasurements(const edm4hep::Measurement2DCollection& coll) {
+  const auto m1 = coll[0];
+  const auto [m1d1, m1d2] = m1.getDimensions<Cartesian>();
+  if (m1d1 != Cartesian::X) return false;
+  if (m1d2 != Cartesian::Z) return false;
+
+  if (m1.getValue(Cartesian::Z) != 3.14f) return false;
+  if (m1.getValue(Cartesian::X) != 1.234f) return false;
+  if (m1.getCov(Cartesian::X, Cartesian::X) != 42.0f) return false;
+  // Check that also transposition works
+  if (m1.getCov(Cartesian::Z, Cartesian::X) != 3.14f) return false;
+  if (m1.getCov(Cartesian::Z, Cartesian::Z) != 2.34f) return false;
+
+
+  const auto m2 = coll[1];
+  // We can also work with the array here directly
+  const auto& m2Dims = m2.getDimensions<Polar>();
+  if (m2Dims[0] != Polar::PHI) return false;
+  if (m2Dims[1] != Polar::THETA) return false;
+
+  if (m2.getValue(Polar::THETA) != 200.f) return false;
+  if (m2.getValue(Polar::PHI) != 100.f) return false;
+
+  if (m2.getCov(Polar::THETA, Polar::THETA) != 30.f) return false;
+  if (m2.getCov(Polar::THETA, Polar::PHI) != 20.f) return false;
+  if (m2.getCov(Polar::PHI, Polar::PHI) != 10.f) return false;
+
+
+  const auto m3 = coll[2];
+
+  // As is the case with all instances we can access this with arbitrary enums
+  const auto m3DimsPolar = m3.getDimensions<Polar>();
+  if (m3DimsPolar[0] != Polar::THETA) return false;
+  if (m3DimsPolar[1] != Polar::PHI) return false;
+
+  const auto m3DimsCart = m3.getDimensions<Cartesian>();
+  if (m3DimsCart[0] != Cartesian::Z) return false;
+  if (m3DimsCart[1] != Cartesian::Y) return false;
+
+  if (m3.getValue(Polar::THETA) != 1.234f) return false;
+  if (m3.getValue(Cartesian::Y) != 5.678f) return false;
+
+  if (m3.getCov(Polar::THETA, Polar::THETA) != 1.1f) return false;
+  if (m3.getCov(Cartesian::Y, Cartesian::Y) != 3.3f) return false;
+  if (m3.getCov(Polar::THETA, Polar::PHI) != m3.getCov(Cartesian::Z, Cartesian::Y)) return false;
+
+  return true;
+}
+
 void processEvent(podio::EventStore& store, bool verboser, unsigned eventNum) {
   auto& mcps   = store.get<edm4hep::MCParticleCollection>("MCParticles");
   auto& sths   = store.get<edm4hep::SimTrackerHitCollection>("SimTrackerHits");
@@ -159,6 +212,11 @@ void processEvent(podio::EventStore& store, bool verboser, unsigned eventNum) {
  //    throw std::runtime_error("Collection 'SimCalorimeterHitContributions' should be present");
  //  }
 
+  auto& measurements = store.get<edm4hep::Measurement2DCollection>("measurement2D");
+  if (!checkMeasurements(measurements)) {
+    throw std::runtime_error("Measurement not as expected");
+  }
+
 }
 
 template<typename ReaderT>

test/write_events.h

@@ -1,11 +1,14 @@
 #ifndef EDM4HEP_TEST_WRITE_EVENTS_H__
 #define EDM4HEP_TEST_WRITE_EVENTS_H__
 
+#include "dimensions.h"
+
 // Data model
 #include "edm4hep/MCParticleCollection.h"
 #include "edm4hep/SimTrackerHitCollection.h"
 #include "edm4hep/CaloHitContributionCollection.h"
 #include "edm4hep/SimCalorimeterHitCollection.h"
+#include "edm4hep/Measurement2DCollection.h"
 
 // STL
 #include <iostream>
@@ -14,6 +17,32 @@
 // podio specific includes
 #include "podio/EventStore.h"
 
+void fillCollection(edm4hep::Measurement2DCollection& coll) {
+  auto m1 = coll.create();
+  // Set things via the tacked on interface
+  m1.setDimensions(Cartesian::X, Cartesian::Z);
+  m1.setValue(1.234f, Cartesian::X);
+  m1.setValue(3.14f, Cartesian::Z);
+  m1.setCov(42.0f, Cartesian::X, Cartesian::X);
+  m1.setCov(3.14f, Cartesian::X, Cartesian::Z);
+  m1.setCov(2.34f, Cartesian::Z, Cartesian::Z);
+
+  // Again but with enum class
+  auto m2 = coll.create();
+  m2.setDimensions(Polar::PHI, Polar::THETA);
+  m2.setValue(100.f, Polar::PHI);
+  m2.setValue(200.f, Polar::THETA);
+  m2.setCov(10.f, Polar::PHI, Polar::PHI);
+  m2.setCov(20.f, Polar::THETA, Polar::PHI);
+  m2.setCov(30.f, Polar::THETA, Polar::THETA);
+
+  // We can also go through the "native" interface
+  auto m3 = coll.create();
+  m3.setDims({2, 1});
+  m3.setVals({1.234f, 5.678f});
+  m3.setCovMatrix({1.1f, 2.2f, 3.3f});
+}
+
 template<class WriterT>
 void write(std::string outfilename) {
   std::cout<<"start processing"<<std::endl;
@@ -34,6 +63,8 @@ void write(std::string outfilename) {
   auto& sccons = store.create<edm4hep::CaloHitContributionCollection>("SimCalorimeterHitContributions");
   writer.registerForWrite("SimCalorimeterHitContributions");
 
+  auto& measurements = store.create<edm4hep::Measurement2DCollection>("measurement2D");
+  writer.registerForWrite("measurement2D");
 
   unsigned nevents = 10 ;
 
@@ -181,6 +212,8 @@ void write(std::string outfilename) {
     std::cout << "\n collection:  " << "SimCalorimeterHits" <<  " of type " <<  schs.getValueTypeName() << "\n\n"
         << schs << std::endl ;
 
+    fillCollection(measurements);
+
     //===============================================================================
 
     writer.writeEvent();