acts/src/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 "GaudiAlg/GaudiAlgorithm.h"
0007 #include "GaudiAlg/GaudiTool.h"
0008 #include "GaudiAlg/Transformer.h"
0009 #include "GaudiKernel/RndmGenerators.h"
0010 #include "GaudiKernel/ToolHandle.h"
0011
0012 #include "JugBase/DataHandle.h"
0013 #include "JugBase/IGeoSvc.h"
0014
0015 #include "DD4hep/DD4hepUnits.h"
0016 #include "DD4hep/Volumes.h"
0017 #include "DDRec/CellIDPositionConverter.h"
0018 #include "DDRec/Surface.h"
0019 #include "DDRec/SurfaceManager.h"
0020
0021 #include "Acts/Definitions/Common.hpp"
0022 #include "Acts/Definitions/Units.hpp"
0023 #include "Acts/Geometry/TrackingGeometry.hpp"
0024 #include "Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp"
0025 #include "Acts/Surfaces/Surface.hpp"
0026
0027 #include "JugTrack/GeometryContainers.hpp"
0028 #include "JugTrack/Index.hpp"
0029 #include "JugTrack/IndexSourceLink.hpp"
0030 #include "JugTrack/Measurement.hpp"
0031 #include "JugTrack/ProtoTrack.hpp"
0032
0033 #include "edm4eic/TrackerHitCollection.h"
0034
0035 namespace Jug::Reco {
0036
0037 /** Single Track source Linker and proto tracks.
0038  *
0039  * This algorithm assumes only single track events.
0040  *
0041  * \ingroup tracking
0042  */
0043 class SingleTrackSourceLinker : public GaudiAlgorithm {
0044 private:
0045   DataHandle<edm4eic::TrackerHitCollection> m_inputHitCollection{"inputHitCollection", Gaudi::DataHandle::Reader, this};
0046   DataHandle<std::list<IndexSourceLink>> m_sourceLinkStorage{"sourceLinkStorage", Gaudi::DataHandle::Writer, this};
0047   DataHandle<IndexSourceLinkContainer> m_outputSourceLinks{"outputSourceLinks", Gaudi::DataHandle::Writer, this};
0048   DataHandle<MeasurementContainer> m_outputMeasurements{"outputMeasurements", Gaudi::DataHandle::Writer, this};
0049   DataHandle<ProtoTrackContainer> m_outputProtoTracks{"outputProtoTracks", Gaudi::DataHandle::Writer, this};
0050   /// Pointer to the geometry service
0051   SmartIF<IGeoSvc> m_geoSvc;
0052
0053 public:
0054   SingleTrackSourceLinker(const std::string& name, ISvcLocator* svcLoc) : GaudiAlgorithm(name, svcLoc) {
0055     declareProperty("inputHitCollection", m_inputHitCollection, "");
0056     declareProperty("sourceLinkStorage", m_sourceLinkStorage, "");
0057     declareProperty("outputSourceLinks", m_outputSourceLinks, "");
0058     declareProperty("outputMeasurements", m_outputMeasurements, "");
0059     declareProperty("outputProtoTracks", m_outputProtoTracks, "");
0060   }
0061
0062   StatusCode initialize() override {
0063     if (GaudiAlgorithm::initialize().isFailure()) {
0064       return StatusCode::FAILURE;
0065     }
0066     m_geoSvc = service("GeoSvc");
0067     if (!m_geoSvc) {
0068       error() << "Unable to locate Geometry Service. "
0069               << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
0070       return StatusCode::FAILURE;
0071     }
0072     return StatusCode::SUCCESS;
0073   }
0074
0075   StatusCode execute() override {
0076     // input collection
0077     const edm4eic::TrackerHitCollection* hits = m_inputHitCollection.get();
0078     // Create output collections
0079     auto* linkStorage  = m_sourceLinkStorage.createAndPut();
0080     auto* sourceLinks  = m_outputSourceLinks.createAndPut();
0081     auto* measurements = m_outputMeasurements.createAndPut();
0082     auto* protoTracks  = m_outputProtoTracks.createAndPut();
0083     // IndexMultimap<ActsFatras::Barcode> hitParticlesMap;
0084     // IndexMultimap<Index> hitSimHitsMap;
0085     sourceLinks->reserve(hits->size());
0086     measurements->reserve(hits->size());
0087
0088     // assume single track --> one ProtoTrack
0089     ProtoTrack track;
0090     track.reserve((*hits).size());
0091
0092     if (msgLevel(MSG::DEBUG)) {
0093       debug() << (*hits).size() << " hits " << endmsg;
0094     }
0095     int ihit = 0;
0096     for (const auto& ahit : *hits) {
0097
0098       track.emplace_back(ihit);
0099
0100       const auto* vol_ctx = m_geoSvc->cellIDPositionConverter()->findContext(ahit.getCellID());
0101       auto vol_id   = vol_ctx->identifier;
0102       const auto is = m_geoSvc->surfaceMap().find(vol_id);
0103       if (is == m_geoSvc->surfaceMap().end()) {
0104         if (msgLevel(MSG::DEBUG)) {
0105           debug() << " vol_id (" << vol_id << ")  not found in m_surfaces!!!!" << endmsg;
0106         }
0107         continue;
0108       }
0109       const Acts::Surface* surface = is->second;
0110
0111       // NOTE
0112       // Here is where the important hit and tracking geometry is connected.
0113       // A "Measurement" is constructed to for each hit which makes the connection to
0114       // the tracking surface and covariance matrix
0115
0116       auto acts_pos = surface
0117                           ->globalToLocal(Acts::GeometryContext(),
0118                                           {ahit.getPosition().x, ahit.getPosition().y, ahit.getPosition().z}, {0, 0, 0})
0119                           .value();
0120
0121       if (msgLevel(MSG::DEBUG)) {
0122         auto volman  = m_geoSvc->detector()->volumeManager();
0123         auto detelem = volman.lookupDetElement(vol_id);
0124         auto local_pos =
0125             detelem.nominal().worldToLocal({ahit.getPosition().x, ahit.getPosition().y, ahit.getPosition().z});
0126         debug() << "===== Debugging hit =====" << endmsg;
0127         debug() << "DD4hep global pos (" << ahit.getPosition().x << "," << ahit.getPosition().y << ","
0128                 << ahit.getPosition().z << ")" << endmsg;
0129         debug() << "DD4hep local  pos (" << local_pos.x() << "," << local_pos.y() << "," << local_pos.z() << ")"
0130                 << endmsg;
0131         debug() << "ACTS local position : (" << acts_pos[0] << "," << acts_pos[1] << ")" << endmsg;
0132         debug() << "ACTS surface center : " << surface->center(Acts::GeometryContext()).transpose() << endmsg;
0133         debug() << "DD4hep DetElement center : "
0134                 << detelem.nominal().localToWorld(detelem.placement().position()) / dd4hep::mm << endmsg;
0135       }
0136       // construct the vector of measured parameters (2d getPosition in this case)
0137       Acts::Vector2 pos(acts_pos.x(), acts_pos.y());
0138
0139       // construct the covariance matrix
0140       Acts::SymMatrix2 cov = Acts::SymMatrix2::Zero();
0141       cov(0, 0)            = ahit.getPositionError().xx * Acts::UnitConstants::mm * Acts::UnitConstants::mm;
0142       cov(1, 1)            = ahit.getPositionError().yy * Acts::UnitConstants::mm * Acts::UnitConstants::mm;
0143
0144       // Above we only consider the two position coordinates the comment below shows how to add time
0145       // which we will probably want to try later.
0146       //
0147       // Acts::SymMatrix3 cov;
0148       // cov << 0.05, 0., 0., 0., 0.05, 0., 0., 0., 900. * Acts::UnitConstants::ps * Acts::UnitConstants::ps;
0149       // Acts::Vector3 par(localX, localY, simHit.time());
0150
0151       linkStorage->emplace_back(surface->geometryId(), ihit);
0152       IndexSourceLink& sourceLink = linkStorage->back();
0153       auto meas =
0154           Acts::makeMeasurement(sourceLink, pos, cov, Acts::eBoundLoc0, Acts::eBoundLoc1); //, Acts::eBoundTime);
0155
0156       // add to output containers. since the input is already geometry-order,
0157       // new elements in geometry containers can just be appended at the end.
0158       sourceLinks->emplace_hint(sourceLinks->end(), sourceLink);
0159       measurements->emplace_back(std::move(meas));
0160
0161       ihit++;
0162     }
0163     // add proto track to the output collection
0164     protoTracks->emplace_back(std::move(track));
0165     return StatusCode::SUCCESS;
0166   }
0167 };
0168 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0169 DECLARE_COMPONENT(SingleTrackSourceLinker)
0170
0171 } // namespace Jug::Reco