Utilities/src/ProtoTracksToParameters.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Utilities/ProtoTracksToParameters.hpp"
0010
0011 #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
0012 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0013 #include "ActsExamples/EventData/ProtoTrack.hpp"
0014 #include "ActsExamples/EventData/SpacePoint.hpp"
0015
0016 #include <algorithm>
0017 #include <tuple>
0018
0019 using namespace Acts;
0020
0021 namespace ActsExamples {
0022
0023 ProtoTracksToParameters::ProtoTracksToParameters(
0024     Config cfg, std::unique_ptr<const Acts::Logger> logger)
0025     : IAlgorithm("ProtoTracksToParameters", std::move(logger)),
0026       m_cfg(std::move(cfg)) {
0027   m_outputSeeds.initialize(m_cfg.outputSeeds);
0028   m_outputProtoTracks.initialize(m_cfg.outputProtoTracks);
0029   m_inputProtoTracks.initialize(m_cfg.inputProtoTracks);
0030   m_inputSpacePoints.initialize(m_cfg.inputSpacePoints);
0031   m_outputParameters.initialize(m_cfg.outputParameters);
0032
0033   if (m_cfg.geometry == nullptr) {
0034     throw std::invalid_argument("No geometry given");
0035   }
0036   if (m_cfg.magneticField == nullptr) {
0037     throw std::invalid_argument("No magnetic field given");
0038   }
0039
0040   // Set up the track parameters covariance (the same for all tracks)
0041   for (std::size_t i = eBoundLoc0; i < eBoundSize; ++i) {
0042     m_covariance(i, i) = m_cfg.initialVarInflation[i] * m_cfg.initialSigmas[i] *
0043                          m_cfg.initialSigmas[i];
0044   }
0045 }
0046
0047 ProtoTracksToParameters::~ProtoTracksToParameters() = default;
0048
0049 ProcessCode ProtoTracksToParameters::execute(
0050     const AlgorithmContext &ctx) const {
0051   auto bCache = m_cfg.magneticField->makeCache(ctx.magFieldContext);
0052   const auto &sps = m_inputSpacePoints(ctx);
0053   auto protoTracks = m_inputProtoTracks(ctx);
0054
0055   // Make some lookup tables and pre-allocate some space
0056   // Note this is a heuristic, since it is not garantueed that each measurement
0057   // is part of a space point
0058   std::vector<SpacePointIndex> indexToSpacePoint(2 * sps.size(),
0059                                                  kSpacePointIndex2Invalid);
0060   std::vector<GeometryIdentifier> indexToGeoId(2 * sps.size(),
0061                                                GeometryIdentifier{0});
0062
0063   for (const auto &sp : sps) {
0064     for (const auto &sl : sp.sourceLinks()) {
0065       const auto &isl = sl.template get<IndexSourceLink>();
0066       if (isl.index() >= indexToSpacePoint.size()) {
0067         indexToSpacePoint.resize(isl.index() + 1, kSpacePointIndex2Invalid);
0068         indexToGeoId.resize(isl.index() + 1, GeometryIdentifier{0});
0069       }
0070       indexToSpacePoint.at(isl.index()) = sp.index();
0071       indexToGeoId.at(isl.index()) = isl.geometryId();
0072     }
0073   }
0074
0075   ProtoTrackContainer seededTracks;
0076   seededTracks.reserve(protoTracks.size());
0077
0078   SeedContainer seeds;
0079   seeds.assignSpacePointContainer(sps);
0080   seeds.reserve(protoTracks.size());
0081
0082   TrackParametersContainer parameters;
0083   parameters.reserve(protoTracks.size());
0084
0085   // Loop over the proto tracks to make seeds
0086   ProtoTrack tmpTrack;
0087   std::vector<SpacePointIndex> tmpSps;
0088   std::size_t skippedTracks = 0;
0089   for (auto &track : protoTracks) {
0090     ACTS_VERBOSE("Try to get seed from proto track with " << track.size()
0091                                                           << " hits");
0092     // Make proto track unique with respect to volume and layer so we don't get
0093     // a seed where we have two space points on the same layer
0094
0095     // Here, we want to create a seed only if the proto track with removed
0096     // unique layer-volume space points has 3 or more hits. However, if this is
0097     // the case, we want to keep the whole proto track. Therefore, we operate on
0098     // a tmpTrack.
0099     std::ranges::sort(track, {}, [&](const auto &t) {
0100       return std::make_tuple(indexToGeoId[t].volume(), indexToGeoId[t].layer());
0101     });
0102
0103     tmpTrack.clear();
0104     std::unique_copy(
0105         track.begin(), track.end(), std::back_inserter(tmpTrack),
0106         [&](auto a, auto b) {
0107           return indexToGeoId[a].volume() == indexToGeoId[b].volume() &&
0108                  indexToGeoId[a].layer() == indexToGeoId[b].layer();
0109         });
0110
0111     // in this case we cannot seed properly
0112     if (tmpTrack.size() < 3) {
0113       ACTS_DEBUG(
0114           "Cannot seed because less then three hits with unique (layer, "
0115           "volume)");
0116       skippedTracks++;
0117       continue;
0118     }
0119
0120     // Make the seed
0121     auto result =
0122         track | std::views::filter([&](auto i) {
0123           return i < indexToSpacePoint.size() &&
0124                  indexToSpacePoint.at(i) != kSpacePointIndex2Invalid;
0125         }) |
0126         std::views::transform([&](auto i) { return indexToSpacePoint.at(i); });
0127     tmpSps.clear();
0128     std::ranges::copy(result, std::back_inserter(tmpSps));
0129
0130     if (tmpSps.size() < 3) {
0131       ACTS_WARNING("Could not find all space points, skip");
0132       skippedTracks++;
0133       continue;
0134     }
0135
0136     std::ranges::sort(
0137         tmpSps, {}, [&sps](const SpacePointIndex &t) { return sps.at(t).r(); });
0138
0139     // Simply use r = m*z + t and solve for r=0 to find z vertex position...
0140     // Probably not the textbook way to do
0141     const float m = (sps.at(tmpSps.back()).r() - sps.at(tmpSps.front()).r()) /
0142                     (sps.at(tmpSps.back()).z() - sps.at(tmpSps.front()).z());
0143     const float t = sps.at(tmpSps.front()).r() - m * sps.at(tmpSps.front()).z();
0144     const float vertexZ = -t / m;
0145     const std::size_t s = tmpSps.size();
0146
0147     auto seed = seeds.createSeed();
0148     if (m_cfg.buildTightSeeds) {
0149       seed.assignSpacePointIndices(
0150           std::array{tmpSps.at(0), tmpSps.at(1), tmpSps.at(2)});
0151     } else {
0152       seed.assignSpacePointIndices(
0153           std::array{tmpSps.at(0), tmpSps.at(s / 2), tmpSps.at(s - 1)});
0154     }
0155     seed.vertexZ() = vertexZ;
0156
0157     // Compute parameters
0158
0159     const ConstSpacePointProxy bottomSp = seed.spacePoints()[0];
0160     const ConstSpacePointProxy middleSp = seed.spacePoints()[1];
0161     const ConstSpacePointProxy topSp = seed.spacePoints()[2];
0162
0163     const Acts::Vector3 bottomSpVec{bottomSp.x(), bottomSp.y(), bottomSp.z()};
0164     const Acts::Vector3 middleSpVec{middleSp.x(), middleSp.y(), middleSp.z()};
0165     const Acts::Vector3 topSpVec{topSp.x(), topSp.y(), topSp.z()};
0166
0167     const auto bottomGeoId =
0168         bottomSp.sourceLinks()[0].template get<IndexSourceLink>().geometryId();
0169     const Surface *bottomSurface = m_cfg.geometry->findSurface(bottomGeoId);
0170     if (bottomSurface == nullptr) {
0171       ACTS_WARNING(
0172           "Surface from source link is not found in the tracking geometry");
0173       continue;
0174     }
0175
0176     // Get the magnetic field at the bottom space point
0177     const auto fieldRes = m_cfg.magneticField->getField(bottomSpVec, bCache);
0178     if (!fieldRes.ok()) {
0179       ACTS_ERROR("Field lookup error: " << fieldRes.error());
0180       return ProcessCode::ABORT;
0181     }
0182     const Acts::Vector3 &field = *fieldRes;
0183
0184     if (field.norm() < m_cfg.bFieldMin) {
0185       ACTS_WARNING("Magnetic field at seed is too small " << field.norm());
0186       continue;
0187     }
0188
0189     // Estimate the track parameters from seed
0190     Acts::Result<Acts::BoundVector> boundParams =
0191         Acts::estimateTrackParamsFromSeed(
0192             ctx.geoContext, *bottomSurface, bottomSpVec,
0193             std::isnan(bottomSp.time()) ? 0.0 : bottomSp.time(), middleSpVec,
0194             topSpVec, field);
0195     if (!boundParams.ok()) {
0196       ACTS_WARNING("Failed to estimate track parameters from seed: "
0197                    << boundParams.error().message());
0198       continue;
0199     }
0200
0201     seededTracks.push_back(track);
0202     parameters.emplace_back(bottomSurface->getSharedPtr(), *boundParams,
0203                             m_covariance, m_cfg.particleHypothesis);
0204   }
0205
0206   if (skippedTracks > 0) {
0207     ACTS_WARNING("Skipped seeding of " << skippedTracks);
0208   }
0209
0210   ACTS_DEBUG("Seeded " << seeds.size() << " out of " << protoTracks.size()
0211                        << " proto tracks");
0212
0213   m_outputSeeds(ctx, std::move(seeds));
0214   m_outputProtoTracks(ctx, std::move(seededTracks));
0215   m_outputParameters(ctx, std::move(parameters));
0216
0217   return ProcessCode::SUCCESS;
0218 }
0219
0220 }  // namespace ActsExamples