src/components/TrackParamACTSSeeding.cpp

0001 #include <cmath>
0002 #include <algorithm>
0003 #include <unordered_map>
0004
0005 #include "Acts/ActsVersion.hpp"
0006 #include "Acts/Definitions/Units.hpp"
0007 #include "Acts/Definitions/Common.hpp"
0008 #include "Acts/Geometry/TrackingGeometry.hpp"
0009
0010 #include "Acts/Seeding/BinFinder.hpp"
0011 #include "Acts/Seeding/BinnedSPGroup.hpp"
0012 #include "Acts/Seeding/SpacePointGrid.hpp"
0013 #include "Acts/Seeding/SeedFilter.hpp"
0014 #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
0015 #include "Acts/Surfaces/PerigeeSurface.hpp"
0016 #include "Acts/Seeding/SeedFinderConfig.hpp"
0017 #include "Acts/Seeding/SeedFinder.hpp"
0018
0019 // Gaudi
0020 #include "GaudiAlg/GaudiAlgorithm.h"
0021 #include "GaudiKernel/ToolHandle.h"
0022 #include "GaudiAlg/Transformer.h"
0023 #include "GaudiAlg/GaudiTool.h"
0024 #include "GaudiKernel/RndmGenerators.h"
0025 #include "Gaudi/Property.h"
0026
0027 #include <k4FWCore/DataHandle.h>
0028 #include <k4Interface/IGeoSvc.h>
0029 #include "JugTrack/IActsGeoSvc.h"
0030 #include "JugTrack/DD4hepBField.h"
0031 #include "ActsExamples/EventData/Index.hpp"
0032 #include "ActsExamples/EventData/Measurement.hpp"
0033 #include "ActsExamples/EventData/Track.hpp"
0034
0035 #include "DDRec/CellIDPositionConverter.h"
0036
0037 #include "edm4eic/TrackerHitCollection.h"
0038 #include "Math/Vector3D.h"
0039
0040
0041   ///// (Reconstructed) track parameters e.g. close to the vertex.
0042   //using TrackParameters = Acts::CurvilinearTrackParameters;
0043
0044   ///// Container of reconstructed track states for multiple tracks.
0045   //using TrackParametersContainer = std::vector<TrackParameters>;
0046
0047   ///// MultiTrajectory definition
0048   //using Trajectory = Acts::MultiTrajectory<SourceLink>;
0049
0050   ///// Container for the truth fitting/finding track(s)
0051   //using TrajectoryContainer = std::vector<SimMultiTrajectory>;
0052
0053 namespace Jug::Reco {
0054
0055     /** Initial Track parameters from MC truth.
0056      *
0057      *  TrackParmetersContainer
0058      *  \ingroup tracking
0059      */
0060     class TrackParamACTSSeeding : public GaudiAlgorithm {
0061     public:
0062         DataHandle<edm4eic::TrackerHitCollection>
0063         m_inputHitCollection { "inputHitCollection",
0064             Gaudi::DataHandle::Reader, this };
0065         DataHandle<ActsExamples::TrackParametersContainer>
0066         m_outputInitialTrackParameters {
0067             "outputInitialTrackParameters",
0068             Gaudi::DataHandle::Writer, this };
0069
0070         SmartIF<IGeoSvc> m_geoSvc;
0071         SmartIF<IActsGeoSvc> m_actsGeoSvc;
0072         std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> m_converter;
0073
0074         // Acts::GeometryContext m_geoContext;
0075         std::shared_ptr<const Jug::BField::DD4hepBField> m_BField =
0076             nullptr;
0077         Acts::MagneticFieldContext m_fieldContext;
0078
0079         /// Index type to reference elements in a container.
0080         ///
0081         /// We do not expect to have more than 2^32 elements in any
0082         /// given container so a fixed sized integer type is
0083         /// sufficient.
0084         //using Index = eic::Index;
0085
0086         /// Space point representation of eic::TrackerHitData suitable
0087         /// for ACTS track seeding.
0088         class SpacePoint
0089         // : edm4eic::TrackerHitData
0090         {
0091         public:
0092             // Acts::Vector3 _dummy[16];
0093             Acts::Vector3 _position;
0094             Acts::Vector3 _positionError;
0095             int32_t _measurementIndex;
0096             const Acts::Surface *_surface;
0097             // Constructor to circumvent the fact that eic::TrackerHit
0098             // and associated classes are all non-polymorphic
0099             SpacePoint(const edm4eic::TrackerHit h,
0100                        const int32_t measurementIndex,
0101                        SmartIF<IActsGeoSvc> m_actsGeoSvc,
0102                        std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> m_converter)
0103                 : _measurementIndex(measurementIndex)
0104             {
0105                 _position[0] = h.getPosition().x;
0106                 _position[1] = h.getPosition().y;
0107                 _position[2] = h.getPosition().z;
0108                 _positionError[0] = h.getPositionError().xx;
0109                 _positionError[1] = h.getPositionError().yy;
0110                 _positionError[2] = h.getPositionError().zz;
0111                 const auto volumeId =
0112                     m_converter->findContext(h.getCellID())->identifier;
0113                 const auto its = m_actsGeoSvc->surfaceMap().find(volumeId);
0114                 if (its == m_actsGeoSvc->surfaceMap().end()) {
0115                     _surface = nullptr;
0116                 }
0117                 else {
0118                     _surface = its->second;
0119                 }
0120             }
0121             SpacePoint(const SpacePoint &sp)
0122                 : _position(sp._position),
0123                   _positionError(sp._positionError),
0124                   _measurementIndex(sp._measurementIndex),
0125                   _surface(sp._surface)
0126             {
0127             }
0128             float x() const { return _position[0]; }
0129             float y() const { return _position[1]; }
0130             float z() const { return _position[2]; }
0131             float r() const { return std::hypot(x(), y()); }
0132             float varianceR() const
0133             {
0134                 return (std::pow(x(), 2) * _positionError[0] +
0135                         std::pow(y(), 2) * _positionError[1]) /
0136                     (std::pow(x(), 2) + std::pow(y(), 2));
0137             }
0138             float varianceZ() const { return _positionError[2]; }
0139             constexpr uint32_t measurementIndex() const {
0140                 return _measurementIndex; }
0141             bool isOnSurface() const {
0142                 if (_surface == nullptr) {
0143                     return false;
0144                 }
0145                 return _surface->isOnSurface(
0146                      Acts::GeometryContext(), {x(), y(), z()},
0147                      {0, 0, 0});
0148             }
0149         };
0150
0151         static bool spCompare(SpacePoint r, SpacePoint s)
0152         {
0153             return
0154                 std::hypot(r.x(), r.y(), r.z()) <
0155                 std::hypot(s.x(), s.y(), s.z());
0156         }
0157
0158         /// Container of sim seed
0159         using SeedContainer = std::vector<Acts::Seed<SpacePoint>>;
0160
0161         /// A proto track is a collection of hits identified by their
0162         /// indices.
0163         using ProtoTrack = std::vector<ActsExamples::Index>;
0164         /// Container of proto tracks. Each proto track is identified by
0165         /// its index.
0166         using ProtoTrackContainer = std::vector<ProtoTrack>;
0167
0168         using SpacePointContainer = std::vector<SpacePoint>;
0169
0170         struct Config {
0171             /// Input space point collections.
0172             ///
0173             /// We allow multiple space point collections to allow
0174             /// different parts of the detector to use different
0175             /// algorithms for space point construction, e.g.
0176             /// single-hit space points for pixel-like detectors or
0177             /// double-hit space points for strip-like detectors.
0178             std::vector<std::string> inputSpacePoints;
0179             /// Output track seed collection.
0180             std::string outputSeeds;
0181             /// Output proto track collection.
0182             std::string outputProtoTracks;
0183
0184             float cotThetaMax = std::sinh(4.01);
0185             float minPt = 100 * Acts::UnitConstants::MeV / cotThetaMax;
0186             float rMax = 440 * Acts::UnitConstants::mm;
0187             float zMin = -1500 * Acts::UnitConstants::mm;
0188             float zMax = 1700 * Acts::UnitConstants::mm;
0189             float deltaRMin = 0 * Acts::UnitConstants::mm;
0190             float deltaRMax = 600 * Acts::UnitConstants::mm;
0191             //
0192             float collisionRegionMin = -250 * Acts::UnitConstants::mm;
0193             float collisionRegionMax = 250 * Acts::UnitConstants::mm;
0194             float maxSeedsPerSpM = 2;
0195             float sigmaScattering = 5;
0196             float radLengthPerSeed = 0.1;
0197             float impactMax = 3 * Acts::UnitConstants::mm;
0198             //
0199             float bFieldInZ = 1.7 * Acts::UnitConstants::T;
0200             float beamPosX = 0 * Acts::UnitConstants::mm;
0201             float beamPosY = 0 * Acts::UnitConstants::mm;
0202
0203             /// The minimum magnetic field to trigger the track
0204             /// parameters estimation
0205             double bFieldMin = 0.1 * Acts::UnitConstants::T;
0206
0207             /// Constant term of the loc0 resolution.
0208             double sigmaLoc0 = 25 * Acts::UnitConstants::um;
0209             /// Constant term of the loc1 resolution.
0210             double sigmaLoc1 = 100 * Acts::UnitConstants::um;
0211             /// Phi angular resolution.
0212             double sigmaPhi = 0.02 * Acts::UnitConstants::degree;
0213             /// Theta angular resolution.
0214             double sigmaTheta = 0.02 * Acts::UnitConstants::degree;
0215             /// q/p resolution.
0216             double sigmaQOverP = 0.1 / Acts::UnitConstants::GeV;
0217             /// Time resolution.
0218             double sigmaT0 = 1400 * Acts::UnitConstants::s;
0219
0220         int numPhiNeighbors = 1;
0221
0222             float deltaRMiddleMinSPRange = 10. * Acts::UnitConstants::mm;
0223             float deltaRMiddleMaxSPRange = 10. * Acts::UnitConstants::mm;
0224
0225             // vector containing the map of z bins in the top and bottom layers
0226             std::vector<std::pair<int, int> > zBinNeighborsTop;
0227             std::vector<std::pair<int, int> > zBinNeighborsBottom;
0228         } m_cfg;
0229         Acts::SpacePointGridConfig m_gridCfg;
0230         Acts::SpacePointGridOptions m_gridOpt;
0231         Acts::SeedFinderConfig<SpacePoint> m_finderCfg;
0232         Acts::SeedFinderOptions m_finderOpt;
0233         /// The track parameters covariance (assumed to be the same
0234         /// for all estimated track parameters for the moment)
0235         Acts::BoundSquareMatrix m_covariance =
0236             Acts::BoundSquareMatrix::Zero();
0237
0238     public:
0239         TrackParamACTSSeeding(const std::string &name,
0240                               ISvcLocator* svcLoc)
0241             : GaudiAlgorithm(name, svcLoc)
0242         {
0243             declareProperty("inputHitCollection",
0244                             m_inputHitCollection, "");
0245             declareProperty("outputInitialTrackParameters",
0246                             m_outputInitialTrackParameters, "");
0247         }
0248
0249         StatusCode initialize() override;
0250
0251         StatusCode execute() override;
0252     };
0253
0254
0255     StatusCode TrackParamACTSSeeding::initialize()
0256     {
0257         if (GaudiAlgorithm::initialize().isFailure()) {
0258             return StatusCode::FAILURE;
0259         }
0260
0261         m_geoSvc = service("GeoSvc");
0262         if (m_geoSvc == nullptr) {
0263             error() << "Unable to locate Geometry Service. " << endmsg;
0264             return StatusCode::FAILURE;
0265         }
0266         m_actsGeoSvc = service("ActsGeoSvc");
0267         if (m_actsGeoSvc == nullptr) {
0268             error() << "Unable to locate ACTS Geometry Service. " << endmsg;
0269             return StatusCode::FAILURE;
0270         }
0271         m_converter = std::make_shared<const dd4hep::rec::CellIDPositionConverter>(*(m_geoSvc->getDetector()));
0272
0273         m_BField = std::dynamic_pointer_cast<
0274             const Jug::BField::DD4hepBField>(
0275                 m_actsGeoSvc->getFieldProvider());
0276         m_fieldContext = Jug::BField::BFieldVariant(m_BField);
0277
0278         m_gridCfg.minPt = m_cfg.minPt;
0279         m_gridCfg.rMax = m_cfg.rMax;
0280         m_gridCfg.zMax = m_cfg.zMax;
0281         m_gridCfg.zMin = m_cfg.zMin;
0282         m_gridCfg.cotThetaMax = m_cfg.cotThetaMax;
0283
0284         m_gridCfg.deltaRMax = m_cfg.deltaRMax;
0285
0286         // Construct seed filter
0287         Acts::SeedFilterConfig filterCfg;
0288         filterCfg.maxSeedsPerSpM = m_cfg.maxSeedsPerSpM;
0289         m_finderCfg.seedFilter =
0290             std::make_unique<Acts::SeedFilter<SpacePoint>>(
0291                 Acts::SeedFilter<SpacePoint>(filterCfg));
0292
0293         m_finderCfg.rMax = m_cfg.rMax;
0294         m_finderCfg.deltaRMin = m_cfg.deltaRMin;
0295         m_finderCfg.deltaRMax = m_cfg.deltaRMax;
0296         m_finderCfg.collisionRegionMin = m_cfg.collisionRegionMin;
0297         m_finderCfg.collisionRegionMax = m_cfg.collisionRegionMax;
0298         m_finderCfg.zMin = m_cfg.zMin;
0299         m_finderCfg.zMax = m_cfg.zMax;
0300         m_finderCfg.maxSeedsPerSpM = m_cfg.maxSeedsPerSpM;
0301         m_finderCfg.cotThetaMax = m_cfg.cotThetaMax;
0302         m_finderCfg.sigmaScattering = m_cfg.sigmaScattering;
0303         m_finderCfg.radLengthPerSeed = m_cfg.radLengthPerSeed;
0304         m_finderCfg.minPt = m_cfg.minPt;
0305         m_finderCfg.impactMax = m_cfg.impactMax;
0306
0307         m_finderOpt.bFieldInZ = m_cfg.bFieldInZ;
0308         m_finderOpt.beamPos =
0309             Acts::Vector2(m_cfg.beamPosX, m_cfg.beamPosY);
0310
0311         m_finderCfg =
0312           m_finderCfg.toInternalUnits().calculateDerivedQuantities();
0313         m_finderOpt =
0314           m_finderOpt.toInternalUnits().calculateDerivedQuantities(
0315             m_finderCfg);
0316
0317         // Set up the track parameters covariance (the same for all
0318         // tracks)
0319         m_covariance(Acts::eBoundLoc0, Acts::eBoundLoc0) =
0320             std::pow(m_cfg.sigmaLoc0, 2);
0321         m_covariance(Acts::eBoundLoc1, Acts::eBoundLoc1) =
0322             std::pow(m_cfg.sigmaLoc1, 2);
0323         m_covariance(Acts::eBoundPhi, Acts::eBoundPhi) =
0324             std::pow(m_cfg.sigmaPhi, 2);
0325         m_covariance(Acts::eBoundTheta, Acts::eBoundTheta) =
0326             std::pow(m_cfg.sigmaTheta, 2);
0327         m_covariance(Acts::eBoundQOverP, Acts::eBoundQOverP) =
0328             std::pow(m_cfg.sigmaQOverP, 2);
0329         m_covariance(Acts::eBoundTime, Acts::eBoundTime) =
0330             std::pow(m_cfg.sigmaT0, 2);
0331
0332         return StatusCode::SUCCESS;
0333     }
0334
0335     StatusCode TrackParamACTSSeeding::execute()
0336     {
0337         const edm4eic::TrackerHitCollection *hits =
0338             m_inputHitCollection.get();
0339         // Create output collections
0340         auto initTrackParameters =
0341             m_outputInitialTrackParameters.createAndPut();
0342
0343         static SeedContainer seeds;
0344         static Acts::SeedFinder<SpacePoint>::SeedingState state;
0345
0346         // Sadly, eic::TrackerHit and eic::TrackerHitData are
0347     // non-polymorphic
0348         std::vector<SpacePoint> spacePoint;
0349         std::vector<const SpacePoint *> spacePointPtrs;
0350
0351         std::shared_ptr<const Acts::TrackingGeometry>
0352             trackingGeometry = m_actsGeoSvc->trackingGeometry();
0353
0354 #ifdef USE_LOCAL_COORD
0355         // Currently broken, possibly geometry issues
0356 #error Do not use, broken
0357         if (msgLevel(MSG::DEBUG)) {
0358             debug() << __FILE__ << ':' << __LINE__ << ": "
0359                     << sourceLinks->size() << ' '
0360                     << hits->size() << endmsg;
0361         }
0362         auto its = sourceLinks->begin();
0363         auto itm = measurements->begin();
0364         for (; its != sourceLinks->end() &&
0365                  itm != measurements->end();
0366              its++, itm++) {
0367             const Acts::Surface *surface =
0368                 trackingGeometry->findSurface(its->get().geometryId());
0369             if (surface != nullptr) {
0370                 Acts::Vector3 v =
0371                     surface->localToGlobal(
0372                         Acts::GeometryContext(),
0373                         {std::get<Acts::Measurement<
0374                          Acts::BoundIndices, 2>>(*itm).
0375                          parameters()[0],
0376                          std::get<Acts::Measurement<
0377                          Acts::BoundIndices, 2>>(*itm).
0378                          parameters()[1]},
0379                         {0, 0, 0});
0380                 if (msgLevel(MSG::DEBUG)) {
0381                     debug() << __FILE__ << ':' << __LINE__ << ": "
0382                             << its - sourceLinks->begin() << ' '
0383                         // << itm - measurements->begin() << ' '
0384                             << v[0] << ' ' << v[1] << ' ' << v[2]
0385                             << endmsg;
0386                 }
0387                 spacePoint.push_back(
0388                     SpacePoint(
0389                         edm4eic::TrackerHit(
0390                             static_cast<uint64_t>(spacePoint.size() + 1),
0391                             edm4eic::Vector3f(v[0], v[1], v[2]),
0392                             // Dummy uncertainties
0393                             edm4eic::CovDiag3f(25.0e-6 / 3.0,
0394                                                25.0e-6 / 3.0, 0.0),
0395                             0.0, 10.0, 0.05, 0.0),
0396                         static_cast<int32_t>(spacePoint.size())));
0397                 spacePointPtrs.push_back(&spacePoint.back());
0398         }
0399 #else // USE_LOCAL_COORD
0400         for(const auto &h : *hits) {
0401             spacePoint.push_back(SpacePoint(
0402                 h, static_cast<int32_t>(spacePoint.size() + 1),
0403                 m_actsGeoSvc, m_converter));
0404             spacePointPtrs.push_back(&spacePoint.back());
0405             if (msgLevel(MSG::DEBUG)) {
0406                 debug() << __FILE__ << ':' << __LINE__ << ": "
0407                         << ' ' << h.getPosition().x
0408                         << ' ' << h.getPosition().y
0409                         << ' ' << h.getPosition().z
0410                         << ' ' << h.getPositionError().xx
0411                         << ' ' << h.getPositionError().yy
0412                         << ' ' << h.getPositionError().zz
0413                         << ' ' << h.getTime()
0414                         << ' ' << h.getTimeError()
0415                         << ' ' << h.getEdep()
0416                         << ' ' << h.getEdepError()
0417                         << ' ' << spacePointPtrs.back()
0418                         << ' ' << spacePointPtrs.back()->measurementIndex()
0419                         << ' ' << spacePointPtrs.back()->isOnSurface()
0420                         << endmsg;
0421             }
0422         }
0423 #endif // USE_LOCAL_COORD
0424         if (msgLevel(MSG::DEBUG)) {
0425             debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg;
0426         }
0427
0428         auto extractGlobalQuantities =
0429             [=](const SpacePoint& sp, float, float, float) ->
0430             std::pair<Acts::Vector3, Acts::Vector2> {
0431                 Acts::Vector3 position { sp.x(), sp.y(), sp.z() };
0432                 Acts::Vector2 covariance {
0433                     sp.varianceR(), sp.varianceZ() };
0434                 return std::make_pair(position, covariance);
0435         };
0436         if (msgLevel(MSG::DEBUG)) {
0437             debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg;
0438         }
0439
0440         // extent used to store r range for middle spacepoint
0441         Acts::Extent rRangeSPExtent;
0442
0443         auto bottomBinFinder =
0444             std::make_shared<Acts::BinFinder<SpacePoint>>(
0445                 Acts::BinFinder<SpacePoint>(m_cfg.zBinNeighborsBottom,
0446                                             m_cfg.numPhiNeighbors));
0447         if (msgLevel(MSG::DEBUG)) {
0448             debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg;
0449         }
0450         auto topBinFinder =
0451             std::make_shared<Acts::BinFinder<SpacePoint>>(
0452                 Acts::BinFinder<SpacePoint>(m_cfg.zBinNeighborsTop,
0453                                             m_cfg.numPhiNeighbors));
0454         if (msgLevel(MSG::DEBUG)) {
0455             debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg;
0456         }
0457
0458         const float bFieldInZSave = m_gridOpt.bFieldInZ;
0459         const float minPtSave = m_gridCfg.minPt;
0460         m_gridCfg.minPt = 400 * Acts::UnitConstants::MeV;
0461         m_gridOpt.bFieldInZ =
0462             (m_gridCfg.minPt / Acts::UnitConstants::MeV) /
0463             (150.0 * (1.0 + FLT_EPSILON) *
0464              (m_gridCfg.rMax / Acts::UnitConstants::mm)) *
0465             1000.0 * Acts::UnitConstants::T;
0466         if (msgLevel(MSG::DEBUG)) {
0467             debug() << "createGrid() with temporary minPt = "
0468                     << m_gridCfg.minPt / Acts::UnitConstants::MeV
0469                     << " MeV, bFieldInZ = "
0470                     << m_gridOpt.bFieldInZ / (1000 * Acts::UnitConstants::T)
0471                     << " kT" << endmsg;
0472         }
0473         auto grid =
0474             Acts::SpacePointGridCreator::createGrid<SpacePoint>(
0475                 m_gridCfg, m_gridOpt);
0476         m_gridOpt.bFieldInZ = bFieldInZSave;
0477         m_gridCfg.minPt = minPtSave;
0478         if (msgLevel(MSG::DEBUG)) {
0479             debug() << "phiBins = "
0480                     << grid->axes().front()->getNBins()
0481                     << ", zBins = "
0482                     << grid->axes().back()->getNBins() << endmsg;
0483         }
0484
0485         auto spacePointsGrouping =
0486             Acts::BinnedSPGroup<SpacePoint>(
0487                 spacePointPtrs.begin(), spacePointPtrs.end(),
0488                 extractGlobalQuantities, bottomBinFinder,
0489                 topBinFinder, std::move(grid),
0490                 rRangeSPExtent,
0491                 m_finderCfg, m_finderOpt);
0492         auto finder = Acts::SeedFinder<SpacePoint>(m_finderCfg);
0493
0494         if (msgLevel(MSG::DEBUG)) {
0495             debug() << __FILE__ << ':' << __LINE__
0496                     << ": spacePointsGrouping.size() = "
0497                     << spacePointsGrouping.size() << endmsg;
0498         }
0499
0500         const Acts::Range1D<float> rMiddleSPRange(
0501             std::floor(rRangeSPExtent.min(Acts::binR) / 2) * 2 +
0502             m_cfg.deltaRMiddleMinSPRange,
0503             std::floor(rRangeSPExtent.max(Acts::binR) / 2) * 2 -
0504             m_cfg.deltaRMiddleMaxSPRange);
0505
0506         // Run the seeding
0507         seeds.clear();
0508
0509         for (const auto [bottom, middle, top] : spacePointsGrouping) {
0510             finder.createSeedsForGroup(
0511                 m_finderOpt, state, spacePointsGrouping.grid(),
0512                 std::back_inserter(seeds), bottom, middle, top, rMiddleSPRange
0513             );
0514         }
0515
0516         if (msgLevel(MSG::DEBUG)) {
0517             debug() << "seeds.size() = " << seeds.size() << endmsg;
0518         }
0519
0520         ActsExamples::TrackParametersContainer trackParameters;
0521         ProtoTrackContainer tracks;
0522         trackParameters.reserve(seeds.size());
0523         tracks.reserve(seeds.size());
0524
0525         std::shared_ptr<const Acts::MagneticFieldProvider>
0526             magneticField = m_actsGeoSvc->getFieldProvider();
0527
0528         // if (msgLevel(MSG::DEBUG)) { debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg; }
0529         auto bCache = magneticField->makeCache(m_fieldContext);
0530
0531         std::unordered_map<size_t, bool> spTaken;
0532
0533         for (size_t iseed = 0; iseed < seeds.size(); iseed++) {
0534             const auto &seed = seeds[iseed];
0535             // Get the bottom space point and its reference surface
0536             const auto bottomSP = seed.sp().front();
0537             auto hitIdx = bottomSP->measurementIndex();
0538             // const Acts::Surface *surface = nullptr;
0539             const Acts::Surface *surface = bottomSP->_surface;
0540             if (surface == nullptr) {
0541                 if (msgLevel(MSG::DEBUG)) {
0542                     debug() << "hit " << hitIdx << " ("
0543                             << bottomSP->x() << ", " << bottomSP->y()
0544                             << ", " << bottomSP->z()
0545                             << ") lost its surface" << endmsg;
0546                 }
0547             }
0548             if (std::find_if(spacePoint.begin(), spacePoint.end(),
0549                              [&surface](const SpacePoint sp) {
0550                                  return surface == sp._surface;
0551                              }) == spacePoint.end()) {
0552                 if (msgLevel(MSG::DEBUG)) {
0553                     debug() << "hit " << hitIdx
0554                             << " has a surface that was never "
0555                             << "associated with a hit" << endmsg;
0556                 }
0557                 surface = nullptr;
0558             }
0559             if (surface == nullptr && msgLevel(MSG::DEBUG)) {
0560                 debug() << "hit " << hitIdx
0561                         << " is not found in the tracking gemetry"
0562                         << endmsg;
0563                 continue;
0564             }
0565
0566             // Get the magnetic field at the bottom space point
0567             auto fieldRes = magneticField->getField(
0568                 {bottomSP->x(), bottomSP->y(), bottomSP->z()},
0569                 bCache);
0570             // if (msgLevel(MSG::DEBUG)) { debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg; }
0571             // Estimate the track parameters from seed
0572             auto optParams = Acts::estimateTrackParamsFromSeed(
0573                 Acts::GeometryContext(),
0574                 seed.sp().begin(), seed.sp().end(),
0575                 *surface, *fieldRes, m_cfg.bFieldMin);
0576             // if (msgLevel(MSG::DEBUG)) { debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg; }
0577             if (not optParams.has_value()) {
0578                 debug() << "Estimation of track parameters for seed "
0579                         << iseed << " failed." << endmsg;
0580                 continue;
0581             }
0582             else if (!(spTaken[seed.sp()[0]->measurementIndex()] ||
0583                        spTaken[seed.sp()[1]->measurementIndex()] ||
0584                        spTaken[seed.sp()[2]->measurementIndex()])) {
0585                 const auto& params = optParams.value();
0586                 initTrackParameters->emplace_back(
0587                     surface->getSharedPtr(), params,
0588                     m_covariance, Acts::ParticleHypothesis::pion());
0589                 // Activate/deactivate for unique seed filtering
0590 #if 0
0591                 spTaken[seed.sp()[0]->measurementIndex()] = true;
0592                 spTaken[seed.sp()[1]->measurementIndex()] = true;
0593                 spTaken[seed.sp()[2]->measurementIndex()] = true;
0594 #endif
0595             }
0596             // if (msgLevel(MSG::DEBUG)) { debug() << __FILE__ << ':' << __LINE__ << ": " << endmsg; }
0597         }
0598
0599         return StatusCode::SUCCESS;
0600     }
0601
0602     DECLARE_COMPONENT(TrackParamACTSSeeding)
0603 } // namespace Jug::reco