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