algorithms/tracking/TrackerMeasurementFromHits.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 - 2025 Shujie Li, Wouter Deconinck
0003
0004 #include "TrackerMeasurementFromHits.h"
0005
0006 #include <Acts/Definitions/Algebra.hpp>
0007 #include <Acts/Definitions/TrackParametrization.hpp>
0008 #include <Acts/Definitions/Units.hpp>
0009 #include <Acts/Geometry/GeometryContext.hpp>
0010 #include <Acts/Geometry/GeometryIdentifier.hpp>
0011 #include <Acts/Surfaces/Surface.hpp>
0012 #include <Acts/Utilities/Result.hpp>
0013 #include <DD4hep/Alignments.h>
0014 #include <DD4hep/DetElement.h>
0015 #include <DD4hep/VolumeManager.h>
0016 #include <DDRec/CellIDPositionConverter.h>
0017 #include <Evaluator/DD4hepUnits.h>
0018 #include <Math/GenVector/Cartesian3D.h>
0019 #include <Math/GenVector/DisplacementVector3D.h>
0020 #include <algorithms/logger.h>
0021 #include <edm4eic/Cov3f.h>
0022 #include <edm4eic/CovDiag3f.h>
0023 #include <edm4hep/Vector2f.h>
0024 #include <edm4hep/Vector3f.h>
0025 #include <fmt/core.h>
0026 #include <Eigen/Core>
0027 #include <exception>
0028 #include <unordered_map>
0029 #include <utility>
0030
0031 #include "ActsGeometryProvider.h"
0032
0033 namespace eicrecon {
0034
0035 void TrackerMeasurementFromHits::init() {
0036   m_detid_b0tracker = m_dd4hepGeo->constant<unsigned long>("B0Tracker_Station_1_ID");
0037 }
0038
0039 void TrackerMeasurementFromHits::process(const Input& input, const Output& output) const {
0040   const auto [trk_hits] = input;
0041   auto [meas2Ds]        = output;
0042
0043   constexpr double mm_acts = Acts::UnitConstants::mm;
0044   constexpr double mm_conv = mm_acts / dd4hep::mm; // = 1/0.1
0045
0046   // output collections
0047   auto const& surfaceMap = m_acts_context->surfaceMap();
0048
0049   // To do: add clustering to allow forming one measurement from several hits.
0050   // For now, one hit = one measurement.
0051   for (const auto& hit : *trk_hits) {
0052
0053     Acts::SquareMatrix2 cov = Acts::SquareMatrix2::Zero();
0054     cov(0, 0)               = hit.getPositionError().xx * mm_acts * mm_acts; // note mm = 1 (Acts)
0055     cov(1, 1)               = hit.getPositionError().yy * mm_acts * mm_acts;
0056     cov(0, 1) = cov(1, 0) = 0.0;
0057
0058     const auto* vol_ctx = m_converter->findContext(hit.getCellID());
0059     auto vol_id         = vol_ctx->identifier;
0060
0061     // trace("Hit preparation information: {}", hit_index);
0062     trace("   System id: {}, Cell id: {}", hit.getCellID() & 0xFF, hit.getCellID());
0063     trace("   cov matrix:      {:>12.2e} {:>12.2e}", cov(0, 0), cov(0, 1));
0064     trace("                    {:>12.2e} {:>12.2e}", cov(1, 0), cov(1, 1));
0065     trace("   surfaceMap size: {}", surfaceMap.size());
0066
0067     const auto is = surfaceMap.find(vol_id);
0068     if (is == surfaceMap.end()) {
0069       warning(" WARNING: vol_id ({})  not found in m_surfaces.", vol_id);
0070       continue;
0071     }
0072     const Acts::Surface* surface = is->second;
0073     // variable surf_center not used anywhere;
0074
0075     const auto& hit_pos = hit.getPosition(); // 3d position
0076
0077     Acts::Vector2 pos;
0078     auto hit_det = hit.getCellID() & 0xFF;
0079     auto onSurfaceTolerance =
0080         0.1 *
0081         Acts::UnitConstants::um; // By default, ACTS uses 0.1 micron as the on surface tolerance
0082     if (hit_det == m_detid_b0tracker) {
0083       onSurfaceTolerance =
0084           1 *
0085           Acts::UnitConstants::
0086               um; // FIXME Ugly hack for testing B0. Should be a way to increase this tolerance in geometry.
0087     }
0088
0089     try {
0090       // transform global position into local coordinates
0091       // geometry context contains nothing here
0092       pos = surface
0093                 ->globalToLocal(Acts::GeometryContext(), {hit_pos.x, hit_pos.y, hit_pos.z},
0094                                 {0, 0, 0}, onSurfaceTolerance)
0095                 .value();
0096
0097     } catch (std::exception& ex) {
0098       warning("Can't convert globalToLocal for hit: vol_id={} det_id={} CellID={} x={} y={} z={}",
0099               vol_id, hit.getCellID() & 0xFF, hit.getCellID(), hit_pos.x, hit_pos.y, hit_pos.z);
0100       continue;
0101     }
0102
0103     if (level() <= algorithms::LogLevel::kTrace) {
0104
0105       Acts::Vector2 loc     = Acts::Vector2::Zero();
0106       loc[Acts::eBoundLoc0] = pos[0];
0107       loc[Acts::eBoundLoc1] = pos[1];
0108
0109       auto volman         = m_acts_context->dd4hepDetector()->volumeManager();
0110       auto alignment      = volman.lookupDetElement(vol_id).nominal();
0111       auto local_position = (alignment.worldToLocal(
0112                                 {hit_pos.x / mm_conv, hit_pos.y / mm_conv, hit_pos.z / mm_conv})) *
0113                             mm_conv;
0114       double surf_center_x = surface->center(Acts::GeometryContext()).transpose()[0];
0115       double surf_center_y = surface->center(Acts::GeometryContext()).transpose()[1];
0116       double surf_center_z = surface->center(Acts::GeometryContext()).transpose()[2];
0117       trace("   hit position     : {:>10.2f} {:>10.2f} {:>10.2f}", hit_pos.x, hit_pos.y, hit_pos.z);
0118       trace("   local position   : {:>10.2f} {:>10.2f} {:>10.2f}", local_position.x(),
0119             local_position.y(), local_position.z());
0120       trace("   surface center   : {:>10.2f} {:>10.2f} {:>10.2f}", surf_center_x, surf_center_y,
0121             surf_center_z);
0122       trace("   acts local center: {:>10.2f} {:>10.2f}", pos.transpose()[0], pos.transpose()[1]);
0123       trace("   acts loc pos     : {:>10.2f} {:>10.2f}", loc[Acts::eBoundLoc0],
0124             loc[Acts::eBoundLoc1]);
0125     }
0126
0127     auto meas2D = meas2Ds->create();
0128     meas2D.setSurface(surface->geometryId().value()); // Surface for bound coordinates (geometryID)
0129     meas2D.setLoc(
0130         {static_cast<float>(pos[0]), static_cast<float>(pos[1])}); // 2D location on surface
0131     meas2D.setTime(hit.getTime());                                 // Measurement time
0132     // fixme: no off-diagonal terms. cov(0,1) = cov(1,0)??
0133     meas2D.setCovariance({cov(0, 0), cov(1, 1), hit.getTimeError() * hit.getTimeError(),
0134                           cov(0, 1)}); // Covariance on location and time
0135     meas2D.addToWeights(1.0);          // Weight for each of the hits, mirrors hits array
0136     meas2D.addToHits(hit);
0137   }
0138
0139   debug("All hits processed. Hits size: {}  measurements->size: {}", trk_hits->size(),
0140         meas2Ds->size());
0141 }
0142
0143 } // namespace eicrecon