EIC code displayed by LXR |
|
|
|
|||
File indexing completed on 2025-12-15 09:42:16
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/EventData/CompositeSpacePoint.hpp" 0012 #include "Acts/Utilities/ArrayHelpers.hpp" 0013 #include "Acts/Utilities/Logger.hpp" 0014 #include "Acts/Utilities/detail/Line3DWithPartialDerivatives.hpp" 0015 0016 #include <cstdint> 0017 0018 namespace Acts { 0019 template <Experimental::CompositeSpacePoint SpacePoint_t> 0020 /// @brief Print the position, the drift radius & the sensor directions of a space point. 0021 /// If the space point is shipped with an ostream operator, this oone is 0022 /// used 0023 /// @param measurement: Reference to the space point to print 0024 std::string toString(const SpacePoint_t& measurement); 0025 } 0026 0027 namespace Acts::Experimental::detail { 0028 /// @brief Helper class to calculate the residual between a straight line and 0029 /// a CompositeSpacePoint measurement as well as the partial derivatives. 0030 /// The residual is expressed as a 3D vector, where first two components 0031 /// describe the spatial residual w.r.t. the precision & non-precision 0032 /// direction, and the last one expresses the residual between the 0033 /// parametrized time of arrival of the track and the measurement's 0034 /// recorded time. 0035 /// 0036 /// For straw type measurements, the precision residual is calculated as 0037 /// the difference between the signed line distance between straw wire & 0038 /// line, and the signed drift radius where the sign simply encodes the 0039 /// left/right amibiguity. If the measurement additionally carries 0040 /// information about the passage along the wire, the non-precision 0041 /// residual is then defined to be the distance along the wire. 0042 /// 0043 /// For strip type measurements, the residual is the distance in the 0044 /// strip-readout plane between the point along the line that intersects 0045 /// the plane spanned by the measurement. 0046 /// 0047 /// If the strip measurements provide also the time of their record, the 0048 /// residual calculation can be extended to the time of arrival parameter 0049 /// The residual is then defined as the strip's recorded time minus the 0050 /// time of flight of a particle on a straight line minus the time 0051 /// offset. 0052 /// 0053 /// Straw measurements are indirectly influenced by the time offset 0054 /// parameter. The primary electrons produced by the traversing ionizing 0055 /// particle drift towards the central straw wire. The drift time can be 0056 /// directly mapped to the drift radius during the calibration procedure 0057 /// where the rt relation and its first two derivatives are known 0058 /// apprxomately. Further, it's assumed that the chip records are 0059 /// composed timestamp that includes the drift time & the time of flight 0060 /// of the ionizing particle. An update of the point of closest approach 0061 /// leads to an indirect update of the drift radius and hence additional 0062 /// terms need to be considered when calculating the residual's 0063 /// derivatives. 0064 0065 class CompSpacePointAuxiliaries { 0066 public: 0067 using Line_t = Acts::detail::Line3DWithPartialDerivatives<double>; 0068 using LineIndex = Line_t::ParIndex; 0069 using Vector = Line_t::Vector; 0070 enum class FitParIndex : std::uint8_t { 0071 x0 = toUnderlying(LineIndex::x0), 0072 y0 = toUnderlying(LineIndex::y0), 0073 theta = toUnderlying(LineIndex::theta), 0074 phi = toUnderlying(LineIndex::phi), 0075 t0 = 4, // time offset 0076 nPars = 5 0077 }; 0078 static constexpr auto s_nPars = toUnderlying(FitParIndex::nPars); 0079 /// @brief Prints a fit parameter as string 0080 static std::string parName(const FitParIndex idx); 0081 /// @brief Assignment of the residual components. 0082 enum class ResidualIdx : std::uint8_t { 0083 nonBending = 0, 0084 bending = 1, 0085 time = 2 0086 }; 0087 /// @brief Configuration object of the residual calculator 0088 struct Config { 0089 /// @brief Transform to place the composite station frame inside the 0090 /// global experiment's frame. Needed for the time residual 0091 /// calculation with time of flight 0092 Acts::Transform3 localToGlobal{Acts::Transform3::Identity()}; 0093 /// @brief Flag toggling whether the hessian of the residual shall be calculated 0094 bool useHessian{false}; 0095 /// @brief Flag toggling whether the along the wire component of straws shall be calculated 0096 /// if provided by the straw measurement. 0097 bool calcAlongStraw{true}; 0098 /// @brief Flag toggling whether the residual along the strip direction 0099 /// shall be calculated if the space point does not measure both 0100 /// spatial coordinates on the plane 0101 bool calcAlongStrip{true}; 0102 /// @brief Include the time of flight assuming that the particle travels with the 0103 /// speed of light in the time residual calculations 0104 bool includeToF{true}; 0105 /// @brief List of fit parameters to which the partial derivative of the 0106 /// residual shall be calculated 0107 std::vector<FitParIndex> parsToUse{FitParIndex::x0, FitParIndex::y0, 0108 FitParIndex::theta, FitParIndex::phi}; 0109 }; 0110 /// @brief Constructor to instantiate a new instance 0111 /// @param cfg: Configuration object to toggle the calculation of the complementary residual components & the full evaluation of the second derivative 0112 /// @param logger: New logging object for debugging 0113 explicit CompSpacePointAuxiliaries( 0114 const Config& cfg, 0115 std::unique_ptr<const Logger> logger = 0116 getDefaultLogger("CompSpacePointAuxiliaries", Logging::Level::INFO)); 0117 0118 /// @brief Returns the config object 0119 const Config& config() const { return m_cfg; } 0120 /// @brief Updates the spatial residual components between the line and the passed 0121 /// measurement. The result is cached internally and can be later 0122 /// fetched by the residual(), gradient() and hessian() methods. If the 0123 /// residual calculation fails due to parallel line & measurement, all 0124 /// components are set to zero. 0125 /// @param line: Reference to the line to which the residual is calculated 0126 /// @param spacePoint: Reference to the space point measurement to which the residual is calculated 0127 template <CompositeSpacePoint Point_t> 0128 void updateSpatialResidual(const Line_t& line, const Point_t& spacePoint); 0129 /// @brief Updates all residual components between the line and the passed measurement 0130 /// First the spatial components are calculated and then if the 0131 /// measurement also provides time information, the time residual & its 0132 /// derivatives are evaluated. 0133 /// @param line: Reference to the line to which the residual is calculated 0134 /// @param timeOffet: Value of the t0 fit parameter. 0135 /// @param spacePoint: Reference to the space point measurement to which the residual is calculated 0136 /// @param driftV: Associated drift velocity given as the derivative of the r-t relation 0137 /// @param driftA: Associated drift acceleration given as the second derivative of the r-t relation 0138 template <CompositeSpacePoint Point_t> 0139 void updateFullResidual(const Line_t& line, const double timeOffset, 0140 const Point_t& spacePoint, const double driftV = 0., 0141 const double driftA = 0.); 0142 0143 /// @brief Helper struct to calculate the overall chi2 from the composite space points 0144 struct ChiSqWithDerivatives { 0145 /// @brief Chi2 squared term 0146 double chi2{0.}; 0147 /// @brief First derivative of the chi2 w.r.t. the fit parameters 0148 Acts::ActsVector<s_nPars> gradient{Acts::ActsVector<s_nPars>::Zero()}; 0149 /// @brief Second derivative of the chi2 w.r.t. the fit parameters 0150 Acts::ActsSquareMatrix<s_nPars> hessian{ 0151 Acts::ActsSquareMatrix<s_nPars>::Zero()}; 0152 /// @brief Set the chi2, the gradient and hessian back to zero 0153 void reset(); 0154 }; 0155 /// @brief Updates the passed chi2 object by adding up the residual contributions 0156 /// from the previous composite space point used to update the 0157 /// Auxiliary class. For the Hessian term only the lower triangle is 0158 /// updated. The other triangle needs to be copied later from the lower 0159 /// one. 0160 /// @param chiSqObj: Chi2 & derivatives to be updated 0161 /// @param cov: The composite space point's covariance values 0162 void updateChiSq(ChiSqWithDerivatives& chiSqObj, 0163 const std::array<double, 3>& cov) const; 0164 /// @brief Fills the upper triangle of the Hessian matrix with the 0165 /// values from the lower triangle 0166 /// @param chiSqObj: Reference to the chiSqObj carrying the Hessian 0167 void symmetrizeHessian(ChiSqWithDerivatives& chiSqObj) const; 0168 /// @brief Returns the previously calculated residual. 0169 const Vector& residual() const; 0170 /// @brief Returns the gradient of the previously calculated residual 0171 /// @param par: Index of the partiald derivative 0172 const Vector& gradient(const FitParIndex param) const; 0173 /// @brief Returns the gradient of the previously calculated residual 0174 /// @param param1: First index of the second partial derivative 0175 /// @param param2: Second index of the second partial derivative 0176 const Vector& hessian(const FitParIndex param1, 0177 const FitParIndex param2) const; 0178 /// @brief Returns whether the passed parameter describes a direction angle 0179 static constexpr bool isDirectionParam(const FitParIndex param) { 0180 return param == FitParIndex::theta || param == FitParIndex::phi; 0181 } 0182 /// @brief Returns whether the passed parameter describes the displacement 0183 /// in the reference plane 0184 static constexpr bool isPositionParam(const FitParIndex param) { 0185 return param == FitParIndex::x0 || param == FitParIndex::y0; 0186 } 0187 /// @brief Extrapolates the straight line onto the space point's strip-plane defined by the 0188 /// `localPosition()` and the `planeNormal()` 0189 /// @param line: Reference to the line of interest 0190 /// @param hit: Reference to the hit of interest 0191 template <CompositeSpacePoint SpacePoint_t> 0192 static Vector extrapolateToPlane(const Line_t& line, const SpacePoint_t& hit); 0193 /// @brief Extrapolates the straight line parametrized by a point and a direction 0194 /// onto the space point's strip-plane defined by the 0195 /// `localPosition()` and the `planeNormal()` 0196 /// @param pos: Point on the line 0197 /// @param dir: Direction of the line 0198 /// @param hit: Reference to the hit of interest 0199 template <CompositeSpacePoint SpacePoint_t> 0200 static Vector extrapolateToPlane(const Vector& pos, const Vector& dir, 0201 const SpacePoint_t& hit); 0202 /// @brief Calculates the spaztial chiSq term of a composite space point to a line 0203 /// @param line: Reference to the line of interest 0204 /// @param hit: Reference to the hit of interest 0205 template <CompositeSpacePoint SpacePoint_t> 0206 static double chi2Term(const Line_t& line, const SpacePoint_t& hit); 0207 /// @brief Calculates the spatial chiSq term of a composite space point to a line parameterized 0208 /// by a position & direction vector 0209 /// @param pos: Point on the line 0210 /// @param dir: Direction of the line 0211 /// @param hit: Reference to the hit of interest 0212 template <CompositeSpacePoint SpacePoint_t> 0213 static double chi2Term(const Vector& pos, const Vector& dir, 0214 const SpacePoint_t& hit); 0215 /// @brief Calculates the chiSq term of a composite space point to a line taking into account 0216 /// the time offset t0 0217 /// @param line: Reference to the line of interest 0218 /// @param t0: Time off set evaluated at the measurement's plane 0219 /// @param hit: Reference to the hit of interest 0220 template <CompositeSpacePoint SpacePoint_t> 0221 static double chi2Term(const Line_t& line, const double t0, 0222 const SpacePoint_t& hit); 0223 /// @brief Calculates the chiSq term of a composite space point to a line taking into account 0224 /// the time offset t0 0225 /// @param pos: Point on the line 0226 /// @param dir: Direction of the line 0227 /// @param t0: Time off set evaluated at the measurement's plane 0228 /// @param hit: Reference to the hit of interest 0229 template <CompositeSpacePoint SpacePoint_t> 0230 static double chi2Term(const Vector& pos, const Vector& dir, const double t0, 0231 const SpacePoint_t& hit); 0232 /// @brief Calculates the chiSq term of a composite space point taking into account 0233 /// the time offset t0 & the time of arrival of the particle assuming 0234 /// the speed of light 0235 /// @param line: Reference to the line of interest 0236 /// @param t0: Time off set evaluated at the measurement's plane 0237 /// @param hit: Reference to the hit of interest 0238 template <CompositeSpacePoint SpacePoint_t> 0239 static double chi2Term(const Line_t& line, 0240 const Acts::Transform3& localToGlobal, const double t0, 0241 const SpacePoint_t& hit); 0242 /// @brief Calculates the chiSq term of a composite space point taking into account 0243 /// the time offset t0 & the time of arrival of the particle assuming 0244 /// the speed of light 0245 /// @param pos: Point on the line 0246 /// @param dir: Direction of the line 0247 /// @param t0: Time off set evaluated at the measurement's plane 0248 /// @param hit: Reference to the hit of interest 0249 template <CompositeSpacePoint SpacePoint_t> 0250 static double chi2Term(const Vector& pos, const Vector& dir, 0251 const Acts::Transform3& localToGlobal, const double t0, 0252 const SpacePoint_t& hit); 0253 /// @brief Calculate whether the track passed on the left (-1) or the right (1) side 0254 /// of the straw wire. Returns 0 for strips 0255 /// @param line: Reference to the line of interest 0256 /// @param strawSp: Straw measurement of interest 0257 template <CompositeSpacePoint Point_t> 0258 static int strawSign(const Line_t& line, const Point_t& strawSp); 0259 /// @brief Calculate whether a generic line is on the left (-1) or right (1) side 0260 /// of the straw wire. Return 0 for strip 0261 /// @param pos: Point on the line 0262 /// @param dir: Direction of the line 0263 /// @param strawSp: Straw measurement of interest 0264 template <CompositeSpacePoint Point_t> 0265 static int strawSign(const Vector& pos, const Vector& dir, 0266 const Point_t& strawSp); 0267 /// @brief Calculate the straw signs for a set of measurements 0268 /// @param line: Reference to the line to which the residual is calculated 0269 /// @param measurements: List of straw measurements to calculate the signs for 0270 template <CompositeSpacePointContainer StrawCont_t> 0271 static std::vector<int> strawSigns(const Line_t& line, 0272 const StrawCont_t& measurements); 0273 /// @brief Calculate the straw signs for a set of measurements 0274 /// @param pos: Point on the line 0275 /// @param dir: Direction of the line 0276 /// @param measurements: List of straw measurements to calculate the signs for 0277 template <CompositeSpacePointContainer StrawCont_t> 0278 static std::vector<int> strawSigns(const Vector& pos, const Vector& dir, 0279 const StrawCont_t& measurements); 0280 0281 private: 0282 /// @brief Reference to the logging object 0283 const Logger& logger() const { return *m_logger; } 0284 /// @brief Update the auxiliary variables needed to calculate the residuals 0285 /// of a straw measurement to the current line. They are the 0286 /// projection of the line direction vector into the straw's wire 0287 /// plane and its derivatives. Returns fals if line and wire are 0288 /// parallel to each other 0289 /// @param line: Reference to the line to project 0290 /// @param wireDir: Reference to the straw wire direction vector 0291 bool updateStrawAuxiliaries(const Line_t& line, const Vector& wireDir); 0292 /// @brief Updates the residuals of a straw measurement to a line ignoring 0293 /// the distance along the straw wire. Returns false if the residual 0294 /// cannot be updated due to parallelism between straw wire and line 0295 /// @param line: Reference to the line to which the residual is calculated 0296 /// @param hitMinSeg: Difference of the line reference point & the straw 0297 /// position 0298 /// @param wireDir: Direction vector of the straw wire 0299 /// @param driftRadius: Measured radius of the straw measurement 0300 bool updateStrawResidual(const Line_t& line, const Vector& hitMinSeg, 0301 const Vector& wireDir, const double driftRadius); 0302 /// @brief Calculates the along-the wire component of the straw 0303 /// measurement's residual to the line 0304 /// @param line: Reference to the line to which the residual is calculated 0305 /// @param hitMinSeg: Difference of the straw position & the line reference point 0306 /// @param wireDir: Direction vector of the straw wire 0307 void updateAlongTheStraw(const Line_t& line, const Vector& hitMinSeg, 0308 const Vector& wireDir); 0309 /// @brief Calculates the residual of a strip measurement w.r.t. the line 0310 /// @param line: Reference to the line to which the residual is calculated 0311 /// @param normal: Reference to the vector normal on the strip measurement plane 0312 /// @param sensorN: Reference to the first basis vector inside the strip measruement plane, 0313 /// which is given by the sensor normal 0314 /// @param sensorD: Reference to the second basis vector inside the strip measruement plane, 0315 /// which is given by the sensor direction 0316 /// @param stripPos: Position of the strip measurement 0317 /// @param isBending: Flag toggling whether the precision direction is constrained 0318 /// @param isNonBending: Flag toggling whether the non-precision direction is constrained 0319 void updateStripResidual(const Line_t& line, const Vector& normal, 0320 const Vector& sensorN, const Vector& sensorD, 0321 const Vector& stripPos, const bool isBending, 0322 const bool isNonBending); 0323 /// @brief Calculates the reidual of a strip measurement w.r.t. the time offset parameter 0324 /// @param sensorN: Reference to the first basis vector inside the strip measruement plane, 0325 /// which is given by the sensor normal 0326 /// @param sensorD: Reference to the second basis vector inside the strip measruement plane, 0327 /// which is given by the sensor direction 0328 /// @param stripPos: Position of the strip measurement 0329 /// @param isBending: Flag toggling whether the precision direction is constrained 0330 /// @param recordTime: Time of the measurement 0331 /// @param timeOffet: Value of the t0 fit parameter. 0332 void updateTimeStripRes(const Vector& sensorN, const Vector& sensorD, 0333 const Vector& stripPos, const bool isBending, 0334 const double recordTime, const double timeOffset); 0335 /// @brief Calculates the residual derivatives of a straw tube measurement when the drift radius is 0336 // an implicit function of the point of closest approach and the time 0337 // offset parameter 0338 /// @param line: Reference to the line to which the residual is calculated 0339 /// @param hitMinSeg: Difference of the straw position & the line reference point 0340 /// @param wireDir: The direction along the wire 0341 /// @param driftR: Current drift radius of the straw measurement 0342 /// @param driftV: Associated drift velocity given as the derivative of the r-t relation 0343 /// @param driftA: Associated drift acceleration given as the second derivative of the r-t relation 0344 void updateTimeStrawRes(const Line_t& line, const Vector& hitMinSeg, 0345 const Vector& wireDir, const double driftR, 0346 const double driftV, const double driftA); 0347 /// @brief Resets the residual and all partial derivatives to zero. 0348 void reset(); 0349 /// @brief Resets the time residual and the partial derivatives 0350 void resetTime(); 0351 Config m_cfg{}; 0352 std::unique_ptr<const Logger> m_logger{}; 0353 0354 /// @brief Cached residual vector calculated from the measurement & the parametrized line 0355 Vector m_residual{Vector::Zero()}; 0356 /// @brief Partial derivatives of the residual w.r.t. the fit parameters parameters 0357 std::array<Vector3, s_nPars> m_gradient{ 0358 filledArray<Vector3, s_nPars>(Vector3::Zero())}; 0359 /// @brief Second partial derivatives of the residual w.r.t. the fit parameters parameters 0360 std::array<Vector3, sumUpToN(s_nPars)> m_hessian{ 0361 filledArray<Vector3, sumUpToN(s_nPars)>(Vector3::Zero())}; 0362 0363 // Auxiliary parameters needed to calculate the residual of the 0364 // point of closest approach and its derivatives 0365 0366 /// @brief Number of spatial line parameters 0367 static constexpr std::uint8_t s_nLinePars = Line_t::s_nPars; 0368 /// @brief projection of the segment direction onto the wire plane 0369 Vector m_projDir{Vector::Zero()}; 0370 /// @brief Partial derivatives of the dir projection w.r.t. line parameters 0371 std::array<Vector, s_nLinePars> m_gradProjDir{ 0372 filledArray<Vector, s_nLinePars>(Vector::Zero())}; 0373 /// @brief Component of the direction vector parallel to the wire 0374 double m_wireProject{0.}; 0375 /// @brief Length squared of the projected direction 0376 double m_invProjDirLenSq{0.}; 0377 /// @brief Inverse of the projected direction length 0378 double m_invProjDirLen{0.}; 0379 /// @brief Partial derivatives of the dir projection lengths w.r.t line parameters 0380 std::array<double, s_nLinePars> m_projDirLenPartial{ 0381 filledArray<double, s_nLinePars>(0.)}; 0382 0383 std::array<Vector, sumUpToN(s_nLinePars)> m_hessianProjDir{ 0384 filledArray<Vector, sumUpToN(s_nLinePars)>(Vector::Zero())}; 0385 0386 /// @brief Gradient vector of the point of closest approach 0387 std::array<Vector, s_nLinePars> m_gradCloseApproach{ 0388 filledArray<Vector, s_nLinePars>(Vector::Zero())}; 0389 /// @brief Partial derivative of the actual distance of the closest approach 0390 std::array<double, s_nLinePars> m_partialApproachDist{ 0391 filledArray<double, s_nLinePars>(0.)}; 0392 /// @brief Tansform matrix to treat stereo angles amongst the strips 0393 ActsSquareMatrix<2> m_stereoTrf{ActsSquareMatrix<2>::Identity()}; 0394 }; 0395 0396 } // namespace Acts::Experimental::detail 0397 #include "Acts/Seeding/detail/CompSpacePointAuxiliaries.ipp"
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|