Acts/TrackFinding/MeasurementSelector.ipp

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 #pragma once
0010
0011 #include "Acts/TrackFinding/MeasurementSelector.hpp"
0012
0013 namespace Acts {
0014
0015 template <typename traj_t>
0016 Result<
0017     std::pair<typename std::vector<typename traj_t::TrackStateProxy>::iterator,
0018               typename std::vector<typename traj_t::TrackStateProxy>::iterator>>
0019 MeasurementSelector::select(
0020     std::vector<typename traj_t::TrackStateProxy>& candidates, bool& isOutlier,
0021     const Logger& logger) const {
0022   using Result = Result<std::pair<
0023       typename std::vector<typename traj_t::TrackStateProxy>::iterator,
0024       typename std::vector<typename traj_t::TrackStateProxy>::iterator>>;
0025
0026   ACTS_VERBOSE("Invoked MeasurementSelector");
0027
0028   // Return if no measurement
0029   if (candidates.empty()) {
0030     return Result::success(std::pair(candidates.begin(), candidates.end()));
0031   }
0032
0033   // Get geoID of this surface
0034   GeometryIdentifier geoID = candidates.front().referenceSurface().geometryId();
0035   // Get the theta of the first track state
0036   const double theta = candidates.front().predicted()[eBoundTheta];
0037   // Find the appropriate cuts
0038   const auto cutsResult = getCuts(geoID, theta);
0039   if (!cutsResult.ok()) {
0040     return cutsResult.error();
0041   }
0042   const Cuts& cuts = *cutsResult;
0043
0044   if (cuts.numMeasurements == 0ul) {
0045     return CombinatorialKalmanFilterError::MeasurementSelectionFailed;
0046   }
0047
0048   double minChi2 = std::numeric_limits<double>::max();
0049   std::size_t minIndex = std::numeric_limits<std::size_t>::max();
0050
0051   isOutlier = false;
0052
0053   // Loop over all measurements to select the compatible measurements
0054   // Sort track states which do not satisfy the chi2 cut to the end.
0055   // When done trackStateIterEnd will point to the first element that
0056   // does not satisfy the chi2 cut.
0057   std::size_t passedCandidates = 0ul;
0058   for (std::size_t i(0ul); i < candidates.size(); ++i) {
0059     auto& trackState = candidates[i];
0060
0061     // We access the dynamic size of the matrix here but use them later
0062     // through a template function which accesses the data pointer
0063     // with compile time size. That way the Eigen math operations are
0064     // still done with compile time size and no dynamic memory allocation
0065     // is needed.
0066     double chi2 = calculateChi2(
0067         trackState.effectiveCalibrated().data(),
0068         trackState.effectiveCalibratedCovariance().data(),
0069         trackState.predicted(), trackState.predictedCovariance(),
0070         trackState.projectorSubspaceIndices(), trackState.calibratedSize());
0071     trackState.chi2() = chi2;
0072
0073     if (chi2 < minChi2) {
0074       minChi2 = chi2;
0075       minIndex = i;
0076     }
0077
0078     // only consider track states which pass the chi2 cut
0079     if (chi2 >= cuts.chi2Measurement) {
0080       continue;
0081     }
0082
0083     if (passedCandidates != i) {
0084       std::swap(candidates[passedCandidates], candidates[i]);
0085       if (minIndex == i) {
0086         minIndex = passedCandidates;
0087       }
0088     }
0089     ++passedCandidates;
0090   }
0091
0092   // Handle if there are no measurements below the chi2 cut off
0093   if (passedCandidates == 0ul) {
0094     if (minChi2 < cuts.chi2Outlier) {
0095       ACTS_VERBOSE(
0096           "No measurement candidate. Return an outlier measurement chi2="
0097           << minChi2);
0098       isOutlier = true;
0099       return Result::success(std::pair(candidates.begin() + minIndex,
0100                                        candidates.begin() + minIndex + 1));
0101     } else {
0102       ACTS_VERBOSE("No measurement candidate. Return empty chi2=" << minChi2);
0103       return Result::success(std::pair(candidates.begin(), candidates.begin()));
0104     }
0105   }
0106
0107   if (passedCandidates <= 1ul) {
0108     // return single item range, no sorting necessary
0109     ACTS_VERBOSE("Returning only 1 element chi2=" << minChi2);
0110     return Result::success(std::pair(candidates.begin() + minIndex,
0111                                      candidates.begin() + minIndex + 1));
0112   }
0113
0114   std::sort(
0115       candidates.begin(), candidates.begin() + passedCandidates,
0116       [](const auto& tsa, const auto& tsb) { return tsa.chi2() < tsb.chi2(); });
0117
0118   ACTS_VERBOSE("Number of selected measurements: "
0119                << passedCandidates << ", max: " << cuts.numMeasurements);
0120
0121   return Result::success(std::pair(
0122       candidates.begin(),
0123       candidates.begin() + std::min(cuts.numMeasurements, passedCandidates)));
0124 }
0125
0126 }  // namespace Acts