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