16 #include "Acts/Utilities/Result.hpp"
17 #include "Acts/Utilities/Delegate.hpp"
18 #include "Acts/EventData/SourceLink.hpp"
19 #include "Acts/TrackFinding/CombinatorialKalmanFilterError.hpp"
20 #include "Acts/Definitions/Algebra.hpp"
21 #include "Acts/Surfaces/Surface.hpp"
22 #include "Acts/Geometry/GeometryHierarchyMap.hpp"
23 #include "Acts/EventData/Types.hpp"
24 #include "Acts/EventData/TrackStatePropMask.hpp"
25 #include "boost/container/small_vector.hpp"
28 #include "Acts/EventData/TrackStateProxy.hpp"
29 #include "Acts/Utilities/CalibrationContext.hpp"
30 #include "Acts/EventData/TrackParameters.hpp"
33 #include "Acts/EventData/MultiTrajectory.hpp"
36 #include <type_traits>
42 template <
typename derived_t>
46 template <std::
size_t N>
50 template <std::
size_t N>
55 template <std::
size_t N>
60 template <std::
size_t N>
64 template <std::
size_t N>
65 using Predicted =
typename Acts::detail_lt::FixedSizeTypes<N>::Coefficients;
68 template <std::
size_t N>
72 template <
typename T_Container>
87 using BoundState = std::tuple<BoundTrackParameters, BoundMatrix, double>;
101 template <
class T_Matrix>
102 static constexpr std::size_t
matrixColumns() {
return T_Matrix::ColsAtCompileTime;}
103 template <
class T_Matrix>
104 static constexpr std::size_t
matrixRows() {
return T_Matrix::RowsAtCompileTime;}
106 template <
typename T_Float,
class T_Matrix>
109 template <
class T_Matrix>
112 template <
class T_Matrix>
126 template <std::
size_t N>
127 using type = std::array<unsigned char, N>;
129 template <std::
size_t N>
132 for(
int i=0;
i<
N; ++
i) {
141 template <std::
size_t N,
class T_ResultType,
class T_Matrix>
144 using MatrixIndexMapType =
unsigned char;
146 using MatrixIndexType =
unsigned int;
153 if constexpr(MeasurementSelectorMatrixTraits::matrixColumns<T_Matrix>() == 1) {
155 for (MatrixIndexType meas_i=0; meas_i<
N; ++meas_i) {
156 assert( meas_i < parameter_map.size() );
157 ret(meas_i,0) =
matrix( parameter_map[meas_i], 0);
163 for (MatrixIndexType meas_i=0; meas_i<
N; ++meas_i) {
164 assert( meas_i < parameter_map.size());
165 MatrixIndexType param_i = parameter_map[meas_i];
166 for (MatrixIndexType meas_j=0; meas_j<
N; ++meas_j) {
167 assert( meas_j < parameter_map.size());
168 ret(meas_i,meas_j) =
matrix(param_i, parameter_map[meas_j]);
179 template <
typename measurement_vector_t,
typename measurement_cov_matrix_t,
180 typename predicted_vector_t,
typename predicted_cov_matrix_t>
182 const measurement_cov_matrix_t &a_cov,
183 const predicted_vector_t &
b,
184 const predicted_cov_matrix_t &b_cov) {
188 static_assert( MeasurementSelectorMatrixTraits::matrixColumns<measurement_vector_t>() == 1);
189 static_assert( MeasurementSelectorMatrixTraits::matrixColumns<predicted_vector_t>() == 1);
190 static_assert( MeasurementSelectorMatrixTraits::matrixRows<measurement_cov_matrix_t>()
191 == MeasurementSelectorMatrixTraits::matrixColumns<measurement_cov_matrix_t>() );
192 static_assert( MeasurementSelectorMatrixTraits::matrixRows<predicted_cov_matrix_t>()
193 == MeasurementSelectorMatrixTraits::matrixColumns<predicted_cov_matrix_t>() );
194 static_assert( MeasurementSelectorMatrixTraits::matrixRows<measurement_cov_matrix_t>()
195 == MeasurementSelectorMatrixTraits::matrixRows<measurement_vector_t>() );
196 static_assert( MeasurementSelectorMatrixTraits::matrixRows<predicted_cov_matrix_t>()
197 == MeasurementSelectorMatrixTraits::matrixRows<predicted_vector_t>() );
198 static_assert( MeasurementSelectorMatrixTraits::matrixRows<measurement_cov_matrix_t>()
199 == MeasurementSelectorMatrixTraits::matrixRows<predicted_cov_matrix_t>() );
211 template <std::
size_t N,
class PayloadType >
250 void acceptAndSort(std::function<
bool(
const PayloadType &
a,
const PayloadType &
b)> comparison) {
271 typename std::array<IndexType, N+1>::const_iterator
begin() {
return m_order.begin(); }
297 template <std::size_t NMeasMax,
302 using Config = Acts::GeometryHierarchyMap<AtlasMeasurementSelectorCuts>;
313 const derived_t &
derived()
const {
return *
static_cast<const derived_t *
>(
this); }
317 template <std::
size_t N>
320 constexpr std::size_t nrows = Acts::MultiTrajectoryTraits::MeasurementSizeMax;
321 constexpr std::size_t ncols = Acts::eBoundSize;
323 std::bitset<nrows * ncols> proj_bitset {};
325 for (
unsigned int col_i=0; col_i<
N; ++col_i) {
326 unsigned int row_i = parameter_map[col_i];
327 unsigned int idx = col_i *nrows + row_i;
328 proj_bitset[ (nrows * ncols - 1) -
idx ] = 1;
330 return proj_bitset.to_ullong();
343 const Acts::BoundSubspaceIndices& subspaceIndices,
344 boost::container::small_vector< typename TrackStateProxy::IndexType, s_maxBranchesPerSurface> &track_states,
345 const Acts::Logger&
logger,
346 bool outlier_states) {
348 const auto &boundParams = derived_t::boundParams(boundState);
349 const auto &pathLength = derived_t::pathLength(boundState);
351 using PM = Acts::TrackStatePropMask;
353 track_states.reserve( n_new_states );
354 std::optional<TrackStateProxy> firstTrackState{
356 for (
unsigned int state_i=0; state_i<n_new_states; ++state_i) {
357 PM
mask = PM::Predicted | PM::Filtered | PM::Jacobian | PM::Calibrated;
359 if (firstTrackState.has_value()) {
362 mask &= ~PM::Predicted & ~PM::Jacobian;
365 if (outlier_states) {
367 mask &= ~PM::Filtered;
371 ACTS_VERBOSE(
"Create SourceLink output track state #"
372 << trackState.index() <<
" with mask: " <<
mask);
374 if (firstTrackState.has_value()) {
375 trackState.shareFrom(*firstTrackState, PM::Predicted);
376 trackState.shareFrom(*firstTrackState, PM::Jacobian);
380 trackState.predicted() = boundParams.parameters();
381 if (boundParams.covariance()) {
382 trackState.predictedCovariance() = *boundParams.covariance();
384 trackState.jacobian() = derived_t::boundJacobiMatrix(boundState);
385 firstTrackState = trackState;
387 trackState.pathLength() = pathLength;
389 trackState.setReferenceSurface(boundParams.referenceSurface().getSharedPtr());
391 trackState.setProjectorSubspaceIndices(subspaceIndices);
394 if (trackState.referenceSurface().surfaceMaterial() !=
nullptr) {
395 typeFlags.set(Acts::TrackStateFlag::MaterialFlag);
397 typeFlags.set(Acts::TrackStateFlag::ParameterFlag);
401 track_states.push_back( trackState.index());
407 template <std::
size_t DIM,
typename T_SourceLink>
424 template <
typename T>
432 template <
typename Iterator>
447 template <std::
size_t DIM,
typename source_link_iterator_t,
typename T_Container>
448 Acts::Result<boost::container::small_vector< typename TrackStateProxy::IndexType, s_maxBranchesPerSurface> >
450 const Acts::CalibrationContext& calibrationContext,
451 const Acts::Surface& surface,
453 const source_link_iterator_t& sourceLinkBegin,
454 const source_link_iterator_t& sourceLinkEnd,
457 const Acts::Logger&
logger,
458 const std::size_t numMeasurementsCut,
459 const std::pair<float,float>& maxChi2Cut,
460 const T_Container &container)
const {
461 Acts::Result<boost::container::small_vector< typename TrackStateProxy::IndexType, s_maxBranchesPerSurface> >
462 result = boost::container::small_vector< typename TrackStateProxy::IndexType, s_maxBranchesPerSurface>{};
464 using BaseElementType = std::remove_cv_t<std::remove_pointer_t< container_value_t > >;
466 using TheMatchingMeasurement = MatchingMeasurement<DIM, container_value_t >;
467 using MeasCovPair = TheMatchingMeasurement::MeasCovPair;
475 Acts::SubspaceIndices<DIM> parameter_map =
derived().template parameterMap<DIM>(geometryContext,
479 = std::make_pair( project<DIM, Predicted>(parameter_map,
481 project<DIM, PredictedCovariance>(parameter_map,
482 derived().boundParams(boundState).covariance().
value()));
485 auto preCalibrator =
derived().template preCalibrator<DIM, BaseElementType>();
488 for (
const auto &measurement : MeasurementRange<T_Container>(container, sourceLinkBegin, sourceLinkEnd) ) {
489 TheMatchingMeasurement &matching_measurement=selected_measurements.
slot();
490 matching_measurement.m_measurement = preCalibrator(geometryContext,
492 derived().
template forwardToCalibrator(measurement),
493 derived().boundParams(boundState));
494 matching_measurement.m_chi2 =
computeChi2(matching_measurement.m_measurement.first,
495 matching_measurement.m_measurement.second,
499 if (matching_measurement.m_chi2<maxChi2Cut.second) {
500 matching_measurement.m_sourceLink=measurement;
501 selected_measurements.
acceptAndSort([](
const TheMatchingMeasurement &
a,
502 const TheMatchingMeasurement &
b) {
503 return a.m_chi2 <
b.m_chi2;
510 auto postCalibrator =
derived().template postCalibrator<DIM, BaseElementType>();
511 using post_calib_meas_cov_pair_t
517 static constexpr
bool pre_and_post_calib_types_agree
523 using Empty =
struct {};
524 typename std::conditional< !pre_and_post_calib_types_agree,
525 std::array< post_calib_meas_cov_pair_t, NMeasMax>,
527 if (postCalibrator) {
532 idx: selected_measurements) {
533 TheMatchingMeasurement &a_selected_measurement = selected_measurements.
getSlot(
idx);
537 post_calib_meas_cov_pair_t &calibrated_measurement
538 = [&calibrated, &a_selected_measurement, calibrated_meas_cov_i]() -> post_calib_meas_cov_pair_t & {
539 if constexpr(pre_and_post_calib_types_agree) {
541 (void) calibrated_meas_cov_i;
542 return a_selected_measurement.m_measurement;
545 (void) a_selected_measurement;
546 assert(calibrated_meas_cov_i < calibrated.size());
547 return calibrated[calibrated_meas_cov_i];
552 calibrated_measurement = postCalibrator(geometryContext,
554 derived().
template forwardToCalibrator(a_selected_measurement.m_sourceLink.value()),
555 derived().boundParams(boundState));
557 a_selected_measurement.m_chi2 =
computeChi2(calibrated_measurement.first,
558 calibrated_measurement.second,
562 a_selected_measurement.m_isOutLier = (a_selected_measurement.m_chi2 >= maxChi2Cut.first);
563 if constexpr(!pre_and_post_calib_types_agree) {
564 ++calibrated_meas_cov_i;
571 idx: selected_measurements) {
572 TheMatchingMeasurement &a_selected_measurement = selected_measurements.
getSlot(
idx);
573 a_selected_measurement.m_isOutLier = (a_selected_measurement.m_chi2 >= maxChi2Cut.first);
583 Acts::BoundSubspaceIndices boundSubspaceIndices;
584 std::copy(parameter_map.begin(), parameter_map.end(), boundSubspaceIndices.begin());
589 boundSubspaceIndices,
592 (!selected_measurements.
empty()
593 ? selected_measurements.
getSlot( *(selected_measurements.
begin())).m_isOutLier
595 assert(
result->size() == selected_measurements.
size() );
598 auto use_calibrated_storage = [&postCalibrator]() ->
bool {
599 if constexpr(pre_and_post_calib_types_agree) {
600 (void) postCalibrator;
606 return postCalibrator;
611 unsigned int state_i=0;
613 idx: selected_measurements) {
614 assert( state_i < result->
size());
616 TheMatchingMeasurement &a_selected_measurement = selected_measurements.
getSlot(
idx);
617 trackState.setUncalibratedSourceLink(
derived().makeSourceLink(std::move(a_selected_measurement.m_sourceLink.value())));
619 trackState.typeFlags().set( a_selected_measurement.m_isOutLier
620 ? Acts::TrackStateFlag::OutlierFlag
621 : Acts::TrackStateFlag::MeasurementFlag );
622 trackState.allocateCalibrated(DIM);
623 if (use_calibrated_storage()) {
626 assert( use_calibrated_storage() == !pre_and_post_calib_types_agree);
627 if constexpr(!pre_and_post_calib_types_agree) {
628 assert( state_i < calibrated.size());
629 trackState.template calibrated<DIM>()
630 = MeasurementSelectorMatrixTraits::matrixTypeCast<typename MeasurementSelectorTraits<derived_t>::MatrixFloatType>(calibrated[state_i].first);
631 trackState.template calibratedCovariance<DIM>()
632 = MeasurementSelectorMatrixTraits::matrixTypeCast<typename MeasurementSelectorTraits<derived_t>::MatrixFloatType>(calibrated[state_i].second);
633 trackState.chi2() = a_selected_measurement.m_chi2;
637 trackState.template calibrated<DIM>()
638 = MeasurementSelectorMatrixTraits::matrixTypeCast<typename MeasurementSelectorTraits<derived_t>::MatrixFloatType>(a_selected_measurement.m_measurement.first);
639 trackState.template calibratedCovariance<DIM>()
640 = MeasurementSelectorMatrixTraits::matrixTypeCast<typename MeasurementSelectorTraits<derived_t>::MatrixFloatType>(a_selected_measurement.m_measurement.second);
641 trackState.chi2() = a_selected_measurement.m_chi2;
648 template <
typename parameters_t>
649 static std::size_t
getEtaBin(
const parameters_t& boundParameters,
650 const std::vector<float> &
etaBins) {
654 const float eta = std::abs(std::atanh(
std::cos(boundParameters.parameters()[Acts::eBoundTheta])));
668 std::tuple< std::size_t, std::pair<float,float> >
getCuts(
const Acts::Surface& surface,
670 const Acts::Logger&
logger)
const {
671 std::tuple< std::size_t, std::pair<float,float> >
result;
672 std::size_t &numMeasurementsCut = std::get<0>(
result);
673 std::pair<float,float> &maxChi2Cut = std::get<1>(
result);
675 auto geoID = surface.geometryId();
680 numMeasurementsCut = 0;
686 numMeasurementsCut = (!
cuts->numMeasurementsCutOff.empty()
689 maxChi2Cut = !
cuts->chi2CutOff.empty()
693 ACTS_VERBOSE(
"Get cut for eta-bin="
695 <<
": chi2 (max,max-outlier) " << maxChi2Cut.first <<
", " << maxChi2Cut.second
696 <<
" max.meas." << numMeasurementsCut);
712 template <std::size_t NMeasMax,
714 typename measurement_container_variant_t >
728 template <
typename source_link_iterator_t>
729 Acts::Result<boost::container::small_vector< typename TrackStateProxy::IndexType, s_maxBranchesPerSurface> >
731 const Acts::CalibrationContext& calibrationContext,
732 const Acts::Surface& surface,
734 source_link_iterator_t sourceLinkBegin,
735 source_link_iterator_t sourceLinkEnd,
738 [[maybe_unused]] std::vector<TrackStateProxy> &trackStateCandidates,
740 const Acts::Logger&
logger)
const {
741 Acts::Result<boost::container::small_vector< typename TrackStateProxy::IndexType, s_maxBranchesPerSurface> >
742 result = Acts::CombinatorialKalmanFilterError::MeasurementSelectionFailed;
743 if (sourceLinkBegin != sourceLinkEnd) {
744 auto [numMeasurementsCut, maxChi2Cut] = this->
getCuts(surface,boundState,
logger);
747 if (numMeasurementsCut>0) {
748 const std::vector<measurement_container_variant_t> &
749 measurementContainer = sourceLinkBegin.m_iterator.measurementContainerList();
751 assert( sourceLinkBegin.m_iterator.containerIndex() == sourceLinkEnd.m_iterator.containerIndex() );
752 const measurement_container_variant_t &a_measurement_container_variant = measurementContainer.at(sourceLinkBegin.m_iterator.containerIndex());
753 result = std::visit( [
this,
764 &maxChi2Cut] (
const auto &measurement_container_with_dimension) {
765 using ArgType = std::remove_cv_t<std::remove_reference_t< decltype(measurement_container_with_dimension) > >;
766 constexpr std::size_t DIM = ArgType::dimension();
767 return this->
template selectMeasurementsCreateTrackStates<DIM>(geometryContext,
778 measurement_container_with_dimension.container());
780 a_measurement_container_variant);
787 template <std::size_t NMeasMax,
789 typename measurement_container_variant_t>
795 return std::get<0>(boundState);
800 return std::get<1>(boundState);
804 return std::get<2>(boundState);
808 template <
typename T_Value>
810 return Acts::SourceLink{
value};
816 template <
typename T>
818 if constexpr( std::is_same<
T, std::remove_pointer_t<T> >::
value ) {
830 template <std::
size_t DIM>
833 const Acts::CalibrationContext&,
834 const Acts::Surface&) {
835 return ParameterMapping::identity<DIM>();
839 template <std::
size_t DIM,
typename measurement_t>
843 (
const Acts::GeometryContext&,
844 const Acts::CalibrationContext&,
845 const measurement_t &,
850 template <std::
size_t DIM,
typename measurement_t>
854 (
const Acts::GeometryContext&,
855 const Acts::CalibrationContext&,
856 const measurement_t &,
862 template <std::
size_t DIM,
typename measurement_t>
870 template <std::
size_t DIM,
typename measurement_t>