TrackFinding/src/SpacePointMaker.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020-2022 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/TrackFinding/SpacePointMaker.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/EventData/SourceLink.hpp"
0014 #include "Acts/SpacePointFormation/SpacePointBuilderConfig.hpp"
0015 #include "Acts/SpacePointFormation/SpacePointBuilderOptions.hpp"
0016 #include "ActsExamples/EventData/GeometryContainers.hpp"
0017 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0018 #include "ActsExamples/EventData/Measurement.hpp"
0019 #include "ActsExamples/EventData/SimSpacePoint.hpp"
0020 #include "ActsExamples/Framework/AlgorithmContext.hpp"
0021 #include "ActsExamples/Utilities/GroupBy.hpp"
0022 #include "ActsExamples/Utilities/Range.hpp"
0023
0024 #include <algorithm>
0025 #include <functional>
0026 #include <iterator>
0027 #include <ostream>
0028 #include <stdexcept>
0029 #include <utility>
0030 #include <variant>
0031
0032 ActsExamples::SpacePointMaker::SpacePointMaker(Config cfg,
0033                                                Acts::Logging::Level lvl)
0034     : IAlgorithm("SpacePointMaker", lvl), m_cfg(std::move(cfg)) {
0035   if (m_cfg.inputSourceLinks.empty()) {
0036     throw std::invalid_argument("Missing source link input collection");
0037   }
0038   if (m_cfg.inputMeasurements.empty()) {
0039     throw std::invalid_argument("Missing measurement input collection");
0040   }
0041   if (m_cfg.outputSpacePoints.empty()) {
0042     throw std::invalid_argument("Missing space point output collection");
0043   }
0044   if (!m_cfg.trackingGeometry) {
0045     throw std::invalid_argument("Missing tracking geometry");
0046   }
0047   if (m_cfg.geometrySelection.empty()) {
0048     throw std::invalid_argument("Missing space point maker geometry selection");
0049   }
0050
0051   m_inputSourceLinks.initialize(m_cfg.inputSourceLinks);
0052   m_inputMeasurements.initialize(m_cfg.inputMeasurements);
0053   m_outputSpacePoints.initialize(m_cfg.outputSpacePoints);
0054
0055   // ensure geometry selection contains only valid inputs
0056   for (const auto& geoId : m_cfg.geometrySelection) {
0057     if ((geoId.approach() != 0u) || (geoId.boundary() != 0u) ||
0058         (geoId.sensitive() != 0u)) {
0059       throw std::invalid_argument(
0060           "Invalid geometry selection: only volume and layer are allowed to be "
0061           "set");
0062     }
0063   }
0064   // remove geometry selection duplicates
0065   //
0066   // the geometry selections must be mutually exclusive, i.e. if we have a
0067   // selection that contains both a volume and a layer within that same volume,
0068   // we would create the space points for the layer twice.
0069   auto isDuplicate = [](Acts::GeometryIdentifier ref,
0070                         Acts::GeometryIdentifier cmp) {
0071     // code assumes ref < cmp and that only volume and layer can be non-zero
0072     // root node always contains everything
0073     if (ref.volume() == 0) {
0074       return true;
0075     }
0076     // unequal volumes always means separate hierarchies
0077     if (ref.volume() != cmp.volume()) {
0078       return false;
0079     }
0080     // within the same volume hierarchy only consider layers
0081     return (ref.layer() == cmp.layer());
0082   };
0083   auto geoSelBeg = m_cfg.geometrySelection.begin();
0084   auto geoSelEnd = m_cfg.geometrySelection.end();
0085   // sort geometry selection so the unique filtering works
0086   std::sort(geoSelBeg, geoSelEnd);
0087   auto geoSelLastUnique = std::unique(geoSelBeg, geoSelEnd, isDuplicate);
0088   if (geoSelLastUnique != geoSelEnd) {
0089     ACTS_WARNING("Removed " << std::distance(geoSelLastUnique, geoSelEnd)
0090                             << " geometry selection duplicates");
0091     m_cfg.geometrySelection.erase(geoSelLastUnique, geoSelEnd);
0092   }
0093   ACTS_INFO("Space point geometry selection:");
0094   for (const auto& geoId : m_cfg.geometrySelection) {
0095     ACTS_INFO("  " << geoId);
0096   }
0097   auto spBuilderConfig = Acts::SpacePointBuilderConfig();
0098   spBuilderConfig.trackingGeometry = m_cfg.trackingGeometry;
0099
0100   m_slSurfaceAccessor.emplace(
0101       IndexSourceLink::SurfaceAccessor{*m_cfg.trackingGeometry});
0102   spBuilderConfig.slSurfaceAccessor
0103       .connect<&IndexSourceLink::SurfaceAccessor::operator()>(
0104           &m_slSurfaceAccessor.value());
0105
0106   auto spConstructor =
0107       [](const Acts::Vector3& pos, std::optional<double> t,
0108          const Acts::Vector2& cov, std::optional<double> varT,
0109          boost::container::static_vector<Acts::SourceLink, 2> slinks)
0110       -> SimSpacePoint {
0111     return SimSpacePoint(pos, t, cov[0], cov[1], varT, std::move(slinks));
0112   };
0113
0114   m_spacePointBuilder = Acts::SpacePointBuilder<SimSpacePoint>(
0115       spBuilderConfig, spConstructor,
0116       Acts::getDefaultLogger("SpacePointBuilder", lvl));
0117 }
0118
0119 ActsExamples::ProcessCode ActsExamples::SpacePointMaker::execute(
0120     const AlgorithmContext& ctx) const {
0121   const auto& sourceLinks = m_inputSourceLinks(ctx);
0122   const auto& measurements = m_inputMeasurements(ctx);
0123
0124   // TODO Support strip measurements
0125   Acts::SpacePointBuilderOptions spOpt;
0126
0127   spOpt.paramCovAccessor = [&measurements](Acts::SourceLink slink) {
0128     const auto islink = slink.get<IndexSourceLink>();
0129     const auto& meas = measurements[islink.index()];
0130
0131     return std::visit(
0132         [](const auto& measurement) {
0133           auto expander = measurement.expander();
0134           Acts::BoundVector par = expander * measurement.parameters();
0135           Acts::BoundSquareMatrix cov =
0136               expander * measurement.covariance() * expander.transpose();
0137           return std::make_pair(par, cov);
0138         },
0139         meas);
0140   };
0141
0142   SimSpacePointContainer spacePoints;
0143   for (Acts::GeometryIdentifier geoId : m_cfg.geometrySelection) {
0144     // select volume/layer depending on what is set in the geometry id
0145     auto range = selectLowestNonZeroGeometryObject(sourceLinks, geoId);
0146     // groupByModule only works with geometry containers, not with an
0147     // arbitrary range. do the equivalent grouping manually
0148     auto groupedByModule = makeGroupBy(range, detail::GeometryIdGetter());
0149
0150     for (auto [moduleGeoId, moduleSourceLinks] : groupedByModule) {
0151       for (auto& sourceLink : moduleSourceLinks) {
0152         m_spacePointBuilder.buildSpacePoint(
0153             ctx.geoContext, {Acts::SourceLink{sourceLink}}, spOpt,
0154             std::back_inserter(spacePoints));
0155       }
0156     }
0157   }
0158
0159   spacePoints.shrink_to_fit();
0160
0161   ACTS_DEBUG("Created " << spacePoints.size() << " space points");
0162   m_outputSpacePoints(ctx, std::move(spacePoints));
0163
0164   return ActsExamples::ProcessCode::SUCCESS;
0165 }