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/GeometryIdentifier.hpp>
0010 #include <Acts/Surfaces/Surface.hpp>
0011 #include <Acts/Utilities/Result.hpp>
0012 #include <DD4hep/Alignments.h>
0013 #include <DD4hep/DetElement.h>
0014 #include <DD4hep/Readout.h>
0015 #include <DD4hep/Segmentations.h>
0016 #include <DD4hep/VolumeManager.h>
0017 #include <DD4hep/detail/SegmentationsInterna.h>
0018 #include <DDRec/CellIDPositionConverter.h>
0019 #include <DDSegmentation/CartesianGridUV.h>
0020 #include <DDSegmentation/MultiSegmentation.h>
0021 #include <DDSegmentation/Segmentation.h>
0022 #include <Evaluator/DD4hepUnits.h>
0023 #include <Math/GenVector/Cartesian3D.h>
0024 #include <Math/GenVector/DisplacementVector3D.h>
0025 // Access "algorithms:GeoSvc"
0026 #include <algorithms/geo.h>
0027 #include <algorithms/logger.h>
0028 #include <edm4eic/Cov3f.h>
0029 #include <edm4eic/CovDiag3f.h>
0030 #include <edm4hep/Vector2f.h>
0031 #include <edm4hep/Vector3f.h>
0032 #include <Eigen/Core>
0033 #include <cmath>
0034 #include <exception>
0035 #include <stdexcept>
0036 #include <tuple>
0037 #include <unordered_map>
0038 #include <utility>
0039 #include <vector>
0040
0041 #include "ActsGeometryProvider.h"
0042
0043 namespace eicrecon {
0044
0045 void TrackerMeasurementFromHits::init() {
0046   // ***** B0Tracker
0047   m_detid_b0tracker = m_dd4hepGeo->constant<unsigned long>("B0Tracker_Station_1_ID");
0048   // ***** OuterMPGDBarrel
0049   std::string readout = "OuterMPGDBarrelHits";
0050   //  If "CartesianGridUV" segmentation, we need to rotate the cov. matrix.
0051   // Particularly when 2D-strip readout, with large uncertainty along strips.
0052   // i) Determine whether "CartesianGridUV", possibly embedded in a
0053   //   MultiSegmentation.
0054   // ii) If indeed:
0055   //   - set "m_detid_OuterMPGD" to single out corresponding hits,
0056   //   - AND set "m_outermpgd_UVsegmentation_mode", to be also required so as
0057   //    to avoid accidentally hitting a coincidence for malformed cellID with
0058   //    systemID = 0.
0059   const dd4hep::Detector* detector = algorithms::GeoSvc::instance().detector();
0060   dd4hep::Segmentation seg;
0061   m_outermpgd_UVsegmentation_mode = true;
0062   try {
0063     seg = detector->readout(readout).segmentation();
0064   } catch (const std::runtime_error&) {
0065     debug(R"(Failed to load Segmentation for "{}" readout.)", readout);
0066     m_outermpgd_UVsegmentation_mode = false;
0067   }
0068   if (m_outermpgd_UVsegmentation_mode) {
0069     debug(R"(Parsing Segmentation for "{}" readout: is it UV?)", readout);
0070     // Segmentation
0071     using Segmentation               = dd4hep::DDSegmentation::Segmentation;
0072     const Segmentation* segmentation = seg->segmentation;
0073     if (segmentation->type() == "MultiSegmentation") {
0074       // MultiSegmentation
0075       // - Parse all subSegmentations
0076       // - Require consistency: either all are UV, w/ same gridAngle, or none.
0077       //  We could invent a sophisticated scheme where cov. matrix would be
0078       //  rotated by a subSegmentation-dependent angle, possibly null. But for
0079       //  the time being, let's keep things simple: throw an exception upon any
0080       //  inconsistency.
0081       using MultiSegmentation    = dd4hep::DDSegmentation::MultiSegmentation;
0082       const auto* multiSeg       = dynamic_cast<const MultiSegmentation*>(segmentation);
0083       unsigned int subSegPattern = 0;
0084       for (const auto entry : multiSeg->subSegmentations()) {
0085         const Segmentation* subSegmentation = entry.segmentation;
0086         if (subSegmentation->type() == "CartesianGridUV") {
0087           const dd4hep::DDSegmentation::CartesianGridUV* gridUV =
0088               dynamic_cast<const dd4hep::DDSegmentation::CartesianGridUV*>(subSegmentation);
0089           double gridAngle = gridUV->gridAngle();
0090           if ((subSegPattern & 0x1) && fabs(gridAngle - m_gridAngle) > 1e-6) {
0091             critical(R"(Inconsistent Segmentation for "{}" readout: UV gridAngle not unique.)",
0092                      readout);
0093             throw std::runtime_error("Inconsistent MultiSegmentation");
0094           }
0095           subSegPattern |= 0x1;
0096           m_gridAngle = gridAngle;
0097         } else {
0098           subSegPattern |= 0x2;
0099         }
0100       }
0101       if (subSegPattern == 0x3) {
0102         critical(R"(Inconsistent Segmentation for "{}" readout: only partially UV.)", readout);
0103         throw std::runtime_error("Inconsistent MultiSegmentation");
0104       } else if (subSegPattern == 0x2) {
0105         m_outermpgd_UVsegmentation_mode = false;
0106       }
0107     } else {
0108       if (segmentation->type() != "CartesianGridUV") {
0109         m_outermpgd_UVsegmentation_mode = false;
0110       }
0111     }
0112     if (m_outermpgd_UVsegmentation_mode) {
0113       m_detid_OuterMPGD = m_dd4hepGeo->constant<unsigned long>("TrackerBarrel_5_ID");
0114     }
0115   }
0116 }
0117
0118 void TrackerMeasurementFromHits::process(const Input& input, const Output& output) const {
0119   const auto [trk_hits] = input;
0120   auto [meas2Ds]        = output;
0121
0122   constexpr double mm_acts = Acts::UnitConstants::mm;
0123   constexpr double mm_conv = mm_acts / dd4hep::mm; // = 1/0.1
0124
0125   // Get run-scoped geometry context from service
0126   const auto& gctx = m_acts_context->getActsGeometryContext();
0127
0128   // output collections
0129   auto const& surfaceMap = m_acts_context->surfaceMap();
0130
0131   // To do: add clustering to allow forming one measurement from several hits.
0132   // For now, one hit = one measurement.
0133   for (const auto& hit : *trk_hits) {
0134
0135     auto hit_det = hit.getCellID() & 0xFF;
0136
0137     Acts::SquareMatrix2 cov = Acts::SquareMatrix2::Zero();
0138     cov(0, 0)               = hit.getPositionError().xx * mm_acts * mm_acts; // note mm = 1 (Acts)
0139     cov(1, 1)               = hit.getPositionError().yy * mm_acts * mm_acts;
0140     cov(0, 1) = cov(1, 0) = 0.0;
0141     if (m_outermpgd_UVsegmentation_mode && hit_det == m_detid_OuterMPGD) {
0142       // Note: I(Y.B.) don't how to initialize an "Acts::RotationMatrix2" w/
0143       // an argument angle. The following turns out to fail:
0144       // const Acts::RotationMatrix2 rot2(m_gridAngle);
0145       const double sA = sin(m_gridAngle), cA = cos(m_gridAngle);
0146       const Acts::RotationMatrix2 rot2{{cA, sA}, {-sA, cA}};
0147       const Acts::RotationMatrix2 inv2{{cA, -sA}, {sA, cA}};
0148       cov = rot2 * cov * inv2;
0149     }
0150
0151     const auto* vol_ctx = m_converter->findContext(hit.getCellID());
0152     auto vol_id         = vol_ctx->identifier;
0153
0154     // trace("Hit preparation information: {}", hit_index);
0155     trace("   System id: {}, Cell id: {}", hit.getCellID() & 0xFF, hit.getCellID());
0156     trace("   cov matrix:      {:>12.2e} {:>12.2e}", cov(0, 0), cov(0, 1));
0157     trace("                    {:>12.2e} {:>12.2e}", cov(1, 0), cov(1, 1));
0158     trace("   surfaceMap size: {}", surfaceMap.size());
0159
0160     const auto is = surfaceMap.find(vol_id);
0161     if (is == surfaceMap.end()) {
0162       warning(" WARNING: vol_id ({})  not found in m_surfaces.", vol_id);
0163       continue;
0164     }
0165     const Acts::Surface* surface = is->second;
0166     // variable surf_center not used anywhere;
0167
0168     const auto& hit_pos = hit.getPosition(); // 3d position
0169
0170     Acts::Vector2 pos;
0171     auto onSurfaceTolerance =
0172         0.1 *
0173         Acts::UnitConstants::um; // By default, ACTS uses 0.1 micron as the on surface tolerance
0174     if (hit_det == m_detid_b0tracker) {
0175       onSurfaceTolerance =
0176           1 *
0177           Acts::UnitConstants::
0178               um; // FIXME Ugly hack for testing B0. Should be a way to increase this tolerance in geometry.
0179     }
0180
0181     try {
0182       // transform global position into local coordinates
0183       pos = surface
0184                 ->globalToLocal(gctx, {hit_pos.x, hit_pos.y, hit_pos.z}, {0, 0, 0},
0185                                 onSurfaceTolerance)
0186                 .value();
0187
0188     } catch (std::exception& ex) {
0189       warning("Can't convert globalToLocal for hit: vol_id={} det_id={} CellID={} x={} y={} z={}",
0190               vol_id, hit.getCellID() & 0xFF, hit.getCellID(), hit_pos.x, hit_pos.y, hit_pos.z);
0191       continue;
0192     }
0193
0194     if (level() <= algorithms::LogLevel::kTrace) {
0195
0196       Acts::Vector2 loc     = Acts::Vector2::Zero();
0197       loc[Acts::eBoundLoc0] = pos[0];
0198       loc[Acts::eBoundLoc1] = pos[1];
0199
0200       auto volman          = m_acts_context->dd4hepDetector()->volumeManager();
0201       auto alignment       = volman.lookupDetElement(vol_id).nominal();
0202       auto local_position  = (alignment.worldToLocal(
0203                                  {hit_pos.x / mm_conv, hit_pos.y / mm_conv, hit_pos.z / mm_conv})) *
0204                              mm_conv;
0205       double surf_center_x = surface->center(gctx).transpose()[0];
0206       double surf_center_y = surface->center(gctx).transpose()[1];
0207       double surf_center_z = surface->center(gctx).transpose()[2];
0208       trace("   hit position     : {:>10.2f} {:>10.2f} {:>10.2f}", hit_pos.x, hit_pos.y, hit_pos.z);
0209       trace("   local position   : {:>10.2f} {:>10.2f} {:>10.2f}", local_position.x(),
0210             local_position.y(), local_position.z());
0211       trace("   surface center   : {:>10.2f} {:>10.2f} {:>10.2f}", surf_center_x, surf_center_y,
0212             surf_center_z);
0213       trace("   acts local center: {:>10.2f} {:>10.2f}", pos.transpose()[0], pos.transpose()[1]);
0214       trace("   acts loc pos     : {:>10.2f} {:>10.2f}", loc[Acts::eBoundLoc0],
0215             loc[Acts::eBoundLoc1]);
0216     }
0217
0218     auto meas2D = meas2Ds->create();
0219     meas2D.setSurface(surface->geometryId().value()); // Surface for bound coordinates (geometryID)
0220     meas2D.setLoc(
0221         {static_cast<float>(pos[0]), static_cast<float>(pos[1])}); // 2D location on surface
0222     meas2D.setTime(hit.getTime());                                 // Measurement time
0223     // fixme: no off-diagonal terms. cov(0,1) = cov(1,0)??
0224     meas2D.setCovariance({cov(0, 0), cov(1, 1), hit.getTimeError() * hit.getTimeError(),
0225                           cov(0, 1)}); // Covariance on location and time
0226     meas2D.addToWeights(1.0);          // Weight for each of the hits, mirrors hits array
0227     meas2D.addToHits(hit);
0228   }
0229
0230   debug("All hits processed. Hits size: {}  measurements->size: {}", trk_hits->size(),
0231         meas2Ds->size());
0232 }
0233
0234 } // namespace eicrecon