From 81eb8551687ced6273944d71b58f7697347282ad Mon Sep 17 00:00:00 2001 From: "se.gorbunov" <se.gorbunov@gsi.de> Date: Wed, 23 Oct 2024 14:53:31 +0200 Subject: [PATCH] Kf: replace KfTrackParamBase with KfTrackParam --- algo/ca/core/pars/CaDefs.h | 2 +- algo/kf/core/KfTrackKalmanFilter.cxx | 17 +- algo/kf/core/KfTrackKalmanFilter.h | 18 +- algo/kf/core/data/KfTrackParam.h | 573 ++++++++++++++------------- algo/kf/core/geo/KfFieldSlice.h | 2 +- reco/KF/CbmKfTrackFitter.cxx | 2 +- reco/KF/CbmKfTrackFitter.h | 4 +- 7 files changed, 316 insertions(+), 302 deletions(-) diff --git a/algo/ca/core/pars/CaDefs.h b/algo/ca/core/pars/CaDefs.h index 99e4c1cdd8..8b4a2280d5 100644 --- a/algo/ca/core/pars/CaDefs.h +++ b/algo/ca/core/pars/CaDefs.h @@ -17,7 +17,7 @@ namespace cbm::algo::ca { - using cbm::algo::kf::TrackParamBase; + using cbm::algo::kf::TrackParam; using cbm::algo::kf::TrackParamV; } // namespace cbm::algo::ca diff --git a/algo/kf/core/KfTrackKalmanFilter.cxx b/algo/kf/core/KfTrackKalmanFilter.cxx index 17801e77bf..1f0966a330 100644 --- a/algo/kf/core/KfTrackKalmanFilter.cxx +++ b/algo/kf/core/KfTrackKalmanFilter.cxx @@ -327,10 +327,9 @@ namespace cbm::algo::kf } template<typename DataT> - void - TrackKalmanFilter<DataT>::FilterExtrapolatedXY(const kf::MeasurementXy<DataT>& m, DataT extrX, DataT extrY, - const std::array<DataT, kf::TrackParamBase<DataT>::kNtrackParam>& Jx, - const std::array<DataT, kf::TrackParamBase<DataT>::kNtrackParam>& Jy) + void TrackKalmanFilter<DataT>::FilterExtrapolatedXY(const kf::MeasurementXy<DataT>& m, DataT extrX, DataT extrY, + const std::array<DataT, kf::TrackParam<DataT>::kNtrackParam>& Jx, + const std::array<DataT, kf::TrackParam<DataT>::kNtrackParam>& Jy) { // add a 2-D measurenent (x,y) at some z, that differs from fTr.GetZ() // extrX, extrY are extrapolated track parameters at z, Jx, Jy are derivatives of the extrapolation @@ -991,7 +990,7 @@ namespace cbm::algo::kf DataT txtx1 = DataT(1.) + txtx; DataT tyty1 = DataT(1.) + tyty; DataT t = sqrt(txtx1 + tyty); - DataT qpt = qp * t; + // DataT qpt = qp * t; DataT lg = DataT(.0136) * (DataT(1.) + DataT(0.038) * log(radThick * t)); lg = kf::utils::iif(lg > DataT(0.), lg, DataT(0.)); @@ -1161,10 +1160,10 @@ namespace cbm::algo::kf template<typename DataT> - void - TrackKalmanFilter<DataT>::GetExtrapolatedXYline(DataT z, const kf::FieldRegion<DataT>& F, DataT& extrX, DataT& extrY, - std::array<DataT, kf::TrackParamBase<DataT>::kNtrackParam>& Jx, - std::array<DataT, kf::TrackParamBase<DataT>::kNtrackParam>& Jy) const + void TrackKalmanFilter<DataT>::GetExtrapolatedXYline(DataT z, const kf::FieldRegion<DataT>& F, DataT& extrX, + DataT& extrY, + std::array<DataT, kf::TrackParam<DataT>::kNtrackParam>& Jx, + std::array<DataT, kf::TrackParam<DataT>::kNtrackParam>& Jy) const { // extrapolate track assuming it is straight (qp==0) // return the extrapolated X, Y and the derivatives of the extrapolated X and Y diff --git a/algo/kf/core/KfTrackKalmanFilter.h b/algo/kf/core/KfTrackKalmanFilter.h index ad18cab106..ccb09396c6 100644 --- a/algo/kf/core/KfTrackKalmanFilter.h +++ b/algo/kf/core/KfTrackKalmanFilter.h @@ -46,12 +46,12 @@ namespace cbm::algo::kf TrackKalmanFilter() = default; - TrackKalmanFilter(const kf::TrackParamBase<DataT>& t) { SetTrack(t); } + TrackKalmanFilter(const kf::TrackParam<DataT>& t) { SetTrack(t); } TrackKalmanFilter(const DataTmask& m, bool fitV) : fMask(m), fDoFitVelocity(fitV) {} template<typename T> - TrackKalmanFilter(const kf::TrackParamBase<T>& t) + TrackKalmanFilter(const kf::TrackParam<T>& t) { SetTrack(t); } @@ -59,7 +59,7 @@ namespace cbm::algo::kf void SetMask(const DataTmask& m) { fMask = m; } template<typename T> - void SetTrack(const kf::TrackParamBase<T>& t) + void SetTrack(const kf::TrackParam<T>& t) { fTr.Set(t); fQp0 = fTr.GetQp(); @@ -67,7 +67,7 @@ namespace cbm::algo::kf void SetQp0(DataT qp0) { fQp0 = qp0; } - kf::TrackParamBase<DataT>& Tr() { return fTr; } + kf::TrackParam<DataT>& Tr() { return fTr; } DataT& Qp0() { return fQp0; } @@ -164,8 +164,8 @@ namespace cbm::algo::kf /// extrapolate track as a line, return the extrapolated X, Y and the Jacobians void GetExtrapolatedXYline(DataT z, const kf::FieldRegion<DataT>& F, DataT& extrX, DataT& extrY, - std::array<DataT, kf::TrackParamBase<DataT>::kNtrackParam>& Jx, - std::array<DataT, kf::TrackParamBase<DataT>::kNtrackParam>& Jy) const; + std::array<DataT, kf::TrackParam<DataT>::kNtrackParam>& Jx, + std::array<DataT, kf::TrackParam<DataT>::kNtrackParam>& Jy) const; /// filter the track with the XY measurement placed at different Z /// \param m - measurement @@ -174,8 +174,8 @@ namespace cbm::algo::kf /// \param Jx - Jacobian of the extrapolated X /// \param Jy - Jacobian of the extrapolated Y void FilterExtrapolatedXY(const kf::MeasurementXy<DataT>& m, DataT extrX, DataT extrY, - const std::array<DataT, kf::TrackParamBase<DataT>::kNtrackParam>& Jx, - const std::array<DataT, kf::TrackParamBase<DataT>::kNtrackParam>& Jy); + const std::array<DataT, kf::TrackParam<DataT>::kNtrackParam>& Jx, + const std::array<DataT, kf::TrackParam<DataT>::kNtrackParam>& Jy); /// extrapolate the track to the given Z using linearization at the straight line, /// \param z_out - Z coordinate to extrapolate to @@ -244,7 +244,7 @@ namespace cbm::algo::kf DataTmask fMask{true}; ///< mask of active elements in simd vectors - kf::TrackParamBase<DataT> fTr{}; ///< track parameters + kf::TrackParam<DataT> fTr{}; ///< track parameters DataT fQp0{0.}; DataT fMass{0.10565800}; ///< particle mass (muon mass by default) diff --git a/algo/kf/core/data/KfTrackParam.h b/algo/kf/core/data/KfTrackParam.h index fb5fd53a77..d439bbd8c8 100644 --- a/algo/kf/core/data/KfTrackParam.h +++ b/algo/kf/core/data/KfTrackParam.h @@ -22,14 +22,15 @@ namespace cbm::algo::kf { /// \class cbm::algo::kf::TrackParamBase - /// \brief Class represents track parameters in the CA tracking algorithm + /// \brief It is a technical base class of kf::TrackParam + /// \brief that represents fitted track parameters at given z /// /// This is a template that can be instantiated for different floating point types. /// Track parameters: variable parameterts {x, y, tx, ty, q/p, t, vi} at fixed z /// Covariation matrix: low-diagopnal C[28] /// - template<typename DataT> - class TrackParamBase { + template<typename T> + class alignas(VcMemAlign) TrackParamBase { public: friend class boost::serialization::access; @@ -37,39 +38,39 @@ namespace cbm::algo::kf static constexpr int kNtrackParam{7}; ///< N of variable track parameters: {x, y, tx, ty, q/p, t, vi} static constexpr int kNcovParam{(kNtrackParam) * (kNtrackParam + 1) / 2}; ///< N of covariance matrix parameters - typedef std::array<DataT, kNcovParam> CovMatrix_t; ///< covariance matrix type + typedef std::array<T, kNcovParam> CovMatrix_t; ///< covariance matrix type TrackParamBase() = default; - template<typename T> - TrackParamBase(const TrackParamBase<T>& tr) + template<typename T1> + TrackParamBase(const TrackParamBase<T1>& tr) { Set(tr); } /// Set all SIMD entries from all SIMD entries of the other class /// It works for scalar and fvec types, - /// except of the case when DataT is scalar and TdataB is fvec. - template<typename TdataB> - void Set(const TrackParamBase<TdataB>& Tb) + /// except of the case when T is scalar and TdataB is fvec. + template<typename T1> + void Set(const TrackParamBase<T1>& Tb) { - CopyBase<TdataB, true, true>(0, Tb, 0); + CopyBase<T1, true, true>(0, Tb, 0); } /// Set all SIMD entries from one SIMD entry of the other class - /// It also works when DataT is scalar + /// It also works when T is scalar void Set(const TrackParamBase<fvec>& Tb, const int ib) { CopyBase<fvec, true, false>(0, Tb, ib); } /// Set one SIMD entry from one SIMD entry of the other class - /// It only works when DataT is fvec, TdataB is scalar - template<typename TdataB> - void SetOneEntry(const int ia, const TrackParamBase<TdataB>& Tb) + /// It only works when T is fvec, TdataB is scalar + template<typename T1> + void SetOneEntry(const int ia, const TrackParamBase<T1>& Tb) { - CopyBase<TdataB, false, true>(ia, Tb, 0); + CopyBase<T1, false, true>(ia, Tb, 0); } /// Set one SIMD entry from one SIMD entry of the other class - /// It only works when DataT is fvec, TdataB is fvec + /// It only works when T is fvec, TdataB is fvec void SetOneEntry(const int ia, const TrackParamBase<fvec>& Tb, const int ib) { CopyBase<fvec, false, false>(ia, Tb, ib); @@ -80,46 +81,46 @@ namespace cbm::algo::kf /// /// \brief Gets z position [cm] - DataT Z() const { return fZ; } + T Z() const { return fZ; } /// \brief Gets x position [cm] - DataT X() const { return fX; } + T X() const { return fX; } /// \brief Gets y position [cm] - DataT Y() const { return fY; } + T Y() const { return fY; } /// \brief Gets slope along x-axis - DataT Tx() const { return fTx; } + T Tx() const { return fTx; } /// \brief Gets slope along y-axis - DataT Ty() const { return fTy; } + T Ty() const { return fTy; } /// \brief Gets charge over momentum [ec/GeV] - DataT Qp() const { return fQp; } + T Qp() const { return fQp; } /// \brief Gets time [ns] - DataT Time() const { return fT; } + T Time() const { return fT; } /// \brief Gets inverse velocity [ns/cm] - DataT Vi() const { return fVi; } + T Vi() const { return fVi; } /// \brief Gets Chi-square of track fit model - DataT ChiSq() const { return fChiSq; } + T ChiSq() const { return fChiSq; } /// \brief Gets NDF of track fit model - DataT Ndf() const { return fNdf; } + T Ndf() const { return fNdf; } /// \brief Gets Chi-square of time measurements - DataT ChiSqTime() const { return fChiSqTime; } + T ChiSqTime() const { return fChiSqTime; } /// \brief Gets NDF of time measurements - DataT NdfTime() const { return fNdfTime; } + T NdfTime() const { return fNdfTime; } /// \brief Get covariance matrix element /// \param i row /// \param j column /// \return covariance matrix element - DataT C(int i, int j) const + T C(int i, int j) const { int ind = (j <= i) ? i * (1 + i) / 2 + j : j * (1 + j) / 2 + i; return fCovMatrix[ind]; @@ -130,7 +131,7 @@ namespace cbm::algo::kf /// \param j column /// \return matrix element template<int i, int j> - DataT C() const + T C() const { constexpr int ind = (j <= i) ? i * (1 + i) / 2 + j : j * (1 + j) / 2 + i; return fCovMatrix[ind]; @@ -139,110 +140,110 @@ namespace cbm::algo::kf /// \brief Individual getters for covariance matrix elements /// \return covariance matrix element /// - DataT C00() const { return C<0, 0>(); } - DataT C01() const { return C<0, 1>(); } - DataT C02() const { return C<0, 2>(); } - DataT C03() const { return C<0, 3>(); } - DataT C04() const { return C<0, 4>(); } - DataT C05() const { return C<0, 5>(); } - DataT C06() const { return C<0, 6>(); } - - DataT C10() const { return C<1, 0>(); } - DataT C11() const { return C<1, 1>(); } - DataT C12() const { return C<1, 2>(); } - DataT C13() const { return C<1, 3>(); } - DataT C14() const { return C<1, 4>(); } - DataT C15() const { return C<1, 5>(); } - DataT C16() const { return C<1, 6>(); } - - DataT C20() const { return C<2, 0>(); } - DataT C21() const { return C<2, 1>(); } - DataT C22() const { return C<2, 2>(); } - DataT C23() const { return C<2, 3>(); } - DataT C24() const { return C<2, 4>(); } - DataT C25() const { return C<2, 5>(); } - DataT C26() const { return C<2, 6>(); } - - DataT C30() const { return C<3, 0>(); } - DataT C31() const { return C<3, 1>(); } - DataT C32() const { return C<3, 2>(); } - DataT C33() const { return C<3, 3>(); } - DataT C34() const { return C<3, 4>(); } - DataT C35() const { return C<3, 5>(); } - DataT C36() const { return C<3, 6>(); } - - DataT C40() const { return C<4, 0>(); } - DataT C41() const { return C<4, 1>(); } - DataT C42() const { return C<4, 2>(); } - DataT C43() const { return C<4, 3>(); } - DataT C44() const { return C<4, 4>(); } - DataT C45() const { return C<4, 5>(); } - DataT C46() const { return C<4, 6>(); } - - DataT C50() const { return C<5, 0>(); } - DataT C51() const { return C<5, 1>(); } - DataT C52() const { return C<5, 2>(); } - DataT C53() const { return C<5, 3>(); } - DataT C54() const { return C<5, 4>(); } - DataT C55() const { return C<5, 5>(); } - DataT C56() const { return C<5, 6>(); } - - DataT C60() const { return C<6, 0>(); } - DataT C61() const { return C<6, 1>(); } - DataT C62() const { return C<6, 2>(); } - DataT C63() const { return C<6, 3>(); } - DataT C64() const { return C<6, 4>(); } - DataT C65() const { return C<6, 5>(); } - DataT C66() const { return C<6, 6>(); } + T C00() const { return C<0, 0>(); } + T C01() const { return C<0, 1>(); } + T C02() const { return C<0, 2>(); } + T C03() const { return C<0, 3>(); } + T C04() const { return C<0, 4>(); } + T C05() const { return C<0, 5>(); } + T C06() const { return C<0, 6>(); } + + T C10() const { return C<1, 0>(); } + T C11() const { return C<1, 1>(); } + T C12() const { return C<1, 2>(); } + T C13() const { return C<1, 3>(); } + T C14() const { return C<1, 4>(); } + T C15() const { return C<1, 5>(); } + T C16() const { return C<1, 6>(); } + + T C20() const { return C<2, 0>(); } + T C21() const { return C<2, 1>(); } + T C22() const { return C<2, 2>(); } + T C23() const { return C<2, 3>(); } + T C24() const { return C<2, 4>(); } + T C25() const { return C<2, 5>(); } + T C26() const { return C<2, 6>(); } + + T C30() const { return C<3, 0>(); } + T C31() const { return C<3, 1>(); } + T C32() const { return C<3, 2>(); } + T C33() const { return C<3, 3>(); } + T C34() const { return C<3, 4>(); } + T C35() const { return C<3, 5>(); } + T C36() const { return C<3, 6>(); } + + T C40() const { return C<4, 0>(); } + T C41() const { return C<4, 1>(); } + T C42() const { return C<4, 2>(); } + T C43() const { return C<4, 3>(); } + T C44() const { return C<4, 4>(); } + T C45() const { return C<4, 5>(); } + T C46() const { return C<4, 6>(); } + + T C50() const { return C<5, 0>(); } + T C51() const { return C<5, 1>(); } + T C52() const { return C<5, 2>(); } + T C53() const { return C<5, 3>(); } + T C54() const { return C<5, 4>(); } + T C55() const { return C<5, 5>(); } + T C56() const { return C<5, 6>(); } + + T C60() const { return C<6, 0>(); } + T C61() const { return C<6, 1>(); } + T C62() const { return C<6, 2>(); } + T C63() const { return C<6, 3>(); } + T C64() const { return C<6, 4>(); } + T C65() const { return C<6, 5>(); } + T C66() const { return C<6, 6>(); } /// --------------------------------------------------------------------------------------------------------------------- /// Getters with 'Get' prefix /// Some of them involve calculations /// \brief Gets z position [cm] - DataT GetZ() const { return fZ; } + T GetZ() const { return fZ; } /// \brief Gets x position [cm] - DataT GetX() const { return fX; } + T GetX() const { return fX; } /// \brief Gets x position error [cm] - DataT GetXError() const { return sqrt(C00()); } + T GetXError() const { return sqrt(C00()); } /// \brief Gets y position [cm] - DataT GetY() const { return fY; } + T GetY() const { return fY; } /// \brief Gets y position error [cm] - DataT GetYError() const { return sqrt(C11()); } + T GetYError() const { return sqrt(C11()); } /// \brief Gets slope along x-axis - DataT GetTx() const { return fTx; } + T GetTx() const { return fTx; } /// \brief Gets error of slope along x-axis - DataT GetTxError() const { return sqrt(C22()); } + T GetTxError() const { return sqrt(C22()); } /// \brief Gets slope along y-axis - DataT GetTy() const { return fTy; } + T GetTy() const { return fTy; } /// \brief Gets error of slope along y-axis - DataT GetTyError() const { return sqrt(C33()); } + T GetTyError() const { return sqrt(C33()); } /// \brief Gets charge over momentum [ec/GeV] - DataT GetQp() const { return fQp; } + T GetQp() const { return fQp; } /// \brief Gets error of charge over momentum [ec/GeV] - DataT GetQpError() const { return sqrt(C44()); } + T GetQpError() const { return sqrt(C44()); } /// \brief Gets time [ns] - DataT GetTime() const { return fT; } + T GetTime() const { return fT; } /// \brief Gets time error [ns] - DataT GetTimeError() const { return sqrt(C55()); } + T GetTimeError() const { return sqrt(C55()); } /// \brief Gets inverse velocity [ns/cm] in downstream direction - DataT GetVi() const { return fVi; } + T GetVi() const { return fVi; } /// \brief Gets inverse velocity error [ns/cm] - DataT GetViError() const { return sqrt(C66()); } + T GetViError() const { return sqrt(C66()); } /// \brief Gets covariance matrix const CovMatrix_t& GetCovMatrix() const { return fCovMatrix; } @@ -251,55 +252,52 @@ namespace cbm::algo::kf /// \param i row /// \param j column /// \return covariance matrix element - DataT GetCovariance(int i, int j) const { return C(i, j); } + T GetCovariance(int i, int j) const { return C(i, j); } /// Gets Chi-square of track fit model - DataT GetChiSq() const { return fChiSq; } + T GetChiSq() const { return fChiSq; } /// Gets NDF of track fit model - DataT GetNdf() const { return fNdf; } + T GetNdf() const { return fNdf; } /// Gets Chi-square of time measurements - DataT GetChiSqTime() const { return fChiSqTime; } + T GetChiSqTime() const { return fChiSqTime; } /// Gets NDF of time measurements - DataT GetNdfTime() const { return fNdfTime; } + T GetNdfTime() const { return fNdfTime; } /// Gets charge - DataT GetCharge() const { return utils::iif(GetQp() > DataT(0.), DataT(1.), DataT(-1.)); } + T GetCharge() const { return utils::iif(GetQp() > T(0.), T(1.), T(-1.)); } /// \brief Gets azimuthal angle [rad] - DataT GetPhi() const { return atan2(GetTy(), GetTx()); } + T GetPhi() const { return atan2(GetTy(), GetTx()); } /// \brief Gets azimuthal angle error [rad] - DataT GetPhiError() const; + T GetPhiError() const; /// \brief Gets polar angle [rad] - DataT GetTheta() const { return atan(sqrt(GetTx() * GetTx() + GetTy() * GetTy())); } + T GetTheta() const { return atan(sqrt(GetTx() * GetTx() + GetTy() * GetTy())); } /// \brief Gets polar angle error [rad] - DataT GetThetaError() const; + T GetThetaError() const; /// Gets momentum [GeV/ec]. For the straight tracks returns 1.e4 [GeV/ec] - DataT GetP() const - { - return utils::iif(utils::fabs(GetQp()) > DataT(1.e-4), DataT(1.) / utils::fabs(GetQp()), DataT(1.e4)); - } + T GetP() const { return utils::iif(utils::fabs(GetQp()) > T(1.e-4), T(1.) / utils::fabs(GetQp()), T(1.e4)); } /// Gets z-component of the momentum [GeV/ec] - DataT GetPz() const { return GetP() / sqrt(DataT(1.) + GetTx() * GetTx() + GetTy() * GetTy()); } + T GetPz() const { return GetP() / sqrt(T(1.) + GetTx() * GetTx() + GetTy() * GetTy()); } /// Gets x-component of the momentum [GeV/ec] - DataT GetPx() const { return GetPz() * GetTx(); } + T GetPx() const { return GetPz() * GetTx(); } /// Gets y-component of the momentum [GeV/ec] - DataT GetPy() const { return GetPz() * GetTy(); } + T GetPy() const { return GetPz() * GetTy(); } /// Gets transverse momentum - DataT GetPt() const + T GetPt() const { - DataT t2 = GetTx() * GetTx() + GetTy() * GetTy(); - return GetP() * sqrt(t2 / (DataT(1.) + t2)); + T t2 = GetTx() * GetTx() + GetTy() * GetTy(); + return GetP() * sqrt(t2 / (T(1.) + t2)); } /// --------------------------------------------------------------------------------------------------------------------- @@ -307,40 +305,40 @@ namespace cbm::algo::kf /// /// \brief Sets z position [cm] - void SetZ(DataT v) { fZ = v; } + void SetZ(T v) { fZ = v; } /// \brief Sets x position [cm] - void SetX(DataT v) { fX = v; } + void SetX(T v) { fX = v; } /// \brief Sets y position [cm] - void SetY(DataT v) { fY = v; } + void SetY(T v) { fY = v; } /// \brief Sets slope along x-axis - void SetTx(DataT v) { fTx = v; } + void SetTx(T v) { fTx = v; } /// \brief Sets slope along y-axis - void SetTy(DataT v) { fTy = v; } + void SetTy(T v) { fTy = v; } /// \brief Sets charge over momentum [ec/GeV] - void SetQp(DataT v) { fQp = v; } + void SetQp(T v) { fQp = v; } /// \brief Sets time [ns] - void SetTime(DataT v) { fT = v; } + void SetTime(T v) { fT = v; } /// \brief Sets inverse velocity [ns/cm] - void SetVi(DataT v) { fVi = v; } + void SetVi(T v) { fVi = v; } /// \brief Sets Chi-square of track fit model - void SetChiSq(DataT v) { fChiSq = v; } + void SetChiSq(T v) { fChiSq = v; } /// \brief Sets NDF of track fit model - void SetNdf(DataT v) { fNdf = v; } + void SetNdf(T v) { fNdf = v; } /// \brief Sets Chi-square of time measurements - void SetChiSqTime(DataT v) { fChiSqTime = v; } + void SetChiSqTime(T v) { fChiSqTime = v; } /// \brief Sets NDF of time measurements - void SetNdfTime(DataT v) { fNdfTime = v; } + void SetNdfTime(T v) { fNdfTime = v; } /// \brief Sets covariance matrix void SetCovMatrix(const CovMatrix_t& val) { fCovMatrix = val; } @@ -349,7 +347,7 @@ namespace cbm::algo::kf /// \param i row /// \param j column /// \param val covariance matrix element - void SetCovariance(int i, int j, DataT val) + void SetCovariance(int i, int j, T val) { int ind = (j <= i) ? i * (1 + i) / 2 + j : j * (1 + j) / 2 + i; fCovMatrix[ind] = val; @@ -360,7 +358,7 @@ namespace cbm::algo::kf /// \param j column /// \param val matrix element template<int i, int j> - void SetCovariance(DataT val) + void SetCovariance(T val) { constexpr int ind = (j <= i) ? i * (1 + i) / 2 + j : j * (1 + j) / 2 + i; fCovMatrix[ind] = val; @@ -368,34 +366,34 @@ namespace cbm::algo::kf /// \brief Individual setters for covariance matrix elements /// - void SetC00(DataT val) { SetCovariance<0, 0>(val); } - void SetC10(DataT val) { SetCovariance<1, 0>(val); } - void SetC11(DataT val) { SetCovariance<1, 1>(val); } - void SetC20(DataT val) { SetCovariance<2, 0>(val); } - void SetC21(DataT val) { SetCovariance<2, 1>(val); } - void SetC22(DataT val) { SetCovariance<2, 2>(val); } - void SetC30(DataT val) { SetCovariance<3, 0>(val); } - void SetC31(DataT val) { SetCovariance<3, 1>(val); } - void SetC32(DataT val) { SetCovariance<3, 2>(val); } - void SetC33(DataT val) { SetCovariance<3, 3>(val); } - void SetC40(DataT val) { SetCovariance<4, 0>(val); } - void SetC41(DataT val) { SetCovariance<4, 1>(val); } - void SetC42(DataT val) { SetCovariance<4, 2>(val); } - void SetC43(DataT val) { SetCovariance<4, 3>(val); } - void SetC44(DataT val) { SetCovariance<4, 4>(val); } - void SetC50(DataT val) { SetCovariance<5, 0>(val); } - void SetC51(DataT val) { SetCovariance<5, 1>(val); } - void SetC52(DataT val) { SetCovariance<5, 2>(val); } - void SetC53(DataT val) { SetCovariance<5, 3>(val); } - void SetC54(DataT val) { SetCovariance<5, 4>(val); } - void SetC55(DataT val) { SetCovariance<5, 5>(val); } - void SetC60(DataT val) { SetCovariance<6, 0>(val); } - void SetC61(DataT val) { SetCovariance<6, 1>(val); } - void SetC62(DataT val) { SetCovariance<6, 2>(val); } - void SetC63(DataT val) { SetCovariance<6, 3>(val); } - void SetC64(DataT val) { SetCovariance<6, 4>(val); } - void SetC65(DataT val) { SetCovariance<6, 5>(val); } - void SetC66(DataT val) { SetCovariance<6, 6>(val); } + void SetC00(T val) { SetCovariance<0, 0>(val); } + void SetC10(T val) { SetCovariance<1, 0>(val); } + void SetC11(T val) { SetCovariance<1, 1>(val); } + void SetC20(T val) { SetCovariance<2, 0>(val); } + void SetC21(T val) { SetCovariance<2, 1>(val); } + void SetC22(T val) { SetCovariance<2, 2>(val); } + void SetC30(T val) { SetCovariance<3, 0>(val); } + void SetC31(T val) { SetCovariance<3, 1>(val); } + void SetC32(T val) { SetCovariance<3, 2>(val); } + void SetC33(T val) { SetCovariance<3, 3>(val); } + void SetC40(T val) { SetCovariance<4, 0>(val); } + void SetC41(T val) { SetCovariance<4, 1>(val); } + void SetC42(T val) { SetCovariance<4, 2>(val); } + void SetC43(T val) { SetCovariance<4, 3>(val); } + void SetC44(T val) { SetCovariance<4, 4>(val); } + void SetC50(T val) { SetCovariance<5, 0>(val); } + void SetC51(T val) { SetCovariance<5, 1>(val); } + void SetC52(T val) { SetCovariance<5, 2>(val); } + void SetC53(T val) { SetCovariance<5, 3>(val); } + void SetC54(T val) { SetCovariance<5, 4>(val); } + void SetC55(T val) { SetCovariance<5, 5>(val); } + void SetC60(T val) { SetCovariance<6, 0>(val); } + void SetC61(T val) { SetCovariance<6, 1>(val); } + void SetC62(T val) { SetCovariance<6, 2>(val); } + void SetC63(T val) { SetCovariance<6, 3>(val); } + void SetC64(T val) { SetCovariance<6, 4>(val); } + void SetC65(T val) { SetCovariance<6, 5>(val); } + void SetC66(T val) { SetCovariance<6, 6>(val); } ///--------------------------------------------------------------------------------------------------------------------- /// References to parameters to ease the math formulae @@ -403,49 +401,49 @@ namespace cbm::algo::kf /// /// \brief Reference to z position [cm] - DataT& Z() { return fZ; } + T& Z() { return fZ; } /// \brief Reference to x position [cm] - DataT& X() { return fX; } + T& X() { return fX; } /// \brief Reference to y position [cm] - DataT& Y() { return fY; } + T& Y() { return fY; } /// \brief Reference to slope along x-axis - DataT& Tx() { return fTx; } + T& Tx() { return fTx; } /// \brief Reference to slope along y-axis - DataT& Ty() { return fTy; } + T& Ty() { return fTy; } /// \brief Reference to charge over momentum [ec/GeV] - DataT& Qp() { return fQp; } + T& Qp() { return fQp; } /// \brief Reference to time [ns] - DataT& Time() { return fT; } + T& Time() { return fT; } /// \brief Reference to inverse velocity [ns/cm] - DataT& Vi() { return fVi; } + T& Vi() { return fVi; } /// \brief Reference to covariance matrix CovMatrix_t& CovMatrix() { return fCovMatrix; } /// \brief Reference to Chi-square of track fit model - DataT& ChiSq() { return fChiSq; } + T& ChiSq() { return fChiSq; } /// \brief Reference to NDF of track fit model - DataT& Ndf() { return fNdf; } + T& Ndf() { return fNdf; } /// \brief Reference to Chi-square of time measurements - DataT& ChiSqTime() { return fChiSqTime; } + T& ChiSqTime() { return fChiSqTime; } /// \brief Reference to NDF of time measurements - DataT& NdfTime() { return fNdfTime; } + T& NdfTime() { return fNdfTime; } /// \brief Get a reference to the covariance matrix element /// \param i row /// \param j column /// \return matrix element - DataT& C(int i, int j) + T& C(int i, int j) { int ind = (j <= i) ? i * (1 + i) / 2 + j : j * (1 + j) / 2 + i; return fCovMatrix[ind]; @@ -456,7 +454,7 @@ namespace cbm::algo::kf /// \param j column /// \param val matrix element template<int i, int j> - DataT& C() + T& C() { constexpr int ind = (j <= i) ? i * (1 + i) / 2 + j : j * (1 + j) / 2 + i; return fCovMatrix[ind]; @@ -464,61 +462,61 @@ namespace cbm::algo::kf /// \brief Individual references to covariance matrix elements /// - DataT& C00() { return C<0, 0>(); } - DataT& C01() { return C<0, 1>(); } - DataT& C02() { return C<0, 2>(); } - DataT& C03() { return C<0, 3>(); } - DataT& C04() { return C<0, 4>(); } - DataT& C05() { return C<0, 5>(); } - DataT& C06() { return C<0, 6>(); } - - DataT& C10() { return C<1, 0>(); } - DataT& C11() { return C<1, 1>(); } - DataT& C12() { return C<1, 2>(); } - DataT& C13() { return C<1, 3>(); } - DataT& C14() { return C<1, 4>(); } - DataT& C15() { return C<1, 5>(); } - DataT& C16() { return C<1, 6>(); } - - DataT& C20() { return C<2, 0>(); } - DataT& C21() { return C<2, 1>(); } - DataT& C22() { return C<2, 2>(); } - DataT& C23() { return C<2, 3>(); } - DataT& C24() { return C<2, 4>(); } - DataT& C25() { return C<2, 5>(); } - DataT& C26() { return C<2, 6>(); } - - DataT& C30() { return C<3, 0>(); } - DataT& C31() { return C<3, 1>(); } - DataT& C32() { return C<3, 2>(); } - DataT& C33() { return C<3, 3>(); } - DataT& C34() { return C<3, 4>(); } - DataT& C35() { return C<3, 5>(); } - DataT& C36() { return C<3, 6>(); } - - DataT& C40() { return C<4, 0>(); } - DataT& C41() { return C<4, 1>(); } - DataT& C42() { return C<4, 2>(); } - DataT& C43() { return C<4, 3>(); } - DataT& C44() { return C<4, 4>(); } - DataT& C45() { return C<4, 5>(); } - DataT& C46() { return C<4, 6>(); } - - DataT& C50() { return C<5, 0>(); } - DataT& C51() { return C<5, 1>(); } - DataT& C52() { return C<5, 2>(); } - DataT& C53() { return C<5, 3>(); } - DataT& C54() { return C<5, 4>(); } - DataT& C55() { return C<5, 5>(); } - DataT& C56() { return C<5, 6>(); } - - DataT& C60() { return C<6, 0>(); } - DataT& C61() { return C<6, 1>(); } - DataT& C62() { return C<6, 2>(); } - DataT& C63() { return C<6, 3>(); } - DataT& C64() { return C<6, 4>(); } - DataT& C65() { return C<6, 5>(); } - DataT& C66() { return C<6, 6>(); } + T& C00() { return C<0, 0>(); } + T& C01() { return C<0, 1>(); } + T& C02() { return C<0, 2>(); } + T& C03() { return C<0, 3>(); } + T& C04() { return C<0, 4>(); } + T& C05() { return C<0, 5>(); } + T& C06() { return C<0, 6>(); } + + T& C10() { return C<1, 0>(); } + T& C11() { return C<1, 1>(); } + T& C12() { return C<1, 2>(); } + T& C13() { return C<1, 3>(); } + T& C14() { return C<1, 4>(); } + T& C15() { return C<1, 5>(); } + T& C16() { return C<1, 6>(); } + + T& C20() { return C<2, 0>(); } + T& C21() { return C<2, 1>(); } + T& C22() { return C<2, 2>(); } + T& C23() { return C<2, 3>(); } + T& C24() { return C<2, 4>(); } + T& C25() { return C<2, 5>(); } + T& C26() { return C<2, 6>(); } + + T& C30() { return C<3, 0>(); } + T& C31() { return C<3, 1>(); } + T& C32() { return C<3, 2>(); } + T& C33() { return C<3, 3>(); } + T& C34() { return C<3, 4>(); } + T& C35() { return C<3, 5>(); } + T& C36() { return C<3, 6>(); } + + T& C40() { return C<4, 0>(); } + T& C41() { return C<4, 1>(); } + T& C42() { return C<4, 2>(); } + T& C43() { return C<4, 3>(); } + T& C44() { return C<4, 4>(); } + T& C45() { return C<4, 5>(); } + T& C46() { return C<4, 6>(); } + + T& C50() { return C<5, 0>(); } + T& C51() { return C<5, 1>(); } + T& C52() { return C<5, 2>(); } + T& C53() { return C<5, 3>(); } + T& C54() { return C<5, 4>(); } + T& C55() { return C<5, 5>(); } + T& C56() { return C<5, 6>(); } + + T& C60() { return C<6, 0>(); } + T& C61() { return C<6, 1>(); } + T& C62() { return C<6, 2>(); } + T& C63() { return C<6, 3>(); } + T& C64() { return C<6, 4>(); } + T& C65() { return C<6, 5>(); } + T& C66() { return C<6, 6>(); } ///--------------------------------------------------------------------------------------------------------------------- @@ -532,7 +530,7 @@ namespace cbm::algo::kf /// \param c44 Variance of charge over momentum [(ec/GeV)2] /// \param c55 Variance of time [ns2] /// \param c66 Variance of inverse velocity [1/c2] - void ResetErrors(DataT c00, DataT c11, DataT c22, DataT c33, DataT c44, DataT c55, DataT c66); + void ResetErrors(T c00, T c11, T c22, T c33, T c44, T c55, T c66); /// \brief Prints parameters to a string /// \param i Index of SIMD vector element (when i== -1, the entire vector is printed) @@ -584,8 +582,8 @@ namespace cbm::algo::kf /// \param ia Index of SIMD vector element of the current class /// \param Tb Other class /// \param ib Index of SIMD vector element of the other class - template<typename TdataB, bool TDoAllA, bool TDoAllB> - void CopyBase(const int ia, const TrackParamBase<TdataB>& Tb, const int ib); + template<typename T1, bool TDoAllA, bool TDoAllB> + void CopyBase(const int ia, const TrackParamBase<T1>& Tb, const int ib); private: @@ -607,41 +605,60 @@ namespace cbm::algo::kf 0., 0., 0., 0., 0., 0., 1.}; ///< covariance matrix // clang-format on - DataT fX{0.}; ///< x-position [cm] - DataT fY{0.}; ///< y-position [cm] - DataT fZ{0.}; ///< z-position [cm] - DataT fTx{0.}; ///< slope along x-axis - DataT fTy{0.}; ///< slope along y-axis - DataT fQp{0.}; ///< charge over momentum [ec/GeV] - DataT fT{0.}; ///< time [ns] - DataT fVi{0.}; ///< inverse velocity in downstream direction [ns/cm] + T fX{0.}; ///< x-position [cm] + T fY{0.}; ///< y-position [cm] + T fZ{0.}; ///< z-position [cm] + T fTx{0.}; ///< slope along x-axis + T fTy{0.}; ///< slope along y-axis + T fQp{0.}; ///< charge over momentum [ec/GeV] + T fT{0.}; ///< time [ns] + T fVi{0.}; ///< inverse velocity in downstream direction [ns/cm] - DataT fChiSq{0.}; ///< chi^2 of track fit, spatial measurements - DataT fNdf{0.}; ///< NDF of track fit, spatial measurements - DataT fChiSqTime{0.}; ///< chi^2 of track fit, time measurements - DataT fNdfTime{0.}; ///< NDF of track fit, time measurements + T fChiSq{0.}; ///< chi^2 of track fit, spatial measurements + T fNdf{0.}; ///< NDF of track fit, spatial measurements + T fChiSqTime{0.}; ///< chi^2 of track fit, time measurements + T fNdfTime{0.}; ///< NDF of track fit, time measurements - } _fvecalignment; // class TrackParamBase + }; // class TrackParamBase /// \class cbm::algo::ca:: TrackParamBaseScalar /// \brief Scalar version of TrackParamBase /// /// It contains extra methods that are difficult to implement in SIMD version - template<typename DataT> - class TrackParamBaseScalar : public TrackParamBase<DataT> { + template<typename T> + class TrackParamBaseScalar : public TrackParamBase<T> { public: + using TrackParamBase<T>::TrackParamBase; + + /// ---------------------------------------------------------------------------------------------------------------------using /// \brief Gets pseudo-rapidity - DataT GetEta() const { return -log(tan(this->GetTheta() * DataT(0.5))); } + T GetEta() const { return -log(tan(this->GetTheta() * T(0.5))); } + }; - } _fvecalignment; // --------------------------------------------------------------------------------------------------------------------- /// TrackParam classes of different types - typedef TrackParamBase<fvec> TrackParamV; - typedef TrackParamBaseScalar<fscal> TrackParamS; - typedef TrackParamBaseScalar<double> TrackParamD; + /// \class cbm::algo::kf::TrackParam + /// \brief Class represents fitted track parameters at given z + /// + template<typename T> + class TrackParam : public TrackParamBaseScalar<T> { + public: + using TrackParamBaseScalar<T>::TrackParamBaseScalar; + }; + + template<> + class TrackParam<fvec> : public TrackParamBase<fvec> { + public: + using TrackParamBase<fvec>::TrackParamBase; + }; + + + typedef TrackParam<fvec> TrackParamV; + typedef TrackParam<fscal> TrackParamS; + typedef TrackParam<double> TrackParamD; /// --------------------------------------------------------------------------------------------------------------------- @@ -650,42 +667,41 @@ namespace cbm::algo::kf /// --------------------------------------------------------------------------------------------------------------------- /// - template<typename DataT> - inline DataT TrackParamBase<DataT>::GetPhiError() const + template<typename T> + inline T TrackParamBase<T>::GetPhiError() const { // phi = atan( tx / ty ); } - DataT phiDdenom = GetTx() * GetTx() + GetTy() * GetTy(); - DataT phiDTx = -GetTy() / phiDdenom; // partial derivative of phi over Tx - DataT phiDTy = +GetTx() / phiDdenom; // partial derivative of phi over Ty + T phiDdenom = GetTx() * GetTx() + GetTy() * GetTy(); + T phiDTx = -GetTy() / phiDdenom; // partial derivative of phi over Tx + T phiDTy = +GetTx() / phiDdenom; // partial derivative of phi over Ty - DataT varTx = C22(); // variance of Tx - DataT varTy = C33(); // variance of Ty - DataT covTxTy = C32(); // covariance of Tx and Ty + T varTx = C22(); // variance of Tx + T varTy = C33(); // variance of Ty + T covTxTy = C32(); // covariance of Tx and Ty - DataT varPhi = phiDTx * phiDTx * varTx + phiDTy * phiDTy * varTy + DataT(2.) * phiDTx * phiDTy * covTxTy; + T varPhi = phiDTx * phiDTx * varTx + phiDTy * phiDTy * varTy + T(2.) * phiDTx * phiDTy * covTxTy; return sqrt(varPhi); } /// --------------------------------------------------------------------------------------------------------------------- /// - template<typename DataT> - inline DataT TrackParamBase<DataT>::GetThetaError() const + template<typename T> + inline T TrackParamBase<T>::GetThetaError() const { // theta = atan(sqrt( tx * tx + ty * ty) ) - DataT sumSqSlopes = GetTx() * GetTx() + GetTy() * GetTy(); - DataT thetaDdenom = sqrt(sumSqSlopes) * (DataT(1.) + sumSqSlopes); - DataT thetaDTx = GetTx() / thetaDdenom; - DataT thetaDTy = GetTy() / thetaDdenom; + T sumSqSlopes = GetTx() * GetTx() + GetTy() * GetTy(); + T thetaDdenom = sqrt(sumSqSlopes) * (T(1.) + sumSqSlopes); + T thetaDTx = GetTx() / thetaDdenom; + T thetaDTy = GetTy() / thetaDdenom; - DataT varTx = C22(); // variance of Tx - DataT varTy = C33(); // variance of Ty - DataT covTxTy = C32(); // covariance of Tx and Ty + T varTx = C22(); // variance of Tx + T varTy = C33(); // variance of Ty + T covTxTy = C32(); // covariance of Tx and Ty - DataT varTheta = - thetaDTx * thetaDTx * varTx + thetaDTy * thetaDTy * varTy + DataT(2.) * thetaDTx * thetaDTy * covTxTy; + T varTheta = thetaDTx * thetaDTx * varTx + thetaDTy * thetaDTy * varTy + T(2.) * thetaDTx * thetaDTy * covTxTy; return sqrt(varTheta); } @@ -724,9 +740,8 @@ namespace cbm::algo::kf // --------------------------------------------------------------------------------------------------------------------- // - template<typename DataT> - inline void TrackParamBase<DataT>::ResetErrors(DataT c00, DataT c11, DataT c22, DataT c33, DataT c44, DataT c55, - DataT c66) + template<typename T> + inline void TrackParamBase<T>::ResetErrors(T c00, T c11, T c22, T c33, T c44, T c55, T c66) { fCovMatrix.fill(0.); @@ -746,8 +761,8 @@ namespace cbm::algo::kf // --------------------------------------------------------------------------------------------------------------------- // - template<typename DataT> - inline void TrackParamBase<DataT>::InitVelocityRange(fscal minP) + template<typename T> + inline void TrackParamBase<T>::InitVelocityRange(fscal minP) { // initialise the velocity range with respect to the minimal momentum minP {GeV/c} using defs::ProtonMass; diff --git a/algo/kf/core/geo/KfFieldSlice.h b/algo/kf/core/geo/KfFieldSlice.h index 34a6e201ab..a626ca4a24 100644 --- a/algo/kf/core/geo/KfFieldSlice.h +++ b/algo/kf/core/geo/KfFieldSlice.h @@ -81,7 +81,7 @@ namespace cbm::algo::kf /// \brief Gets field value for the intersection with a straight track /// \param trackParam Track parameter set /// \return Field value [kG] - constexpr FieldValue<T> GetFieldValueForLine(const TrackParamBase<T>& trkPar) const + constexpr FieldValue<T> GetFieldValueForLine(const TrackParam<T>& trkPar) const { T dz = fZref - trkPar.GetZ(); return GetFieldValue(trkPar.GetX() + trkPar.GetTx() * dz, trkPar.GetY() + trkPar.GetTy() * dz); diff --git a/reco/KF/CbmKfTrackFitter.cxx b/reco/KF/CbmKfTrackFitter.cxx index b29c924b31..805a1cd82c 100644 --- a/reco/KF/CbmKfTrackFitter.cxx +++ b/reco/KF/CbmKfTrackFitter.cxx @@ -733,7 +733,7 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t) return ok; } -bool CbmKfTrackFitter::Smooth(kf::TrackParamBase<double>& t1, const kf::TrackParamBase<double>& t2) +bool CbmKfTrackFitter::Smooth(kf::TrackParamD& t1, const kf::TrackParamD& t2) { // TODO: move to the CaTrackFit class diff --git a/reco/KF/CbmKfTrackFitter.h b/reco/KF/CbmKfTrackFitter.h index 50f6182c3c..7449b49423 100644 --- a/reco/KF/CbmKfTrackFitter.h +++ b/reco/KF/CbmKfTrackFitter.h @@ -59,7 +59,7 @@ class CbmKfTrackFitter { double fZ{0.}; ///< Z coordinate - kf::TrackParamBase<double> fTrack{}; ///< fitted track + kf::TrackParamD fTrack{}; ///< fitted track /// == Material information (if present) @@ -152,7 +152,7 @@ class CbmKfTrackFitter { void AddMaterialEffects(FitNode& n, kf::FitDirection direction); // combine two tracks - bool Smooth(kf::TrackParamBase<double>& t1, const kf::TrackParamBase<double>& t2); + bool Smooth(kf::TrackParamD& t1, const kf::TrackParamD& t2); private: // input data arrays -- GitLab