From b7ab8942f910206f3707622fd5001f142765a62c Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Sun, 15 Oct 2023 16:14:13 +0000
Subject: [PATCH] Ca: rewrite L1*Measurement classes, move them to /algo
---
algo/ca/core/CMakeLists.txt | 6 +
algo/ca/core/data/CaMeasurementTime.cxx | 37 +++
algo/ca/core/data/CaMeasurementTime.h | 99 ++++++++
algo/ca/core/data/CaMeasurementU.cxx | 40 ++++
algo/ca/core/data/CaMeasurementU.h | 124 ++++++++++
algo/ca/core/data/CaMeasurementXy.cxx | 42 ++++
algo/ca/core/data/CaMeasurementXy.h | 213 ++++++++++++++++++
algo/ca/core/data/CaTrackParam.h | 24 +-
algo/ca/core/utils/CaUtils.h | 58 ++++-
reco/KF/ParticleFitter/CbmL1PFFitter.cxx | 68 +++---
reco/L1/CMakeLists.txt | 16 +-
reco/L1/CbmL1.cxx | 1 -
reco/L1/CbmL1.h | 7 -
reco/L1/L1Algo/L1Algo.h | 2 +-
reco/L1/L1Algo/L1CaTrackFinderSlice.cxx | 10 +-
reco/L1/L1Algo/L1Fit.cxx | 176 ++++++++-------
reco/L1/L1Algo/L1Fit.h | 25 +-
reco/L1/L1Algo/L1TrackFitter.cxx | 106 +++++----
reco/L1/L1Algo/L1TripletConstructor.cxx | 48 ++--
reco/L1/L1Algo/L1UMeasurementInfo.cxx | 39 ----
reco/L1/L1Algo/L1UMeasurementInfo.h | 42 ----
reco/L1/L1Algo/L1XYMeasurementInfo.cxx | 30 ---
reco/L1/L1Algo/L1XYMeasurementInfo.h | 51 -----
reco/L1/L1Algo/{ => inactive}/L1Event.cxx | 0
reco/L1/L1Algo/{ => inactive}/L1Event.h | 0
.../{ => inactive}/L1EventEfficiencies.h | 0
.../L1/L1Algo/{ => inactive}/L1EventMatch.cxx | 0
reco/L1/L1Algo/{ => inactive}/L1EventMatch.h | 0
reco/L1/L1Algo/{ => inactive}/L1MCEvent.cxx | 0
reco/L1/L1Algo/{ => inactive}/L1MCEvent.h | 0
30 files changed, 857 insertions(+), 407 deletions(-)
create mode 100644 algo/ca/core/data/CaMeasurementTime.cxx
create mode 100644 algo/ca/core/data/CaMeasurementTime.h
create mode 100644 algo/ca/core/data/CaMeasurementU.cxx
create mode 100644 algo/ca/core/data/CaMeasurementU.h
create mode 100644 algo/ca/core/data/CaMeasurementXy.cxx
create mode 100644 algo/ca/core/data/CaMeasurementXy.h
delete mode 100644 reco/L1/L1Algo/L1UMeasurementInfo.cxx
delete mode 100644 reco/L1/L1Algo/L1UMeasurementInfo.h
delete mode 100644 reco/L1/L1Algo/L1XYMeasurementInfo.cxx
delete mode 100644 reco/L1/L1Algo/L1XYMeasurementInfo.h
rename reco/L1/L1Algo/{ => inactive}/L1Event.cxx (100%)
rename reco/L1/L1Algo/{ => inactive}/L1Event.h (100%)
rename reco/L1/L1Algo/{ => inactive}/L1EventEfficiencies.h (100%)
rename reco/L1/L1Algo/{ => inactive}/L1EventMatch.cxx (100%)
rename reco/L1/L1Algo/{ => inactive}/L1EventMatch.h (100%)
rename reco/L1/L1Algo/{ => inactive}/L1MCEvent.cxx (100%)
rename reco/L1/L1Algo/{ => inactive}/L1MCEvent.h (100%)
diff --git a/algo/ca/core/CMakeLists.txt b/algo/ca/core/CMakeLists.txt
index 9f0a9e5372..5660908bf5 100644
--- a/algo/ca/core/CMakeLists.txt
+++ b/algo/ca/core/CMakeLists.txt
@@ -12,6 +12,9 @@ set(SRCS
${CMAKE_CURRENT_SOURCE_DIR}/data/CaTrackParam.cxx
${CMAKE_CURRENT_SOURCE_DIR}/data/CaGrid.cxx
${CMAKE_CURRENT_SOURCE_DIR}/data/CaTriplet.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/data/CaMeasurementU.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/data/CaMeasurementXy.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/data/CaMeasurementTime.cxx
${CMAKE_CURRENT_SOURCE_DIR}/pars/CaConfigReader.cxx
${CMAKE_CURRENT_SOURCE_DIR}/pars/CaField.cxx
${CMAKE_CURRENT_SOURCE_DIR}/pars/CaInitManager.cxx
@@ -55,6 +58,9 @@ install(
data/CaInputData.h
data/CaTrackParam.h
data/CaTrack.h
+ data/CaMeasurementU.h
+ data/CaMeasurementXy.h
+ data/CaMeasurementTime.h
pars/CaConfigReader.h
data/CaGridEntry.h
data/CaGrid.h
diff --git a/algo/ca/core/data/CaMeasurementTime.cxx b/algo/ca/core/data/CaMeasurementTime.cxx
new file mode 100644
index 0000000000..05e0bcc266
--- /dev/null
+++ b/algo/ca/core/data/CaMeasurementTime.cxx
@@ -0,0 +1,37 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+/// \file CaMeasurementTime.cxx
+/// \brief Implementation of the CaMeasurementTime class
+
+#include "CaMeasurementTime.h"
+
+#include <iomanip>
+#include <sstream> // for stringstream
+
+namespace cbm::algo::ca
+{
+ //----------------------------------------------------------------------------------------------------------------------
+ //
+ template<typename DataT>
+ std::string MeasurementTime<DataT>::ToString(int indentLevel) const
+ {
+ std::stringstream aStream {};
+ // TODO: possibly it is better to place the indentChar into ca::Parameters (S.Zharko)
+ constexpr char indentChar = '\t';
+ std::string indent(indentLevel, indentChar);
+ aStream << indent << "t: " << std::setw(12) << std::setfill(' ') << T() << '\n';
+ aStream << indent << "dt2: " << std::setw(12) << std::setfill(' ') << Dt2() << '\n';
+ aStream << indent << "ndfT: " << std::setw(12) << std::setfill(' ') << NdfT() << '\n';
+ return aStream.str();
+ }
+
+ ///----------------------------------------------------------------------------------------------------------------------
+ /// All necessary instantiations of the template class
+
+ template class cbm::algo::ca::MeasurementTime<fvec>;
+ template class cbm::algo::ca::MeasurementTime<float>;
+ template class cbm::algo::ca::MeasurementTime<double>;
+
+} // namespace cbm::algo::ca
diff --git a/algo/ca/core/data/CaMeasurementTime.h b/algo/ca/core/data/CaMeasurementTime.h
new file mode 100644
index 0000000000..f74597788e
--- /dev/null
+++ b/algo/ca/core/data/CaMeasurementTime.h
@@ -0,0 +1,99 @@
+/* Copyright (C) 2007-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov [committer], Sergei Zharko */
+
+/// \file CaMeasurementTime
+/// \brief Definition of the CaMeasurementTime class
+
+#pragma once // include this header only once per compilation unit
+
+
+#include <boost/serialization/access.hpp>
+
+#include <string>
+
+#include "CaConstants.h"
+#include "CaSimd.h"
+#include "CaUtils.h"
+
+namespace cbm::algo::ca
+{
+
+ /// \brief The class describes a time measurement
+ ///
+ /// The measurement has a finite resolution, i.e. the measurement is not a value, but a distribution
+ /// with a certain rms.
+ /// The measurement may be used in the chi2 calculation or not
+ /// The measurement may be a SIMD vector of values, when DataT is fvec type
+ ///
+ template<typename DataT>
+ class MeasurementTime {
+
+ public:
+ friend class boost::serialization::access;
+
+ /// default constructor
+ MeasurementTime() = default;
+
+ /// constructor
+ /// \param t time coordinate of the measurement
+ /// \param dt2 rms^2 of the time coordinate measurement
+ /// \param ndfT number of degrees of freedom for the time coordinate measurement
+ /// if ndfT == 1, the measurement is used in fit and in the chi2 calculation
+ /// if ndfT == 0, the measurement is used neither in fit nor in the chi2 calculation
+ MeasurementTime(DataT t, DataT dt2, DataT ndfT) : fT(t), fDt2(dt2), fNdfT(ndfT) {}
+
+ ///------------------------------
+ /// Setters and getters
+
+ void SetT(DataT t) { fT = t; }
+ void SetDt2(DataT dt2) { fDt2 = dt2; }
+ void SetNdfT(DataT ndfT) { fNdfT = ndfT; }
+
+ DataT T() const { return fT; }
+ DataT Dt2() const { return fDt2; }
+ DataT NdfT() const { return fNdfT; }
+
+
+ ///------------------------------
+ /// Methods for debugging
+
+ /// String representation of class contents
+ /// \param indentLevel number of indent characters in the output
+ std::string ToString(int indentLevel = 0) const;
+
+ /// Checks, if all fields are finite
+ bool IsFinite() const { return (utils::IsFinite(T()) && utils::IsFinite(Dt2()) && utils::IsFinite(NdfT())); }
+
+ /// Checks, if some fields are undefined
+ bool IsUndefined() const
+ {
+ return (utils::IsUndefined(T()) || utils::IsUndefined(Dt2()) || utils::IsUndefined(NdfT()));
+ }
+
+ /// Serialization function
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int)
+ {
+ ar& fT;
+ ar& fDt2;
+ ar& fNdfT;
+ }
+
+ private:
+ ///------------------------------
+ /// Data members
+
+ DataT fT {constants::Undef<DataT>}; ///< time coordinate of the measurement
+ DataT fDt2 {constants::Undef<DataT>}; ///< rms^2 of the time coordinate measurement
+
+ /// number of degrees of freedom (used for chi2 calculation)
+ /// if ndf == 1, the measurement is used in fit and in the chi2 calculation
+ /// if ndf == 0, the measurement is used neither in fit nor in the chi2 calculation
+
+ DataT fNdfT = constants::Undef<DataT>; ///< ndf for the time coordinate measurement
+
+ } _fvecalignment;
+
+} // namespace cbm::algo::ca
diff --git a/algo/ca/core/data/CaMeasurementU.cxx b/algo/ca/core/data/CaMeasurementU.cxx
new file mode 100644
index 0000000000..2d02b120d6
--- /dev/null
+++ b/algo/ca/core/data/CaMeasurementU.cxx
@@ -0,0 +1,40 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+/// \file CaMeasurementU.cxx
+/// \brief Implementation of the CaMeasurementU class
+
+#include "CaMeasurementU.h"
+
+#include <iomanip>
+#include <sstream>
+
+namespace cbm::algo::ca
+{
+
+ //----------------------------------------------------------------------------------------------------------------------
+ //
+ template<typename DataT>
+ std::string MeasurementU<DataT>::ToString(int indentLevel) const
+ {
+ std::stringstream aStream {};
+ // TODO: possibly it is better to place the indentChar into ca::Parameters (S.Zharko)
+ constexpr char indentChar = '\t';
+ std::string indent(indentLevel, indentChar);
+ aStream << indent << "cos(phi): " << std::setw(12) << std::setfill(' ') << CosPhi() << '\n';
+ aStream << indent << "sin(phi): " << std::setw(12) << std::setfill(' ') << SinPhi() << '\n';
+ aStream << indent << "u: " << std::setw(12) << std::setfill(' ') << U() << '\n';
+ aStream << indent << "du2: " << std::setw(12) << std::setfill(' ') << Du2() << '\n';
+ aStream << indent << "ndf: " << std::setw(12) << std::setfill(' ') << Ndf() << '\n';
+ return aStream.str();
+ }
+
+ ///----------------------------------------------------------------------------------------------------------------------
+ /// All necessary instantiations of the template class
+
+ template class MeasurementU<fvec>;
+ template class MeasurementU<float>;
+ template class MeasurementU<double>;
+
+} // namespace cbm::algo::ca
diff --git a/algo/ca/core/data/CaMeasurementU.h b/algo/ca/core/data/CaMeasurementU.h
new file mode 100644
index 0000000000..1a62925148
--- /dev/null
+++ b/algo/ca/core/data/CaMeasurementU.h
@@ -0,0 +1,124 @@
+/* Copyright (C) 2007-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov [committer], Igor Kulakov, Sergei Zharko */
+
+/// \file CaMeasurementU
+/// \brief Definition of the CaMeasurementU class
+
+#pragma once // include this header only once per compilation unit
+
+
+#include <boost/serialization/access.hpp>
+
+#include <string>
+
+#include "CaConstants.h"
+#include "CaSimd.h"
+#include "CaUtils.h"
+
+namespace cbm::algo::ca
+{
+ /// \brief The class describes a 1D - measurement U in XY coordinate system
+ ///
+ /// The measurement is defined as
+ ///
+ /// u = x * cos(phi) + y * sin(phi)
+ ///
+ /// where phi is the azimuthal angle of the measurement direction
+ /// and x, y are the coordinates of the measurement
+ /// The measurement may be scaled, i.e. cos^2 + sin^2 must not necessarily be equal to 1
+ /// The measurement has a finite resolution, i.e. the measurement is not a point, but a distribution
+ /// with a certain rms.
+ /// The measurement may be used in the chi2 calculation or not
+ /// The measurement may be a SIMD vector of values, when DataT is fvec type
+ ///
+ template<typename DataT>
+ class MeasurementU {
+
+ public:
+ friend class boost::serialization::access;
+
+ /// default constructor
+ MeasurementU() = default;
+
+ /// constructor
+ /// \param cosPhi cos(phi)
+ /// \param sinPhi sin(phi)
+ /// \param u measurement, u = x * cos(phi) + y * sin(phi)
+ /// \param du2 rms^2 of the measurement
+ /// \param ndf number of degrees of freedom (used for chi2 calculation)
+ /// if ndf == 1, the measurement is used in the chi2 calculation
+ /// if ndf == 0, the measurement is not used in the chi2 calculation
+ MeasurementU(DataT cosPhi, DataT sinPhi, DataT u, DataT du2, DataT ndf)
+ : fCosPhi(cosPhi)
+ , fSinPhi(sinPhi)
+ , fU(u)
+ , fDu2(du2)
+ , fNdf(ndf)
+ {
+ }
+
+ ///------------------------------
+ /// Setters and getters
+
+ void SetCosPhi(DataT cosPhi) { fCosPhi = cosPhi; }
+ void SetSinPhi(DataT sinPhi) { fSinPhi = sinPhi; }
+ void SetU(DataT u) { fU = u; }
+ void SetDu2(DataT du2) { fDu2 = du2; }
+ void SetNdf(DataT ndf) { fNdf = ndf; }
+
+ DataT CosPhi() const { return fCosPhi; }
+ DataT SinPhi() const { return fSinPhi; }
+ DataT U() const { return fU; }
+ DataT Du2() const { return fDu2; }
+ DataT Ndf() const { return fNdf; }
+
+ ///------------------------------
+ /// Methods for debugging
+
+ /// String representation of class contents
+ /// \param indentLevel number of indent characters in the output
+ std::string ToString(int indentLevel = 0) const;
+
+ /// Checks, if all fields are finite
+ bool IsFinite() const
+ {
+ return (utils::IsFinite(CosPhi()) && utils::IsFinite(SinPhi()) && utils::IsFinite(U()) && utils::IsFinite(Du2())
+ && utils::IsFinite(Ndf()));
+ }
+
+ /// Checks, if some fields are undefined
+ bool IsUndefined() const
+ {
+ return (utils::IsUndefined(CosPhi()) || utils::IsUndefined(SinPhi()) || utils::IsUndefined(U())
+ || utils::IsUndefined(Du2()) || utils::IsUndefined(Ndf()));
+ }
+
+ /// Serialization function
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int)
+ {
+ ar& fCosPhi;
+ ar& fSinPhi;
+ ar& fU;
+ ar& fDu2;
+ ar& fNdf;
+ }
+
+ private:
+ ///------------------------------
+ /// Data members
+
+ DataT fCosPhi = constants::Undef<DataT>; ///< cos(phi)
+ DataT fSinPhi = constants::Undef<DataT>; ///< sin(phi)
+ DataT fU = constants::Undef<DataT>; ///< measurement, u = x * cos(phi) + y * sin(phi)
+ DataT fDu2 = constants::Undef<DataT>; ///< rms^2 of the measurement
+
+ /// \brief number of degrees of freedom (used for chi2 calculation)
+ /// if ndf == 1, the measurement is used in the chi2 calculation
+ /// if ndf == 0, the measurement is not used in the chi2 calculation
+ DataT fNdf = constants::Undef<DataT>;
+ } _fvecalignment;
+
+} // namespace cbm::algo::ca
diff --git a/algo/ca/core/data/CaMeasurementXy.cxx b/algo/ca/core/data/CaMeasurementXy.cxx
new file mode 100644
index 0000000000..fea779e4e7
--- /dev/null
+++ b/algo/ca/core/data/CaMeasurementXy.cxx
@@ -0,0 +1,42 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+/// \file CaMeasurementXy.cxx
+/// \brief Implementation of the CaMeasurementXy class
+
+#include "CaMeasurementXy.h"
+
+#include <iomanip>
+#include <sstream> // for stringstream
+
+namespace cbm::algo::ca
+{
+
+ //----------------------------------------------------------------------------------------------------------------------
+ //
+ template<typename DataT>
+ std::string MeasurementXy<DataT>::ToString(int indentLevel) const
+ {
+ std::stringstream aStream {};
+ // TODO: possibly it is better to place the indentChar into ca::Parameters (S.Zharko)
+ constexpr char indentChar = '\t';
+ std::string indent(indentLevel, indentChar);
+ aStream << indent << "x: " << std::setw(12) << std::setfill(' ') << X() << '\n';
+ aStream << indent << "y: " << std::setw(12) << std::setfill(' ') << Y() << '\n';
+ aStream << indent << "dx2: " << std::setw(12) << std::setfill(' ') << Dx2() << '\n';
+ aStream << indent << "dy2: " << std::setw(12) << std::setfill(' ') << Dy2() << '\n';
+ aStream << indent << "dxy: " << std::setw(12) << std::setfill(' ') << Dxy() << '\n';
+ aStream << indent << "ndfX: " << std::setw(12) << std::setfill(' ') << NdfX() << '\n';
+ aStream << indent << "ndfY: " << std::setw(12) << std::setfill(' ') << NdfY() << '\n';
+ return aStream.str();
+ }
+
+ ///----------------------------------------------------------------------------------------------------------------------
+ /// All necessary instantiations of the template class
+
+ template class MeasurementXy<fvec>;
+ template class MeasurementXy<float>;
+ template class MeasurementXy<double>;
+
+} // namespace cbm::algo::ca
diff --git a/algo/ca/core/data/CaMeasurementXy.h b/algo/ca/core/data/CaMeasurementXy.h
new file mode 100644
index 0000000000..546c47d115
--- /dev/null
+++ b/algo/ca/core/data/CaMeasurementXy.h
@@ -0,0 +1,213 @@
+/* Copyright (C) 2007-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov [committer], Sergei Zharko */
+
+/// \file CaMeasurementXy
+/// \brief Definition of the CaMeasurementXy class
+
+#pragma once // include this header only once per compilation unit
+
+
+#include <boost/serialization/access.hpp>
+
+#include <string>
+
+#include "CaConstants.h"
+#include "CaSimd.h"
+#include "CaUtils.h"
+
+namespace cbm::algo::ca
+{
+
+ /// \brief The class describes a 2D - measurement (x, y) in XY coordinate system
+ ///
+ /// The measurement has a finite resolution, i.e. the measurement is not a point, but a distribution
+ /// with a certain rms.
+ /// The measurement components may be used in the chi2 calculation or not
+ /// The measurement may be a SIMD vector of values, when DataT is fvec type
+ ///
+ template<typename DataT>
+ class MeasurementXy {
+
+ public:
+ friend class boost::serialization::access;
+
+ /// default constructor
+ MeasurementXy() = default;
+
+ /// constructor
+ /// \param x x coordinate of the measurement
+ /// \param y y coordinate of the measurement
+ /// \param dx2 rms^2 of the x coordinate measurement
+ /// \param dy2 rms^2 of the y coordinate measurement
+ /// \param dxy covariance of the x and y coordinate measurements
+ /// \param ndfX number of degrees of freedom for the x coordinate measurement
+ /// if ndfX == 1, the measurement is used in fit and in the chi2 calculation
+ /// if ndfX == 0, the measurement is used in fit, but not used in the chi2 calculation
+ /// \param ndfY number of degrees of freedom for the y coordinate measurement
+ /// if ndfY == 1, the measurement is used in fit and in the chi2 calculation
+ /// if ndfY == 0, the measurement is used in fit, but not used in the chi2 calculation
+ MeasurementXy(DataT x, DataT y, DataT dx2, DataT dy2, DataT dxy, DataT ndfX, DataT ndfY)
+ : fX(x)
+ , fY(y)
+ , fDx2(dx2)
+ , fDy2(dy2)
+ , fDxy(dxy)
+ , fNdfX(ndfX)
+ , fNdfY(ndfY)
+ {
+ }
+
+
+ /// 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 DataTb>
+ void Set(const MeasurementXy<DataTb>& m)
+ {
+ CopyBase<DataTb, true, true>(0, m, 0);
+ }
+
+ /// Set all SIMD entries from one SIMD entry of the other class
+ /// It also works when DataT is scalar
+ void Set(const MeasurementXy<fvec>& m, const int im) { CopyBase<fvec, true, false>(0, m, im); }
+
+ /// Set one SIMD entry from one SIMD entry of the other class
+ /// It only works when DataT is fvec, TdataB is scalar
+ template<typename DataTb>
+ void SetOneEntry(const int i, const MeasurementXy<DataTb>& m)
+ {
+ CopyBase<DataTb, false, true>(i, m, 0);
+ }
+
+ /// Set one SIMD entry from one SIMD entry of the other class
+ /// It only works when DataT is fvec, TdataB is fvec
+ void SetOneEntry(const int i, const MeasurementXy<fvec>& m, const int im)
+ {
+ CopyBase<fvec, false, false>(i, m, im);
+ }
+
+
+ ///------------------------------
+ /// Setters and getters
+
+ void SetX(DataT x) { fX = x; }
+ void SetY(DataT y) { fY = y; }
+ void SetDx2(DataT dx2) { fDx2 = dx2; }
+ void SetDy2(DataT dy2) { fDy2 = dy2; }
+ void SetDxy(DataT dxy) { fDxy = dxy; }
+ void SetNdfX(DataT ndfX) { fNdfX = ndfX; }
+ void SetNdfY(DataT ndfY) { fNdfY = ndfY; }
+ void SetCov(DataT dx2, DataT dxy, DataT dy2)
+ {
+ fDx2 = dx2;
+ fDxy = dxy;
+ fDy2 = dy2;
+ }
+
+ DataT X() const { return fX; }
+ DataT Y() const { return fY; }
+ DataT Dx2() const { return fDx2; }
+ DataT Dy2() const { return fDy2; }
+ DataT Dxy() const { return fDxy; }
+ DataT NdfX() const { return fNdfX; }
+ DataT NdfY() const { return fNdfY; }
+
+ ///------------------------------
+ /// references, to ease assignment to SIMD vector components when DataT has fvec type
+
+ DataT& X() { return fX; }
+ DataT& Y() { return fY; }
+ DataT& Dx2() { return fDx2; }
+ DataT& Dy2() { return fDy2; }
+ DataT& Dxy() { return fDxy; }
+ DataT& NdfX() { return fNdfX; }
+ DataT& NdfY() { return fNdfY; }
+
+ ///------------------------------
+ /// Methods for debugging
+
+ /// String representation of class contents
+ /// \param indentLevel number of indent characters in the output
+ std::string ToString(int indentLevel = 0) const;
+
+ /// Checks, if all fields are finite
+ bool IsFinite() const
+ {
+ return (utils::IsFinite(X()) && utils::IsFinite(Y()) && utils::IsFinite(Dx2()) && utils::IsFinite(Dy2())
+ && utils::IsFinite(Dxy()) && utils::IsFinite(NdfX()) && utils::IsFinite(NdfY()));
+ }
+
+ /// Checks, if some fields are undefined
+ bool IsUndefined() const
+ {
+ return (utils::IsUndefined(X()) || utils::IsUndefined(Y()) || utils::IsUndefined(Dx2())
+ || utils::IsUndefined(Dy2()) || utils::IsUndefined(Dxy()) || utils::IsUndefined(NdfX())
+ || utils::IsUndefined(NdfY()));
+ }
+
+ /// Serialization function
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int)
+ {
+ ar& fX;
+ ar& fY;
+ ar& fDx2;
+ ar& fDy2;
+ ar& fDxy;
+ ar& fNdfX;
+ ar& fNdfY;
+ }
+
+ private:
+ /// \brief Copies all/one entries from the other class
+ /// \tparam TdataB Type of the other class
+ /// \tparam TDoAllA If true, all entries of the current class must be set
+ /// \tparam TDoAllB If true, all entries of the other class must be used
+ /// \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 MeasurementXy<TdataB>& Tb, const int ib);
+
+ private:
+ ///------------------------------
+ /// Data members
+
+ DataT fX {constants::Undef<DataT>}; ///< x coordinate of the measurement
+ DataT fY {constants::Undef<DataT>}; ///< y coordinate of the measurement
+ DataT fDx2 {constants::Undef<DataT>}; ///< rms^2 of the x coordinate measurement
+ DataT fDy2 {constants::Undef<DataT>}; ///< rms^2 of the y coordinate measurement
+ DataT fDxy {constants::Undef<DataT>}; ///< covariance of the x and y coordinate measurements
+
+ /// number of degrees of freedom (used for chi2 calculation)
+ /// if ndf == 1, the measurement is used in the chi2 calculation
+ /// if ndf == 0, the measurement is not used in the chi2 calculation
+
+ DataT fNdfX = constants::Undef<DataT>; ///< ndf for the x coordinate measurement
+ DataT fNdfY = constants::Undef<DataT>; ///< ndf for the y coordinate measurement
+
+ } _fvecalignment;
+
+
+ // ---------------------------------------------------------------------------------------------------------------------
+ //
+ template<typename TdataA>
+ template<typename TdataB, bool TDoAllA, bool TDoAllB>
+ inline void MeasurementXy<TdataA>::CopyBase(const int ia, const MeasurementXy<TdataB>& Tb, const int ib)
+ {
+ auto copy = [&](TdataA& a, const TdataB& b) {
+ utils::VecCopy<TdataA, TdataB, TDoAllA, TDoAllB>::CopyEntries(a, ia, b, ib);
+ };
+
+ copy(fX, Tb.X());
+ copy(fY, Tb.Y());
+ copy(fDx2, Tb.Dx2());
+ copy(fDy2, Tb.Dy2());
+ copy(fDxy, Tb.Dxy());
+ copy(fNdfX, Tb.NdfX());
+ copy(fNdfY, Tb.NdfY());
+ } // CopyBase
+
+} // namespace cbm::algo::ca
diff --git a/algo/ca/core/data/CaTrackParam.h b/algo/ca/core/data/CaTrackParam.h
index eb2b4e33b2..8f1cf3a425 100644
--- a/algo/ca/core/data/CaTrackParam.h
+++ b/algo/ca/core/data/CaTrackParam.h
@@ -635,28 +635,7 @@ namespace cbm::algo::ca
inline void TrackParamBase<TdataA>::CopyBase(const int ia, const TrackParamBase<TdataB>& Tb, const int ib)
{
auto copy = [&](TdataA& a, const TdataB& b) {
- if constexpr (TDoAllA || !std::is_same<TdataA, fvec>::value) {
- // set all entries or TdataA is a scalar
- if constexpr (TDoAllB || !std::is_same<TdataB, fvec>::value) {
- // set all entries or TdataB is a scalar
- a = b;
- }
- else {
- // TdataB is a vector
- a = b[ib];
- }
- }
- else {
- // set only one entry and TdataA is a vector
- if constexpr (TDoAllB || !std::is_same<TdataB, fvec>::value) {
- // set all entries or TdataB is a scalar
- a[ia] = b;
- }
- else {
- // TdataB is a vector
- a[ia] = b[ib];
- }
- }
+ utils::VecCopy<TdataA, TdataB, TDoAllA, TDoAllB>::CopyEntries(a, ia, b, ib);
};
copy(fX, Tb.GetX());
@@ -718,4 +697,3 @@ namespace cbm::algo::ca
}
} // namespace cbm::algo::ca
-
diff --git a/algo/ca/core/utils/CaUtils.h b/algo/ca/core/utils/CaUtils.h
index 6a0fa0ee70..bd2d6eb175 100644
--- a/algo/ca/core/utils/CaUtils.h
+++ b/algo/ca/core/utils/CaUtils.h
@@ -100,13 +100,69 @@ namespace cbm::algo::ca::utils
}
}
+
+ template<typename TdataA, typename TdataB, bool TDoAllA, bool TDoAllB>
+ class VecCopySpec {
+ public:
+ [[gnu::always_inline]] static void CopyEntries(TdataA& a, int ia, const TdataB& b, int ib);
+ };
+
+ template<typename TdataA, typename TdataB>
+ class VecCopySpec<TdataA, TdataB, true, true> {
+ public:
+ [[gnu::always_inline]] static void CopyEntries(TdataA& a, int, const TdataB& b, int) { a = b; }
+ };
+
+ template<typename TdataA, typename TdataB>
+ class VecCopySpec<TdataA, TdataB, true, false> {
+ public:
+ [[gnu::always_inline]] static void CopyEntries(TdataA& a, int, const TdataB& b, int ib) { a = b[ib]; }
+ };
+
+ template<typename TdataA, typename TdataB>
+ class VecCopySpec<TdataA, TdataB, false, true> {
+ public:
+ [[gnu::always_inline]] static void CopyEntries(TdataA& a, int ia, const TdataB& b, int) { a[ia] = b; }
+ };
+
+ template<typename TdataA, typename TdataB>
+ class VecCopySpec<TdataA, TdataB, false, false> {
+ public:
+ [[gnu::always_inline]] static void CopyEntries(TdataA& a, int ia, const TdataB& b, int ib) { a[ia] = b[ib]; }
+ };
+
+
+ /// \brief Copies all/one SIMD entries from one class to the other class
+ /// \details It helps to work with templates when classes might be SIMD or scalar
+ /// \details If one of the classes is scalar, its SIMD element index is ignored
+ ///
+ /// \tparam TdataA Type of the first class
+ /// \tparam TdataB Type of the second class
+ /// \tparam TDoAllA If true, all entries of the first class must be set
+ /// \tparam TDoAllB If true, all entries of the second class must be used
+ /// \param a First class
+ /// \param ia Index of SIMD vector element of the first class
+ /// \param b Second class
+ /// \param ib Index of SIMD vector element of the second class
+ template<typename TdataA, typename TdataB, bool TDoAllA, bool TDoAllB>
+ class VecCopy {
+ public:
+ static void CopyEntries(TdataA& a, int ia, const TdataB& b, int ib)
+ {
+ constexpr bool allA {TDoAllA || !std::is_same<TdataA, fvec>::value};
+ constexpr bool allB {TDoAllB || !std::is_same<TdataB, fvec>::value};
+ VecCopySpec<TdataA, TdataB, allA, allB>::CopyEntries(a, ia, b, ib);
+ }
+ };
+
+
/// \brief Checks, if a particular value lies within selected limits
/// \param name Name of parameters
/// \param value Value of parameter
/// \param lLimit Lower limit of parameter
/// \param uLimit Upper limit of parameter
template<typename T>
- bool CheckValueLimits(const std::string& name, T value, T lLimit, T uLimit)
+ inline bool CheckValueLimits(const std::string& name, T value, T lLimit, T uLimit)
{
if (value < lLimit || value > uLimit) {
LOG(error) << "parameter\033[1;32m" << name << "\033[0m (" << value << ") runs out of range: [" << lLimit << ','
diff --git a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
index 8b13196a4d..e8a0973bd1 100644
--- a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
+++ b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
@@ -139,16 +139,16 @@ inline int CbmL1PFFitter::GetStsStationIndex(const CbmStsHit* hit)
}
-void FilterFirst(L1Fit& fit, fvec& x, fvec& y, fvec& t, L1XYMeasurementInfo& cov, fvec& dt2)
+void FilterFirst(L1Fit& fit, ca::MeasurementXy<fvec>& mxy, fvec& t, fvec& dt2)
{
TrackParamV& tr = fit.Tr();
- tr.ResetErrors(cov.C00, cov.C11, 1., 1., 1., dt2, 1.e2);
- tr.C10() = cov.C10;
- tr.X() = x;
- tr.Y() = y;
+ tr.ResetErrors(mxy.Dx2(), mxy.Dy2(), 1., 1., 1., dt2, 1.e2);
+ tr.C10() = mxy.Dxy();
+ tr.X() = mxy.X();
+ tr.Y() = mxy.Y();
tr.Time() = t;
tr.Vi() = 0.;
- tr.Ndf() = -3.;
+ tr.Ndf() = mxy.NdfX() + mxy.NdfY() - fvec(5.);
tr.NdfTime() = -2.;
}
@@ -179,13 +179,13 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
fvec y[nHits];
fvec z[nHits];
fvec t[nHits];
- L1XYMeasurementInfo cov[nHits];
+ ca::MeasurementXy<fvec> mxy[nHits];
fvec dt2[nHits];
fmask w[nHits];
// fvec y_temp;
fvec x_first, y_first, t_first, x_last, y_last, t_last;
- L1XYMeasurementInfo cov_first, cov_last;
+ ca::MeasurementXy<fvec> mxy_first, mxy_last;
fvec dt2_first, dt2_last;
fvec z0, z1, z2, dz, z_start, z_end;
@@ -266,10 +266,14 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
z[ista][iVec] = hit->GetZ();
t[ista][iVec] = hit->GetTime();
- cov[ista].C00[iVec] = hit->GetDx() * hit->GetDx();
- cov[ista].C10[iVec] = hit->GetDxy();
- cov[ista].C11[iVec] = hit->GetDy() * hit->GetDy();
- dt2[ista][iVec] = hit->GetTimeError() * hit->GetTimeError();
+ mxy[ista].X()[iVec] = hit->GetX();
+ mxy[ista].Y()[iVec] = hit->GetY();
+ mxy[ista].Dx2()[iVec] = hit->GetDx() * hit->GetDx();
+ mxy[ista].Dy2()[iVec] = hit->GetDy() * hit->GetDy();
+ mxy[ista].Dxy()[iVec] = hit->GetDxy();
+ mxy[ista].NdfX()[iVec] = 1.;
+ mxy[ista].NdfY()[iVec] = 1.;
+ dt2[ista][iVec] = hit->GetTimeError() * hit->GetTimeError();
w[ista][iVec] = true;
@@ -279,23 +283,31 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
fB[ista].y[iVec] = fB_temp.y[iVec];
fB[ista].z[iVec] = fB_temp.z[iVec];
if (i == 0) {
- z_start[iVec] = z[ista][iVec];
- x_first[iVec] = x[ista][iVec];
- y_first[iVec] = y[ista][iVec];
- t_first[iVec] = t[ista][iVec];
- cov_first.C00[iVec] = cov[ista].C00[iVec];
- cov_first.C10[iVec] = cov[ista].C10[iVec];
- cov_first.C11[iVec] = cov[ista].C11[iVec];
- dt2_first[iVec] = dt2[ista][iVec];
+ z_start[iVec] = z[ista][iVec];
+ x_first[iVec] = x[ista][iVec];
+ y_first[iVec] = y[ista][iVec];
+ t_first[iVec] = t[ista][iVec];
+ mxy_first.X()[iVec] = mxy[ista].X()[iVec];
+ mxy_first.Y()[iVec] = mxy[ista].Y()[iVec];
+ mxy_first.Dx2()[iVec] = mxy[ista].Dx2()[iVec];
+ mxy_first.Dy2()[iVec] = mxy[ista].Dy2()[iVec];
+ mxy_first.Dxy()[iVec] = mxy[ista].Dxy()[iVec];
+ mxy_first.NdfX()[iVec] = mxy[ista].NdfX()[iVec];
+ mxy_first.NdfY()[iVec] = mxy[ista].NdfY()[iVec];
+ dt2_first[iVec] = dt2[ista][iVec];
}
if (i == nHitsTrack - 1) {
z_end[iVec] = z[ista][iVec];
x_last[iVec] = x[ista][iVec];
y_last[iVec] = y[ista][iVec];
t_last[iVec] = t[ista][iVec];
- cov_last.C00[iVec] = cov[ista].C00[iVec];
- cov_last.C10[iVec] = cov[ista].C10[iVec];
- cov_last.C11[iVec] = cov[ista].C11[iVec];
+ mxy_last.X()[iVec] = mxy[ista].X()[iVec];
+ mxy_last.Y()[iVec] = mxy[ista].Y()[iVec];
+ mxy_last.Dx2()[iVec] = mxy[ista].Dx2()[iVec];
+ mxy_last.Dy2()[iVec] = mxy[ista].Dy2()[iVec];
+ mxy_last.Dxy()[iVec] = mxy[ista].Dxy()[iVec];
+ mxy_last.NdfX()[iVec] = mxy[ista].NdfX()[iVec];
+ mxy_last.NdfY()[iVec] = mxy[ista].NdfY()[iVec];
dt2_last[iVec] = dt2[ista][iVec];
}
}
@@ -305,7 +317,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
int i = 0;
- FilterFirst(fit, x_first, y_first, t_first, cov_first, dt2_first);
+ FilterFirst(fit, mxy_first, t_first, dt2_first);
fit.SetQp0(fit.Tr().Qp());
z1 = z[i];
@@ -332,7 +344,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
fit.EnergyLossCorrection(radThick, -fvec::One());
fit.SetMask(initialised && w[i]);
- fit.FilterXY(cov[i], x[i], y[i]);
+ fit.FilterXY(mxy[i]);
fit.FilterTime(t[i], dt2[i], sta[i].timeInfo);
b2 = b1;
@@ -366,7 +378,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
i = nHits - 1;
- FilterFirst(fit, x_last, y_last, t_last, cov_last, dt2_last);
+ FilterFirst(fit, mxy_last, t_last, dt2_last);
z1 = z[i];
b1 = sta[i].fieldSlice.GetFieldValue(T.X(), T.Y());
@@ -393,7 +405,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
fit.EnergyLossCorrection(radThick, fvec::One());
fit.SetMask(initialised && w[i]);
- fit.FilterXY(cov[i], x[i], y[i]);
+ fit.FilterXY(mxy[i]);
fit.FilterTime(t[i], dt2[i], sta[i].timeInfo);
b2 = b1;
@@ -466,7 +478,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
int nStations = fNmvdStationsActive + fNstsStationsActive;
const ca::Station* sta = CbmL1::Instance()->fpAlgo->GetParameters()->GetStations().begin();
- fvec* zSta = new fvec[nStations];
+ fvec* zSta = new fvec[nStations];
for (int iSta = 0; iSta < nStations; iSta++) {
zSta[iSta] = sta[iSta].fZ;
}
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 43ce0dd84b..7c9fce58a5 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -42,18 +42,15 @@ set(SRCS
CbmCaMCModule.cxx
CbmCaTimeSliceReader.cxx
L1Algo/L1Fit.cxx
- L1Algo/L1MCEvent.cxx
CbmL1MCTrack.cxx
CbmL1Track.cxx
L1Algo/L1MaterialMonitor.cxx
- L1Algo/L1UMeasurementInfo.cxx
- L1Algo/L1XYMeasurementInfo.cxx
L1Algo/L1CloneMerger.cxx
L1Algo/utils/L1AlgoDraw.cxx
L1Algo/utils/L1AlgoEfficiencyPerformance.cxx
L1Algo/utils/L1AlgoPulls.cxx
L1Algo/utils/CaUvConverter.cxx
-
+
catools/CaToolsHitRecord.cxx
catools/CaToolsMCData.cxx
catools/CaToolsMCPoint.cxx
@@ -81,8 +78,9 @@ set(SRCS
)
set(NO_DICT_SRCS
- L1Algo/L1Event.cxx
- L1Algo/L1EventMatch.cxx
+ L1Algo/inactive/L1Event.cxx
+ L1Algo/inactive/L1MCEvent.cxx
+ L1Algo/inactive/L1EventMatch.cxx
)
set(HEADERS
@@ -185,14 +183,14 @@ add_dependencies(G__L1 KFPARTICLE)
install(FILES CbmL1Counters.h
L1Algo/L1Assert.h
- L1Algo/L1EventEfficiencies.h
- L1Algo/L1Event.h
- L1Algo/L1EventMatch.h
L1Algo/L1Utils.h
utils/CbmCaIdealHitProducer.h
utils/CbmCaIdealHitProducerDet.h
L1Algo/utils/CaUvConverter.h
qa/CbmCaInputQaBase.h
+# L1Algo/inactive/L1EventEfficiencies.h
+# L1Algo/inactive/L1Event.h
+# L1Algo/inactive/L1EventMatch.h
DESTINATION include
)
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index afbc212022..df22d57610 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -68,7 +68,6 @@
#include "CaHit.h"
#include "L1Algo/L1Algo.h"
#include "L1Algo/L1Assert.h"
-#include "L1Event.h"
using namespace cbm::algo; // TODO: remove this line
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 2d5bf5321d..b86e6567ce 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -58,14 +58,10 @@
#include "CaMonitor.h"
#include "CaVector.h"
#include "L1Algo/L1Algo.h"
-#include "L1EventEfficiencies.h"
#include "L1MaterialMonitor.h"
class L1Algo;
-class L1Event;
-class CbmL1ParticlesFinder;
class CbmL1MCTrack;
-class KFTopoPerformance;
class CbmMCDataObject;
class CbmEvent;
class TProfile2D;
@@ -685,9 +681,6 @@ private:
bool fExtrapolateToTheEndOfSTS {false};
- KFTopoPerformance* fTopoPerformance {nullptr};
- L1EventEfficiencies fEventEfficiency {}; // average efficiencies
-
std::vector<L1MaterialMonitor> fMaterialMonitor {}; ///< material monitoring
ClassDef(CbmL1, 0);
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index b3df979331..c4e6509aff 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -380,7 +380,7 @@ private:
bool fIsTargetField {false}; ///< is the magnetic field present at the target
ca::FieldValue fTargB _fvecalignment {}; // field in the target point (modifiable, do not touch!!)
- L1XYMeasurementInfo TargetXYInfo _fvecalignment {}; // target constraint [cm]
+ ca::MeasurementXy<fvec> TargetMeasurement _fvecalignment {}; // target constraint [cm]
int fGhostSuppression {1}; // NOTE: Should be equal to 0 in TRACKS_FROM_TRIPLETS mode!
diff --git a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
index 27e29a8b2c..7fd6185839 100644
--- a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
+++ b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
@@ -335,9 +335,13 @@ void L1Algo::CaTrackFinderSlice()
fIsTargetField = (fabs(fTargB.y[0]) > 0.001);
- TargetXYInfo.C00 = SigmaTargetX * SigmaTargetX;
- TargetXYInfo.C10 = 0;
- TargetXYInfo.C11 = SigmaTargetY * SigmaTargetY;
+ TargetMeasurement.SetX(fTargX);
+ TargetMeasurement.SetY(fTargY);
+ TargetMeasurement.SetDx2(SigmaTargetX * SigmaTargetX);
+ TargetMeasurement.SetDxy(0);
+ TargetMeasurement.SetDy2(SigmaTargetY * SigmaTargetY);
+ TargetMeasurement.SetNdfX(1);
+ TargetMeasurement.SetNdfY(1);
/// Set correction in order to take into account overlaping and iff z.
/// The reason is that low momentum tracks are too curved and goes not from target direction. That's why sort by hit_y/hit_z is not work idealy
diff --git a/reco/L1/L1Algo/L1Fit.cxx b/reco/L1/L1Algo/L1Fit.cxx
index 70d81b0a00..be576ff0c9 100644
--- a/reco/L1/L1Algo/L1Fit.cxx
+++ b/reco/L1/L1Algo/L1Fit.cxx
@@ -5,45 +5,47 @@
#include "L1Fit.h"
#include "CaHit.h"
+#include "CaMeasurementU.h"
+#include "CaMeasurementXy.h"
#include "CaStation.h"
#define cnst const fvec
-void L1Fit::FilterU(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2)
+void L1Fit::Filter1d(const ca::MeasurementU<fvec>& m)
{
fvec zeta, HCH;
fvec F0, F1, F2, F3, F4, F5, F6;
fvec K1, K2, K3, K4, K5, K6;
- zeta = info.cos_phi * fTr.X() + info.sin_phi * fTr.Y() - u;
+ zeta = m.CosPhi() * fTr.X() + m.SinPhi() * fTr.Y() - m.U();
// F = CH'
- F0 = info.cos_phi * fTr.C00() + info.sin_phi * fTr.C10();
- F1 = info.cos_phi * fTr.C10() + info.sin_phi * fTr.C11();
+ F0 = m.CosPhi() * fTr.C00() + m.SinPhi() * fTr.C10();
+ F1 = m.CosPhi() * fTr.C10() + m.SinPhi() * fTr.C11();
- HCH = (F0 * info.cos_phi + F1 * info.sin_phi);
+ HCH = (F0 * m.CosPhi() + F1 * m.SinPhi());
- F2 = info.cos_phi * fTr.C20() + info.sin_phi * fTr.C21();
- F3 = info.cos_phi * fTr.C30() + info.sin_phi * fTr.C31();
- F4 = info.cos_phi * fTr.C40() + info.sin_phi * fTr.C41();
- F5 = info.cos_phi * fTr.C50() + info.sin_phi * fTr.C51();
- F6 = info.cos_phi * fTr.C60() + info.sin_phi * fTr.C61();
+ F2 = m.CosPhi() * fTr.C20() + m.SinPhi() * fTr.C21();
+ F3 = m.CosPhi() * fTr.C30() + m.SinPhi() * fTr.C31();
+ F4 = m.CosPhi() * fTr.C40() + m.SinPhi() * fTr.C41();
+ F5 = m.CosPhi() * fTr.C50() + m.SinPhi() * fTr.C51();
+ F6 = m.CosPhi() * fTr.C60() + m.SinPhi() * fTr.C61();
- const fmask maskDoFilter = (HCH < sigma2 * 16.f);
+ const fmask maskDoFilter = (HCH < m.Du2() * 16.f);
//cnst maskDoFilter = _f32vec4_true;
// correction to HCH is needed for the case when sigma2 is so small
// with respect to HCH that it disappears due to the roundoff error
//
- fvec wi = fMaskF / (sigma2 + fvec(1.0000001) * HCH);
- fvec zetawi = fMaskF * zeta / (iif(maskDoFilter, sigma2, fvec::Zero()) + HCH);
+ fvec wi = fMaskF / (m.Du2() + fvec(1.0000001) * HCH);
+ fvec zetawi = fMaskF * zeta / (iif(maskDoFilter, m.Du2(), fvec::Zero()) + HCH);
//wi = iif(wi > fvec::Zero(), wi, fvec::Zero());
- wi = iif(sigma2 > fvec::Zero(), wi, fvec::Zero());
+ wi = iif(m.Du2() > fvec::Zero(), wi, fvec::Zero());
// fTr.ChiSq() += iif( maskDoFilter, zeta * zetawi, fvec::Zero() );
- fTr.ChiSq() += zeta * zeta * wi;
- fTr.Ndf() += fMaskF;
+ fTr.ChiSq() += m.Ndf() * zeta * zeta * wi;
+ fTr.Ndf() += m.Ndf() * fMaskF;
K1 = F1 * wi;
K2 = F2 * wi;
@@ -181,24 +183,31 @@ void L1Fit::FilterTime(cnst& t, cnst& dt2, cnst& timeInfo)
}
-void L1Fit::FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y)
+void L1Fit::FilterXY(const ca::MeasurementXy<fvec>& mxy)
{
{
- L1UMeasurementInfo xInfo;
- xInfo.cos_phi = fvec::One();
- xInfo.sin_phi = fvec::Zero();
-
- L1UMeasurementInfo uInfo;
- uInfo.cos_phi = -info.C10 / info.C00;
- uInfo.sin_phi = fvec::One();
- fvec u = uInfo.cos_phi * x + y;
- fvec u2 = info.C11 - info.C10 * info.C10 / info.C00;
-
- FilterU(xInfo, x, info.C00);
- FilterU(uInfo, u, u2);
+ ca::MeasurementU<fvec> mx;
+ mx.SetCosPhi(fvec::One());
+ mx.SetSinPhi(fvec::Zero());
+ mx.SetU(mxy.X());
+ mx.SetDu2(mxy.Dx2());
+ mx.SetNdf(fvec::One());
+
+ ca::MeasurementU<fvec> mu;
+ mu.SetCosPhi(-mxy.Dxy() / mxy.Dx2());
+ mu.SetSinPhi(fvec::One());
+ mu.SetU(mu.CosPhi() * mxy.X() + mxy.Y());
+ mu.SetDu2(mxy.Dy2() - mxy.Dxy() * mxy.Dxy() / mxy.Dx2());
+ mu.SetNdf(fvec::One());
+
+ Filter1d(mx);
+ Filter1d(mu);
return;
}
+ //----------------------------------------------------------------------------------------------
+ // the other way: filter 2 dimensions at once
+
cnst TWO(2.);
fvec zeta0, zeta1, S00, S10, S11, si;
@@ -207,8 +216,8 @@ void L1Fit::FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y)
fvec K00, K10, K20, K30, K40, K50, K60;
fvec K01, K11, K21, K31, K41, K51, K61;
- zeta0 = fTr.X() - x;
- zeta1 = fTr.Y() - y;
+ zeta0 = fTr.X() - mxy.X();
+ zeta1 = fTr.Y() - mxy.Y();
// F = CH'
F00 = fTr.C00();
@@ -227,9 +236,9 @@ void L1Fit::FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y)
F51 = fTr.C51();
F61 = fTr.C61();
- S00 = F00 + info.C00;
- S10 = F10 + info.C10;
- S11 = F11 + info.C11;
+ S00 = F00 + mxy.Dx2();
+ S10 = F10 + mxy.Dxy();
+ S11 = F11 + mxy.Dy2();
si = 1.f / (S00 * S11 - S10 * S10);
fvec S00tmp = S00;
@@ -307,16 +316,19 @@ void L1Fit::FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y)
void L1Fit::FilterHit(const ca::Station& sta, const ca::Hit& hit)
{
- L1XYMeasurementInfo info;
- info.C00 = hit.dX2();
- info.C10 = hit.dXY();
- info.C11 = hit.dY2();
- FilterXY(info, hit.X(), hit.Y());
+ ca::MeasurementXy<fvec> m;
+ m.SetDx2(hit.dX2());
+ m.SetDy2(hit.dY2());
+ m.SetDxy(hit.dXY());
+ m.SetX(hit.X());
+ m.SetY(hit.Y());
+ m.SetNdfX(fvec::One());
+ m.SetNdfY(fvec::One());
+ FilterXY(m);
FilterTime(hit.T(), hit.dT2(), sta.timeInfo);
}
-void L1Fit::FilterExtrapolatedXY(cnst& x, cnst& y, const L1XYMeasurementInfo& info, cnst& extrX, cnst& extrY,
- cnst Jx[6], cnst Jy[6])
+void L1Fit::FilterExtrapolatedXY(const ca::MeasurementXy<fvec>& m, cnst& extrX, cnst& extrY, cnst Jx[6], cnst Jy[6])
{
// 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
@@ -329,8 +341,8 @@ void L1Fit::FilterExtrapolatedXY(cnst& x, cnst& y, const L1XYMeasurementInfo& in
//zeta0 = T.x + Jx[2]*T.Tx() + Jx[3]*T.Ty() + Jx[4]*T.qp - x;
//zeta1 = T.y + Jy[2]*T.Tx() + Jy[3]*T.Ty() + Jy[4]*T.qp - y;
- fvec zeta0 = extrX - x;
- fvec zeta1 = extrY - y;
+ fvec zeta0 = extrX - m.X();
+ fvec zeta1 = extrY - m.Y();
// H = 1 0 Jx[2] Jx[3] Jx[4] 0
// 0 1 Jy[2] Jy[3] Jy[4] 0
@@ -347,9 +359,9 @@ void L1Fit::FilterExtrapolatedXY(cnst& x, cnst& y, const L1XYMeasurementInfo& in
fvec F40 = Jx[4] * T.C44();
fvec F41 = Jy[4] * T.C44();
- fvec S00 = info.C00 + F00 + Jx[2] * F20 + Jx[3] * F30 + Jx[4] * F40;
- fvec S10 = info.C10 + F10 + Jy[2] * F20 + Jy[3] * F30 + Jy[4] * F40;
- fvec S11 = info.C11 + F11 + Jy[2] * F21 + Jy[3] * F31 + Jy[4] * F41;
+ fvec S00 = m.Dx2() + F00 + Jx[2] * F20 + Jx[3] * F30 + Jx[4] * F40;
+ fvec S10 = m.Dxy() + F10 + Jy[2] * F20 + Jy[3] * F30 + Jy[4] * F40;
+ fvec S11 = m.Dy2() + F11 + Jy[2] * F21 + Jy[3] * F31 + Jy[4] * F41;
fvec si = fvec(1.) / (S00 * S11 - S10 * S10);
@@ -359,7 +371,7 @@ void L1Fit::FilterExtrapolatedXY(cnst& x, cnst& y, const L1XYMeasurementInfo& in
S11 = si * S00tmp;
T.ChiSq() += zeta0 * zeta0 * S00 + fvec(2.) * zeta0 * zeta1 * S10 + zeta1 * zeta1 * S11;
- T.Ndf() += fvec(2.);
+ T.Ndf() += m.NdfX() + m.NdfY();
fvec K00 = F00 * S00 + F01 * S10;
fvec K01 = F00 * S10 + F01 * S11;
@@ -1700,14 +1712,13 @@ void L1Fit::GetExtrapolatedXYline(cnst& z, const ca::FieldRegion& F, fvec& extrX
}
-void L1Fit::AddTargetToLine(cnst& targX, cnst& targY, cnst& targZ, const L1XYMeasurementInfo& targXYInfo,
- const ca::FieldRegion& F)
+void L1Fit::AddTargetToLine(cnst& targZ, const ca::MeasurementXy<fvec>& targXY, const ca::FieldRegion& F)
{
// Add the target constraint to a straight line track
fvec eX, eY, Jx[6], Jy[6];
GetExtrapolatedXYline(targZ, F, eX, eY, Jx, Jy);
- FilterExtrapolatedXY(targX, targY, targXYInfo, eX, eY, Jx, Jy);
+ FilterExtrapolatedXY(targXY, eX, eY, Jx, Jy);
}
fvec L1Fit::ApproximateBetheBloch(cnst& bg2)
@@ -1803,26 +1814,22 @@ fvec L1Fit::ApproximateBetheBloch(cnst& bg2, cnst& kp0, cnst& kp1, cnst& kp2, cn
}
-void L1Fit::FilterChi2XYC00C10C11(const L1UMeasurementInfo& info, fvec& x, fvec& y, fvec& C00, fvec& C10, fvec& C11,
- fvec& chi2, cnst& u, cnst& du2)
+void L1Fit::FilterChi2XYC00C10C11(const ca::MeasurementU<fvec>& m, fvec& x, fvec& y, fvec& C00, fvec& C10, fvec& C11,
+ fvec& chi2)
{
- fvec wi, zeta, zetawi, HCH;
- fvec F0, F1;
- fvec K1;
-
- zeta = info.cos_phi * x + info.sin_phi * y - u;
+ fvec zeta = m.CosPhi() * x + m.SinPhi() * y - m.U();
// F = CH'
- F0 = info.cos_phi * C00 + info.sin_phi * C10;
- F1 = info.cos_phi * C10 + info.sin_phi * C11;
+ fvec F0 = m.CosPhi() * C00 + m.SinPhi() * C10;
+ fvec F1 = m.CosPhi() * C10 + m.SinPhi() * C11;
- HCH = (F0 * info.cos_phi + F1 * info.sin_phi);
+ fvec HCH = (F0 * m.CosPhi() + F1 * m.SinPhi());
- wi = fvec(1.) / (du2 + HCH);
- zetawi = zeta * wi;
- chi2 += zeta * zetawi;
+ fvec wi = fvec(1.) / (m.Du2() + HCH);
+ fvec zetawi = zeta * wi;
+ chi2 += m.Ndf() * zeta * zetawi;
- K1 = F1 * wi;
+ fvec K1 = F1 * wi;
x -= F0 * zetawi;
y -= F1 * zetawi;
@@ -1833,42 +1840,43 @@ void L1Fit::FilterChi2XYC00C10C11(const L1UMeasurementInfo& info, fvec& x, fvec&
}
-void L1Fit::FilterChi2(const L1UMeasurementInfo& info, cnst& x, cnst& y, cnst& C00, cnst& C10, cnst& C11, fvec& chi2,
- cnst& u, cnst& du2)
+void L1Fit::FilterChi2(const ca::MeasurementU<fvec>& m, cnst& x, cnst& y, cnst& C00, cnst& C10, cnst& C11, fvec& chi2)
{
- fvec zeta, HCH;
- fvec F0, F1;
-
- zeta = info.cos_phi * x + info.sin_phi * y - u;
+ fvec zeta = m.CosPhi() * x + m.SinPhi() * y - m.U();
// F = CH'
- F0 = info.cos_phi * C00 + info.sin_phi * C10;
- F1 = info.cos_phi * C10 + info.sin_phi * C11;
+ fvec F0 = m.CosPhi() * C00 + m.SinPhi() * C10;
+ fvec F1 = m.CosPhi() * C10 + m.SinPhi() * C11;
- HCH = (F0 * info.cos_phi + F1 * info.sin_phi);
+ fvec HCH = (F0 * m.CosPhi() + F1 * m.SinPhi());
- chi2 += zeta * zeta / (du2 + HCH);
+ chi2 += m.Ndf() * zeta * zeta / (m.Du2() + HCH);
}
+
std::tuple<fvec, fvec> L1Fit::GetChi2XChi2U(fvec xm, fvec ym, fvec V00, fvec V10, fvec V11, fvec xt, fvec yt, fvec C00,
fvec C10, fvec C11)
{
- L1UMeasurementInfo xInfo;
- xInfo.cos_phi = fvec::One();
- xInfo.sin_phi = fvec::Zero();
+ ca::MeasurementU<fvec> measX;
+ measX.SetCosPhi(fvec::One());
+ measX.SetSinPhi(fvec::Zero());
+ measX.SetU(xm);
+ measX.SetDu2(V00);
+ measX.SetNdf(fvec::One());
- L1UMeasurementInfo uInfo;
- uInfo.cos_phi = -V10 / V00;
- uInfo.sin_phi = fvec::One();
- fvec u = uInfo.cos_phi * xm + ym;
- fvec du2 = V11 - V10 * V10 / V00;
+ ca::MeasurementU<fvec> measU;
+ measU.SetCosPhi(-V10 / V00);
+ measU.SetSinPhi(fvec::One());
+ measU.SetU(measU.CosPhi() * xm + ym);
+ measU.SetDu2(V11 - V10 * V10 / V00);
+ measU.SetNdf(fvec::One());
fvec chi2x = 0.f;
fvec chi2u = 0.f;
- FilterChi2XYC00C10C11(xInfo, xt, yt, C00, C10, C11, chi2x, xm, V00);
- FilterChi2(uInfo, xt, yt, C00, C10, C11, chi2u, u, du2);
+ FilterChi2XYC00C10C11(measX, xt, yt, C00, C10, C11, chi2x);
+ FilterChi2(measU, xt, yt, C00, C10, C11, chi2u);
return std::tuple<fvec, fvec>(chi2x, chi2u);
}
diff --git a/reco/L1/L1Algo/L1Fit.h b/reco/L1/L1Algo/L1Fit.h
index e0386ee991..94ca661969 100644
--- a/reco/L1/L1Algo/L1Fit.h
+++ b/reco/L1/L1Algo/L1Fit.h
@@ -12,14 +12,15 @@
#define L1Fit_h
#include "CaField.h"
+#include "CaMeasurementU.h"
+#include "CaMeasurementXy.h"
+#include "CaSimd.h"
#include "CaTrackParam.h"
-#include "L1Def.h"
-#include "L1UMeasurementInfo.h"
-#include "L1XYMeasurementInfo.h"
namespace
{
using namespace cbm::algo;
+ using namespace cbm::algo::ca;
} // namespace
namespace cbm::algo::ca
@@ -95,16 +96,15 @@ public:
/// get the particle mass
fvec GetMaxExtrapolationStep() const { return fMaxExtraplationStep; }
- void FilterU(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2);
- void FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y);
+ void Filter1d(const ca::MeasurementU<fvec>& m);
+ void FilterXY(const ca::MeasurementXy<fvec>& m);
void FilterTime(cnst& t, cnst& dt2, cnst& timeInfo);
void FilterHit(const ca::Station& s, const ca::Hit& h);
void FilterVi(fvec vi);
- void FilterExtrapolatedXY(cnst& x, cnst& y, const L1XYMeasurementInfo& info, cnst& extrX, cnst& extrY, cnst Jx[6],
- cnst Jy[6]);
+ void FilterExtrapolatedXY(const ca::MeasurementXy<fvec>& m, cnst& extrX, cnst& extrY, cnst Jx[6], cnst Jy[6]);
void MeasureVelocityWithQp();
@@ -136,18 +136,17 @@ public:
void ExtrapolateC10Line(cnst& z_out, fvec& C10) const;
- void AddTargetToLine(cnst& targX, cnst& targY, cnst& targZ, const L1XYMeasurementInfo& targXYInfo,
- const ca::FieldRegion& F);
+ void AddTargetToLine(cnst& targZ, const ca::MeasurementXy<fvec>& targXYInfo, const ca::FieldRegion& F);
static fvec ApproximateBetheBloch(cnst& bg2);
static fvec ApproximateBetheBloch(cnst& bg2, cnst& kp0, cnst& kp1, cnst& kp2, cnst& kp3, cnst& kp4);
- static void FilterChi2XYC00C10C11(const L1UMeasurementInfo& info, fvec& x, fvec& y, fvec& C00, fvec& C10, fvec& C11,
- fvec& chi2, cnst& u, cnst& du2);
+ static void FilterChi2XYC00C10C11(const ca::MeasurementU<fvec>& m, fvec& x, fvec& y, fvec& C00, fvec& C10, fvec& C11,
+ fvec& chi2);
- static void FilterChi2(const L1UMeasurementInfo& info, cnst& x, cnst& y, cnst& C00, cnst& C10, cnst& C11, fvec& chi2,
- cnst& u, cnst& du2);
+ static void FilterChi2(const ca::MeasurementU<fvec>& m, cnst& x, cnst& y, cnst& C00, cnst& C10, cnst& C11,
+ fvec& chi2);
static std::tuple<fvec, fvec> GetChi2XChi2U(fvec xm, fvec ym, fvec V00, fvec V10, fvec V11, fvec xt, fvec yt,
fvec C00, fvec C10, fvec C11);
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 82917095cc..ecbd80131d 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -49,9 +49,9 @@ void L1Algo::L1KFTrackFitter()
// Spatial-time position of a hit vs. station and track in the portion
- fvec x[constants::size::MaxNstations]; // Hit position along the x-axis [cm]
- fvec y[constants::size::MaxNstations]; // Hit position along the y-axis [cm]
- L1XYMeasurementInfo cov_xy[constants::size::MaxNstations]; // Covariance matrix for x,y
+ fvec x[constants::size::MaxNstations]; // Hit position along the x-axis [cm]
+ fvec y[constants::size::MaxNstations]; // Hit position along the y-axis [cm]
+ ca::MeasurementXy<fvec> mxy[constants::size::MaxNstations]; // Covariance matrix for x,y
fvec z[constants::size::MaxNstations]; // Hit position along the z-axis (precised) [cm]
@@ -60,7 +60,7 @@ void L1Algo::L1KFTrackFitter()
fvec x_first;
fvec y_first;
- L1XYMeasurementInfo cov_xy_fst;
+ ca::MeasurementXy<fvec> mxy_first;
fvec time_first;
fvec wtime_first;
@@ -68,7 +68,7 @@ void L1Algo::L1KFTrackFitter()
fvec x_last;
fvec y_last;
- L1XYMeasurementInfo cov_xy_lst;
+ ca::MeasurementXy<fvec> mxy_last;
fvec time_last;
fvec wtime_last;
@@ -91,10 +91,8 @@ void L1Algo::L1KFTrackFitter()
fvec ZSta[constants::size::MaxNstations];
for (int ista = 0; ista < nStations; ista++) {
- ZSta[ista] = sta[ista].fZ;
- cov_xy[ista].C00 = 1.f;
- cov_xy[ista].C10 = 0.f;
- cov_xy[ista].C11 = 1.f;
+ ZSta[ista] = sta[ista].fZ;
+ mxy[ista].SetCov(1., 0., 1.);
}
unsigned short N_vTracks = fSliceRecoTracks.size(); // number of tracks processed per one SSE register
@@ -144,41 +142,53 @@ void L1Algo::L1KFTrackFitter()
if (!sta[ista].timeInfo) { dt2[ista][iVec] = 1.e4; }
z[ista][iVec] = hit.Z();
fB_temp = sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista]);
- cov_xy[ista].C00[iVec] = hit.dX2();
- cov_xy[ista].C10[iVec] = hit.dXY();
- cov_xy[ista].C11[iVec] = hit.dY2();
+ mxy[ista].X()[iVec] = hit.X();
+ mxy[ista].Y()[iVec] = hit.Y();
+ mxy[ista].Dx2()[iVec] = hit.dX2();
+ mxy[ista].Dy2()[iVec] = hit.dY2();
+ mxy[ista].Dxy()[iVec] = hit.dXY();
+ mxy[ista].NdfX()[iVec] = 1.;
+ mxy[ista].NdfY()[iVec] = 1.;
fB[ista].x[iVec] = fB_temp.x[iVec];
fB[ista].y[iVec] = fB_temp.y[iVec];
fB[ista].z[iVec] = fB_temp.z[iVec];
if (ih == 0) {
- z_start[iVec] = z[ista][iVec];
- x_first[iVec] = x[ista][iVec];
- y_first[iVec] = y[ista][iVec];
- time_first[iVec] = time[ista][iVec];
- wtime_first[iVec] = sta[ista].timeInfo ? 1. : 0.;
- dt2_first[iVec] = dt2[ista][iVec];
- cov_xy_fst.C00[iVec] = cov_xy[ista].C00[iVec];
- cov_xy_fst.C10[iVec] = cov_xy[ista].C10[iVec];
- cov_xy_fst.C11[iVec] = cov_xy[ista].C11[iVec];
+ z_start[iVec] = z[ista][iVec];
+ x_first[iVec] = x[ista][iVec];
+ y_first[iVec] = y[ista][iVec];
+ time_first[iVec] = time[ista][iVec];
+ wtime_first[iVec] = sta[ista].timeInfo ? 1. : 0.;
+ dt2_first[iVec] = dt2[ista][iVec];
+ mxy_first.X()[iVec] = mxy[ista].X()[iVec];
+ mxy_first.Y()[iVec] = mxy[ista].Y()[iVec];
+ mxy_first.Dx2()[iVec] = mxy[ista].Dx2()[iVec];
+ mxy_first.Dy2()[iVec] = mxy[ista].Dy2()[iVec];
+ mxy_first.Dxy()[iVec] = mxy[ista].Dxy()[iVec];
+ mxy_first.NdfX()[iVec] = mxy[ista].NdfX()[iVec];
+ mxy_first.NdfY()[iVec] = mxy[ista].NdfY()[iVec];
}
else if (ih == nHitsTrack - 1) {
- z_end[iVec] = z[ista][iVec];
- x_last[iVec] = x[ista][iVec];
- y_last[iVec] = y[ista][iVec];
- cov_xy_lst.C00[iVec] = cov_xy[ista].C00[iVec];
- cov_xy_lst.C10[iVec] = cov_xy[ista].C10[iVec];
- cov_xy_lst.C11[iVec] = cov_xy[ista].C11[iVec];
- time_last[iVec] = time[ista][iVec];
- dt2_last[iVec] = dt2[ista][iVec];
- wtime_last[iVec] = sta[ista].timeInfo ? 1. : 0.;
+ z_end[iVec] = z[ista][iVec];
+ x_last[iVec] = x[ista][iVec];
+ y_last[iVec] = y[ista][iVec];
+ mxy_last.X()[iVec] = mxy[ista].X()[iVec];
+ mxy_last.Y()[iVec] = mxy[ista].Y()[iVec];
+ mxy_last.Dx2()[iVec] = mxy[ista].Dx2()[iVec];
+ mxy_last.Dy2()[iVec] = mxy[ista].Dy2()[iVec];
+ mxy_last.Dxy()[iVec] = mxy[ista].Dxy()[iVec];
+ mxy_last.NdfX()[iVec] = mxy[ista].NdfX()[iVec];
+ mxy_last.NdfY()[iVec] = mxy[ista].NdfY()[iVec];
+ time_last[iVec] = time[ista][iVec];
+ dt2_last[iVec] = dt2[ista][iVec];
+ wtime_last[iVec] = sta[ista].timeInfo ? 1. : 0.;
}
}
for (int ih = nHitsTrack - 1; ih >= 0; ih--) {
- const int ista = iSta[ih];
+ const int ista = iSta[ih];
const ca::Station& st = fParameters.GetStation(ista);
- By[ista][iVec] = st.fieldSlice.cy[0][0];
+ By[ista][iVec] = st.fieldSlice.cy[0][0];
}
}
@@ -199,10 +209,10 @@ void L1Algo::L1KFTrackFitter()
time_last = iif(w_time[ista], time_last, fvec::Zero());
dt2_last = iif(w_time[ista], dt2_last, fvec(1.e6));
- tr.ResetErrors(cov_xy_lst.C00, cov_xy_lst.C11, 0.1, 0.1, 1.0, dt2_last, 1.e-2);
- tr.C10() = cov_xy_lst.C10;
- tr.X() = x_last;
- tr.Y() = y_last;
+ tr.ResetErrors(mxy_last.Dx2(), mxy_last.Dy2(), 0.1, 0.1, 1.0, dt2_last, 1.e-2);
+ tr.C10() = mxy_last.Dxy();
+ tr.X() = mxy_last.X();
+ tr.Y() = mxy_last.Y();
tr.Time() = time_last;
tr.Vi() = constants::phys::SpeedOfLightInv;
tr.InitVelocityRange(0.5);
@@ -239,7 +249,7 @@ void L1Algo::L1KFTrackFitter()
fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()), fvec(1.f));
fit.SetMask(initialised && w[ista]);
- fit.FilterXY(cov_xy[ista], x[ista], y[ista]);
+ fit.FilterXY(mxy[ista]);
fit.FilterTime(time[ista], dt2[ista], sta[ista].timeInfo);
@@ -255,18 +265,18 @@ void L1Algo::L1KFTrackFitter()
{
fitpv.SetMask(fmask::One());
- L1XYMeasurementInfo vtxInfo;
- vtxInfo.C00 = fvec(1.e-8);
- vtxInfo.C10 = fvec::Zero();
- vtxInfo.C11 = fvec(1.e-8);
+ ca::MeasurementXy<fvec> vtxInfo = TargetMeasurement;
+ vtxInfo.SetDx2(1.e-8);
+ vtxInfo.SetDxy(0.);
+ vtxInfo.SetDy2(1.e-8);
if (kGlobal == fTrackingMode) {
for (int vtxIter = 0; vtxIter < 2; vtxIter++) {
fitpv.SetQp0(fitpv.Tr().Qp());
- fitpv.Tr() = fit.Tr();
+ fitpv.Tr() = fit.Tr();
fitpv.Tr().Qp() = fitpv.Qp0();
fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fld);
- fitpv.FilterXY(vtxInfo, fParameters.GetTargetPositionX(), fParameters.GetTargetPositionY());
+ fitpv.FilterXY(vtxInfo);
}
}
else {
@@ -297,11 +307,11 @@ void L1Algo::L1KFTrackFitter()
ista = 0;
- tr.ResetErrors(cov_xy_fst.C00, cov_xy_fst.C11, 0.1, 0.1, 1., dt2_first, 1.e-2);
- tr.C10() = cov_xy_fst.C10;
+ tr.ResetErrors(mxy_first.Dx2(), mxy_first.Dy2(), 0.1, 0.1, 1., dt2_first, 1.e-2);
+ tr.C10() = mxy_first.Dxy();
- tr.X() = x_first;
- tr.Y() = y_first;
+ tr.X() = mxy_first.X();
+ tr.Y() = mxy_first.Y();
tr.Time() = time_first;
tr.Vi() = constants::phys::SpeedOfLightInv;
tr.InitVelocityRange(0.5);
@@ -335,7 +345,7 @@ void L1Algo::L1KFTrackFitter()
fit.AddMsInMaterial(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()));
fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()), fvec(-1.f));
fit.SetMask(initialised && w[ista]);
- fit.FilterXY(cov_xy[ista], x[ista], y[ista]);
+ fit.FilterXY(mxy[ista]);
fit.FilterTime(time[ista], dt2[ista], sta[ista].timeInfo);
fldB2 = fldB1;
diff --git a/reco/L1/L1Algo/L1TripletConstructor.cxx b/reco/L1/L1Algo/L1TripletConstructor.cxx
index 68fd9c1dbf..a29e269242 100644
--- a/reco/L1/L1Algo/L1TripletConstructor.cxx
+++ b/reco/L1/L1Algo/L1TripletConstructor.cxx
@@ -152,7 +152,7 @@ void L1TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar,
// add the target constraint
- fFit.AddTargetToLine(fAlgo->fTargX, fAlgo->fTargY, fAlgo->fTargZ, fAlgo->TargetXYInfo, fld0);
+ fFit.AddTargetToLine(fAlgo->fTargZ, fAlgo->TargetMeasurement, fld0);
fFit.AddMsInMaterial(fAlgo->fParameters.GetMaterialThickness(fIstaL, T.GetX(), T.GetY()));
@@ -207,21 +207,16 @@ void L1TripletConstructor::FitDoublets()
T2 = fTrL;
fFit.SetQp0(fvec(0.f));
- fvec x_2 = hitm.X();
- fvec y_2 = hitm.Y();
fvec z_2 = hitm.Z();
- fvec t_2 = hitm.T();
- L1XYMeasurementInfo cov_2;
- cov_2.C00 = hitm.dX2();
- cov_2.C10 = hitm.dXY();
- cov_2.C11 = hitm.dY2();
+ ca::MeasurementXy<fvec> m_2(hitm.X(), hitm.Y(), hitm.dX2(), hitm.dY2(), hitm.dXY(), fvec::One(), fvec::One());
+ fvec t_2 = hitm.T();
fvec dt2_2 = hitm.dT2();
// add the middle hit
fFit.ExtrapolateLineNoField(z_2);
- fFit.FilterXY(cov_2, x_2, y_2);
+ fFit.FilterXY(m_2);
fFit.FilterTime(t_2, dt2_2, staM().timeInfo);
@@ -404,18 +399,17 @@ void L1TripletConstructor::FitTriplets()
fscal x[NHits], y[NHits], z[NHits], t[NHits];
fscal dt2[NHits];
- L1XYMeasurementInfo cov[NHits];
+ ca::MeasurementXy<fvec> mxy[NHits];
for (int ih = 0; ih < NHits; ++ih) {
const ca::Hit& hit = fAlgo->fInputData.GetHit(ihit[ih]);
- x[ih] = hit.X();
- y[ih] = hit.Y();
- z[ih] = hit.Z();
- t[ih] = hit.T();
- cov[ih].C00 = hit.dX2();
- cov[ih].C10 = hit.dXY();
- cov[ih].C11 = hit.dY2();
- dt2[ih] = hit.dT2();
+ mxy[ih] = ca::MeasurementXy<fvec>(hit.X(), hit.Y(), hit.dX2(), hit.dY2(), hit.dXY(), fvec::One(), fvec::One());
+
+ x[ih] = hit.X();
+ y[ih] = hit.Y();
+ z[ih] = hit.Z();
+ t[ih] = hit.T();
+ dt2[ih] = hit.dT2();
};
// find the field along the track
@@ -462,22 +456,22 @@ void L1TripletConstructor::FitTriplets()
T.Time() = t[ih0];
T.ResetErrors(1., 1., 1., 1., 100., (sta[ih0].timeInfo ? dt2[ih0] : 1.e6), 1.e2);
- T.C00() = cov[ih0].C00;
- T.C10() = cov[ih0].C10;
- T.C11() = cov[ih0].C11;
+ T.C00() = mxy[ih0].Dx2();
+ T.C10() = mxy[ih0].Dxy();
+ T.C11() = mxy[ih0].Dy2();
T.Ndf() = -ndfTrackModel + 2;
T.NdfTime() = sta[ih0].timeInfo ? -1 : -2;
// add the target constraint
- fit.AddTargetToLine(fAlgo->fTargX, fAlgo->fTargY, fAlgo->fTargZ, fAlgo->TargetXYInfo, fldTarget);
+ fit.AddTargetToLine(fAlgo->fTargZ, fAlgo->TargetMeasurement, fldTarget);
for (int ih = 1; ih < NHits; ++ih) {
fit.Extrapolate(z[ih], fld);
auto radThick = fAlgo->fParameters.GetMaterialThickness(ista[ih], T.X(), T.Y());
fit.AddMsInMaterial(radThick);
fit.EnergyLossCorrection(radThick, fvec(-1.f));
- fit.FilterXY(cov[ih], x[ih], y[ih]);
+ fit.FilterXY(mxy[ih]);
fit.FilterTime(t[ih], dt2[ih], sta[ih].timeInfo);
}
}
@@ -499,9 +493,9 @@ void L1TripletConstructor::FitTriplets()
T.Vi() = 0.;
T.ResetErrors(1., 1., 1., 1., 100., (sta[ih0].timeInfo ? dt2[ih0] : 1.e6), 1.e2);
- T.C00() = cov[ih0].C00;
- T.C10() = cov[ih0].C10;
- T.C11() = cov[ih0].C11;
+ T.C00() = mxy[ih0].Dx2();
+ T.C10() = mxy[ih0].Dxy();
+ T.C11() = mxy[ih0].Dy2();
T.Ndf() = -ndfTrackModel + 2;
T.NdfTime() = sta[ih0].timeInfo ? -1 : -2;
@@ -511,7 +505,7 @@ void L1TripletConstructor::FitTriplets()
auto radThick = fAlgo->fParameters.GetMaterialThickness(ista[ih], T.X(), T.Y());
fit.AddMsInMaterial(radThick);
fit.EnergyLossCorrection(radThick, fvec(1.f));
- fit.FilterXY(cov[ih], x[ih], y[ih]);
+ fit.FilterXY(mxy[ih]);
fit.FilterTime(t[ih], dt2[ih], sta[ih].timeInfo);
}
}
diff --git a/reco/L1/L1Algo/L1UMeasurementInfo.cxx b/reco/L1/L1Algo/L1UMeasurementInfo.cxx
deleted file mode 100644
index 3fc58dd95a..0000000000
--- a/reco/L1/L1Algo/L1UMeasurementInfo.cxx
+++ /dev/null
@@ -1,39 +0,0 @@
-/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov, Sergei Zharko [committer] */
-
-#include "L1UMeasurementInfo.h"
-
-#include <iomanip>
-#include <sstream>
-
-using namespace cbm::algo::ca;
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-std::string L1UMeasurementInfo::ToString(int indentLevel) const
-{
- std::stringstream aStream {};
- // TODO: possibly it is better to place the indentChar into L1Parameters (S.Zharko)
- constexpr char indentChar = '\t';
- std::string indent(indentLevel, indentChar);
- aStream << indent << "cos(phi): " << std::setw(12) << std::setfill(' ') << cos_phi[0] << '\n';
- aStream << indent << "sin(phi): " << std::setw(12) << std::setfill(' ') << sin_phi[0] << '\n';
- return aStream.str();
-}
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-void L1UMeasurementInfo::CheckConsistency() const
-{
- // (i) Checks for the horizontal equality of SIMD vector elements
- L1Utils::CheckSimdVectorEquality(cos_phi, "L1UMeasurementsInfo::cos_phi");
- L1Utils::CheckSimdVectorEquality(sin_phi, "L1UMeasurementsInfo::sin_phi");
- /*
- if (fabs( cos_phi[0]*cos_phi[0]+sin_phi[0]*sin_phi[0]-1.) >1.e-8) {
- std::stringstream msg;
- msg << "L1UMeasurementsInfo: cos_phi and sin_phi do not match: " << cos_phi[0] << " "<<sin_phi[0];
- throw std::logic_error(msg.str());
- }
- */
-}
diff --git a/reco/L1/L1Algo/L1UMeasurementInfo.h b/reco/L1/L1Algo/L1UMeasurementInfo.h
deleted file mode 100644
index 2daa40fafe..0000000000
--- a/reco/L1/L1Algo/L1UMeasurementInfo.h
+++ /dev/null
@@ -1,42 +0,0 @@
-/* Copyright (C) 2007-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer], Igor Kulakov, Sergei Zharko */
-
-#ifndef L1UMeasurementInfo_h
-#define L1UMeasurementInfo_h 1
-
-#include <string>
-
-#include "CaConstants.h"
-#include "CaSimd.h"
-#include "L1Def.h"
-#include "L1Utils.h"
-
-class L1UMeasurementInfo {
-
-public:
- cbm::algo::ca::fvec cos_phi = cbm::algo::ca::constants::Undef<cbm::algo::ca::fvec>;
- cbm::algo::ca::fvec sin_phi = cbm::algo::ca::constants::Undef<cbm::algo::ca::fvec>;
-
- /// String representation of class contents
- /// \param indentLevel number of indent characters in the output
- std::string ToString(int indentLevel = 0) const;
-
- /// Check consistency
- void CheckConsistency() const;
-
- /// Checks, if the fields are NaN
- bool IsNaN() const { return isnan(cos_phi).isNotEmpty() || isnan(sin_phi).isNotEmpty(); }
-
- /// Serialization function
- friend class boost::serialization::access;
- template<class Archive>
- void serialize(Archive& ar, const unsigned int)
- {
- ar& cos_phi;
- ar& sin_phi;
- }
-} _fvecalignment;
-
-
-#endif
diff --git a/reco/L1/L1Algo/L1XYMeasurementInfo.cxx b/reco/L1/L1Algo/L1XYMeasurementInfo.cxx
deleted file mode 100644
index 1eb9a7ef0d..0000000000
--- a/reco/L1/L1Algo/L1XYMeasurementInfo.cxx
+++ /dev/null
@@ -1,30 +0,0 @@
-/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov, Sergei Zharko [committer] */
-
-#include "L1XYMeasurementInfo.h"
-
-#include <iomanip>
-#include <sstream> // for stringstream
-
-#include "L1Utils.h"
-
-void L1XYMeasurementInfo::CheckConsistency() const
-{
- /* (i) Checks for the horizontal equality of SIMD vector elements */
- L1Utils::CheckSimdVectorEquality(C00, "L1XYMeasurementsInfo::C00");
- L1Utils::CheckSimdVectorEquality(C10, "L1XYMeasurementsInfo::C10");
- L1Utils::CheckSimdVectorEquality(C11, "L1XYMeasurementsInfo::C11");
-}
-
-std::string L1XYMeasurementInfo::ToString(int indentLevel) const
-{
- std::stringstream aStream {};
- // TODO: possibly is better to be place indentChar into L1Parameters (S.Zharko)
- constexpr char indentChar = '\t';
- std::string indent(indentLevel, indentChar);
- aStream << indent << "C00: " << std::setw(12) << std::setfill(' ') << C00[0] << '\n';
- aStream << indent << "C10: " << std::setw(12) << std::setfill(' ') << C10[0] << '\n';
- aStream << indent << "C11: " << std::setw(12) << std::setfill(' ') << C11[0];
- return aStream.str();
-}
diff --git a/reco/L1/L1Algo/L1XYMeasurementInfo.h b/reco/L1/L1Algo/L1XYMeasurementInfo.h
deleted file mode 100644
index a3a3fe7f77..0000000000
--- a/reco/L1/L1Algo/L1XYMeasurementInfo.h
+++ /dev/null
@@ -1,51 +0,0 @@
-/* Copyright (C) 2007-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer], Sergei Zharko */
-
-#ifndef L1XYMeasurementInfo_h
-#define L1XYMeasurementInfo_h 1
-
-#include <string>
-
-#include "CaConstants.h"
-#include "CaSimd.h"
-#include "CaUtils.h"
-#include "L1Def.h"
-
-using namespace cbm::algo::ca;
-
-class L1XYMeasurementInfo {
-public:
- /// Consistency checker
- void CheckConsistency() const;
-
- /// String representation of class contents
- /// \param indentLevel number of indent characters in the output
- std::string ToString(int indentLevel = 0) const;
-
- /// Checks, if the fields are NaN
- bool IsNaN() const
- {
- return cbm::algo::ca::utils::IsUndefined(C00) || cbm::algo::ca::utils::IsUndefined(C10)
- || cbm::algo::ca::utils::IsUndefined(C11);
- }
-
- /// Serialization function
- friend class boost::serialization::access;
- template<class Archive>
- void serialize(Archive& ar, const unsigned int)
- {
- ar& C00;
- ar& C10;
- ar& C11;
- }
-
-public:
- cbm::algo::ca::fvec C00 {cbm::algo::ca::constants::Undef<cbm::algo::ca::fvec>};
- cbm::algo::ca::fvec C10 {cbm::algo::ca::constants::Undef<cbm::algo::ca::fvec>};
- cbm::algo::ca::fvec C11 {cbm::algo::ca::constants::Undef<cbm::algo::ca::fvec>};
-
-} _fvecalignment;
-
-
-#endif
diff --git a/reco/L1/L1Algo/L1Event.cxx b/reco/L1/L1Algo/inactive/L1Event.cxx
similarity index 100%
rename from reco/L1/L1Algo/L1Event.cxx
rename to reco/L1/L1Algo/inactive/L1Event.cxx
diff --git a/reco/L1/L1Algo/L1Event.h b/reco/L1/L1Algo/inactive/L1Event.h
similarity index 100%
rename from reco/L1/L1Algo/L1Event.h
rename to reco/L1/L1Algo/inactive/L1Event.h
diff --git a/reco/L1/L1Algo/L1EventEfficiencies.h b/reco/L1/L1Algo/inactive/L1EventEfficiencies.h
similarity index 100%
rename from reco/L1/L1Algo/L1EventEfficiencies.h
rename to reco/L1/L1Algo/inactive/L1EventEfficiencies.h
diff --git a/reco/L1/L1Algo/L1EventMatch.cxx b/reco/L1/L1Algo/inactive/L1EventMatch.cxx
similarity index 100%
rename from reco/L1/L1Algo/L1EventMatch.cxx
rename to reco/L1/L1Algo/inactive/L1EventMatch.cxx
diff --git a/reco/L1/L1Algo/L1EventMatch.h b/reco/L1/L1Algo/inactive/L1EventMatch.h
similarity index 100%
rename from reco/L1/L1Algo/L1EventMatch.h
rename to reco/L1/L1Algo/inactive/L1EventMatch.h
diff --git a/reco/L1/L1Algo/L1MCEvent.cxx b/reco/L1/L1Algo/inactive/L1MCEvent.cxx
similarity index 100%
rename from reco/L1/L1Algo/L1MCEvent.cxx
rename to reco/L1/L1Algo/inactive/L1MCEvent.cxx
diff --git a/reco/L1/L1Algo/L1MCEvent.h b/reco/L1/L1Algo/inactive/L1MCEvent.h
similarity index 100%
rename from reco/L1/L1Algo/L1MCEvent.h
rename to reco/L1/L1Algo/inactive/L1MCEvent.h
--
GitLab