TrackFindingGnn/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/TrackFindingGnn/PrototracksToParameters.hpp"
0010
0011 #include "Acts/Seeding/BinnedGroup.hpp"
0012 #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
0013 #include "Acts/Seeding/SeedFilter.hpp"
0014 #include "Acts/Seeding/SeedFinder.hpp"
0015 #include "Acts/Seeding/SeedFinderConfig.hpp"
0016 #include "Acts/Utilities/Zip.hpp"
0017 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0018 #include "ActsExamples/EventData/ProtoTrack.hpp"
0019 #include "ActsExamples/EventData/SimSeed.hpp"
0020 #include "ActsExamples/Framework/WhiteBoard.hpp"
0021 #include "ActsExamples/Utilities/EventDataTransforms.hpp"
0022
0023 #include <algorithm>
0024 #include <tuple>
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 and pre-allocate some space
0063   // Note this is a heuristic, since it is not garantueed that each measurement
0064   // is part of a spacepoint
0065   std::vector<const SimSpacePoint *> indexToSpacepoint(2 * sps.size(), nullptr);
0066   std::vector<Acts::GeometryIdentifier> indexToGeoId(
0067       2 * sps.size(), Acts::GeometryIdentifier{0});
0068
0069   for (const auto &sp : sps) {
0070     for (const auto &sl : sp.sourceLinks()) {
0071       const auto &isl = sl.template get<IndexSourceLink>();
0072       if (isl.index() >= indexToSpacepoint.size()) {
0073         indexToSpacepoint.resize(isl.index() + 1, nullptr);
0074         indexToGeoId.resize(isl.index() + 1, Acts::GeometryIdentifier{0});
0075       }
0076       indexToSpacepoint.at(isl.index()) = &sp;
0077       indexToGeoId.at(isl.index()) = isl.geometryId();
0078     }
0079   }
0080
0081   ProtoTrackContainer seededTracks;
0082   seededTracks.reserve(prototracks.size());
0083
0084   SimSeedContainer seeds;
0085   seeds.reserve(prototracks.size());
0086
0087   TrackParametersContainer parameters;
0088   parameters.reserve(prototracks.size());
0089
0090   // Loop over the prototracks to make seeds
0091   ProtoTrack tmpTrack;
0092   std::vector<const SimSpacePoint *> tmpSps;
0093   std::size_t skippedTracks = 0;
0094   for (auto &track : prototracks) {
0095     ACTS_VERBOSE("Try to get seed from prototrack with " << track.size()
0096                                                          << " hits");
0097     // Make prototrack unique with respect to volume and layer
0098     // so we don't get a seed where we have two spacepoints on the same layer
0099
0100     // Here, we want to create a seed only if the prototrack with removed unique
0101     // layer-volume spacepoints has 3 or more hits. However, if this is the
0102     // case, we want to keep the whole prototrack. Therefore, we operate on a
0103     // tmpTrack.
0104     std::ranges::sort(track, {}, [&](const auto &t) {
0105       return std::make_tuple(indexToGeoId[t].volume(), indexToGeoId[t].layer());
0106     });
0107
0108     tmpTrack.clear();
0109     std::unique_copy(
0110         track.begin(), track.end(), std::back_inserter(tmpTrack),
0111         [&](auto a, auto b) {
0112           return indexToGeoId[a].volume() == indexToGeoId[b].volume() &&
0113                  indexToGeoId[a].layer() == indexToGeoId[b].layer();
0114         });
0115
0116     // in this case we cannot seed properly
0117     if (tmpTrack.size() < 3) {
0118       ACTS_DEBUG(
0119           "Cannot seed because less then three hits with unique (layer, "
0120           "volume)");
0121       skippedTracks++;
0122       continue;
0123     }
0124
0125     // Make the seed
0126     auto result =
0127         track | std::views::filter([&](auto i) {
0128           return i < indexToSpacepoint.size() &&
0129                  indexToSpacepoint.at(i) != nullptr;
0130         }) |
0131         std::views::transform([&](auto i) { return indexToSpacepoint.at(i); });
0132     tmpSps.clear();
0133     std::ranges::copy(result, std::back_inserter(tmpSps));
0134
0135     if (tmpSps.size() < 3) {
0136       ACTS_WARNING("Could not find all spacepoints, skip");
0137       skippedTracks++;
0138       continue;
0139     }
0140
0141     std::ranges::sort(tmpSps, {}, [](const auto &t) { return t->r(); });
0142
0143     // Simply use r = m*z + t and solve for r=0 to find z vertex position...
0144     // Probably not the textbook way to do
0145     const auto m = (tmpSps.back()->r() - tmpSps.front()->r()) /
0146                    (tmpSps.back()->z() - tmpSps.front()->z());
0147     const auto t = tmpSps.front()->r() - m * tmpSps.front()->z();
0148     const auto z_vertex = -t / m;
0149     const auto s = tmpSps.size();
0150
0151     SimSeed seed =
0152         m_cfg.buildTightSeeds
0153             ? SimSeed(*tmpSps.at(0), *tmpSps.at(1), *tmpSps.at(2))
0154             : SimSeed(*tmpSps.at(0), *tmpSps.at(s / 2), *tmpSps.at(s - 1));
0155     seed.setVertexZ(z_vertex);
0156
0157     // Compute parameters
0158     const auto &bottomSP = seed.sp().front();
0159     const auto geoId = bottomSP->sourceLinks()
0160                            .front()
0161                            .template get<IndexSourceLink>()
0162                            .geometryId();
0163     const auto &surface = *m_cfg.geometry->findSurface(geoId);
0164
0165     auto fieldRes = m_cfg.magneticField->getField(
0166         {bottomSP->x(), bottomSP->y(), bottomSP->z()}, bCache);
0167     if (!fieldRes.ok()) {
0168       ACTS_ERROR("Field lookup error: " << fieldRes.error());
0169       return ProcessCode::ABORT;
0170     }
0171     Acts::Vector3 field = *fieldRes;
0172
0173     if (field.norm() < m_cfg.bFieldMin) {
0174       ACTS_WARNING("Magnetic field at seed is too small " << field.norm());
0175       continue;
0176     }
0177
0178     auto parsResult = Acts::estimateTrackParamsFromSeed(
0179         ctx.geoContext, seed.sp(), surface, field);
0180     if (!parsResult.ok()) {
0181       ACTS_WARNING("Skip track because of bad params");
0182     }
0183     const auto &pars = *parsResult;
0184
0185     seededTracks.push_back(track);
0186     seeds.emplace_back(std::move(seed));
0187     parameters.push_back(Acts::BoundTrackParameters(
0188         surface.getSharedPtr(), pars, m_covariance, m_cfg.particleHypothesis));
0189   }
0190
0191   if (skippedTracks > 0) {
0192     ACTS_WARNING("Skipped seeding of " << skippedTracks);
0193   }
0194
0195   ACTS_DEBUG("Seeded " << seeds.size() << " out of " << prototracks.size()
0196                        << " prototracks");
0197
0198   m_outputSeeds(ctx, std::move(seeds));
0199   m_outputProtoTracks(ctx, std::move(seededTracks));
0200   m_outputParameters(ctx, std::move(parameters));
0201
0202   return ProcessCode::SUCCESS;
0203 }
0204
0205 }  // namespace ActsExamples