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