src/TrackFinding/MeasurementSelector.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 "Acts/TrackFinding/MeasurementSelector.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/EventData/MeasurementHelpers.hpp"
0013 #include "Acts/EventData/SubspaceHelpers.hpp"
0014 #include "Acts/EventData/Types.hpp"
0015 #include "Acts/Geometry/GeometryIdentifier.hpp"
0016
0017 #include <algorithm>
0018 #include <cstddef>
0019 #include <limits>
0020 #include <numbers>
0021 #include <span>
0022 #include <stdexcept>
0023
0024 namespace Acts {
0025
0026 MeasurementSelector::MeasurementSelector()
0027     : MeasurementSelector({{GeometryIdentifier(), MeasurementSelectorCuts{}}}) {
0028 }
0029
0030 MeasurementSelector::MeasurementSelector(const MeasurementSelectorCuts& cuts)
0031     : MeasurementSelector({{GeometryIdentifier(), cuts}}) {}
0032
0033 MeasurementSelector::MeasurementSelector(const Config& config) {
0034   if (config.empty()) {
0035     throw std::invalid_argument(
0036         "MeasurementSelector: Configuration must not be empty");
0037   }
0038
0039   std::vector<InternalConfig::InputElement> tmp;
0040   tmp.reserve(config.size());
0041   for (std::size_t i = 0; i < config.size(); ++i) {
0042     GeometryIdentifier geoID = config.idAt(i);
0043     MeasurementSelectorCuts cuts = config.valueAt(i);
0044
0045     InternalCutBins internalCuts = convertCutBins(cuts);
0046     tmp.emplace_back(geoID, std::move(internalCuts));
0047   }
0048   m_config = InternalConfig(std::move(tmp));
0049 }
0050
0051 MeasurementSelector::InternalCutBins MeasurementSelector::convertCutBins(
0052     const MeasurementSelectorCuts& config) {
0053   InternalCutBins cutBins;
0054
0055   auto getEtaOrInf = [](const auto& vec, std::size_t bin) {
0056     if (bin >= vec.size()) {
0057       return std::numeric_limits<double>::infinity();
0058     }
0059     assert(vec[bin] >= 0 && "Eta bins must be positive");
0060     return vec[bin];
0061   };
0062
0063   auto getBinOrBackOrMax = [](const auto& vec, std::size_t bin) {
0064     using Value = std::remove_reference_t<decltype(vec[0])>;
0065     static constexpr Value max = std::numeric_limits<Value>::max();
0066     return vec.empty() ? max : (bin < vec.size() ? vec[bin] : vec.back());
0067   };
0068
0069   for (std::size_t bin = 0; bin < config.etaBins.size() + 1; ++bin) {
0070     InternalCutBin cuts;
0071     cuts.maxTheta = getEtaOrInf(config.etaBins, bin);
0072     cuts.maxNumMeasurements =
0073         getBinOrBackOrMax(config.numMeasurementsCutOff, bin);
0074     cuts.maxChi2Measurement = getBinOrBackOrMax(config.chi2CutOff, bin);
0075     cuts.maxChi2Outlier = getBinOrBackOrMax(config.chi2CutOffOutlier, bin);
0076     cutBins.push_back(cuts);
0077   }
0078
0079   return cutBins;
0080 }
0081
0082 double MeasurementSelector::calculateChi2(
0083     const double* fullCalibrated, const double* fullCalibratedCovariance,
0084     TrackStateTraits<kMeasurementSizeMax, false>::Parameters predicted,
0085     TrackStateTraits<kMeasurementSizeMax, false>::Covariance
0086         predictedCovariance,
0087     BoundSubspaceIndices projector, unsigned int calibratedSize) const {
0088   return visit_measurement(
0089       calibratedSize,
0090       [&fullCalibrated, &fullCalibratedCovariance, &predicted,
0091        &predictedCovariance, &projector](auto N) -> double {
0092         constexpr std::size_t kMeasurementSize = decltype(N)::value;
0093
0094         typename TrackStateTraits<kMeasurementSize, true>::Calibrated
0095             calibrated{fullCalibrated};
0096
0097         typename TrackStateTraits<kMeasurementSize, true>::CalibratedCovariance
0098             calibratedCovariance{fullCalibratedCovariance};
0099
0100         using ParametersVector = Vector<kMeasurementSize>;
0101
0102         std::span<const std::uint8_t, kMeasurementSize> validSubspaceIndices(
0103             projector.begin(), projector.begin() + kMeasurementSize);
0104         FixedBoundSubspaceHelper<kMeasurementSize> subspaceHelper(
0105             validSubspaceIndices);
0106
0107         // Get the residuals
0108         ParametersVector res =
0109             calibrated - subspaceHelper.projectVector(predicted);
0110
0111         // Get the chi2
0112         return (res.transpose() *
0113                 (calibratedCovariance +
0114                  subspaceHelper.projectMatrix(predictedCovariance))
0115                     .inverse() *
0116                 res)
0117             .eval()(0, 0);
0118       });
0119 }
0120
0121 MeasurementSelector::Cuts MeasurementSelector::getCutsByTheta(
0122     const InternalCutBins& config, double theta) {
0123   // since theta is in [0, pi] and we have a symmetric cut in eta, we can just
0124   // look at the positive half of the Z axis
0125   const double constrainedTheta = std::min(theta, std::numbers::pi - theta);
0126
0127   auto it = std::ranges::find_if(
0128       config, [constrainedTheta](const InternalCutBin& cuts) {
0129         return constrainedTheta < cuts.maxTheta;
0130       });
0131   assert(it != config.end());
0132   return {it->maxNumMeasurements, it->maxChi2Measurement, it->maxChi2Outlier};
0133 }
0134
0135 Result<MeasurementSelector::Cuts> MeasurementSelector::getCuts(
0136     const GeometryIdentifier& geoID, double theta) const {
0137   // Find the appropriate cuts
0138   auto cuts = m_config.find(geoID);
0139   if (cuts == m_config.end()) {
0140     // for now we consider missing cuts an unrecoverable error
0141     // TODO consider other options e.g. do not add measurements at all (not
0142     // even as outliers)
0143     return CombinatorialKalmanFilterError::MeasurementSelectionFailed;
0144   }
0145   return getCutsByTheta(*cuts, theta);
0146 }
0147
0148 }  // namespace Acts