diff --git a/Core/include/Acts/EventData/MultiTrajectory.hpp b/Core/include/Acts/EventData/MultiTrajectory.hpp index d0f366946273d9c203fa2e49f41d1ec3216a4a40..f290c1029191bea386a1bedce8dda194eaf124f1 100644 --- a/Core/include/Acts/EventData/MultiTrajectory.hpp +++ b/Core/include/Acts/EventData/MultiTrajectory.hpp @@ -108,6 +108,7 @@ struct IndexData { double chi2; double pathLength; + TrackStateType typeFlags; IndexType iuncalibrated = kInvalid; IndexType icalibrated = kInvalid; @@ -465,6 +466,19 @@ class TrackStateProxy { /// @return The path length of this track state double pathLength() const { return data().pathLength; } + /// Getter for the type flags associated with the track state. + /// This overloaded is only enabled if not read-only, and returns a mutable + /// reference. + /// @return reference to the type flags. + template <bool RO = ReadOnly, typename = std::enable_if_t<!RO>> + TrackStateType& typeFlags() { + return data().typeFlags; + } + + /// Getter for the type flags. Returns a copy of the type flags value. + /// @return The type flags of this track state + TrackStateType typeFlags() const { return data().typeFlags; } + private: // Private since it can only be created by the trajectory. TrackStateProxy(ConstIf<MultiTrajectory<SourceLink>, ReadOnly>& trajectory, diff --git a/Core/include/Acts/EventData/MultiTrajectory.ipp b/Core/include/Acts/EventData/MultiTrajectory.ipp index dec9dcd12b697af522b522856d5edcfa3888a70c..3d68920c061b1ffb4c18348754b678af0358d514 100644 --- a/Core/include/Acts/EventData/MultiTrajectory.ipp +++ b/Core/include/Acts/EventData/MultiTrajectory.ipp @@ -222,6 +222,7 @@ inline size_t MultiTrajectory<SL>::addTrackState( nts.chi2() = ts.parameter.chi2; nts.pathLength() = ts.parameter.pathLength; + nts.typeFlags() = ts.type(); return index; } diff --git a/Core/include/Acts/EventData/TrackState.hpp b/Core/include/Acts/EventData/TrackState.hpp index 3df7097951858e3109e3c8b4f87d8b328ceba473..d9bd0b9e0bb306a5dd7bf4bfe0b71ddf3a1d538c 100644 --- a/Core/include/Acts/EventData/TrackState.hpp +++ b/Core/include/Acts/EventData/TrackState.hpp @@ -17,6 +17,20 @@ namespace Acts { +/// @enum TrackStateFlag +/// +/// This enum describes the type of TrackState +enum TrackStateFlag { + MeasurementFlag = 0, + ParameterFlag = 1, + OutlierFlag = 2, + HoleFlag = 3, + MaterialFlag = 4, + NumTrackStateFlags = 5 +}; + +using TrackStateType = std::bitset<TrackStateFlag::NumTrackStateFlags>; + class Surface; /// @class TrackState @@ -44,6 +58,7 @@ class TrackState { /// @param m The measurement object TrackState(SourceLink m) : m_surface(&m.referenceSurface()) { measurement.uncalibrated = std::move(m); + m_typeFlags.set(MeasurementFlag); } /// Constructor from parameters @@ -53,6 +68,7 @@ class TrackState { TrackState(parameters_t p) { m_surface = &p.referenceSurface(); parameter.predicted = std::move(p); + m_typeFlags.set(ParameterFlag); } /// Virtual destructor @@ -64,7 +80,8 @@ class TrackState { TrackState(const TrackState& rhs) : parameter(rhs.parameter), measurement(rhs.measurement), - m_surface(rhs.m_surface) {} + m_surface(rhs.m_surface), + m_typeFlags(rhs.m_typeFlags) {} /// Copy move constructor /// @@ -72,7 +89,8 @@ class TrackState { TrackState(TrackState&& rhs) : parameter(std::move(rhs.parameter)), measurement(std::move(rhs.measurement)), - m_surface(std::move(rhs.m_surface)) {} + m_surface(std::move(rhs.m_surface)), + m_typeFlags(std::move(rhs.m_typeFlags)) {} /// Assignment operator /// @@ -81,6 +99,7 @@ class TrackState { parameter = rhs.parameter; measurement = rhs.measurement; m_surface = rhs.m_surface; + m_typeFlags = rhs.m_typeFlags; return (*this); } @@ -91,12 +110,27 @@ class TrackState { parameter = std::move(rhs.parameter); measurement = std::move(rhs.measurement); m_surface = std::move(rhs.m_surface); + m_typeFlags = std::move(rhs.m_typeFlags); return (*this); } /// @brief return method for the surface const Surface& referenceSurface() const { return (*m_surface); } + /// @brief set the type flag + void setType(const TrackStateFlag& flag, bool status = true) { + m_typeFlags.set(flag, status); + } + + /// @brief test if the tracks state is flagged as a given type + bool isType(const TrackStateFlag& flag) const { + assert(flag < NumTrackStateFlags); + return m_typeFlags.test(flag); + } + + /// @brief return method for the type flags + TrackStateType type() const { return m_typeFlags; } + /// @brief number of Measured parameters, forwarded /// @note This only returns a value if there is a calibrated measurement /// set. If not, this returns boost::none @@ -141,5 +175,7 @@ class TrackState { private: /// The surface of this TrackState const Surface* m_surface = nullptr; + /// The type flag of this TrackState + TrackStateType m_typeFlags; }; } // namespace Acts diff --git a/Core/include/Acts/EventData/detail/covariance_helper.hpp b/Core/include/Acts/EventData/detail/covariance_helper.hpp new file mode 100644 index 0000000000000000000000000000000000000000..87e58c3dd83f1b600f8a22f4e0afe8e28911208d --- /dev/null +++ b/Core/include/Acts/EventData/detail/covariance_helper.hpp @@ -0,0 +1,72 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2019 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once +#include <type_traits> +#include "Acts/Utilities/detail/DefaultParameterDefinitions.hpp" + +namespace Acts { +namespace detail { +/// @brief check and correct covariance matrix +/// +/// @tparam CovMatrix_t The type of covariance matrix +/// @tparam NumIter The number of iterations to run the correction +/// +/// Invocation: +/// - covariance_helper<CovMatrix_t, numIter>::validate(covariance) +/// The 'covariance' is checked against semi-positivedefiniteness +/// and limited number of iterations to replace it with the +/// closest semi-positivedefinite are made if it's not +/// +/// @return The (corrected) covariance is semi-positivedefinite or not +template <typename CovMatrix_t, signed int NumIter = 1> +struct covariance_helper { + /// check if the covariance is semi-positive and correction is attempted + static bool validate(CovMatrix_t& covariance) { + if (covariance.hasNaN()) { + return false; + } + size_t nIteration = 0; + while (nIteration < NumIter) { + if (isSemiPositive(covariance)) { + return true; + } else { + Eigen::JacobiSVD<CovMatrix_t> svdCov( + covariance, Eigen::ComputeFullU | Eigen::ComputeFullV); + CovMatrix_t S = svdCov.singularValues().asDiagonal(); + CovMatrix_t V = svdCov.matrixV(); + CovMatrix_t H = V * S * V.transpose(); + covariance = (covariance + H) / 2; + nIteration++; + } + } + /// check again after the iterations + return isSemiPositive(covariance); + } + + /// check if the covariance is semi-positive + static bool isSemiPositive(const CovMatrix_t& covariance) { + if (covariance.hasNaN()) { + return false; + } + Eigen::LDLT<CovMatrix_t> ldltCov(covariance); + return ldltCov.isPositive(); + } + + /// check if the covariance is positive + static bool isPositive(const CovMatrix_t& covariance) { + if (covariance.hasNaN()) { + return false; + } + Eigen::LLT<CovMatrix_t> lltCov(covariance); + return lltCov.info() == Eigen::Success ? true : false; + } +}; + +} // namespace detail +} // namespace Acts diff --git a/Core/include/Acts/Fitter/GainMatrixSmoother.hpp b/Core/include/Acts/Fitter/GainMatrixSmoother.hpp index 16398be3495259d0fa1a58d8067673b7364e4d27..d672326b02d4fc12dc30e05f8356eec4a684480a 100644 --- a/Core/include/Acts/Fitter/GainMatrixSmoother.hpp +++ b/Core/include/Acts/Fitter/GainMatrixSmoother.hpp @@ -12,6 +12,7 @@ #include <memory> #include "Acts/EventData/MultiTrajectory.hpp" #include "Acts/EventData/TrackParameters.hpp" +#include "Acts/EventData/detail/covariance_helper.hpp" #include "Acts/Fitter/KalmanFitterError.hpp" #include "Acts/Utilities/Logger.hpp" #include "Acts/Utilities/Result.hpp" @@ -117,6 +118,18 @@ class GainMatrixSmoother { G * (prev_ts.predictedCovariance() - prev_ts.smoothedCovariance()) * G.transpose(); + // Check if the covariance matrix is semi-positive definite. + // If not, make one (could do more) attempt to replace it with the + // nearest semi-positive def matrix, + // but it could still be non semi-positive + CovMatrix_t smoothedCov = ts.smoothedCovariance(); + if (not detail::covariance_helper<CovMatrix_t>::validate(smoothedCov)) { + ACTS_DEBUG( + "Smoothed covariance is not positive definite. Could result in " + "negative covariance!"); + } + // Reset smoothed covariance + ts.smoothedCovariance() = smoothedCov; ACTS_VERBOSE("Smoothed covariance is: \n" << ts.smoothedCovariance()); prev_ts = ts; diff --git a/Core/include/Acts/Fitter/KalmanFitter.hpp b/Core/include/Acts/Fitter/KalmanFitter.hpp index 59ed4e52688bd8ce29faa6c368d23bcdf15b08f2..b7f5c0e6110ce0a2dae5480011402e2a6c54609e 100644 --- a/Core/include/Acts/Fitter/KalmanFitter.hpp +++ b/Core/include/Acts/Fitter/KalmanFitter.hpp @@ -18,11 +18,13 @@ #include "Acts/Fitter/detail/VoidKalmanComponents.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/MagneticField/MagneticFieldContext.hpp" +#include "Acts/Material/MaterialProperties.hpp" #include "Acts/Propagator/AbortList.hpp" #include "Acts/Propagator/ActionList.hpp" #include "Acts/Propagator/ConstrainedStep.hpp" #include "Acts/Propagator/DirectNavigator.hpp" #include "Acts/Propagator/Propagator.hpp" +#include "Acts/Propagator/detail/PointwiseMaterialInteraction.hpp" #include "Acts/Propagator/detail/StandardAborters.hpp" #include "Acts/Utilities/CalibrationContext.hpp" #include "Acts/Utilities/Definitions.hpp" @@ -54,11 +56,14 @@ struct KalmanFitterOptions { KalmanFitterOptions(std::reference_wrapper<const GeometryContext> gctx, std::reference_wrapper<const MagneticFieldContext> mctx, std::reference_wrapper<const CalibrationContext> cctx, - const Surface* rSurface = nullptr) + const Surface* rSurface = nullptr, + bool mScattering = true, bool eLoss = true) : geoContext(gctx), magFieldContext(mctx), calibrationContext(cctx), - referenceSurface(rSurface) {} + referenceSurface(rSurface), + multipleScattering(mScattering), + energyLoss(eLoss) {} /// Context object for the geometry std::reference_wrapper<const GeometryContext> geoContext; @@ -69,6 +74,12 @@ struct KalmanFitterOptions { /// The reference Surface const Surface* referenceSurface = nullptr; + + /// Whether to consider multiple scattering. + bool multipleScattering = true; + + /// Whether to consider energy loss. + bool energyLoss = true; }; template <typename source_link_t> @@ -84,6 +95,9 @@ struct KalmanFitterResult { // The optional Parameters at the provided surface boost::optional<BoundParameters> fittedParameters; + // Counter for states with measurements + size_t measurementStates = 0; + // Counter for handled states size_t processedStates = 0; @@ -210,6 +224,12 @@ class KalmanFitter { /// Allows retrieving measurements for a surface std::map<const Surface*, source_link_t> inputMeasurements; + /// Whether to consider multiple scattering. + bool multipleScattering = true; + + /// Whether to consider energy loss. + bool energyLoss = true; + /// @brief Kalman actor operation /// /// @tparam propagator_state_t is the type of Propagagor state @@ -233,9 +253,10 @@ class KalmanFitter { } // Update: - // - Waiting for a current surface that appears in the measurement list + // - Waiting for a current surface that has material + // -> a trackState will be created on surface with material auto surface = state.navigation.currentSurface; - if (surface and not result.smoothed) { + if ((surface and surface->surfaceMaterial()) and not result.smoothed) { // Check if the surface is in the measurement map // -> Get the measurement / calibrate // -> Create the predicted state @@ -252,9 +273,9 @@ class KalmanFitter { // Finalization: // When all track states have been handled or the navigation is breaked - if ((result.processedStates == inputMeasurements.size() or - (state.navigation.navigationBreak and - result.processedStates > 0)) and + if ((result.measurementStates == inputMeasurements.size() or + (result.measurementStates > 0 and + state.navigation.navigationBreak)) and not result.smoothed) { // -> Sort the track states (as now the path length is set) // -> Call the smoothing @@ -311,6 +332,14 @@ class KalmanFitter { // Screen output message ACTS_VERBOSE("Measurement surface " << surface->geoID() << " detected."); + + // Update state and stepper with pre material effects + materialInteractor(surface, state, stepper, preUpdate); + + // Transport & bind the state to the current surface + auto [boundParams, jacobian, pathLength] = + stepper.boundState(state.stepping, *surface, true); + // add a full TrackState entry multi trajectory // (this allocates storage for all components, we will set them later) result.trackTip = result.fittedStates.addTrackState( @@ -323,10 +352,6 @@ class KalmanFitter { // assign the source link to the track state trackStateProxy.uncalibrated() = sourcelink_it->second; - // Transport & bind the state to the current surface - auto [boundParams, jacobian, pathLength] = - stepper.boundState(state.stepping, *surface, true); - // Fill the track state trackStateProxy.predicted() = boundParams.parameters(); trackStateProxy.predictedCovariance() = *boundParams.covariance(); @@ -342,13 +367,19 @@ class KalmanFitter { m_calibrator(trackStateProxy.uncalibrated(), trackStateProxy.predicted())); + // Get and set the type flags + auto& typeFlags = trackStateProxy.typeFlags(); + typeFlags.set(TrackStateFlag::MaterialFlag); + typeFlags.set(TrackStateFlag::MeasurementFlag); + typeFlags.set(TrackStateFlag::ParameterFlag); + // If the update is successful, set covariance and auto updateRes = m_updater(state.geoContext, trackStateProxy); if (!updateRes.ok()) { ACTS_ERROR("Update step failed: " << updateRes.error()); return updateRes.error(); } else { - // Update the stepping state + // Update the stepping state with filtered parameters ACTS_VERBOSE("Filtering step successful, updated parameters are : \n" << trackStateProxy.filtered().transpose()); // update stepping state using filtered parameters after kalman update @@ -356,18 +387,124 @@ class KalmanFitter { // a bit awkward. stepper.update(state.stepping, trackStateProxy.filteredParameters( state.options.geoContext)); + + // Update state and stepper with post material effects + materialInteractor(surface, state, stepper, postUpdate); } + // We count the state with measurement + ++result.measurementStates; // We count the processed state ++result.processedStates; - } else if (surface->associatedDetectorElement() != nullptr) { - // Count the missed surface - ACTS_VERBOSE("Detected hole on " << surface->geoID()); - result.missedActiveSurfaces.push_back(surface); - } + } else { + // add a non-measurement TrackState entry multi trajectory + // (this allocates storage for components except measurements, we will + // set them later) + result.trackTip = result.fittedStates.addTrackState( + ~(TrackStatePropMask::Uncalibrated | + TrackStatePropMask::Calibrated), + result.trackTip); + + // now get track state proxy back + auto trackStateProxy = + result.fittedStates.getTrackState(result.trackTip); + + // Set the surface + trackStateProxy.setReferenceSurface(surface->getSharedPtr()); + + // Set the track state flags + auto& typeFlags = trackStateProxy.typeFlags(); + typeFlags.set(TrackStateFlag::MaterialFlag); + typeFlags.set(TrackStateFlag::ParameterFlag); + + if (surface->associatedDetectorElement() != nullptr) { + ACTS_VERBOSE("Detected hole on " << surface->geoID()); + // If the surface is sensitive, set the hole type flag + typeFlags.set(TrackStateFlag::HoleFlag); + // Count the missed surface + result.missedActiveSurfaces.push_back(surface); + + // Transport & bind the state to the current surface + auto [boundParams, jacobian, pathLength] = + stepper.boundState(state.stepping, *surface, true); + + // Fill the track state + trackStateProxy.predicted() = boundParams.parameters(); + trackStateProxy.predictedCovariance() = *boundParams.covariance(); + trackStateProxy.jacobian() = jacobian; + trackStateProxy.pathLength() = pathLength; + } else { + ACTS_VERBOSE("Detected in-sensitive surface " << surface->geoID()); + + // Transport & get curvilinear state instead of bound state + auto [curvilinearParams, jacobian, pathLength] = + stepper.curvilinearState(state.stepping, true); + + // Fill the track state + trackStateProxy.predicted() = curvilinearParams.parameters(); + trackStateProxy.predictedCovariance() = + *curvilinearParams.covariance(); + trackStateProxy.jacobian() = jacobian; + trackStateProxy.pathLength() = pathLength; + } + + // Update state and stepper with material effects + materialInteractor(surface, state, stepper, fullUpdate); + + // Set the filtered parameter to be the same with predicted parameter + // @Todo: shall we update the filterd parameter with material effects? + // But it seems that the smoothing does not like this + trackStateProxy.filtered() = trackStateProxy.predicted(); + trackStateProxy.filteredCovariance() = + trackStateProxy.predictedCovariance(); + + // We count the processed state + ++result.processedStates; + } return Result<void>::success(); } + /// @brief Kalman actor operation : material interaction + /// + /// @tparam propagator_state_t is the type of Propagagor state + /// @tparam stepper_t Type of the stepper + /// + /// @param surface The surface where the material interaction happens + /// @param state The mutable propagator state object + /// @param stepper The stepper in use + /// @param updateStage The materal update stage + /// + template <typename propagator_state_t, typename stepper_t> + void materialInteractor(const Surface* surface, propagator_state_t& state, + stepper_t& stepper, + const MaterialUpdateStage& updateStage) const { + // Prepare relevant input particle properties + detail::PointwiseMaterialInteraction interaction(surface, state, stepper); + + // Evaluate the material properties + if (interaction.evaluateMaterialProperties(state, updateStage)) { + // Evaluate the material effects + interaction.evaluatePointwiseMaterialInteraction(multipleScattering, + energyLoss); + + ACTS_VERBOSE("Material effects on surface: " + << surface->geoID() << " at update stage: " << updateStage + << " are :"); + ACTS_VERBOSE("eLoss = " + << interaction.Eloss << ", " + << "variancePhi = " << interaction.variancePhi << ", " + << "varianceTheta = " << interaction.varianceTheta << ", " + << "varianceQoverP = " << interaction.varianceQoverP); + + // Update the state and stepper with material effects + interaction.updateState(state, stepper); + } else { + ACTS_VERBOSE("No material effects on surface: " << surface->geoID() + << " at update stage: " + << updateStage); + } + } + /// @brief Kalman actor operation : finalize /// /// @tparam propagator_state_t is the type of Propagagor state @@ -500,6 +637,8 @@ class KalmanFitter { kalmanActor.m_logger = m_logger.get(); kalmanActor.inputMeasurements = std::move(inputMeasurements); kalmanActor.targetSurface = kfOptions.referenceSurface; + kalmanActor.multipleScattering = kfOptions.multipleScattering; + kalmanActor.energyLoss = kfOptions.energyLoss; // also set logger on updater and smoother kalmanActor.m_updater.m_logger = m_logger; @@ -517,6 +656,12 @@ class KalmanFitter { /// Get the result of the fit auto kalmanResult = propRes.template get<KalmanResult>(); + /// It could happen that the fit ends in zero processed states. + /// The result gets meaningless so such case is regarded as fit failure. + if (kalmanResult.result.ok() and not kalmanResult.processedStates) { + kalmanResult.result = Result<void>(KalmanFitterError::PropagationInVain); + } + if (!kalmanResult.result.ok()) { return kalmanResult.result.error(); } @@ -600,6 +745,12 @@ class KalmanFitter { /// Get the result of the fit auto kalmanResult = propRes.template get<KalmanResult>(); + /// It could happen that the fit ends in zero processed states. + /// The result gets meaningless so such case is regarded as fit failure. + if (kalmanResult.result.ok() and not kalmanResult.processedStates) { + kalmanResult.result = Result<void>(KalmanFitterError::PropagationInVain); + } + if (!kalmanResult.result.ok()) { return kalmanResult.result.error(); } diff --git a/Core/include/Acts/Fitter/KalmanFitterError.hpp b/Core/include/Acts/Fitter/KalmanFitterError.hpp index 99793635723f0602d4e58582696ae65e94dc1a01..56be3facb94a0e65c058f797182c7d0d7c30f1c4 100644 --- a/Core/include/Acts/Fitter/KalmanFitterError.hpp +++ b/Core/include/Acts/Fitter/KalmanFitterError.hpp @@ -17,7 +17,8 @@ namespace Acts { enum class KalmanFitterError { UpdateFailed = 1, SmoothFailed = 2, - OutputConversionFailed = 3 + OutputConversionFailed = 3, + PropagationInVain = 4 }; namespace detail { @@ -35,6 +36,8 @@ class KalmanFitterErrorCategory : public std::error_category { return "Kalman smooth failed"; case KalmanFitterError::OutputConversionFailed: return "Kalman output conversion failed"; + case KalmanFitterError::PropagationInVain: + return "No detector observed during the propagation"; default: return "unknown"; } diff --git a/Core/include/Acts/Propagator/detail/PointwiseMaterialInteraction.hpp b/Core/include/Acts/Propagator/detail/PointwiseMaterialInteraction.hpp index c10cda3f37fdb6d9606588e6675995b4f7059eaf..76990252591b314153655a9cf865feb3ef8ea0e6 100644 --- a/Core/include/Acts/Propagator/detail/PointwiseMaterialInteraction.hpp +++ b/Core/include/Acts/Propagator/detail/PointwiseMaterialInteraction.hpp @@ -8,6 +8,7 @@ #pragma once +#include "Acts/Material/ISurfaceMaterial.hpp" #include "Acts/Material/MaterialProperties.hpp" #include "Acts/Surfaces/Surface.hpp" @@ -148,4 +149,4 @@ struct PointwiseMaterialInteraction { double updateVariance(double variance, double change) const; }; } // namespace detail -} // end of namespace Acts \ No newline at end of file +} // end of namespace Acts diff --git a/Core/src/Propagator/detail/PointwiseMaterialInteraction.cpp b/Core/src/Propagator/detail/PointwiseMaterialInteraction.cpp index 75a589917430f78cb0e776ba430ca14a77488620..efb940d8d6627ed0213232786931c121d88e25ca 100644 --- a/Core/src/Propagator/detail/PointwiseMaterialInteraction.cpp +++ b/Core/src/Propagator/detail/PointwiseMaterialInteraction.cpp @@ -51,4 +51,4 @@ double PointwiseMaterialInteraction::updateVariance(double variance, return std::max(0., variance + std::copysign(change, nav)); } } // namespace detail -} // end of namespace Acts \ No newline at end of file +} // end of namespace Acts