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<MultiTrajectoryTraits::MeasurementSizeMax,
0085                      false>::Parameters predicted,
0086     TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
0087                      false>::Covariance predictedCovariance,
0088     BoundSubspaceIndices projector, unsigned int calibratedSize) const {
0089   return visit_measurement(
0090       calibratedSize,
0091       [&fullCalibrated, &fullCalibratedCovariance, &predicted,
0092        &predictedCovariance, &projector](auto N) -> double {
0093         constexpr std::size_t kMeasurementSize = decltype(N)::value;
0094
0095         typename TrackStateTraits<kMeasurementSize, true>::Calibrated
0096             calibrated{fullCalibrated};
0097
0098         typename TrackStateTraits<kMeasurementSize, true>::CalibratedCovariance
0099             calibratedCovariance{fullCalibratedCovariance};
0100
0101         using ParametersVector = ActsVector<kMeasurementSize>;
0102
0103         std::span<const std::uint8_t, kMeasurementSize> validSubspaceIndices(
0104             projector.begin(), projector.begin() + kMeasurementSize);
0105         FixedBoundSubspaceHelper<kMeasurementSize> subspaceHelper(
0106             validSubspaceIndices);
0107
0108         // Get the residuals
0109         ParametersVector res =
0110             calibrated - subspaceHelper.projectVector(predicted);
0111
0112         // Get the chi2
0113         return (res.transpose() *
0114                 (calibratedCovariance +
0115                  subspaceHelper.projectMatrix(predictedCovariance))
0116                     .inverse() *
0117                 res)
0118             .eval()(0, 0);
0119       });
0120 }
0121
0122 MeasurementSelector::Cuts MeasurementSelector::getCutsByTheta(
0123     const InternalCutBins& config, double theta) {
0124   // since theta is in [0, pi] and we have a symmetric cut in eta, we can just
0125   // look at the positive half of the Z axis
0126   const double constrainedTheta = std::min(theta, std::numbers::pi - theta);
0127
0128   auto it = std::ranges::find_if(
0129       config, [constrainedTheta](const InternalCutBin& cuts) {
0130         return constrainedTheta < cuts.maxTheta;
0131       });
0132   assert(it != config.end());
0133   return {it->maxNumMeasurements, it->maxChi2Measurement, it->maxChi2Outlier};
0134 }
0135
0136 Result<MeasurementSelector::Cuts> MeasurementSelector::getCuts(
0137     const GeometryIdentifier& geoID, double theta) const {
0138   // Find the appropriate cuts
0139   auto cuts = m_config.find(geoID);
0140   if (cuts == m_config.end()) {
0141     // for now we consider missing cuts an unrecoverable error
0142     // TODO consider other options e.g. do not add measurements at all (not
0143     // even as outliers)
0144     return CombinatorialKalmanFilterError::MeasurementSelectionFailed;
0145   }
0146   return getCutsByTheta(*cuts, theta);
0147 }
0148
0149 }  // namespace Acts