src/components/SingleTrackSourceLinker.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Whitney Armstrong, Sylvester Joosten, Wouter Deconinck
0003
0004 // Gaudi
0005 #include "Gaudi/Property.h"
0006 #include "Gaudi/Algorithm.h"
0007 #include "GaudiKernel/RndmGenerators.h"
0008 #include "GaudiKernel/ToolHandle.h"
0009
0010 #include <k4FWCore/DataHandle.h>
0011 #include <k4Interface/IGeoSvc.h>
0012 #include "JugTrack/IActsGeoSvc.h"
0013
0014 #include "DD4hep/DD4hepUnits.h"
0015 #include "DD4hep/Volumes.h"
0016 #include "DDRec/CellIDPositionConverter.h"
0017 #include "DDRec/Surface.h"
0018 #include "DDRec/SurfaceManager.h"
0019
0020 #include "Acts/Definitions/Common.hpp"
0021 #include "Acts/Definitions/Units.hpp"
0022 #if Acts_VERSION_MAJOR < 36
0023 #include "Acts/EventData/Measurement.hpp"
0024 #endif
0025 #include "Acts/EventData/MeasurementHelpers.hpp"
0026 #include "Acts/Geometry/TrackingGeometry.hpp"
0027 #include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
0028 #include "Acts/Surfaces/Surface.hpp"
0029
0030 #include "ActsExamples/EventData/GeometryContainers.hpp"
0031 #include "ActsExamples/EventData/Index.hpp"
0032 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0033 #include "ActsExamples/EventData/Measurement.hpp"
0034 #include "ActsExamples/EventData/ProtoTrack.hpp"
0035
0036 #include "edm4eic/TrackerHitCollection.h"
0037
0038 namespace Jug::Reco {
0039
0040 /** Single Track source Linker and proto tracks.
0041  *
0042  * This algorithm assumes only single track events.
0043  *
0044  * \ingroup tracking
0045  */
0046 class SingleTrackSourceLinker : public Gaudi::Algorithm {
0047 private:
0048   mutable DataHandle<edm4eic::TrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
0049   mutable DataHandle<std::list<ActsExamples::IndexSourceLink>> m_sourceLinkStorage{"sourceLinkStorage", Gaudi::DataHandle::Writer, this};
0050 #if Acts_VERSION_MAJOR < 37 || (Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR == 0)
0051   mutable DataHandle<ActsExamples::IndexSourceLinkContainer> m_outputSourceLinks{"outputSourceLinks", Gaudi::DataHandle::Writer, this};
0052 #endif
0053   mutable DataHandle<ActsExamples::MeasurementContainer> m_outputMeasurements{"outputMeasurements", Gaudi::DataHandle::Writer, this};
0054   mutable DataHandle<ActsExamples::ProtoTrackContainer> m_outputProtoTracks{"outputProtoTracks", Gaudi::DataHandle::Writer, this};
0055   /// Pointer to the geometry service
0056   SmartIF<IGeoSvc> m_geoSvc;
0057   SmartIF<IActsGeoSvc> m_actsGeoSvc;
0058   std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> m_converter;
0059
0060 public:
0061   SingleTrackSourceLinker(const std::string& name, ISvcLocator* svcLoc) : Gaudi::Algorithm(name, svcLoc) {
0062     declareProperty("inputHitCollection", m_inputHitCollection, "");
0063     declareProperty("sourceLinkStorage", m_sourceLinkStorage, "");
0064 #if Acts_VERSION_MAJOR < 37 || (Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR == 0)
0065     declareProperty("outputSourceLinks", m_outputSourceLinks, "");
0066 #endif
0067     declareProperty("outputMeasurements", m_outputMeasurements, "");
0068     declareProperty("outputProtoTracks", m_outputProtoTracks, "");
0069   }
0070
0071   StatusCode initialize() override {
0072     if (Gaudi::Algorithm::initialize().isFailure()) {
0073       return StatusCode::FAILURE;
0074     }
0075     m_geoSvc = service("GeoSvc");
0076     if (!m_geoSvc) {
0077       error() << "Unable to locate Geometry Service. "
0078               << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
0079       return StatusCode::FAILURE;
0080     }
0081     m_actsGeoSvc = service("ActsGeoSvc");
0082     if (!m_actsGeoSvc) {
0083       error() << "Unable to locate ACTS Geometry Service. "
0084               << "Make sure you have ActsGeoSvc and SimSvc in the right order in the configuration." << endmsg;
0085       return StatusCode::FAILURE;
0086     }
0087     m_converter = std::make_shared<const dd4hep::rec::CellIDPositionConverter>(*(m_geoSvc->getDetector()));
0088     return StatusCode::SUCCESS;
0089   }
0090
0091   StatusCode execute(const EventContext&) const override {
0092     // input collection
0093     const edm4eic::TrackerHitCollection* hits = m_inputHitCollection.get();
0094     // Create output collections
0095     auto* linkStorage  = m_sourceLinkStorage.createAndPut();
0096 #if Acts_VERSION_MAJOR < 37 || (Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR == 0)
0097     auto* sourceLinks  = m_outputSourceLinks.createAndPut();
0098     sourceLinks->reserve(hits->size());
0099 #endif
0100     auto* measurements = m_outputMeasurements.createAndPut();
0101     measurements->reserve(hits->size());
0102     auto* protoTracks  = m_outputProtoTracks.createAndPut();
0103     // IndexMultimap<ActsFatras::Barcode> hitParticlesMap;
0104     // IndexMultimap<Index> hitSimHitsMap;
0105
0106     // assume single track --> one ProtoTrack
0107     ActsExamples::ProtoTrack track;
0108     track.reserve((*hits).size());
0109
0110     if (msgLevel(MSG::DEBUG)) {
0111       debug() << (*hits).size() << " hits " << endmsg;
0112     }
0113     int ihit = 0;
0114     for (const auto& ahit : *hits) {
0115
0116       track.emplace_back(ihit);
0117
0118       const auto* vol_ctx = m_converter->findContext(ahit.getCellID());
0119       auto vol_id   = vol_ctx->identifier;
0120       const auto is = m_actsGeoSvc->surfaceMap().find(vol_id);
0121       if (is == m_actsGeoSvc->surfaceMap().end()) {
0122         if (msgLevel(MSG::DEBUG)) {
0123           debug() << " vol_id (" << vol_id << ")  not found in m_surfaces!!!!" << endmsg;
0124         }
0125         continue;
0126       }
0127       const Acts::Surface* surface = is->second;
0128
0129       // NOTE
0130       // Here is where the important hit and tracking geometry is connected.
0131       // A "Measurement" is constructed to for each hit which makes the connection to
0132       // the tracking surface and covariance matrix
0133
0134       auto acts_pos = surface
0135                           ->globalToLocal(Acts::GeometryContext(),
0136                                           {ahit.getPosition().x, ahit.getPosition().y, ahit.getPosition().z}, {0, 0, 0})
0137                           .value();
0138
0139       if (msgLevel(MSG::DEBUG)) {
0140         auto volman  = m_geoSvc->getDetector()->volumeManager();
0141         auto detelem = volman.lookupDetElement(vol_id);
0142         auto local_pos =
0143             detelem.nominal().worldToLocal({ahit.getPosition().x, ahit.getPosition().y, ahit.getPosition().z});
0144         debug() << "===== Debugging hit =====" << endmsg;
0145         debug() << "DD4hep global pos (" << ahit.getPosition().x << "," << ahit.getPosition().y << ","
0146                 << ahit.getPosition().z << ")" << endmsg;
0147         debug() << "DD4hep local  pos (" << local_pos.x() << "," << local_pos.y() << "," << local_pos.z() << ")"
0148                 << endmsg;
0149         debug() << "ACTS local position : (" << acts_pos[0] << "," << acts_pos[1] << ")" << endmsg;
0150         debug() << "ACTS surface center : " << surface->center(Acts::GeometryContext()).transpose() << endmsg;
0151         debug() << "DD4hep DetElement center : "
0152                 << detelem.nominal().localToWorld(detelem.placement().position()) / dd4hep::mm << endmsg;
0153       }
0154       // construct the vector of measured parameters (2d getPosition in this case)
0155       Acts::Vector2 pos(acts_pos.x(), acts_pos.y());
0156
0157       // construct the covariance matrix
0158       Acts::SquareMatrix2 cov = Acts::SquareMatrix2::Zero();
0159       cov(0, 0)            = ahit.getPositionError().xx * Acts::UnitConstants::mm * Acts::UnitConstants::mm;
0160       cov(1, 1)            = ahit.getPositionError().yy * Acts::UnitConstants::mm * Acts::UnitConstants::mm;
0161
0162       // Above we only consider the two position coordinates the comment below shows how to add time
0163       // which we will probably want to try later.
0164       //
0165       // Acts::SquareMatrix3 cov;
0166       // cov << 0.05, 0., 0., 0., 0.05, 0., 0., 0., 900. * Acts::UnitConstants::ps * Acts::UnitConstants::ps;
0167       // Acts::Vector3 par(localX, localY, simHit.time());
0168
0169       auto geoId = surface->geometryId();
0170
0171       linkStorage->emplace_back(geoId, ihit);
0172 #if Acts_VERSION_MAJOR < 37 || (Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR == 0)
0173       ActsExamples::IndexSourceLink& sourceLink = linkStorage->back();
0174       sourceLinks->emplace_hint(sourceLinks->end(), sourceLink);
0175 #endif
0176
0177 #if Acts_VERSION_MAJOR > 37 || (Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR >= 1)
0178       std::array<Acts::BoundIndices, 2> indices = {Acts::eBoundLoc0, Acts::eBoundLoc1};
0179       Acts::visit_measurement(
0180         indices.size(), [&](auto dim) -> ActsExamples::VariableBoundMeasurementProxy {
0181           if constexpr (dim == indices.size()) {
0182             return ActsExamples::VariableBoundMeasurementProxy{
0183               measurements->emplaceMeasurement<dim>(geoId, indices, pos, cov)
0184             };
0185           } else {
0186             throw std::runtime_error("Dimension not supported in measurement creation");
0187           }
0188         }
0189       );
0190 #elif Acts_VERSION_MAJOR == 37 && Acts_VERSION_MINOR == 0
0191       std::array<Acts::BoundIndices, 2> indices = {Acts::eBoundLoc0, Acts::eBoundLoc1};
0192       Acts::visit_measurement(
0193         indices.size(), [&](auto dim) -> ActsExamples::VariableBoundMeasurementProxy {
0194           if constexpr (dim == indices.size()) {
0195             return ActsExamples::VariableBoundMeasurementProxy{
0196               measurements->emplaceMeasurement<dim>(Acts::SourceLink{sourceLink}, indices, pos, cov)
0197             };
0198           } else {
0199             throw std::runtime_error("Dimension not supported in measurement creation");
0200           }
0201         }
0202       );
0203 #elif Acts_VERSION_MAJOR == 36 && Acts_VERSION_MINOR >= 1
0204       auto measurement = ActsExamples::makeVariableSizeMeasurement(
0205           Acts::SourceLink{sourceLink}, pos, cov, Acts::eBoundLoc0, Acts::eBoundLoc1);
0206       measurements->emplace_back(std::move(measurement));
0207 #elif Acts_VERSION_MAJOR == 36 && Acts_VERSION_MINOR == 0
0208       auto measurement = ActsExamples::makeFixedSizeMeasurement(
0209           Acts::SourceLink{sourceLink}, pos, cov, Acts::eBoundLoc0, Acts::eBoundLoc1);
0210       measurements->emplace_back(std::move(measurement));
0211 #else
0212       auto measurement =
0213           Acts::makeMeasurement(Acts::SourceLink{sourceLink}, pos, cov, Acts::eBoundLoc0, Acts::eBoundLoc1);
0214       measurements->emplace_back(std::move(measurement));
0215 #endif
0216
0217       ihit++;
0218     }
0219     // add proto track to the output collection
0220     protoTracks->emplace_back(std::move(track));
0221     return StatusCode::SUCCESS;
0222   }
0223 };
0224 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0225 DECLARE_COMPONENT(SingleTrackSourceLinker)
0226
0227 } // namespace Jug::Reco