TrackFindingExaTrkX/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/TrackFindingExaTrkX/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. 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::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     tmpSps.clear();
0122     std::transform(track.begin(), track.end(), std::back_inserter(tmpSps),
0123                    [&](auto i) { return indexToSpacepoint[i]; });
0124     tmpSps.erase(std::remove_if(tmpSps.begin(), tmpSps.end(),
0125                                 [](auto sp) { return sp == nullptr; }),
0126                  tmpSps.end());
0127
0128     if (tmpSps.size() < 3) {
0129       ACTS_WARNING("Could not find all spacepoints, skip");
0130       skippedTracks++;
0131       continue;
0132     }
0133
0134     std::ranges::sort(tmpSps, {}, [](const auto &t) { return t->r(); });
0135
0136     // Simply use r = m*z + t and solve for r=0 to find z vertex position...
0137     // Probably not the textbook way to do
0138     const auto m = (tmpSps.back()->r() - tmpSps.front()->r()) /
0139                    (tmpSps.back()->z() - tmpSps.front()->z());
0140     const auto t = tmpSps.front()->r() - m * tmpSps.front()->z();
0141     const auto z_vertex = -t / m;
0142     const auto s = tmpSps.size();
0143
0144     SimSeed seed = m_cfg.buildTightSeeds
0145                        ? SimSeed(*tmpSps[0], *tmpSps[1], *tmpSps[2])
0146                        : SimSeed(*tmpSps[0], *tmpSps[s / 2], *tmpSps[s - 1]);
0147     seed.setVertexZ(z_vertex);
0148
0149     // Compute parameters
0150     const auto &bottomSP = seed.sp().front();
0151     const auto geoId = bottomSP->sourceLinks()
0152                            .front()
0153                            .template get<IndexSourceLink>()
0154                            .geometryId();
0155     const auto &surface = *m_cfg.geometry->findSurface(geoId);
0156
0157     auto fieldRes = m_cfg.magneticField->getField(
0158         {bottomSP->x(), bottomSP->y(), bottomSP->z()}, bCache);
0159     if (!fieldRes.ok()) {
0160       ACTS_ERROR("Field lookup error: " << fieldRes.error());
0161       return ProcessCode::ABORT;
0162     }
0163     Acts::Vector3 field = *fieldRes;
0164
0165     if (field.norm() < m_cfg.bFieldMin) {
0166       ACTS_WARNING("Magnetic field at seed is too small " << field.norm());
0167       continue;
0168     }
0169
0170     auto parsResult = Acts::estimateTrackParamsFromSeed(
0171         ctx.geoContext, seed.sp(), surface, field);
0172     if (!parsResult.ok()) {
0173       ACTS_WARNING("Skip track because of bad params");
0174     }
0175     const auto &pars = *parsResult;
0176
0177     seededTracks.push_back(track);
0178     seeds.emplace_back(std::move(seed));
0179     parameters.push_back(Acts::BoundTrackParameters(
0180         surface.getSharedPtr(), pars, m_covariance, m_cfg.particleHypothesis));
0181   }
0182
0183   if (skippedTracks > 0) {
0184     ACTS_WARNING("Skipped seeding of " << skippedTracks);
0185   }
0186
0187   ACTS_DEBUG("Seeded " << seeds.size() << " out of " << prototracks.size()
0188                        << " prototracks");
0189
0190   m_outputSeeds(ctx, std::move(seeds));
0191   m_outputProtoTracks(ctx, std::move(seededTracks));
0192   m_outputParameters(ctx, std::move(parameters));
0193
0194   return ProcessCode::SUCCESS;
0195 }
0196
0197 }  // namespace ActsExamples