algorithms/tracking/TrackerMeasurementFromHits.cc

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