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