src/components/TrackerSourceLinker.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Whitney Armstrong, Sylvester Joosten, Wouter Deconinck
0003
0004 #include "ActsExamples/EventData/GeometryContainers.hpp"
0005
0006 // Gaudi
0007 #include "Gaudi/Property.h"
0008 #include "Gaudi/Algorithm.h"
0009 #include "GaudiKernel/RndmGenerators.h"
0010 #include "GaudiKernel/ToolHandle.h"
0011
0012 #include <k4FWCore/DataHandle.h>
0013 #include <k4Interface/IGeoSvc.h>
0014 #include "JugTrack/IActsGeoSvc.h"
0015
0016 #include "DD4hep/DD4hepUnits.h"
0017 #include "DD4hep/Volumes.h"
0018 #include "DDRec/CellIDPositionConverter.h"
0019 #include "DDRec/Surface.h"
0020 #include "DDRec/SurfaceManager.h"
0021
0022 #include "Acts/Definitions/Common.hpp"
0023 #include "Acts/Definitions/Units.hpp"
0024 #if Acts_VERSION_MAJOR < 36
0025 #include "Acts/EventData/Measurement.hpp"
0026 #endif
0027 #include "Acts/EventData/MeasurementHelpers.hpp"
0028 #include "Acts/Geometry/TrackingGeometry.hpp"
0029 #include "Acts/Surfaces/Surface.hpp"
0030
0031 #include "ActsExamples/EventData/Index.hpp"
0032 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0033 #include "ActsExamples/EventData/Measurement.hpp"
0034
0035 #if Acts_VERSION_MAJOR >= 44
0036 #include "ActsPlugins/DD4hep/DD4hepDetectorElement.hpp"
0037 #else
0038 #include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
0039 #endif
0040
0041 #include "edm4eic/TrackerHitCollection.h"
0042
0043 namespace Jug::Reco {
0044
0045 /** Source source Linker.
0046  *
0047  * The source linker creates "source links" which map the hit to the tracking surface.
0048  * It also creates "measurements" which take the hit information and creates a corresponding
0049  * "measurement" which contains the covariance matrix and other geometry related hit information.
0050  *
0051  * \ingroup tracking
0052  */
0053 class TrackerSourceLinker : public Gaudi::Algorithm {
0054 private:
0055   mutable DataHandle<edm4eic::TrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
0056   mutable DataHandle<std::list<ActsExamples::IndexSourceLink>> m_sourceLinkStorage{"sourceLinkStorage", Gaudi::DataHandle::Writer, this};
0057 #if Acts_VERSION_MAJOR < 37 || (Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR == 0)
0058   mutable DataHandle<ActsExamples::IndexSourceLinkContainer> m_outputSourceLinks{"outputSourceLinks", Gaudi::DataHandle::Writer, this};
0059 #endif
0060   mutable DataHandle<ActsExamples::MeasurementContainer> m_outputMeasurements{"outputMeasurements", Gaudi::DataHandle::Writer, this};
0061   /// Pointer to the geometry services
0062   SmartIF<IGeoSvc> m_geoSvc;
0063   SmartIF<IActsGeoSvc> m_actsGeoSvc;
0064   std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> m_converter;
0065
0066 public:
0067   TrackerSourceLinker(const std::string& name, ISvcLocator* svcLoc) : Gaudi::Algorithm(name, svcLoc) {
0068     declareProperty("inputHitCollection", m_inputHitCollection, "");
0069     declareProperty("sourceLinkStorage", m_sourceLinkStorage, "");
0070 #if Acts_VERSION_MAJOR < 37 || (Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR == 0)
0071     declareProperty("outputSourceLinks", m_outputSourceLinks, "");
0072 #endif
0073     declareProperty("outputMeasurements", m_outputMeasurements, "");
0074   }
0075
0076   StatusCode initialize() override {
0077     if (Gaudi::Algorithm::initialize().isFailure()) {
0078       return StatusCode::FAILURE;
0079     }
0080     m_geoSvc = service("GeoSvc");
0081     if (!m_geoSvc) {
0082       error() << "Unable to locate Geometry Service. "
0083               << "Make sure you have GeoSvc in the right place in the configuration." << endmsg;
0084       return StatusCode::FAILURE;
0085     }
0086     m_actsGeoSvc = service("ActsGeoSvc");
0087     if (!m_actsGeoSvc) {
0088       error() << "Unable to locate ACTS Geometry Service. "
0089               << "Make sure you have ActsGeoSvc in the right place in the configuration." << endmsg;
0090       return StatusCode::FAILURE;
0091     }
0092     m_converter = std::make_shared<const dd4hep::rec::CellIDPositionConverter>(*(m_geoSvc->getDetector()));
0093
0094     return StatusCode::SUCCESS;
0095   }
0096
0097   StatusCode execute(const EventContext&) const override {
0098     constexpr double mm_acts = Acts::UnitConstants::mm;
0099     constexpr double mm_conv = mm_acts / dd4hep::mm; // = 1/0.1
0100
0101     // input collection
0102     const edm4eic::TrackerHitCollection* hits = m_inputHitCollection.get();
0103     // Create output collections
0104     auto* linkStorage  = m_sourceLinkStorage.createAndPut();
0105 #if Acts_VERSION_MAJOR < 37 || (Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR == 0)
0106     auto* sourceLinks  = m_outputSourceLinks.createAndPut();
0107     sourceLinks->reserve(hits->size());
0108 #endif
0109     auto* measurements = m_outputMeasurements.createAndPut();
0110     measurements->reserve(hits->size());
0111
0112     if (msgLevel(MSG::DEBUG)) {
0113       debug() << (*hits).size() << " hits " << endmsg;
0114     }
0115     int ihit = 0;
0116     for (const auto& ahit : *hits) {
0117
0118       Acts::SquareMatrix2 cov = Acts::SquareMatrix2::Zero();
0119       cov(0, 0)            = ahit.getPositionError().xx * mm_acts * mm_acts; // note mm = 1 (Acts)
0120       cov(1, 1)            = ahit.getPositionError().yy * mm_acts * mm_acts;
0121       if (msgLevel(MSG::DEBUG)) {
0122         debug() << "cov matrix:\n" << cov << endmsg;
0123       }
0124
0125       const auto* vol_ctx = m_converter->findContext(ahit.getCellID());
0126       auto vol_id = vol_ctx->identifier;
0127
0128       const auto is = m_actsGeoSvc->surfaceMap().find(vol_id);
0129       if (is == m_actsGeoSvc->surfaceMap().end()) {
0130         error() << " vol_id (" << vol_id << ")  not found in m_surfaces." << endmsg;
0131         continue;
0132       }
0133       const Acts::Surface* surface = is->second;
0134       // variable surf_center not used anywhere;
0135       // auto surf_center = surface->center(Acts::GeometryContext());
0136
0137       // transform global position into local coordinates
0138       // geometry context contains nothing here
0139       Acts::Vector2 pos =
0140           surface
0141               ->globalToLocal(Acts::GeometryContext(),
0142                               {ahit.getPosition().x, ahit.getPosition().y, ahit.getPosition().z}, {0, 0, 0})
0143               .value();
0144
0145       Acts::Vector2 loc     = Acts::Vector2::Zero();
0146       loc[Acts::eBoundLoc0] = pos[0];
0147       loc[Acts::eBoundLoc1] = pos[1];
0148
0149       if (msgLevel(MSG::DEBUG)) {
0150         auto volman         = m_geoSvc->getDetector()->volumeManager();
0151         auto alignment      = volman.lookupDetElement(vol_id).nominal();
0152         auto local_position = (alignment.worldToLocal({ahit.getPosition().x / mm_conv, ahit.getPosition().y / mm_conv,
0153                                                        ahit.getPosition().z / mm_conv})) *
0154                               mm_conv;
0155         debug() << " hit position     : " << ahit.getPosition().x << " " << ahit.getPosition().y << " "
0156                 << ahit.getPosition().z << endmsg;
0157         debug() << " dd4hep loc pos   : " << local_position.x() << " " << local_position.y() << " "
0158                 << local_position.z() << endmsg;
0159         debug() << " surface center   :" << surface->center(Acts::GeometryContext()).transpose() << endmsg;
0160         debug() << " acts local center:" << pos.transpose() << endmsg;
0161         debug() << " acts loc pos     : " << loc[Acts::eBoundLoc0] << ", " << loc[Acts::eBoundLoc1] << endmsg;
0162       }
0163
0164       // the measurement container is unordered and the index under which the
0165       // measurement will be stored is known before adding it.
0166       //
0167       // variable hitIdx not used anywhere
0168       // Index hitIdx = measurements->size();
0169
0170       auto geoId = surface->geometryId();
0171
0172       linkStorage->emplace_back(geoId, ihit);
0173 #if Acts_VERSION_MAJOR < 37 || (Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR == 0)
0174       ActsExamples::IndexSourceLink& sourceLink = linkStorage->back();
0175       sourceLinks->emplace_hint(sourceLinks->end(), sourceLink);
0176 #endif
0177
0178 #if Acts_VERSION_MAJOR > 37 || (Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR >= 1)
0179       std::array<Acts::BoundIndices, 2> indices = {Acts::eBoundLoc0, Acts::eBoundLoc1};
0180       Acts::visit_measurement(
0181         indices.size(), [&](auto dim) -> ActsExamples::VariableBoundMeasurementProxy {
0182           if constexpr (dim == indices.size()) {
0183             return ActsExamples::VariableBoundMeasurementProxy{
0184               measurements->emplaceMeasurement<dim>(geoId, indices, loc, cov)
0185             };
0186           } else {
0187             throw std::runtime_error("Dimension not supported in measurement creation");
0188           }
0189         }
0190       );
0191 #elif Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR == 0
0192       std::array<Acts::BoundIndices, 2> indices = {Acts::eBoundLoc0, Acts::eBoundLoc1};
0193       Acts::visit_measurement(
0194         indices.size(), [&](auto dim) -> ActsExamples::VariableBoundMeasurementProxy {
0195           if constexpr (dim == indices.size()) {
0196             return ActsExamples::VariableBoundMeasurementProxy{
0197               measurements->emplaceMeasurement<dim>(Acts::SourceLink{sourceLink}, indices, loc, cov)
0198             };
0199           } else {
0200             throw std::runtime_error("Dimension not supported in measurement creation");
0201           }
0202         }
0203       );
0204 #elif Acts_VERSION_MAJOR == 36 && Acts_VERSION_MINOR >= 1
0205       auto measurement = ActsExamples::makeVariableSizeMeasurement(
0206         Acts::SourceLink{sourceLink}, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1);
0207       measurements->emplace_back(std::move(measurement));
0208 #elif Acts_VERSION_MAJOR == 36 && Acts_VERSION_MINOR == 0
0209       auto measurement = ActsExamples::makeFixedSizeMeasurement(
0210         Acts::SourceLink{sourceLink}, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1);
0211       measurements->emplace_back(std::move(measurement));
0212 #else
0213       auto measurement = Acts::makeMeasurement(
0214         Acts::SourceLink{sourceLink}, loc, cov, Acts::eBoundLoc0, Acts::eBoundLoc1);
0215       measurements->emplace_back(std::move(measurement));
0216 #endif
0217
0218       ihit++;
0219     }
0220     return StatusCode::SUCCESS;
0221   }
0222 };
0223 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0224 DECLARE_COMPONENT(TrackerSourceLinker)
0225
0226 } // namespace Jug::Reco