diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementInputConverter.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementInputConverter.hpp index 947bd561bd3..372f2a18a4f 100644 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementInputConverter.hpp +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementInputConverter.hpp @@ -12,33 +12,33 @@ #include "ActsExamples/EventData/Cluster.hpp" #include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/Framework/DataHandle.hpp" +#include "ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp" #include "ActsExamples/Io/Podio/PodioInputConverter.hpp" +#include #include namespace ActsExamples { -/// Read in a measurement cluster collection as EDM4hep from a @c podio::Frame. -/// -/// Inpersistent information: -/// - hit index -/// - 1D local coords? -/// - segment path -/// -/// Known issues: -/// - cluster channels are read from inappropriate fields -/// - local 2D coordinates and time are read from position +class DD4hepDetector; + +/// Read in a measurement collection from EDM4hep TrackerHitLocal format. class EDM4hepMeasurementInputConverter final : public PodioInputConverter { public: struct Config { /// Where to read the input frame from. std::string inputFrame; + /// Name of the input tracker hit local collection. + std::string inputTrackerHitsLocal; /// Output measurement collection. std::string outputMeasurements; /// Output measurement to sim hit collection. std::string outputMeasurementSimHitsMap; /// Output cluster collection (optional). std::string outputClusters; + + /// DD4hep detector for cellID to geometry identifier resolution. + std::shared_ptr dd4hepDetector; }; /// Construct the cluster reader. @@ -59,6 +59,8 @@ class EDM4hepMeasurementInputConverter final : public PodioInputConverter { private: Config m_cfg; + EDM4hepUtil::MapGeometryIdFrom m_geometryMapper; + WriteDataHandle m_outputMeasurements{ this, "OutputMeasurements"}; diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementOutputConverter.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementOutputConverter.hpp index 164c0791cfe..fc5a66c3e14 100644 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementOutputConverter.hpp +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementOutputConverter.hpp @@ -8,12 +8,14 @@ #pragma once +#include "Acts/Geometry/TrackingGeometry.hpp" #include "ActsExamples/EventData/Cluster.hpp" #include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/Framework/DataHandle.hpp" #include "ActsExamples/Io/Podio/CollectionBaseWriteHandle.hpp" #include "ActsExamples/Io/Podio/PodioOutputConverter.hpp" +#include #include namespace ActsExamples { @@ -33,12 +35,11 @@ class EDM4hepMeasurementOutputConverter final : public PodioOutputConverter { struct Config { /// Which measurement collection to write. std::string inputMeasurements; - /// Which cluster collection to write (optional) - std::string inputClusters; - /// Name of the output tracker hit plane collection. - std::string outputTrackerHitsPlane = "ActsTrackerHitsPlane"; /// Name of the output tracker hit raw collection. - std::string outputTrackerHitsRaw = "ActsTrackerHitsRaw"; + std::string outputTrackerHitsLocal; + + /// Tracking geometry for surface lookup (local-to-global transform). + std::shared_ptr trackingGeometry; }; /// Constructor with @@ -64,12 +65,8 @@ class EDM4hepMeasurementOutputConverter final : public PodioOutputConverter { ReadDataHandle m_inputMeasurements{this, "InputMeasurements"}; - ReadDataHandle m_inputClusters{this, "InputClusters"}; - - CollectionBaseWriteHandle m_outputTrackerHitsPlane{this, - "OutputTrackerHitsPlane"}; - CollectionBaseWriteHandle m_outputTrackerHitsRaw{this, - "OutputTrackerHitsRaw"}; + CollectionBaseWriteHandle m_outputTrackerHitsLocal{this, + "OutputTrackerHitsLocal"}; }; } // namespace ActsExamples diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp index ed3217f947d..544b39014c8 100644 --- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp +++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp @@ -14,7 +14,8 @@ #include "ActsExamples/EventData/SimParticle.hpp" #include "ActsExamples/EventData/Trajectories.hpp" #include "ActsFatras/EventData/Hit.hpp" -#include "ActsPlugins/EDM4hep/TrackerHitCompatibility.hpp" +#include +#include #include @@ -95,21 +96,36 @@ VariableBoundMeasurementProxy readMeasurement( const edm4hep::TrackerHit3DCollection* fromClusters, Cluster* toCluster, const MapGeometryIdFrom& geometryMapper); -/// Writes a measurement cluster to EDM4hep. +/// Writes an ACTS measurement to EDM4hep. /// -/// Inpersistent information: -/// - hit index -/// - 1D local coords? -/// - segment path +/// Converts bound parameters (local coordinates, time) and covariance from an +/// ACTS measurement to a MutableTrackerHitLocal. The surface is used +/// to obtain the detector cell ID from the DD4hepDetectorElement and to +/// transform local coordinates to global position. /// -/// Known issues: -/// - cluster channels are written to inappropriate fields -/// - local 2D coordinates and time are written to position -void writeMeasurement(const ConstVariableBoundMeasurementProxy& from, - edm4hep::MutableTrackerHitPlane to, - const Cluster* fromCluster, - edm4hep::TrackerHit3DCollection& toClusters, - const MapGeometryIdTo& geometryMapper); +/// @param gctx The geometry context for coordinate transformations. +/// @param from The ACTS measurement to convert. +/// @param to The EDM4hep tracker hit to write to. +/// @param surface The surface associated with the measurement (must have a +/// DD4hepDetectorElement placement). +void writeMeasurement(const Acts::GeometryContext& gctx, + const ConstVariableBoundMeasurementProxy& from, + ActsPodioEdm::MutableTrackerHitLocal& to, + const Acts::Surface& surface); + +/// Reads a measurement from an EDM4hep TrackerHitLocal. +/// +/// Converts bound parameters, covariance, and indices from an EDM4hep +/// TrackerHitLocal into an ACTS measurement. The geometry mapper is used to +/// convert the cell ID to a geometry identifier. +/// +/// @param container The measurement container to insert into. +/// @param from The EDM4hep tracker hit to read from. +/// @param geometryMapper Function to map cell ID to geometry identifier. +/// @return Proxy to the created measurement. +VariableBoundMeasurementProxy readMeasurement( + MeasurementContainer& container, const ActsPodioEdm::TrackerHitLocal& from, + const MapGeometryIdFrom& geometryMapper); /// Writes a trajectory to EDM4hep. /// diff --git a/Examples/Io/EDM4hep/src/EDM4hepMeasurementInputConverter.cpp b/Examples/Io/EDM4hep/src/EDM4hepMeasurementInputConverter.cpp index 6de7228db0f..59d60093076 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepMeasurementInputConverter.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepMeasurementInputConverter.cpp @@ -8,15 +8,16 @@ #include "ActsExamples/Io/EDM4hep/EDM4hepMeasurementInputConverter.hpp" +#include "ActsExamples/DD4hepDetector/DD4hepDetector.hpp" #include "ActsExamples/EventData/Cluster.hpp" #include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp" -#include "ActsPlugins/EDM4hep/PodioUtil.hpp" +#include "ActsPlugins/DD4hep/DD4hepDetectorElement.hpp" +#include #include -#include -#include +#include #include namespace ActsExamples { @@ -30,34 +31,46 @@ EDM4hepMeasurementInputConverter::EDM4hepMeasurementInputConverter( if (m_cfg.outputMeasurements.empty()) { throw std::invalid_argument("Missing measurement output collection"); } + if (m_cfg.dd4hepDetector == nullptr) { + throw std::invalid_argument("Missing DD4hep detector"); + } m_outputMeasurements.initialize(m_cfg.outputMeasurements); m_outputMeasurementSimHitsMap.initialize(m_cfg.outputMeasurementSimHitsMap); m_outputClusters.maybeInitialize(m_cfg.outputClusters); + + m_geometryMapper = [detector = m_cfg.dd4hepDetector](std::uint64_t cellId) { + const auto& vm = detector->dd4hepDetector().volumeManager(); + const auto detElement = vm.lookupDetElement(cellId); + + const auto* ext = + detElement.extension( + false); + if (ext == nullptr) { + throw std::runtime_error( + "EDM4hepMeasurementInputConverter: DetElement has no " + "DD4hepDetectorElementExtension for cellId " + + std::to_string(cellId)); + } + return ext->detectorElement().surface().geometryId(); + }; } ProcessCode EDM4hepMeasurementInputConverter::convert( const AlgorithmContext& ctx, const podio::Frame& frame) const { MeasurementContainer measurements; ClusterContainer clusters; - // TODO what about those? IndexMultimap measurementSimHitsMap; - const auto& trackerHitPlaneCollection = - frame.get("ActsTrackerHitsPlane"); - const auto& trackerHitRawCollection = - frame.get("ActsTrackerHitsRaw"); - - for (const auto& trackerHitPlane : trackerHitPlaneCollection) { - Cluster cluster; - EDM4hepUtil::readMeasurement( - measurements, trackerHitPlane, &trackerHitRawCollection, &cluster, - [](std::uint64_t cellId) { return Acts::GeometryIdentifier(cellId); }); + const auto& trackerHitLocalCollection = + frame.get( + m_cfg.inputTrackerHitsLocal); - clusters.push_back(std::move(cluster)); + for (const auto& trackerHitLocal : trackerHitLocalCollection) { + EDM4hepUtil::readMeasurement(measurements, trackerHitLocal, + m_geometryMapper); } - // Write the data to the EventStore m_outputMeasurements(ctx, std::move(measurements)); m_outputMeasurementSimHitsMap(ctx, std::move(measurementSimHitsMap)); if (!m_cfg.outputClusters.empty()) { diff --git a/Examples/Io/EDM4hep/src/EDM4hepMeasurementOutputConverter.cpp b/Examples/Io/EDM4hep/src/EDM4hepMeasurementOutputConverter.cpp index 754a61027e9..39b7906f793 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepMeasurementOutputConverter.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepMeasurementOutputConverter.cpp @@ -11,6 +11,10 @@ #include "ActsExamples/EventData/Cluster.hpp" #include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/Io/EDM4hep/EDM4hepUtil.hpp" +#include "ActsPlugins/EDM4hep/TrackerHitCompatibility.hpp" +#include + +#include #include #include @@ -24,23 +28,18 @@ EDM4hepMeasurementOutputConverter::EDM4hepMeasurementOutputConverter( : PodioOutputConverter("EDM4hepMeasurementOutputConverter", std::move(logger)), m_cfg(config) { + if (m_cfg.trackingGeometry == nullptr) { + throw std::runtime_error( + "EDM4hepMeasurementOutputConverter: trackingGeometry is null"); + } + m_inputMeasurements.initialize(m_cfg.inputMeasurements); - m_inputClusters.maybeInitialize(m_cfg.inputClusters); - m_outputTrackerHitsPlane.initialize(m_cfg.outputTrackerHitsPlane); - m_outputTrackerHitsRaw.initialize(m_cfg.outputTrackerHitsRaw); + m_outputTrackerHitsLocal.initialize(m_cfg.outputTrackerHitsLocal); } ProcessCode EDM4hepMeasurementOutputConverter::execute( const AlgorithmContext& ctx) const { - ClusterContainer clusters; - - edm4hep::TrackerHitPlaneCollection hitsPlane; - edm4hep::TrackerHit3DCollection hits; - - if (!m_cfg.inputClusters.empty()) { - ACTS_VERBOSE("Fetch clusters for writing: " << m_cfg.inputClusters); - clusters = m_inputClusters(ctx); - } + ActsPodioEdm::TrackerHitLocalCollection hits; const auto measurements = m_inputMeasurements(ctx); @@ -50,23 +49,28 @@ ProcessCode EDM4hepMeasurementOutputConverter::execute( for (Index hitIdx = 0u; hitIdx < measurements.size(); ++hitIdx) { ConstVariableBoundMeasurementProxy from = measurements.getMeasurement(hitIdx); - const Cluster* fromCluster = clusters.empty() ? nullptr : &clusters[hitIdx]; - auto to = hitsPlane.create(); - EDM4hepUtil::writeMeasurement( - from, to, fromCluster, hits, - [](Acts::GeometryIdentifier id) { return id.value(); }); + const Acts::Surface* surface = + m_cfg.trackingGeometry->findSurface(from.geometryId()); + if (surface == nullptr) { + throw std::runtime_error( + "EDM4hepMeasurementOutputConverter: surface not found for geometry " + "id " + + std::to_string(from.geometryId().value())); + } + + auto to = hits.create(); + EDM4hepUtil::writeMeasurement(ctx.geoContext, from, to, *surface); } - m_outputTrackerHitsPlane(ctx, std::move(hitsPlane)); - m_outputTrackerHitsRaw(ctx, std::move(hits)); + m_outputTrackerHitsLocal(ctx, std::move(hits)); return ProcessCode::SUCCESS; } std::vector EDM4hepMeasurementOutputConverter::collections() const { - return {m_cfg.outputTrackerHitsPlane, m_cfg.outputTrackerHitsRaw}; + return {m_cfg.outputTrackerHitsLocal}; } } // namespace ActsExamples diff --git a/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp b/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp index acf2e701953..6b18c8bb943 100644 --- a/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp +++ b/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp @@ -16,7 +16,11 @@ #include "ActsExamples/EventData/Index.hpp" #include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/Validation/TrackClassification.hpp" +#include "ActsPlugins/DD4hep/DD4hepDetectorElement.hpp" #include "ActsPlugins/EDM4hep/EDM4hepUtil.hpp" +#include + +#include #include "edm4hep/TrackState.h" @@ -176,37 +180,49 @@ VariableBoundMeasurementProxy EDM4hepUtil::readMeasurement( } void EDM4hepUtil::writeMeasurement( + const Acts::GeometryContext& gctx, const ConstVariableBoundMeasurementProxy& from, - edm4hep::MutableTrackerHitPlane to, const Cluster* /*fromCluster*/, - edm4hep::TrackerHit3DCollection& /*toClusters*/, - const MapGeometryIdTo& geometryMapper) { - Acts::GeometryIdentifier geoId = from.geometryId(); + ActsPodioEdm::MutableTrackerHitLocal& to, const Acts::Surface& surface) { + long dim = from.size(); - if (geometryMapper) { - // no need for digitization as we only want to identify the sensor - to.setCellID(geometryMapper(geoId)); + const auto* placement = surface.surfacePlacement(); + if (placement == nullptr) { + throw std::runtime_error("Surface placement not found"); + } + const auto* dd4hepDetectorElement = + dynamic_cast(placement); + if (dd4hepDetectorElement == nullptr) { + throw std::runtime_error( + "Surface placement is not a DD4hepDetectorElement"); } - const auto& parameters = from.fullParameters(); - const auto& covariance = from.fullCovariance(); - - to.setTime(parameters[Acts::eBoundTime] / Acts::UnitConstants::ns); + std::uint64_t cellId = dd4hepDetectorElement->sourceElement().volumeID(); - to.setType(ActsPlugins::EDM4hepUtil::EDM4HEP_ACTS_POSITION_TYPE); - // TODO set uv (which are in global spherical coordinates with r=1) - to.setPosition({parameters[Acts::eBoundLoc0], parameters[Acts::eBoundLoc1], - parameters[Acts::eBoundTime]}); + ActsPlugins::EDM4hepUtil::writeMeasurement( + gctx, {from.parameters().data(), dim}, + {from.covariance().data(), dim, dim}, from.subspaceHelper().indices(), + cellId, surface, to); +} - to.setCovMatrix({ - static_cast(covariance(Acts::eBoundLoc0, Acts::eBoundLoc0)), - static_cast(covariance(Acts::eBoundLoc1, Acts::eBoundLoc0)), - static_cast(covariance(Acts::eBoundLoc1, Acts::eBoundLoc1)), - 0, - 0, - 0, - }); +VariableBoundMeasurementProxy EDM4hepUtil::readMeasurement( + MeasurementContainer& container, const ActsPodioEdm::TrackerHitLocal& from, + const MapGeometryIdFrom& geometryMapper) { + auto data = ActsPlugins::EDM4hepUtil::readMeasurement(from); + Acts::GeometryIdentifier geometryId = geometryMapper(data.cellId); + + const auto dim = data.indices.size(); + Acts::DynamicVector parameters(dim); + Acts::DynamicMatrix covariance(dim, dim); + for (std::size_t i = 0; i < dim; ++i) { + parameters(i) = data.parameters(data.indices[i]); + for (std::size_t j = 0; j < dim; ++j) { + covariance(i, j) = data.covariance(data.indices[i], data.indices[j]); + } + } - // @TODO: Check if we can write cell info + return container.emplaceMeasurement(static_cast(dim), + geometryId, data.indices, parameters, + covariance); } void EDM4hepUtil::writeTrajectory( @@ -256,8 +272,8 @@ void EDM4hepUtil::writeTrajectory( trackState.omega = converted.values[4]; trackState.time = converted.values[5]; - // Converted parameters are relative to an ad-hoc perigee surface created at - // the hit location + // Converted parameters are relative to an ad-hoc perigee surface created + // at the hit location auto center = converted.surface->center(gctx); trackState.referencePoint.x = center.x(); trackState.referencePoint.y = center.y(); diff --git a/Plugins/EDM4hep/edm.yml b/Plugins/EDM4hep/edm.yml index a9d1bd6bcdf..c3e0a9f7ee5 100644 --- a/Plugins/EDM4hep/edm.yml +++ b/Plugins/EDM4hep/edm.yml @@ -121,3 +121,22 @@ datatypes: std::size_t size() const { return indices_size(); } + + ActsPodioEdm::TrackerHitLocal: + Description: "Sensor-local tracker hit" + Author: "Paul Gessinger, CERN" + Members: + # These are chosen to be compatible with the `edm4hep::TrackerHit` interface + # (https://github.com/key4hep/EDM4hep/blob/c877c7fadf412db3c1126e7b410d8074a78014f2/edm4hep.yaml#L652-L666V) + # to allow for potential upstreaming + - uint64_t cellID // ID of the sensor that created this hit + - int32_t type // type of raw data hit, use depends on writer + - int32_t quality // quality bit flag of the hit + - float time [ns] // time of the hit (depends on measurement type) + - float eDep [GeV] // energy deposited on the hit (optional) + - float eDepError [GeV] // error measured on EDep (optional) + - edm4hep::Vector3d position [mm] // global hit position (depends on measurement type) + + VectorMembers: + - float measurement // Entries in the measurement value vector + - float covariance // Entries in the measurement covariance vector. Order is column major diff --git a/Plugins/EDM4hep/include/ActsPlugins/EDM4hep/EDM4hepUtil.hpp b/Plugins/EDM4hep/include/ActsPlugins/EDM4hep/EDM4hepUtil.hpp index 65e72afa678..814a628f0f7 100644 --- a/Plugins/EDM4hep/include/ActsPlugins/EDM4hep/EDM4hepUtil.hpp +++ b/Plugins/EDM4hep/include/ActsPlugins/EDM4hep/EDM4hepUtil.hpp @@ -15,14 +15,19 @@ #include "Acts/EventData/TrackParameters.hpp" #include "Acts/EventData/TrackProxyConcept.hpp" #include "Acts/EventData/TrackStatePropMask.hpp" +#include "Acts/EventData/Types.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/Logger.hpp" -#include +#include "Acts/Vertexing/Vertex.hpp" +#include "ActsPodioEdm/MutableTrackerHitLocal.h" +#include "ActsPodioEdm/TrackerHitLocal.h" +#include #include +#include #include #include #include @@ -117,12 +122,12 @@ void writeTrack(const Acts::GeometryContext& gctx, track_proxy_t track, auto setParameters = [](edm4hep::TrackState& trackState, const detail::Parameters& params) { - trackState.D0 = params.values[0]; - trackState.Z0 = params.values[1]; - trackState.phi = params.values[2]; - trackState.tanLambda = params.values[3]; - trackState.omega = params.values[4]; - trackState.time = params.values[5]; + trackState.D0 = static_cast(params.values[0]); + trackState.Z0 = static_cast(params.values[1]); + trackState.phi = static_cast(params.values[2]); + trackState.tanLambda = static_cast(params.values[3]); + trackState.omega = static_cast(params.values[4]); + trackState.time = static_cast(params.values[5]); if (params.covariance) { detail::packCovariance(params.covariance.value(), @@ -133,8 +138,7 @@ void writeTrack(const Acts::GeometryContext& gctx, track_proxy_t track, ACTS_VERBOSE("Converting " << track.nTrackStates() << " track states"); for (const auto& state : track.trackStatesReversed()) { - auto typeFlags = state.typeFlags(); - if (!typeFlags.isMeasurement()) { + if (!state.typeFlags().isMeasurement()) { continue; } @@ -160,9 +164,9 @@ void writeTrack(const Acts::GeometryContext& gctx, track_proxy_t track, // Converted parameters are relative to an ad-hoc perigee surface created at // the hit location auto center = converted.surface->center(gctx); - trackState.referencePoint.x = center.x(); - trackState.referencePoint.y = center.y(); - trackState.referencePoint.z = center.z(); + trackState.referencePoint.x = static_cast(center.x()); + trackState.referencePoint.y = static_cast(center.y()); + trackState.referencePoint.z = static_cast(center.z()); ACTS_VERBOSE("- ref surface ctr: " << center.transpose()); } outTrackStates.front().location = edm4hep::TrackState::AtLastHit; @@ -193,13 +197,13 @@ void writeTrack(const Acts::GeometryContext& gctx, track_proxy_t track, // track itself, but if that's not a perigee surface, another ad-hoc perigee // at the position will be created. auto center = converted.surface->center(gctx); - ipState.referencePoint.x = center.x(); - ipState.referencePoint.y = center.y(); - ipState.referencePoint.z = center.z(); + ipState.referencePoint.x = static_cast(center.x()); + ipState.referencePoint.y = static_cast(center.y()); + ipState.referencePoint.z = static_cast(center.z()); ACTS_VERBOSE("- ref surface ctr: " << center.transpose()); - for (auto& trackState : outTrackStates) { + for (const auto& trackState : outTrackStates) { to.addToTrackStates(trackState); } } @@ -218,8 +222,7 @@ void readTrack(const edm4hep::Track& from, track_proxy_t& track, double Bz, std::optional ipState; - auto unpack = - [](const edm4hep::TrackState& trackState) -> detail::Parameters { + auto unpack = [](const edm4hep::TrackState& trackState) { detail::Parameters params; params.covariance = BoundMatrix::Zero(); params.values = BoundVector::Zero(); @@ -336,4 +339,64 @@ constexpr bool kEdm4hepVertexHasTime = void writeVertex(const Acts::Vertex& vertex, edm4hep::MutableVertex to); +namespace detail { +// These functions are exposed here so they can be used from the unit tests +std::uint32_t encodeIndices(std::span indices); +boost::container::static_vector +decodeIndices(std::uint32_t type); +} // namespace detail + +/// Write a measurement to an EDM4hep tracker hit +/// +/// This function converts an ACTS measurement into the EDM4hep format. It +/// handles: +/// - Position conversion from local to global coordinates (in mm) +/// - Time storage (in ns) +/// - Measurement values and covariance matrix storage +/// - Encoding of measurement indices into a 32-bit integer: +/// - First 4 bits: number of indices (max 6) +/// - Next 4 bits per index: which parameter is being measured (0-6) +/// +/// The function will throw if: +/// - The number of indices exceeds 6 +/// - Any index is larger than 6 +/// - There's a size mismatch between parameters and covariance matrix +/// +/// @param gctx The geometry context +/// @param parameters The parameters of the measurement +/// @param covariance The covariance of the measurement +/// @param indices The indices of the measurement +/// @param cellId The cell ID of the measurement +/// @param surface The surface of the measurement +/// @param to The EDM4hep tracker hit to write to +void writeMeasurement(const Acts::GeometryContext& gctx, + const Eigen::Map& parameters, + const Eigen::Map& covariance, + std::span indices, + std::uint64_t cellId, const Acts::Surface& surface, + ActsPodioEdm::MutableTrackerHitLocal& to); + +/// Data extracted when reading a measurement from EDM4hep +struct MeasurementData { + /// Measurement parameters (local coordinates, full bound space) + Acts::BoundVector parameters{Acts::BoundVector::Zero()}; + /// Covariance matrix of the measurement (full bound space) + Acts::BoundMatrix covariance{Acts::BoundMatrix::Zero()}; + /// Indices of the measured parameters (subspace) + boost::container::static_vector + indices; + /// Cell ID of the measurement + std::uint64_t cellId{0}; +}; + +/// Read a measurement from an EDM4hep tracker hit +/// +/// This function extracts measurement parameters, covariance, and indices from +/// an EDM4hep TrackerHitLocal. It is the inverse of writeMeasurement. +/// +/// @param from The EDM4hep tracker hit to read from +/// @return The extracted measurement data (parameters, covariance, indices, +/// cellId) +MeasurementData readMeasurement(const ActsPodioEdm::TrackerHitLocal& from); + } // namespace ActsPlugins::EDM4hepUtil diff --git a/Plugins/EDM4hep/src/EDM4hepUtil.cpp b/Plugins/EDM4hep/src/EDM4hepUtil.cpp index 1cb7caa53e1..8ef05088e19 100644 --- a/Plugins/EDM4hep/src/EDM4hepUtil.cpp +++ b/Plugins/EDM4hep/src/EDM4hepUtil.cpp @@ -105,7 +105,7 @@ void packCovariance(const SquareMatrix<6>& from, float* to) { for (int i = 0; i < from.rows(); i++) { for (int j = 0; j <= i; j++) { std::size_t k = (i + 1) * i / 2 + j; - to[k] = from(i, j); + to[k] = static_cast(from(i, j)); } } } @@ -294,4 +294,148 @@ void writeVertex(const Vertex& vertex, edm4hep::MutableVertex to) { writeVertex(vertex, to); } +namespace detail { +/// Encode a list of bound parameter indices into a 32-bit integer. +/// +/// This function bit-packs up to 6 parameter indices (each 0-5 corresponding to +/// eBoundLoc0, eBoundLoc1, eBoundPhi, eBoundTheta, eBoundQOverP, eBoundTime) +/// into a single 32-bit unsigned integer for storage in EDM4hep format. +/// +/// Bit layout: +/// - Bits 0-3: Number of indices (size) +/// - Bits 4-7: First index +/// - Bits 8-11: Second index +/// - Bits 12-15: Third index +/// - (and so on, up to 6 indices total) +/// +/// @param indices Span of parameter indices to encode (max 6 elements) +/// @return Packed 32-bit unsigned integer containing all indices +std::uint32_t encodeIndices(std::span indices) { + if (indices.size() > eBoundSize) { + throw std::runtime_error( + "Number of indices exceeds maximum of 6 for EDM4hep"); + } + std::uint32_t result = 0; + + std::uint8_t shift = 0; + result |= (indices.size() << 0); + shift += 4; + + for (std::uint8_t index : indices) { + if (index > eBoundSize) { + throw std::runtime_error( + "Index out of range: can only encode indices up to 4 bits (0-15)"); + } + result |= (static_cast(index) << shift); + shift += 4; + } + return result; +} + +/// Decode a 32-bit integer back into a list of bound parameter indices. +/// +/// This function unpacks a bit-packed integer (created by encodeIndices) back +/// into the original list of parameter indices. See encodeIndices for the bit +/// layout specification. +/// +/// @param type Packed 32-bit unsigned integer containing encoded indices +/// @return Vector of decoded parameter indices (0-5 for each bound parameter) +boost::container::static_vector decodeIndices( + std::uint32_t type) { + boost::container::static_vector result; + std::uint8_t size = type & 0xF; + if (size > eBoundSize) { + throw std::runtime_error( + "Number of indices exceeds maximum of 6 for EDM4hep"); + } + result.resize(size); + for (std::size_t i = 0; i < result.size(); ++i) { + result[i] = (type >> ((i + 1) * 4)) & 0xF; + if (result[i] > eBoundSize) { + throw std::runtime_error( + "Index out of range: can only encode indices up to 4 bits (0-15)"); + } + } + return result; +} +} // namespace detail + +void writeMeasurement(const GeometryContext& gctx, + const Eigen::Map& parameters, + const Eigen::Map& covariance, + std::span indices, + std::uint64_t cellId, const Acts::Surface& surface, + ActsPodioEdm::MutableTrackerHitLocal& to) { + if (parameters.size() != covariance.rows() || + covariance.rows() != covariance.cols() || parameters.size() < 0 || + indices.size() != static_cast(parameters.size())) { + throw std::runtime_error( + "Size mismatch between parameters and covariance matrix"); + } + + auto dim = static_cast(parameters.size()); + + if (cellId != 0) { + to.setCellID(cellId); + } + + to.setType(detail::encodeIndices(indices)); + + auto loc0 = std::ranges::find(indices, eBoundLoc0); + auto loc1 = std::ranges::find(indices, eBoundLoc1); + auto time = std::ranges::find(indices, eBoundTime); + + if (loc0 != indices.end() && loc1 != indices.end()) { + Vector2 loc{parameters[std::distance(indices.begin(), loc0)], + parameters[std::distance(indices.begin(), loc1)]}; + Vector3 global = surface.localToGlobal(gctx, loc, Vector3::UnitZ()); + global /= Acts::UnitConstants::mm; + to.setPosition({global.x(), global.y(), global.z()}); + } + + if (time != indices.end()) { + std::size_t timeOffset = std::distance(indices.begin(), time); + to.setTime( + static_cast(parameters[timeOffset] / Acts::UnitConstants::ns)); + } + + for (double value : std::span{parameters.data(), dim}) { + to.addToMeasurement(static_cast(value)); + } + + for (double value : std::span{covariance.data(), dim * dim}) { + to.addToCovariance(static_cast(value)); + } +} + +MeasurementData readMeasurement(const ActsPodioEdm::TrackerHitLocal& from) { + auto indices = detail::decodeIndices(from.getType()); + auto meas = from.getMeasurement(); + auto cov = from.getCovariance(); + + const auto dim = indices.size(); + if (meas.size() != dim || cov.size() != dim * dim) { + throw std::runtime_error( + std::format("Size mismatch in EDM4hep tracker hit: measurement size " + "{}, covariance size {}, indices size {}", + meas.size(), cov.size(), dim)); + } + + MeasurementData result; + result.cellId = from.getCellID(); + result.indices = std::move(indices); + + for (std::size_t i = 0; i < dim; ++i) { + result.parameters(result.indices[i]) = meas[i]; + } + for (std::size_t i = 0; i < dim; ++i) { + for (std::size_t j = 0; j < dim; ++j) { + result.covariance(result.indices[i], result.indices[j]) = + cov[i * dim + j]; // row-major + } + } + + return result; +} + } // namespace ActsPlugins::EDM4hepUtil diff --git a/Python/Examples/src/plugins/EDM4hep.cpp b/Python/Examples/src/plugins/EDM4hep.cpp index 4d99df8e042..67a981f3680 100644 --- a/Python/Examples/src/plugins/EDM4hep.cpp +++ b/Python/Examples/src/plugins/EDM4hep.cpp @@ -100,19 +100,20 @@ PYBIND11_MODULE(ActsExamplesPythonBindingsEDM4hep, m) { } { - auto [alg, c] = + auto [alg, config] = declareAlgorithm( m, "EDM4hepMeasurementInputConverter"); - ACTS_PYTHON_STRUCT(c, inputFrame, outputMeasurements, - outputMeasurementSimHitsMap, outputClusters); + ACTS_PYTHON_STRUCT(config, inputFrame, inputTrackerHitsLocal, + outputMeasurements, outputMeasurementSimHitsMap, + outputClusters, dd4hepDetector); } { auto [alg, config] = declareAlgorithm( m, "EDM4hepMeasurementOutputConverter"); - ACTS_PYTHON_STRUCT(config, inputMeasurements, inputClusters, - outputTrackerHitsPlane, outputTrackerHitsRaw); + ACTS_PYTHON_STRUCT(config, inputMeasurements, outputTrackerHitsLocal, + trackingGeometry); } { diff --git a/Python/Examples/tests/test_edm4hep.py b/Python/Examples/tests/test_edm4hep.py index 92753cfd5ca..efd9705679a 100644 --- a/Python/Examples/tests/test_edm4hep.py +++ b/Python/Examples/tests/test_edm4hep.py @@ -8,6 +8,7 @@ import filelock from helpers import ( + dd4hepEnabled, edm4hepEnabled, podioEnabled, AssertCollectionExistsAlg, @@ -20,6 +21,7 @@ Sequencer, GenericDetector, ) +import acts.examples.json from acts.examples.odd import getOpenDataDetector, getOpenDataDetectorDirectory @@ -53,34 +55,69 @@ def assert_podio( @pytest.mark.edm4hep +@pytest.mark.odd @pytest.mark.skipif(not edm4hepEnabled, reason="EDM4hep is not set up") -def test_edm4hep_measurement_writer(tmp_path, fatras): +@pytest.mark.skipif(not dd4hepEnabled, reason="DD4hep not set up") +def test_edm4hep_measurement_writer(tmp_path, ptcl_gun, rng): from acts.examples.edm4hep import EDM4hepMeasurementOutputConverter, PodioWriter s = Sequencer(numThreads=1, events=10) - _, simAlg, digiAlg = fatras(s) + evGen, h3conv = ptcl_gun(s) - out = tmp_path / "measurements_edm4hep.root" + with getOpenDataDetector() as detector: + trackingGeometry = detector.trackingGeometry() + field = acts.ConstantBField(acts.Vector3(0, 0, 2 * u.T)) - converter = EDM4hepMeasurementOutputConverter( - level=acts.logging.VERBOSE, - inputMeasurements=digiAlg.config.outputMeasurements, - inputClusters=digiAlg.config.outputClusters, - outputTrackerHitsPlane="tracker_hits_plane", - outputTrackerHitsRaw="tracker_hits_raw", - ) - s.addAlgorithm(converter) + simAlg = acts.examples.FatrasSimulation( + level=acts.logging.INFO, + inputParticles=h3conv.config.outputParticles, + outputParticles="particles_simulated", + outputSimHits="simhits", + randomNumbers=rng, + trackingGeometry=trackingGeometry, + magneticField=field, + generateHitsOnSensitive=True, + emScattering=False, + emEnergyLossIonisation=False, + emEnergyLossRadiation=False, + emPhotonConversion=False, + ) + s.addAlgorithm(simAlg) + + digiCfg = acts.examples.DigitizationAlgorithm.Config( + digitizationConfigs=acts.examples.json.readDigiConfigFromJson( + str( + Path(__file__).parent.parent.parent.parent + / "Examples/Configs/odd-digi-smearing-config.json" + ) + ), + surfaceByIdentifier=trackingGeometry.geoIdSurfaceMap(), + randomNumbers=rng, + inputSimHits=simAlg.config.outputSimHits, + ) + digiAlg = acts.examples.DigitizationAlgorithm(digiCfg, acts.logging.INFO) + s.addAlgorithm(digiAlg) - s.addWriter( - PodioWriter( + out = tmp_path / "measurements_edm4hep.root" + + converter = EDM4hepMeasurementOutputConverter( level=acts.logging.VERBOSE, - outputPath=str(out), - category="events", - collections=converter.collections, + inputMeasurements=digiAlg.config.outputMeasurements, + outputTrackerHitsLocal="tracker_hits_local", + trackingGeometry=trackingGeometry, ) - ) + s.addAlgorithm(converter) - s.run() + s.addWriter( + PodioWriter( + level=acts.logging.VERBOSE, + outputPath=str(out), + category="events", + collections=converter.collections, + ) + ) + + s.run() assert os.path.isfile(out) assert os.stat(out).st_size > 10 @@ -88,7 +125,7 @@ def test_edm4hep_measurement_writer(tmp_path, fatras): assert_podio( out, "events", - collections=set(["tracker_hits_plane", "tracker_hits_raw"]), + collections=set(["tracker_hits_local"]), nevents=10, ) @@ -489,8 +526,10 @@ def test_edm4hep_simhit_particle_reader(tmp_path, ddsim_input): @pytest.mark.edm4hep +@pytest.mark.odd @pytest.mark.skipif(not edm4hepEnabled, reason="EDM4hep is not set up") -def test_edm4hep_measurement_reader(tmp_path, fatras): +@pytest.mark.skipif(not dd4hepEnabled, reason="DD4hep not set up") +def test_edm4hep_measurement_reader(tmp_path, ptcl_gun, rng): from acts.examples.edm4hep import ( EDM4hepMeasurementOutputConverter, EDM4hepMeasurementInputConverter, @@ -498,54 +537,91 @@ def test_edm4hep_measurement_reader(tmp_path, fatras): from acts.examples.edm4hep import PodioWriter, PodioReader s = Sequencer(numThreads=1, events=10) - _, simAlg, digiAlg = fatras(s) + evGen, h3conv = ptcl_gun(s) - out = tmp_path / "measurements_edm4hep.root" + with getOpenDataDetector() as detector: + trackingGeometry = detector.trackingGeometry() + field = acts.ConstantBField(acts.Vector3(0, 0, 2 * u.T)) - converter = EDM4hepMeasurementOutputConverter( - level=acts.logging.INFO, - inputMeasurements=digiAlg.config.outputMeasurements, - inputClusters=digiAlg.config.outputClusters, - ) - s.addAlgorithm(converter) - s.addWriter( - PodioWriter( + simAlg = acts.examples.FatrasSimulation( level=acts.logging.INFO, - outputPath=str(out), - category="events", - collections=converter.collections, + inputParticles=h3conv.config.outputParticles, + outputParticles="particles_simulated", + outputSimHits="simhits", + randomNumbers=rng, + trackingGeometry=trackingGeometry, + magneticField=field, + generateHitsOnSensitive=True, + emScattering=False, + emEnergyLossIonisation=False, + emEnergyLossRadiation=False, + emPhotonConversion=False, ) - ) - s.run() + s.addAlgorithm(simAlg) + + digiCfg = acts.examples.DigitizationAlgorithm.Config( + digitizationConfigs=acts.examples.json.readDigiConfigFromJson( + str( + Path(__file__).parent.parent.parent.parent + / "Examples/Configs/odd-digi-smearing-config.json" + ) + ), + surfaceByIdentifier=trackingGeometry.geoIdSurfaceMap(), + randomNumbers=rng, + inputSimHits=simAlg.config.outputSimHits, + ) + digiAlg = acts.examples.DigitizationAlgorithm(digiCfg, acts.logging.INFO) + s.addAlgorithm(digiAlg) - # read back in - s = Sequencer(numThreads=1) + out = tmp_path / "measurements_edm4hep.root" - s.addReader( - PodioReader( - level=acts.logging.WARNING, - inputPath=str(out), - outputFrame="events", - category="events", + converter = EDM4hepMeasurementOutputConverter( + level=acts.logging.INFO, + inputMeasurements=digiAlg.config.outputMeasurements, + outputTrackerHitsLocal="ActsTrackerHitsLocal", + trackingGeometry=trackingGeometry, ) - ) - s.addAlgorithm( - EDM4hepMeasurementInputConverter( - level=acts.logging.WARNING, - inputFrame="events", - outputMeasurements="measurements", - outputMeasurementSimHitsMap="simhitsmap", + s.addAlgorithm(converter) + s.addWriter( + PodioWriter( + level=acts.logging.INFO, + outputPath=str(out), + category="events", + collections=converter.collections, + ) ) - ) + s.run() - alg = AssertCollectionExistsAlg( - ["measurements", "simhitsmap"], "check_alg", acts.logging.WARNING - ) - s.addAlgorithm(alg) + # read back in + s = Sequencer(numThreads=1) - s.run() + s.addReader( + PodioReader( + level=acts.logging.WARNING, + inputPath=str(out), + outputFrame="events", + category="events", + ) + ) + s.addAlgorithm( + EDM4hepMeasurementInputConverter( + level=acts.logging.WARNING, + inputFrame="events", + inputTrackerHitsLocal="ActsTrackerHitsLocal", + outputMeasurements="measurements", + outputMeasurementSimHitsMap="simhitsmap", + dd4hepDetector=detector, + ) + ) - assert alg.events_seen == 10 + alg = AssertCollectionExistsAlg( + ["measurements", "simhitsmap"], "check_alg", acts.logging.WARNING + ) + s.addAlgorithm(alg) + + s.run() + + assert alg.events_seen == 10 @pytest.mark.edm4hep diff --git a/Tests/UnitTests/Plugins/EDM4hep/CMakeLists.txt b/Tests/UnitTests/Plugins/EDM4hep/CMakeLists.txt index 9ab6e9e91ba..62c27d4dc0a 100644 --- a/Tests/UnitTests/Plugins/EDM4hep/CMakeLists.txt +++ b/Tests/UnitTests/Plugins/EDM4hep/CMakeLists.txt @@ -8,9 +8,15 @@ target_link_libraries( set(unittest_extra_libraries Acts::PluginEDM4hep) add_unittest(PodioTrackContainer PodioTrackContainerTest.cpp) add_unittest(PodioTrackStateContainer PodioTrackStateContainerTest.cpp) - add_unittest(EDM4HepVertexWriteTest EDM4HepVertexWriteTests.cpp) + target_link_libraries( ActsUnitTestEDM4HepVertexWriteTest PUBLIC ActsPluginEDM4hep ) + +add_unittest(EDM4HepMeasurementTest EDM4HepMeasurementTests.cpp) +target_link_libraries( + ActsUnitTestEDM4HepMeasurementTest + PUBLIC Acts::PluginEDM4hep +) diff --git a/Tests/UnitTests/Plugins/EDM4hep/EDM4HepMeasurementTests.cpp b/Tests/UnitTests/Plugins/EDM4hep/EDM4HepMeasurementTests.cpp new file mode 100644 index 00000000000..ebbe7f14aef --- /dev/null +++ b/Tests/UnitTests/Plugins/EDM4hep/EDM4HepMeasurementTests.cpp @@ -0,0 +1,274 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#include +#include + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/TrackParametrization.hpp" +#include "Acts/EventData/SubspaceHelpers.hpp" +#include "Acts/EventData/detail/GenerateParameters.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Surfaces/DiscSurface.hpp" +#include "ActsPlugins/EDM4hep/EDM4hepUtil.hpp" +#include "ActsTests/CommonHelpers/FloatComparisons.hpp" +#include + +using namespace Acts; +using namespace ActsPlugins; +using namespace Acts::UnitLiterals; +using namespace Acts::detail::Test; + +namespace { +std::default_random_engine rng(123); +auto gctx = GeometryContext::dangerouslyDefaultConstruct(); +} // namespace + +BOOST_AUTO_TEST_SUITE(EDM4hepMeasurementTest) + +BOOST_AUTO_TEST_CASE(WriteMeasurement) { + auto [parameters, covariance] = + generateParametersCovariance(rng); + + std::uint64_t cellId = 1234; + + ActsPodioEdm::TrackerHitLocalCollection hits; + auto to = hits.create(); + + std::vector indices = {eBoundLoc0, eBoundLoc1}; + FixedSubspaceHelper helper(indices); + + Vector2 measPos = helper.projectVector(parameters); + SquareMatrix<2> measCov = helper.projectMatrix(covariance); + + auto surface = Surface::makeShared( + Transform3::Identity() * Translation3(Vector3{1, 2, 3}), 1_mm, 1_mm); + + Vector3 global = surface->localToGlobal(gctx, measPos, Vector3::UnitZ()); + + EDM4hepUtil::writeMeasurement(gctx, {measPos.data(), 2}, + {measCov.data(), 2, 2}, indices, cellId, + *surface, to); + + BOOST_CHECK_EQUAL(to.getCellID(), cellId); + BOOST_CHECK_EQUAL(to.getPosition()[0] * 1_mm, global.x()); + BOOST_CHECK_EQUAL(to.getPosition()[1] * 1_mm, global.y()); + BOOST_CHECK_EQUAL(to.getPosition()[2] * 1_mm, global.z()); + + // Time should be zero since we don't have time measurement + BOOST_CHECK_EQUAL(to.getTime() * 1_ns, 0.0); + + auto meas = to.getMeasurement(); + BOOST_CHECK_EQUAL(meas.size(), 2); + CHECK_CLOSE_REL(meas[0], measPos.x(), 1e-6); + CHECK_CLOSE_REL(meas[1], measPos.y(), 1e-6); + + auto cov = to.getCovariance(); + BOOST_CHECK_EQUAL(cov.size(), 4); + CHECK_CLOSE_REL(cov[0], measCov(0, 0), 1e-6); + CHECK_CLOSE_REL(cov[1], measCov(0, 1), 1e-6); + CHECK_CLOSE_REL(cov[2], measCov(1, 0), 1e-6); + CHECK_CLOSE_REL(cov[3], measCov(1, 1), 1e-6); + + auto unpackedIndices = EDM4hepUtil::detail::decodeIndices(to.getType()); + BOOST_CHECK_EQUAL(unpackedIndices.size(), 2); + BOOST_CHECK_EQUAL(unpackedIndices[0], eBoundLoc0); + BOOST_CHECK_EQUAL(unpackedIndices[1], eBoundLoc1); + + // Round-trip: read back and verify + auto read = EDM4hepUtil::readMeasurement(to); + BOOST_CHECK_EQUAL(read.cellId, cellId); + BOOST_CHECK_EQUAL(read.indices.size(), 2); + BOOST_CHECK_EQUAL(read.indices[0], eBoundLoc0); + BOOST_CHECK_EQUAL(read.indices[1], eBoundLoc1); + CHECK_CLOSE_REL(read.parameters(eBoundLoc0), measPos.x(), 1e-6); + CHECK_CLOSE_REL(read.parameters(eBoundLoc1), measPos.y(), 1e-6); + CHECK_CLOSE_REL(read.covariance(eBoundLoc0, eBoundLoc0), measCov(0, 0), 1e-6); + CHECK_CLOSE_REL(read.covariance(eBoundLoc0, eBoundLoc1), measCov(0, 1), 1e-6); + CHECK_CLOSE_REL(read.covariance(eBoundLoc1, eBoundLoc0), measCov(1, 0), 1e-6); + CHECK_CLOSE_REL(read.covariance(eBoundLoc1, eBoundLoc1), measCov(1, 1), 1e-6); +} + +BOOST_AUTO_TEST_CASE(WriteMeasurementNoPosition) { + auto [parameters, covariance] = + generateParametersCovariance(rng); + + std::uint64_t cellId = 1234; + + ActsPodioEdm::TrackerHitLocalCollection hits; + auto to = hits.create(); + + // Only measure phi and theta + std::vector indices = {eBoundPhi, eBoundTheta}; + FixedSubspaceHelper helper(indices); + + Vector2 measPos = helper.projectVector(parameters); + SquareMatrix<2> measCov = helper.projectMatrix(covariance); + + auto surface = Surface::makeShared( + Transform3::Identity() * Translation3(Vector3{1, 2, 3}), 1_mm, 1_mm); + + EDM4hepUtil::writeMeasurement(gctx, {measPos.data(), 2}, + {measCov.data(), 2, 2}, indices, cellId, + *surface, to); + + BOOST_CHECK_EQUAL(to.getCellID(), cellId); + // Position should be zero since we don't have loc0/loc1 + BOOST_CHECK_EQUAL(to.getPosition()[0] * 1_mm, 0.0); + BOOST_CHECK_EQUAL(to.getPosition()[1] * 1_mm, 0.0); + BOOST_CHECK_EQUAL(to.getPosition()[2] * 1_mm, 0.0); + + auto meas = to.getMeasurement(); + BOOST_CHECK_EQUAL(meas.size(), 2); + CHECK_CLOSE_REL(meas[0], measPos.x(), 1e-6); + CHECK_CLOSE_REL(meas[1], measPos.y(), 1e-6); + + auto cov = to.getCovariance(); + BOOST_CHECK_EQUAL(cov.size(), 4); + CHECK_CLOSE_REL(cov[0], measCov(0, 0), 1e-6); + CHECK_CLOSE_REL(cov[1], measCov(0, 1), 1e-6); + CHECK_CLOSE_REL(cov[2], measCov(1, 0), 1e-6); + CHECK_CLOSE_REL(cov[3], measCov(1, 1), 1e-6); + + auto unpackedIndices = EDM4hepUtil::detail::decodeIndices(to.getType()); + BOOST_CHECK_EQUAL(unpackedIndices.size(), 2); + BOOST_CHECK_EQUAL(unpackedIndices[0], eBoundPhi); + BOOST_CHECK_EQUAL(unpackedIndices[1], eBoundTheta); +} + +BOOST_AUTO_TEST_CASE(WriteMeasurementWithTime) { + auto [parameters, covariance] = + generateParametersCovariance(rng); + + std::uint64_t cellId = 1234; + + ActsPodioEdm::TrackerHitLocalCollection hits; + auto to = hits.create(); + + // Measure loc0, loc1, and time + std::vector indices = {eBoundLoc0, eBoundLoc1, eBoundTime}; + FixedSubspaceHelper helper(indices); + + Vector3 measPos = helper.projectVector(parameters); + SquareMatrix<3> measCov = helper.projectMatrix(covariance); + + auto surface = Surface::makeShared( + Transform3::Identity() * Translation3(Vector3{1, 2, 3}), 1_mm, 1_mm); + + Vector3 global = + surface->localToGlobal(gctx, measPos.head<2>(), Vector3::UnitZ()); + + EDM4hepUtil::writeMeasurement(gctx, {measPos.data(), 3}, + {measCov.data(), 3, 3}, indices, cellId, + *surface, to); + + BOOST_CHECK_EQUAL(to.getCellID(), cellId); + BOOST_CHECK_EQUAL(to.getPosition()[0] * 1_mm, global.x()); + BOOST_CHECK_EQUAL(to.getPosition()[1] * 1_mm, global.y()); + BOOST_CHECK_EQUAL(to.getPosition()[2] * 1_mm, global.z()); + // Time should be set since we have time measurement + CHECK_CLOSE_REL(to.getTime() * 1_ns, parameters[eBoundTime], 1e-6); + + auto meas = to.getMeasurement(); + BOOST_CHECK_EQUAL(meas.size(), 3); + CHECK_CLOSE_REL(meas[0], measPos.x(), 1e-6); + CHECK_CLOSE_REL(meas[1], measPos.y(), 1e-6); + CHECK_CLOSE_REL(meas[2], measPos.z(), 1e-6); + + auto cov = to.getCovariance(); + BOOST_CHECK_EQUAL(cov.size(), 9); + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + CHECK_CLOSE_REL(cov[i * 3 + j], measCov(i, j), 1e-6); + } + } + + auto unpackedIndices = EDM4hepUtil::detail::decodeIndices(to.getType()); + BOOST_CHECK_EQUAL(unpackedIndices.size(), 3); + BOOST_CHECK_EQUAL(unpackedIndices[0], eBoundLoc0); + BOOST_CHECK_EQUAL(unpackedIndices[1], eBoundLoc1); + BOOST_CHECK_EQUAL(unpackedIndices[2], eBoundTime); +} + +BOOST_AUTO_TEST_CASE(EncodeDecodeIndices) { + // Test empty span + { + std::vector indices = {}; + std::uint32_t encoded = EDM4hepUtil::detail::encodeIndices(indices); + auto decoded = EDM4hepUtil::detail::decodeIndices(encoded); + BOOST_CHECK_EQUAL(decoded.size(), 0); + } + + // Test single value + { + std::vector indices = {3}; + std::uint32_t encoded = EDM4hepUtil::detail::encodeIndices(indices); + auto decoded = EDM4hepUtil::detail::decodeIndices(encoded); + BOOST_CHECK_EQUAL(decoded.size(), 1); + BOOST_CHECK_EQUAL(decoded.at(0), 3); + } + + // Test maximum length (6) + { + std::vector indices = {0, 1, 2, 3, 4, 5}; + auto encoded = EDM4hepUtil::detail::encodeIndices(indices); + auto decoded = EDM4hepUtil::detail::decodeIndices(encoded); + BOOST_CHECK_EQUAL(decoded.size(), 6); + for (std::size_t i = 0; i < 6; ++i) { + BOOST_CHECK_EQUAL(decoded.at(i), i); + } + } + + // Test maximum value (6) + { + std::vector indices = {6, 6, 6}; + auto encoded = EDM4hepUtil::detail::encodeIndices(indices); + auto decoded = EDM4hepUtil::detail::decodeIndices(encoded); + BOOST_CHECK_EQUAL(decoded.size(), 3); + for (std::size_t i = 0; i < 3; ++i) { + BOOST_CHECK_EQUAL(decoded.at(i), 6); + } + } + + // Test mixed values + { + std::vector indices = {2, 5, 1, 4}; + auto encoded = EDM4hepUtil::detail::encodeIndices(indices); + auto decoded = EDM4hepUtil::detail::decodeIndices(encoded); + BOOST_CHECK_EQUAL(decoded.size(), 4); + BOOST_CHECK_EQUAL(decoded.at(0), 2); + BOOST_CHECK_EQUAL(decoded.at(1), 5); + BOOST_CHECK_EQUAL(decoded.at(2), 1); + BOOST_CHECK_EQUAL(decoded.at(3), 4); + } +} + +BOOST_AUTO_TEST_CASE(EncodeDecodeIndicesErrors) { + // Test exceeding maximum length (7 values) + { + std::vector indices = {0, 1, 2, 3, 4, 5, 6}; + BOOST_CHECK_THROW(EDM4hepUtil::detail::encodeIndices(indices), + std::runtime_error); + } + + // Test exceeding maximum value (7) + { + std::vector indices = {7}; + BOOST_CHECK_THROW(EDM4hepUtil::detail::encodeIndices(indices), + std::runtime_error); + } + + // Test mixed valid/invalid values + { + std::vector indices = {2, 7, 1}; + BOOST_CHECK_THROW(EDM4hepUtil::detail::encodeIndices(indices), + std::runtime_error); + } +} + +BOOST_AUTO_TEST_SUITE_END()