TrackFindingExaTrkX/src/PrototracksToParameters.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2024 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/TrackFindingExaTrkX/PrototracksToParameters.hpp"
0010
0011 #include "Acts/Seeding/BinnedGroup.hpp"
0012 #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
0013 #include "Acts/Seeding/InternalSpacePoint.hpp"
0014 #include "Acts/Seeding/SeedFilter.hpp"
0015 #include "Acts/Seeding/SeedFinder.hpp"
0016 #include "Acts/Seeding/SeedFinderConfig.hpp"
0017 #include "Acts/Utilities/Zip.hpp"
0018 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0019 #include "ActsExamples/EventData/ProtoTrack.hpp"
0020 #include "ActsExamples/EventData/SimSeed.hpp"
0021 #include "ActsExamples/Framework/WhiteBoard.hpp"
0022 #include "ActsExamples/Utilities/EventDataTransforms.hpp"
0023
0024 #include <algorithm>
0025
0026 using namespace ActsExamples;
0027 using namespace Acts::UnitLiterals;
0028
0029 namespace ActsExamples {
0030
0031 PrototracksToParameters::PrototracksToParameters(Config cfg,
0032                                                  Acts::Logging::Level lvl)
0033     : IAlgorithm("PrototracksToParsAndSeeds", lvl), m_cfg(std::move(cfg)) {
0034   m_outputSeeds.initialize(m_cfg.outputSeeds);
0035   m_outputProtoTracks.initialize(m_cfg.outputProtoTracks);
0036   m_inputProtoTracks.initialize(m_cfg.inputProtoTracks);
0037   m_inputSpacePoints.initialize(m_cfg.inputSpacePoints);
0038   m_outputParameters.initialize(m_cfg.outputParameters);
0039
0040   if (m_cfg.geometry == nullptr) {
0041     throw std::invalid_argument("No geometry given");
0042   }
0043   if (m_cfg.magneticField == nullptr) {
0044     throw std::invalid_argument("No magnetic field given");
0045   }
0046
0047   // Set up the track parameters covariance (the same for all tracks)
0048   for (std::size_t i = Acts::eBoundLoc0; i < Acts::eBoundSize; ++i) {
0049     m_covariance(i, i) = m_cfg.initialVarInflation[i] * m_cfg.initialSigmas[i] *
0050                          m_cfg.initialSigmas[i];
0051   }
0052 }
0053
0054 PrototracksToParameters::~PrototracksToParameters() {}
0055
0056 ProcessCode PrototracksToParameters::execute(
0057     const AlgorithmContext &ctx) const {
0058   auto bCache = m_cfg.magneticField->makeCache(ctx.magFieldContext);
0059   const auto &sps = m_inputSpacePoints(ctx);
0060   auto prototracks = m_inputProtoTracks(ctx);
0061
0062   // Make some lookup tables. Allocate space for the maximum number of indices
0063   // (max 2 source links per spacepoint)
0064   std::vector<const SimSpacePoint *> indexToSpacepoint(2 * sps.size(), nullptr);
0065   std::vector<Acts::GeometryIdentifier> indexToGeoId(
0066       2 * sps.size(), Acts::GeometryIdentifier{0});
0067
0068   for (const auto &sp : sps) {
0069     for (const auto &sl : sp.sourceLinks()) {
0070       const auto &isl = sl.template get<IndexSourceLink>();
0071       indexToSpacepoint[isl.index()] = &sp;
0072       indexToGeoId[isl.index()] = isl.geometryId();
0073     }
0074   }
0075
0076   ProtoTrackContainer seededTracks;
0077   seededTracks.reserve(prototracks.size());
0078
0079   SimSeedContainer seeds;
0080   seeds.reserve(prototracks.size());
0081
0082   TrackParametersContainer parameters;
0083   parameters.reserve(prototracks.size());
0084
0085   // Loop over the prototracks to make seeds
0086   ProtoTrack tmpTrack;
0087   std::vector<const SimSpacePoint *> tmpSps;
0088   std::size_t skippedTracks = 0;
0089   for (auto &track : prototracks) {
0090     ACTS_VERBOSE("Try to get seed from prototrack with " << track.size()
0091                                                          << " hits");
0092     // Make prototrack unique with respect to volume and layer
0093     // so we don't get a seed where we have two spacepoints on the same layer
0094
0095     // Here, we want to create a seed only if the prototrack with removed unique
0096     // layer-volume spacepoints has 3 or more hits. However, if this is the
0097     // case, we want to keep the whole prototrack. Therefore, we operate on a
0098     // tmpTrack.
0099     std::sort(track.begin(), track.end(), [&](auto a, auto b) {
0100       if (indexToGeoId[a].volume() != indexToGeoId[b].volume()) {
0101         return indexToGeoId[a].volume() < indexToGeoId[b].volume();
0102       }
0103       return indexToGeoId[a].layer() < indexToGeoId[b].layer();
0104     });
0105
0106     tmpTrack.clear();
0107     std::unique_copy(
0108         track.begin(), track.end(), std::back_inserter(tmpTrack),
0109         [&](auto a, auto b) {
0110           return indexToGeoId[a].volume() == indexToGeoId[b].volume() &&
0111                  indexToGeoId[a].layer() == indexToGeoId[b].layer();
0112         });
0113
0114     // in this case we cannot seed properly
0115     if (tmpTrack.size() < 3) {
0116       ACTS_DEBUG(
0117           "Cannot seed because less then three hits with unique (layer, "
0118           "volume)");
0119       skippedTracks++;
0120       continue;
0121     }
0122
0123     // Make the seed
0124     tmpSps.clear();
0125     std::transform(track.begin(), track.end(), std::back_inserter(tmpSps),
0126                    [&](auto i) { return indexToSpacepoint[i]; });
0127     tmpSps.erase(std::remove_if(tmpSps.begin(), tmpSps.end(),
0128                                 [](auto sp) { return sp == nullptr; }),
0129                  tmpSps.end());
0130
0131     if (tmpSps.size() < 3) {
0132       ACTS_WARNING("Could not find all spacepoints, skip");
0133       skippedTracks++;
0134       continue;
0135     }
0136
0137     std::sort(tmpSps.begin(), tmpSps.end(),
0138               [](const auto &a, const auto &b) { return a->r() < b->r(); });
0139
0140     // Simply use r = m*z + t and solve for r=0 to find z vertex position...
0141     // Probably not the textbook way to do
0142     const auto m = (tmpSps.back()->r() - tmpSps.front()->r()) /
0143                    (tmpSps.back()->z() - tmpSps.front()->z());
0144     const auto t = tmpSps.front()->r() - m * tmpSps.front()->z();
0145     const auto z_vertex = -t / m;
0146     const auto s = tmpSps.size();
0147
0148     SimSeed seed =
0149         m_cfg.buildTightSeeds
0150             ? SimSeed(*tmpSps[0], *tmpSps[1], *tmpSps[2], z_vertex)
0151             : SimSeed(*tmpSps[0], *tmpSps[s / 2], *tmpSps[s - 1], z_vertex);
0152
0153     // Compute parameters
0154     const auto &bottomSP = seed.sp().front();
0155     const auto geoId = bottomSP->sourceLinks()
0156                            .front()
0157                            .template get<IndexSourceLink>()
0158                            .geometryId();
0159     const auto &surface = *m_cfg.geometry->findSurface(geoId);
0160
0161     auto field = m_cfg.magneticField->getField(
0162         {bottomSP->x(), bottomSP->y(), bottomSP->z()}, bCache);
0163     if (!field.ok()) {
0164       ACTS_ERROR("Field lookup error: " << field.error());
0165       return ProcessCode::ABORT;
0166     }
0167
0168     auto pars = Acts::estimateTrackParamsFromSeed(
0169         ctx.geoContext, seed.sp().begin(), seed.sp().end(), surface, *field,
0170         m_cfg.bFieldMin);
0171
0172     if (not pars) {
0173       ACTS_WARNING("Skip track because of bad params");
0174     }
0175
0176     seededTracks.push_back(track);
0177     seeds.emplace_back(std::move(seed));
0178     parameters.push_back(Acts::BoundTrackParameters(
0179         surface.getSharedPtr(), *pars, m_covariance, m_cfg.particleHypothesis));
0180   }
0181
0182   if (skippedTracks > 0) {
0183     ACTS_WARNING("Skipped seeding of " << skippedTracks);
0184   }
0185
0186   ACTS_DEBUG("Seeded " << seeds.size() << " out of " << prototracks.size()
0187                        << " prototracks");
0188
0189   m_outputSeeds(ctx, std::move(seeds));
0190   m_outputProtoTracks(ctx, std::move(seededTracks));
0191   m_outputParameters(ctx, std::move(parameters));
0192
0193   return ProcessCode::SUCCESS;
0194 }
0195
0196 }  // namespace ActsExamples