diff --git a/algo/ca/core/CMakeLists.txt b/algo/ca/core/CMakeLists.txt
index 83170ef902fca88d0e4f205cbab2aca9de3a9868..62e882c040da3ab6d2397b70048d85769ba20386 100644
--- a/algo/ca/core/CMakeLists.txt
+++ b/algo/ca/core/CMakeLists.txt
@@ -8,6 +8,7 @@ set(INCLUDE_DIRECTORIES
set(SRCS
${CMAKE_CURRENT_SOURCE_DIR}/data/CaTrackParam.cxx
${CMAKE_CURRENT_SOURCE_DIR}/data/CaTrack.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/data/CaField.cxx
)
add_library(CaCore SHARED ${SRCS})
@@ -23,13 +24,9 @@ target_compile_definitions(CaCore PUBLIC NO_ROOT)
target_link_libraries(CaCore
Vc::Vc
Boost::serialization
- Boost::program_options
- Boost::filesystem
- Boost::headers
- OnlineDataLog
- external::fles_logging
- external::fles_ipc
-# external::fles_monitoring # in test
+ OnlineDataLog # needed for the logger
+ external::fles_logging # needed for the logger
+ external::fles_ipc # needed for the logger
)
install(TARGETS CaCore DESTINATION lib)
@@ -42,6 +39,7 @@ install(
data/CaTrackParam.h
data/CaTrack.h
data/CaHit.h
+ data/CaField.h
pars/CaConstants.h
utils/CaSimd.h
utils/CaSimdVc.h
diff --git a/reco/L1/L1Algo/L1Field.cxx b/algo/ca/core/data/CaField.cxx
similarity index 70%
rename from reco/L1/L1Algo/L1Field.cxx
rename to algo/ca/core/data/CaField.cxx
index f661795944f8656bc210e40a8564e8bacf991451..bb1f2051b1755d4749452cd68392ab7be8996b07 100644
--- a/reco/L1/L1Algo/L1Field.cxx
+++ b/algo/ca/core/data/CaField.cxx
@@ -2,41 +2,43 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Sergey Gorbunov, Sergei Zharko [committer] */
-#include "L1Field.h"
-
-#include "FairField.h"
-#include "FairRunAna.h"
+#include "CaField.h"
#include <iomanip>
#include <iostream>
#include <sstream>
#include "CaTrackParam.h"
-#include "L1Utils.h"
+#include "CaUtils.h"
using namespace cbm::algo::ca;
//
-// L1FieldValue methods
+// FieldValue methods
//
-bool L1FieldRegion::fgkForceUseOfOriginalField = false;
+bool FieldRegion::fgForceUseOfOriginalField = false;
+std::function<std::tuple<double, double, double>(double x, double y, double z)>
+ FieldRegion::fgOdiginalField([](double, double, double) {
+ return std::tuple(constants::Undef<double>, constants::Undef<double>, constants::Undef<double>);
+ });
+
//----------------------------------------------------------------------------------------------------------------------
//
-void L1FieldValue::CheckConsistency() const
+void FieldValue::CheckConsistency() const
{
/* Check SIMD data vectors for consistent initialization */
- L1Utils::CheckSimdVectorEquality(x, "L1FieldValue::x");
- L1Utils::CheckSimdVectorEquality(y, "L1FieldValue::y");
- L1Utils::CheckSimdVectorEquality(z, "L1FieldValue::z");
+ utils::CheckSimdVectorEquality(x, "FieldValue::x");
+ utils::CheckSimdVectorEquality(y, "FieldValue::y");
+ utils::CheckSimdVectorEquality(z, "FieldValue::z");
// TODO: Any other checks? (S.Zharko)
}
//----------------------------------------------------------------------------------------------------------------------
// TODO:
-std::string L1FieldValue::ToString(int indentLevel) const
+std::string FieldValue::ToString(int indentLevel) const
{
std::stringstream aStream {};
constexpr char indentChar = '\t';
@@ -49,18 +51,15 @@ std::string L1FieldValue::ToString(int indentLevel) const
//----------------------------------------------------------------------------------------------------------------------
//
-std::ostream& operator<<(std::ostream& out, const L1FieldValue& B)
-{
- return out << B.x << " | " << B.y << " | " << B.z;
-};
+std::ostream& operator<<(std::ostream& out, const FieldValue& B) { return out << B.x << " | " << B.y << " | " << B.z; };
//
-// L1FieldSlice methods
+// FieldSlice methods
//
//----------------------------------------------------------------------------------------------------------------------
//
-L1FieldSlice::L1FieldSlice()
+FieldSlice::FieldSlice()
{
for (int i = 0; i < constants::size::MaxNFieldApproxCoefficients; ++i) {
cx[i] = constants::Undef<fscal>;
@@ -72,20 +71,20 @@ L1FieldSlice::L1FieldSlice()
//----------------------------------------------------------------------------------------------------------------------
//
-void L1FieldSlice::CheckConsistency() const
+void FieldSlice::CheckConsistency() const
{
/* Check SIMD data vectors for consistent initialization */
for (int i = 0; i < constants::size::MaxNFieldApproxCoefficients; ++i) {
- L1Utils::CheckSimdVectorEquality(cx[i], "L1FieldSlice: cx");
- L1Utils::CheckSimdVectorEquality(cy[i], "L1FieldSlice: cy");
- L1Utils::CheckSimdVectorEquality(cz[i], "L1FieldSlice: cz");
+ utils::CheckSimdVectorEquality(cx[i], "FieldSlice: cx");
+ utils::CheckSimdVectorEquality(cy[i], "FieldSlice: cy");
+ utils::CheckSimdVectorEquality(cz[i], "FieldSlice: cz");
}
- L1Utils::CheckSimdVectorEquality(z, "L1FieldSlice: z");
+ utils::CheckSimdVectorEquality(z, "FieldSlice: z");
}
//----------------------------------------------------------------------------------------------------------------------
// TODO: Should it be inline? (S.Zharko)
-void L1FieldSlice::GetFieldValue(const fvec& x, const fvec& y, L1FieldValue& B) const
+FieldValue FieldSlice::GetFieldValue(const fvec& x, const fvec& y) const
{
fvec x2 = x * x;
fvec y2 = y * y;
@@ -109,6 +108,8 @@ void L1FieldSlice::GetFieldValue(const fvec& x, const fvec& y, L1FieldValue& B)
fvec x3y2 = x3 * y2;
fvec x4y = x4 * y;
+ FieldValue B;
+
B.x = cx[0] + cx[1] * x + cx[2] * y + cx[3] * x2 + cx[4] * xy + cx[5] * y2 + cx[6] * x3 + cx[7] * x2y + cx[8] * xy2
+ cx[9] * y3 + cx[10] * x4 + cx[11] * x3y + cx[12] * x2y2 + cx[13] * xy3 + cx[14] * y4 + cx[15] * x5
+ cx[16] * x4y + cx[17] * x3y2 + cx[18] * x2y3 + cx[19] * xy4 + cx[20] * y5;
@@ -120,17 +121,18 @@ void L1FieldSlice::GetFieldValue(const fvec& x, const fvec& y, L1FieldValue& B)
B.z = cz[0] + cz[1] * x + cz[2] * y + cz[3] * x2 + cz[4] * xy + cz[5] * y2 + cz[6] * x3 + cz[7] * x2y + cz[8] * xy2
+ cz[9] * y3 + cz[10] * x4 + cz[11] * x3y + cz[12] * x2y2 + cz[13] * xy3 + cz[14] * y4 + cz[15] * x5
+ cz[16] * x4y + cz[17] * x3y2 + cz[18] * x2y3 + cz[19] * xy4 + cz[20] * y5;
+ return B;
}
-void L1FieldSlice::GetFieldValueForLine(const cbm::algo::ca::TrackParamV& t, L1FieldValue& B) const
+FieldValue FieldSlice::GetFieldValueForLine(const TrackParamV& t) const
{
fvec dz = z - t.GetZ();
- GetFieldValue(t.GetX() + t.GetTx() * dz, t.GetY() + t.GetTy() * dz, B);
+ return GetFieldValue(t.GetX() + t.GetTx() * dz, t.GetY() + t.GetTy() * dz);
}
//----------------------------------------------------------------------------------------------------------------------
//
-std::string L1FieldSlice::ToString(int indentLevel) const
+std::string FieldSlice::ToString(int indentLevel) const
{
std::stringstream aStream {};
// TODO: possibly it is better to place the indentChar into L1Parameters (S.Zharko)
@@ -148,12 +150,12 @@ std::string L1FieldSlice::ToString(int indentLevel) const
}
//
-// L1FieldRegion methdos
+// FieldRegion methdos
//
//----------------------------------------------------------------------------------------------------------------------
//
-L1FieldRegion::L1FieldRegion(float reg[10]) noexcept
+FieldRegion::FieldRegion(float reg[10]) noexcept
: cx0(reg[0])
, cx1(reg[1])
, cx2(reg[2])
@@ -169,53 +171,51 @@ L1FieldRegion::L1FieldRegion(float reg[10]) noexcept
//----------------------------------------------------------------------------------------------------------------------
//
-void L1FieldRegion::CheckConsistency() const
+void FieldRegion::CheckConsistency() const
{
/* Check SIMD data vectors for consistent initialization */
- L1Utils::CheckSimdVectorEquality(cx0, "L1FieldRegion::cx0");
- L1Utils::CheckSimdVectorEquality(cx1, "L1FieldRegion::cx1");
- L1Utils::CheckSimdVectorEquality(cx2, "L1FieldRegion::cx2");
- L1Utils::CheckSimdVectorEquality(cy0, "L1FieldRegion::cy0");
- L1Utils::CheckSimdVectorEquality(cy1, "L1FieldRegion::cy1");
- L1Utils::CheckSimdVectorEquality(cy2, "L1FieldRegion::cy2");
- L1Utils::CheckSimdVectorEquality(cz0, "L1FieldRegion::cz0");
- L1Utils::CheckSimdVectorEquality(cz1, "L1FieldRegion::cz1");
- L1Utils::CheckSimdVectorEquality(cz2, "L1FieldRegion::cz2");
- L1Utils::CheckSimdVectorEquality(z0, "L1FieldRegion::z0");
+ utils::CheckSimdVectorEquality(cx0, "FieldRegion::cx0");
+ utils::CheckSimdVectorEquality(cx1, "FieldRegion::cx1");
+ utils::CheckSimdVectorEquality(cx2, "FieldRegion::cx2");
+ utils::CheckSimdVectorEquality(cy0, "FieldRegion::cy0");
+ utils::CheckSimdVectorEquality(cy1, "FieldRegion::cy1");
+ utils::CheckSimdVectorEquality(cy2, "FieldRegion::cy2");
+ utils::CheckSimdVectorEquality(cz0, "FieldRegion::cz0");
+ utils::CheckSimdVectorEquality(cz1, "FieldRegion::cz1");
+ utils::CheckSimdVectorEquality(cz2, "FieldRegion::cz2");
+ utils::CheckSimdVectorEquality(z0, "FieldRegion::z0");
// TODO: Any other checks? (S.Zharko)
}
//----------------------------------------------------------------------------------------------------------------------
// TODO: Should it be inline? (S.Zharko)
-void L1FieldRegion::Get(const fvec& x, const fvec& y, const fvec& z, fvec* B) const
+FieldValue FieldRegion::Get(const fvec& x, const fvec& y, const fvec& z) const
{
- if (fgkForceUseOfOriginalField || fUseOriginalField) {
- // TODO: replace with functional object
- assert(FairRunAna::Instance());
+ FieldValue B;
+ if (fgForceUseOfOriginalField || fUseOriginalField) {
for (size_t i = 0; i < fvec::size(); i++) {
- double inPos[3] = {x[i], y[i], z[i]};
- double outB[3] = {0., 0., 0.};
- if (FairRunAna::Instance()->GetField()) { FairRunAna::Instance()->GetField()->GetFieldValue(inPos, outB); }
- B[0][i] = outB[0];
- B[1][i] = outB[1];
- B[2][i] = outB[2];
+ auto [bx, by, bz] = fgOdiginalField(x[i], y[i], z[i]);
+ B.x[i] = bx;
+ B.y[i] = by;
+ B.z[i] = bz;
}
}
else {
fvec dz = z - z0;
fvec dz2 = dz * dz;
- B[0] = cx0 + cx1 * dz + cx2 * dz2;
- B[1] = cy0 + cy1 * dz + cy2 * dz2;
- B[2] = cz0 + cz1 * dz + cz2 * dz2;
+ B.x = cx0 + cx1 * dz + cx2 * dz2;
+ B.y = cy0 + cy1 * dz + cy2 * dz2;
+ B.z = cz0 + cz1 * dz + cz2 * dz2;
}
+ return B;
}
//----------------------------------------------------------------------------------------------------------------------
// TODO: Should it be inline? (S.Zharko)
-void L1FieldRegion::Set(const L1FieldValue& b0, const fvec b0z, const L1FieldValue& b1, const fvec b1z,
- const L1FieldValue& b2, const fvec b2z)
+void FieldRegion::Set(const FieldValue& b0, const fvec b0z, const FieldValue& b1, const fvec b1z, const FieldValue& b2,
+ const fvec b2z)
{
z0 = b0z;
fvec dz1 = b1z - b0z;
@@ -248,7 +248,7 @@ void L1FieldRegion::Set(const L1FieldValue& b0, const fvec b0z, const L1FieldVal
//----------------------------------------------------------------------------------------------------------------------
// TODO: Should it be inline? (S.Zharko)
-void L1FieldRegion::Set(const L1FieldValue& b0, const fvec b0z, const L1FieldValue& b1, const fvec b1z)
+void FieldRegion::Set(const FieldValue& b0, const fvec b0z, const FieldValue& b1, const fvec b1z)
{
z0 = b0z;
fvec dzi = fvec::One() / (b1z - b0z);
@@ -265,7 +265,7 @@ void L1FieldRegion::Set(const L1FieldValue& b0, const fvec b0z, const L1FieldVal
//----------------------------------------------------------------------------------------------------------------------
// TODO: Should it be inline? (S.Zharko)
-void L1FieldRegion::Shift(fvec z)
+void FieldRegion::Shift(fvec z)
{
fvec dz = z - z0;
fvec cx2dz = cx2 * dz;
@@ -282,7 +282,7 @@ void L1FieldRegion::Shift(fvec z)
//----------------------------------------------------------------------------------------------------------------------
// TODO: Should it be inline? (S.Zharko)
-void L1FieldRegion::SetOneEntry(const int i0, const L1FieldRegion& f1, const int i1)
+void FieldRegion::SetOneEntry(const int i0, const FieldRegion& f1, const int i1)
{
cx0[i0] = f1.cx0[i1];
cx1[i0] = f1.cx1[i1];
@@ -298,7 +298,7 @@ void L1FieldRegion::SetOneEntry(const int i0, const L1FieldRegion& f1, const int
//----------------------------------------------------------------------------------------------------------------------
// TODO: Should it be inline? (S.Zharko)
-void L1FieldRegion::SetOneEntry(const L1FieldRegion& f1, const int i1)
+void FieldRegion::SetOneEntry(const FieldRegion& f1, const int i1)
{
cx0 = f1.cx0[i1];
cx1 = f1.cx1[i1];
@@ -312,25 +312,9 @@ void L1FieldRegion::SetOneEntry(const L1FieldRegion& f1, const int i1)
z0 = f1.z0[i1];
}
-//----------------------------------------------------------------------------------------------------------------------
-// TODO: Should it be inline? (S.Zharko)
-void L1FieldRegion::GetOneEntry(float reg[10], const int iVec)
-{
- reg[0] = cx0[iVec];
- reg[1] = cx1[iVec];
- reg[2] = cx2[iVec];
- reg[3] = cy0[iVec];
- reg[4] = cy1[iVec];
- reg[5] = cy2[iVec];
- reg[6] = cz0[iVec];
- reg[7] = cz1[iVec];
- reg[8] = cz2[iVec];
- reg[9] = z0[iVec];
-}
-
//----------------------------------------------------------------------------------------------------------------------
//
-std::string L1FieldRegion::ToString(int indentLevel) const
+std::string FieldRegion::ToString(int indentLevel) const
{
std::stringstream aStream {};
// TODO: possibly it is better to place the indentChar into L1Parameters (S.Zharko)
diff --git a/algo/ca/core/data/CaField.h b/algo/ca/core/data/CaField.h
new file mode 100644
index 0000000000000000000000000000000000000000..0e56f96e03e1e102b1ea7f07b2de88a98613d76a
--- /dev/null
+++ b/algo/ca/core/data/CaField.h
@@ -0,0 +1,232 @@
+/* Copyright (C) 2007-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov [committer], Igor Kulakov, Maksym Zyzak, Sergei Zharko */
+
+#ifndef CaField_h
+#define CaField_h 1
+
+#include <boost/serialization/access.hpp>
+
+#include <functional>
+#include <string>
+
+#include "CaConstants.h"
+#include "CaSimd.h"
+#include "CaTrackParam.h"
+
+namespace cbm::algo::ca
+{
+
+ class FieldValue {
+
+ public:
+ fvec x {0.f}; //< x-component of the field
+ fvec y {0.f}; //< y-component of the field
+ fvec z {0.f}; //< z-component of the field
+
+ /// Combines the magnetic field with another field value object using weight
+ /// \param B other field value to combine with
+ /// \param w weight from 0 to 1
+ void Combine(const FieldValue& B, const fvec& w);
+
+ void Combine(const FieldValue& B, const fmask& w);
+
+ /// Consistency checker
+ void CheckConsistency() const;
+
+ /// Operator << overloading
+ friend std::ostream& operator<<(std::ostream& out, const FieldValue& B);
+
+ /// String representation of class contents
+ /// \param indentLevel number of indent characters in the output
+ std::string ToString(int indentLevel) const;
+
+ /// Serialization function
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int)
+ {
+ ar& x;
+ ar& y;
+ ar& z;
+ }
+ } _fvecalignment;
+
+ [[gnu::always_inline]] inline void FieldValue::Combine(const FieldValue& B, const fvec& w)
+ {
+ x += w * (B.x - x);
+ y += w * (B.y - y);
+ z += w * (B.z - z);
+ }
+
+ [[gnu::always_inline]] inline void FieldValue::Combine(const FieldValue& B, const fmask& w)
+ {
+ x(w) = B.x;
+ y(w) = B.y;
+ z(w) = B.z;
+ }
+
+ /// Class represents a set of magnetic field approximation coefficients
+ ///
+ // TODO: Crosscheck the default content (S.Zharko)
+ class FieldSlice {
+ public:
+ /// Default constructor
+ FieldSlice();
+
+ /// Consistency checker
+ void CheckConsistency() const;
+
+ /// Gets field value from (x, y) fvec point
+ /// \param x x-coordinate of input
+ /// \param y y-coordinate of input
+ /// \return B the FieldValue output
+ FieldValue GetFieldValue(const fvec& x, const fvec& y) const;
+
+ /// Gets field value for the intersection with a straight track
+ /// \param t straight track
+ /// \return B the FieldValue output
+ FieldValue GetFieldValueForLine(const TrackParamV& t) const;
+
+ /// String representation of class contents
+ /// \param indentLevel number of indent characters in the output
+ std::string ToString(int indentLevel = 0) const;
+
+ public:
+ // NOTE: We don't use an initialization of arrays here because we cannot be sure
+ // if the underlying type (fvec) has a default constructor, but
+ // we are sure, that it can be initialized with a float. (S.Zharko)
+ static constexpr auto kMaxNFieldApproxCoefficients = constants::size::MaxNFieldApproxCoefficients;
+ fvec cx[kMaxNFieldApproxCoefficients]; ///< Polynomial coefficients for x-component of the field value
+ fvec cy[kMaxNFieldApproxCoefficients]; ///< Polynomial coefficients for y-component of the field value
+ fvec cz[kMaxNFieldApproxCoefficients]; ///< Polynomial coefficients for z-component of the field value
+
+ fvec z {constants::Undef<fscal>}; ///< z coordinate of the slice
+
+ /// Serialization function
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int)
+ {
+ ar& cx;
+ ar& cy;
+ ar& cz;
+ ar& z;
+ }
+ } _fvecalignment;
+
+
+ class FieldRegion {
+ public:
+ // NOTE: When a custom constructor is defined, default constructor also should be provided (S.Zharko)
+ FieldRegion() = default;
+
+ FieldRegion(float reg[10]) noexcept;
+
+ /// Consistency checker
+ void CheckConsistency() const;
+
+ /// Gets the field vector and writes it into B pointer
+ /// \param x x-coordinate of the point to calculate the field
+ /// \param y y-coordinate of the point to calculate the field
+ /// \param z z-coordinate of the point to calculate the field
+ /// \return the magnetic field value
+ FieldValue Get(const fvec& x, const fvec& y, const fvec& z) const;
+
+ /// Interpolates the magnetic field by three nodal points and sets the result to this FieldRegion object
+ /// The result is a quadratic interpolation of the field as a function of z
+ /// \param b0 field value in the first nodal point
+ /// \param b0z z-coordinate of the first nodal point
+ /// \param b1 field value in the second nodal point
+ /// \param b1z z-coordinate of the second nodal point
+ /// \param b2 field value in the third nodal point
+ /// \param b2z z-coordinate of the third nodal point
+ /// TODO: does the sequence of b0z, b1z and b2z matter? If yes, probalby we need an assertion (S.Zharko)
+ void Set(const FieldValue& b0, const fvec b0z, const FieldValue& b1, const fvec b1z, const FieldValue& b2,
+ const fvec b2z);
+
+ /// Interpolates the magnetic field by thwo nodal points and sets the result to this FieldRegion object.
+ /// The result is a linear interpolation of the field as a function of z
+ /// \param b0 field value in the first nodal point
+ /// \param b0z z-coordinate of the first nodal point
+ /// \param b1 field value in the second nodal point
+ /// \param b1z z-coordinate of the second nodal point
+ /// TODO: does the sequence of b0z and b1z matter? If yes, probalby we need an assertion (S.Zharko)
+ void Set(const FieldValue& b0, const fvec b0z, const FieldValue& b1, const fvec b1z);
+
+ /// Shifts the coefficients to new central point
+ /// \param z z-coordinate of the new central point
+ void Shift(fvec z);
+
+ /// Replaces the selected layer of the coefficients with one from another
+ /// FieldRegion object
+ /// \param i0 the index of the fvect layer in this FieldRegion object
+ /// \param fl the other FieldRegion object, which layer is going to be replaced
+ /// \param i1 the index of the source fvect layer to copy
+ void SetOneEntry(const int i0, const FieldRegion& f1, const int i1);
+
+ /// Replaces all the layers of the coefficients with one selected layer from another
+ /// FieldRegion object
+ /// \param fl the other FieldRegion object, which layer is going to be replaced
+ /// \param i1 the index of the source fvect layer to copy
+ void SetOneEntry(const FieldRegion& f1, const int i1);
+
+ /// String representation of class contents
+ /// \param indentLevel number of indent characters in the output
+ std::string ToString(int indentLevel = 0) const;
+
+ /// Use original field instead of the field approximation
+ void SetUseOriginalField(bool v = true) { fUseOriginalField = v; }
+
+ /// Force using the original field
+ static void ForceUseOfOriginalField(bool v = true) { fgForceUseOfOriginalField = v; }
+
+ /// Set the original field
+ static void SetOriginalField(std::function<std::tuple<double, double, double>(double x, double y, double z)> f)
+ {
+ fgOdiginalField = f;
+ }
+
+ public:
+ static bool fgForceUseOfOriginalField;
+ static std::function<std::tuple<double, double, double>(double x, double y, double z)> fgOdiginalField;
+
+ // TODO: Probably it's better to have arrays instead of separate fvec values? (S.Zharko)
+ // Bx(z) = cx0 + cx1*(z-z0) + cx2*(z-z0)^2
+ fvec cx0 {0.f};
+ fvec cx1 {0.f};
+ fvec cx2 {0.f};
+
+ // By(z) = cy0 + cy1*(z-z0) + cy2*(z-z0)^2
+ fvec cy0 {0.f};
+ fvec cy1 {0.f};
+ fvec cy2 {0.f};
+
+ // Bz(z) = cz0 + cz1*(z-z0) + cz2*(z-z0)^2
+ fvec cz0 {0.f};
+ fvec cz1 {0.f};
+ fvec cz2 {0.f};
+
+ fvec z0 {0.f}; ///< z-coordinate of the field region central point
+
+ bool fUseOriginalField {false};
+
+ /// Serialization function
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int)
+ {
+ ar& cx0;
+ ar& cx1;
+ ar& cx2;
+ ar& cy0;
+ ar& cy1;
+ ar& cy2;
+ ar& cz0;
+ ar& cz1;
+ ar& cz2;
+ }
+ } _fvecalignment;
+
+} // namespace cbm::algo::ca
+#endif
diff --git a/algo/ca/core/utils/CaUtils.h b/algo/ca/core/utils/CaUtils.h
index a75a067e54410d6305b1fd8adf52ca62b7308367..a448de531b868ff2a99612d0ed83f94dafee9342 100644
--- a/algo/ca/core/utils/CaUtils.h
+++ b/algo/ca/core/utils/CaUtils.h
@@ -10,6 +10,8 @@
#ifndef CaUtils_h
#define CaUtils_h 1
+#include <sstream>
+
#include "CaConstants.h"
#include "CaSimd.h"
@@ -71,6 +73,19 @@ namespace cbm::algo::ca::utils
}
}
+ [[gnu::always_inline]] inline void CheckSimdVectorEquality(fvec v, const char* name)
+ {
+ bool ok = true;
+ for (size_t i = 1; i < fvec::size(); i++) {
+ ok = ok && (v[i] == v[0]);
+ }
+ if (!ok) {
+ std::stringstream msg;
+ msg << name << " SIMD vector is inconsistent, not all of the words are equal each other: " << v;
+ throw std::logic_error(msg.str());
+ }
+ }
+
} // namespace cbm::algo::ca::utils
diff --git a/reco/KF/KFQA/CbmKFTrackQa.cxx b/reco/KF/KFQA/CbmKFTrackQa.cxx
index 29805b8633d036f54b419e47f0e9814654c56131..61d655ed4763b6293bf1bc14de5ec1ad6051fee8 100644
--- a/reco/KF/KFQA/CbmKFTrackQa.cxx
+++ b/reco/KF/KFQA/CbmKFTrackQa.cxx
@@ -18,11 +18,11 @@
#include "CbmKFVertex.h"
#include "CbmL1PFFitter.h"
+#include "CaField.h"
#include "KFMCTrack.h"
#include "KFParticleMatch.h"
#include "KFParticleTopoReconstructor.h"
#include "KFTopoPerformance.h"
-#include "L1Field.h"
//ROOT headers
#include "TClonesArray.h"
diff --git a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
index 8d2dadc30520684a0c67bdd71a656976ab501125..fcc622372b58e8ea5ebeee24da2e5fe80f5e894f 100644
--- a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
+++ b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
@@ -34,11 +34,11 @@
#include "TDatabasePDG.h"
+#include "CaField.h"
#include "CaTrackParam.h"
#include "KFParticleDatabase.h"
#include "L1Algo.h" // Also defines L1Parameters
#include "L1Def.h"
-#include "L1Field.h"
#include "L1Fit.h"
#include "L1Station.h"
@@ -57,7 +57,7 @@ namespace NS_L1TrackFitter
using namespace NS_L1TrackFitter;
-inline void CbmL1PFFitter::PFFieldRegion::setFromL1FieldRegion(const L1FieldRegion& fld, int i)
+inline void CbmL1PFFitter::PFFieldRegion::setFromL1FieldRegion(const ca::FieldRegion& fld, int i)
{
fP[0] = fld.cx0[i];
fP[1] = fld.cx1[i];
@@ -74,7 +74,7 @@ inline void CbmL1PFFitter::PFFieldRegion::setFromL1FieldRegion(const L1FieldRegi
fP[9] = fld.z0[i];
}
-inline void CbmL1PFFitter::PFFieldRegion::getL1FieldRegion(L1FieldRegion& fld, int i)
+inline void CbmL1PFFitter::PFFieldRegion::getL1FieldRegion(ca::FieldRegion& fld, int i)
{
fld.cx0[i] = fP[0];
fld.cx1[i] = fP[1];
@@ -91,7 +91,7 @@ inline void CbmL1PFFitter::PFFieldRegion::getL1FieldRegion(L1FieldRegion& fld, i
fld.z0[i] = fP[9];
}
-inline CbmL1PFFitter::PFFieldRegion::PFFieldRegion(const L1FieldRegion& fld, int i) { setFromL1FieldRegion(fld, i); }
+inline CbmL1PFFitter::PFFieldRegion::PFFieldRegion(const ca::FieldRegion& fld, int i) { setFromL1FieldRegion(fld, i); }
CbmL1PFFitter::CbmL1PFFitter() {}
@@ -155,8 +155,8 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
Initialize();
- L1FieldValue b0, b1, b2 _fvecalignment;
- L1FieldRegion fld _fvecalignment;
+ ca::FieldValue b0, b1, b2 _fvecalignment;
+ ca::FieldRegion fld _fvecalignment;
// fld.SetUseOriginalField();
static int nHits = CbmL1::Instance()->fpAlgo->GetParameters()->GetNstationsActive();
@@ -186,8 +186,8 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
fvec dt2_first, dt2_last;
fvec z0, z1, z2, dz, z_start, z_end;
- L1FieldValue fB[nHits];
- L1FieldValue fB_temp _fvecalignment;
+ ca::FieldValue fB[nHits];
+ ca::FieldValue fB_temp _fvecalignment;
unsigned short N_vTracks = Tracks.size();
@@ -271,7 +271,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
w[ista][iVec] = true;
- sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista], fB_temp);
+ fB_temp = sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista]);
fB[ista].x[iVec] = fB_temp.x[iVec];
fB[ista].y[iVec] = fB_temp.y[iVec];
fB[ista].z[iVec] = fB_temp.z[iVec];
@@ -306,17 +306,17 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
fit.SetQp0(fit.Tr().Qp());
z1 = z[i];
- sta[i].fieldSlice.GetFieldValue(T.X(), T.Y(), b1);
+ b1 = sta[i].fieldSlice.GetFieldValue(T.X(), T.Y());
b1.Combine(fB[i], w[i]);
z2 = z[i + 2];
dz = z2 - z1;
- sta[i].fieldSlice.GetFieldValue(T.X() + T.Tx() * dz, T.Y() + T.Ty() * dz, b2);
+ b2 = sta[i].fieldSlice.GetFieldValue(T.X() + T.Tx() * dz, T.Y() + T.Ty() * dz);
b2.Combine(fB[i + 2], w[i + 2]);
fld.Set(b2, z2, b1, z1, b0, z0);
for (++i; i < nHits; i++) {
z0 = z[i];
dz = (z1 - z0);
- sta[i].fieldSlice.GetFieldValue(T.X() - T.Tx() * dz, T.Y() - T.Ty() * dz, b0);
+ b0 = sta[i].fieldSlice.GetFieldValue(T.X() - T.Tx() * dz, T.Y() - T.Ty() * dz);
b0.Combine(fB[i], w[i]);
fld.Set(b0, z0, b1, z1, b2, z2);
@@ -366,18 +366,18 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
FilterFirst(fit, x_last, y_last, t_last, cov_last, dt2_last);
z1 = z[i];
- sta[i].fieldSlice.GetFieldValue(T.X(), T.Y(), b1);
+ b1 = sta[i].fieldSlice.GetFieldValue(T.X(), T.Y());
b1.Combine(fB[i], w[i]);
z2 = z[i - 2];
dz = z2 - z1;
- sta[i].fieldSlice.GetFieldValue(T.X() + T.Tx() * dz, T.Y() + T.Ty() * dz, b2);
+ b2 = sta[i].fieldSlice.GetFieldValue(T.X() + T.Tx() * dz, T.Y() + T.Ty() * dz);
b2.Combine(fB[i - 2], w[i - 2]);
fld.Set(b2, z2, b1, z1, b0, z0);
for (--i; i >= 0; i--) {
z0 = z[i];
dz = (z1 - z0);
- sta[i].fieldSlice.GetFieldValue(T.X() - T.Tx() * dz, T.Y() - T.Ty() * dz, b0);
+ b0 = sta[i].fieldSlice.GetFieldValue(T.X() - T.Tx() * dz, T.Y() - T.Ty() * dz);
b0.Combine(fB[i], w[i]);
fld.Set(b0, z0, b1, z1, b2, z2);
@@ -470,8 +470,8 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
field.reserve(Tracks.size());
- L1FieldRegion fld _fvecalignment;
- L1FieldValue fB[3], fB_temp _fvecalignment;
+ ca::FieldRegion fld _fvecalignment;
+ ca::FieldValue fB[3], fB_temp _fvecalignment;
fvec zField[3];
unsigned short N_vTracks = Tracks.size();
@@ -526,7 +526,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
if (ista < 0) continue;
}
- sta[ista].fieldSlice.GetFieldValue(posx, posy, fB_temp);
+ fB_temp = sta[ista].fieldSlice.GetFieldValue(posx, posy);
fB[iH + 1].x[iVec] = fB_temp.x[iVec];
fB[iH + 1].y[iVec] = fB_temp.y[iVec];
fB[iH + 1].z[iVec] = fB_temp.z[iVec];
@@ -602,7 +602,7 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<PFF
field.reserve(Tracks.size());
- L1FieldRegion fld _fvecalignment;
+ ca::FieldRegion fld _fvecalignment;
int nTracks_SIMD = fvec::size();
TrackParamV T; // fitting parametr coresponding to current track
@@ -611,7 +611,7 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<PFF
int ista;
const L1Station* sta = CbmL1::Instance()->fpAlgo->GetParameters()->GetStations().begin();
- L1FieldValue fB[3], fB_temp _fvecalignment;
+ ca::FieldValue fB[3], fB_temp _fvecalignment;
fvec zField[3];
unsigned short N_vTracks = Tracks.size();
@@ -650,7 +650,7 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<PFF
if (ista < 0) continue;
}
- sta[ista].fieldSlice.GetFieldValue(posx, posy, fB_temp);
+ fB_temp = sta[ista].fieldSlice.GetFieldValue(posx, posy);
fB[iH + 1].x[iVec] = fB_temp.x[iVec];
fB[iH + 1].y[iVec] = fB_temp.y[iVec];
fB[iH + 1].z[iVec] = fB_temp.z[iVec];
@@ -673,7 +673,7 @@ void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks,
field.reserve(Tracks.size());
- L1FieldRegion fld _fvecalignment;
+ ca::FieldRegion fld _fvecalignment;
int nTracks_SIMD = fvec::size();
TrackParamV T; // fitting parametr coresponding to current track
@@ -682,7 +682,7 @@ void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks,
int ista;
const L1Station* sta = CbmL1::Instance()->fpAlgo->GetParameters()->GetStations().begin();
- L1FieldValue fB[3], fB_temp _fvecalignment;
+ ca::FieldValue fB[3], fB_temp _fvecalignment;
fvec zField[3];
unsigned short N_vTracks = Tracks.size();
@@ -724,7 +724,7 @@ void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks,
if (ista < 0) continue;
}
- sta[ista].fieldSlice.GetFieldValue(posx, posy, fB_temp);
+ fB_temp = sta[ista].fieldSlice.GetFieldValue(posx, posy);
fB[iH].x[iVec] = fB_temp.x[iVec];
fB[iH].y[iVec] = fB_temp.y[iVec];
diff --git a/reco/KF/ParticleFitter/CbmL1PFFitter.h b/reco/KF/ParticleFitter/CbmL1PFFitter.h
index 8fe2bfb3ce2f7d4635d11f190974484f01d67614..784cfa5df3591a35b63afd7679b10a191861a92a 100644
--- a/reco/KF/ParticleFitter/CbmL1PFFitter.h
+++ b/reco/KF/ParticleFitter/CbmL1PFFitter.h
@@ -25,18 +25,24 @@
class CbmMvdHit;
class CbmStsHit;
class CbmStsTrack;
-class L1FieldRegion;
+namespace cbm::algo::ca
+{
+ class FieldRegion;
+}
+
+using namespace cbm::algo;
+
class CbmKFVertex;
class TClonesArray;
class CbmL1PFFitter {
public:
- // A container for parameters of L1FieldRegion
+ // A container for parameters of ca::FieldRegion
struct PFFieldRegion {
PFFieldRegion() {}
- PFFieldRegion(const L1FieldRegion&, int i);
- void setFromL1FieldRegion(const L1FieldRegion&, int i);
- void getL1FieldRegion(L1FieldRegion&, int i);
+ PFFieldRegion(const ca::FieldRegion&, int i);
+ void setFromL1FieldRegion(const ca::FieldRegion&, int i);
+ void getL1FieldRegion(ca::FieldRegion&, int i);
float fP[10] {0.};
};
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index b72aa3c63ee12045d2f0c8bc9804f59f038644c5..885b80f460a737f74d6b0690fcbe4dade1f598ba 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -51,7 +51,6 @@ set(SRCS
L1Algo/L1MaterialMonitor.cxx
L1Algo/L1UMeasurementInfo.cxx
L1Algo/L1XYMeasurementInfo.cxx
- L1Algo/L1Field.cxx
L1Algo/L1CAIteration.cxx
L1Algo/L1BaseStationInfo.cxx
L1Algo/L1InitManager.cxx
@@ -218,7 +217,6 @@ install(FILES CbmL1Counters.h
install(FILES L1Algo/L1Algo.h
L1Algo/L1Branch.h
- L1Algo/L1Field.h
L1Algo/L1EArray.h
DESTINATION include/L1Algo
)
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index c9e50cf9b70b52b99d149efc8fcdce5d72744a3c..b5893002d8235d75cf3a3a345c1b6f24eb357eb0 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -64,13 +64,15 @@
#include <iostream>
#include <list>
+#include "CaField.h"
#include "CaHit.h"
#include "L1Algo/L1Algo.h"
#include "L1Algo/L1Assert.h"
#include "L1Algo/L1Branch.h"
-#include "L1Algo/L1Field.h"
#include "L1Event.h"
+using namespace cbm::algo; // TODO: remove this line
+
using std::cout;
using std::endl;
using std::ios;
@@ -622,6 +624,21 @@ InitStatus CbmL1::Init()
fpAlgo = &gAlgo;
fpAlgo->Init(fTrackingMode);
+
+ /// pass the original magnetic field to L1Algo
+
+ {
+ static auto fld = [&](double x, double y, double z) {
+ assert(FairRunAna::Instance());
+ double pos[3] = {x, y, z};
+ double B[3] = {0., 0., 0.};
+ if (FairRunAna::Instance()->GetField()) { FairRunAna::Instance()->GetField()->GetFieldValue(pos, B); }
+ return std::tuple(B[0], B[1], B[2]);
+ };
+ ca::FieldRegion::SetOriginalField(fld);
+ }
+
+
//
// ** Send formed parameters object to L1Algo instance **
//
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 55df90294d27021eb8c7b21eb7ef24cbe1b4d21c..378aee59928eaa558b5ad330b3e39b5c9e935435 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -64,7 +64,6 @@
class L1Algo;
class L1Event;
class CbmL1ParticlesFinder;
-class L1FieldSlice;
class CbmL1MCTrack;
class KFTopoPerformance;
class CbmMCDataObject;
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 016f136db6e1e30534efe60fa1efc538c827c9f8..9657daa9ff25c85ba3a405ea76cf526fc249226a 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -1378,7 +1378,7 @@ void CbmL1::TrackFitPerformance()
if (!mc.IsPrimary()) { // secondary
{ // extrapolate to vertex
- L1FieldRegion fld _fvecalignment;
+ ca::FieldRegion fld _fvecalignment;
fld.SetUseOriginalField();
fit.Extrapolate(mc.z, fld);
// add material
@@ -1441,7 +1441,7 @@ void CbmL1::TrackFitPerformance()
else { // primary
{ // extrapolate to vertex
- L1FieldRegion fld _fvecalignment;
+ ca::FieldRegion fld _fvecalignment;
fld.SetUseOriginalField();
// add material
@@ -1586,7 +1586,7 @@ void CbmL1::FillFitHistos(TrackParamV& track, const CbmL1MCPoint& mcP, bool isTi
//fit.SetMaxExtrapolationStep(10.);
fit.SetDoFitVelocity(true);
fit.SetTrack(track);
- L1FieldRegion fld _fvecalignment;
+ ca::FieldRegion fld _fvecalignment;
fld.SetUseOriginalField();
fit.Extrapolate(mcP.zOut, fld);
track = fit.Tr();
@@ -1662,8 +1662,8 @@ void CbmL1::FieldApproxCheck()
static_cast<int>(NbinsY + 1), -(Ymax + ddy / 2.), (Ymax + ddy / 2.));
Double_t r[3], B[3];
- L1FieldSlice FSl;
- L1FieldValue B_L1;
+ ca::FieldSlice FSl;
+ ca::FieldValue B_L1;
Double_t bbb, bbb_L1;
const int M = 5; // polinom order
@@ -1694,7 +1694,7 @@ void CbmL1::FieldApproxCheck()
MF->GetFieldValue(r, B);
bbb = sqrt(B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
- FSl.GetFieldValue(x, y, B_L1);
+ B_L1 = FSl.GetFieldValue(x, y);
bbb_L1 = sqrt(B_L1.x[0] * B_L1.x[0] + B_L1.y[0] * B_L1.y[0] + B_L1.z[0] * B_L1.z[0]);
stB->SetBinContent(ii, jj, (bbb - bbb_L1));
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index ae6c81db779ed1b8ccbbb53b5b3c4710695d1b34..1c127c2e2d96482bd4f6bee8e8e89258d11a7755 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -90,7 +90,7 @@ void L1Algo::ReceiveParameters(L1Parameters&& parameters)
fGhostSuppression = fParameters.GetGhostSuppression();
- L1FieldRegion::ForceUseOfOriginalField(fParameters.DevIsUseOfOriginalField());
+ ca::FieldRegion::ForceUseOfOriginalField(fParameters.DevIsUseOfOriginalField());
}
std::pair<fscal, fscal> L1Algo::GetHitCoorOnGrid(const ca::Hit& h)
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index f733cb4a745fc06cfa367362ff0a540bf76510b5..48e32fc97dd7a269576341fb995c815e631ccbcf 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -29,13 +29,13 @@ class L1AlgoDraw;
#include <map>
#include "CaConstants.h"
+#include "CaField.h"
#include "CaHit.h"
#include "CaTrack.h"
#include "CaTrackParam.h"
#include "CaVector.h"
#include "L1Branch.h"
#include "L1CloneMerger.h"
-#include "L1Field.h"
#include "L1Fit.h"
#include "L1Grid.h"
#include "L1HitPoint.h"
@@ -382,7 +382,7 @@ private:
bool fIsTargetField {false}; ///< is the magnetic field present at the target
- L1FieldValue fTargB _fvecalignment {}; // field in the target point (modifiable, do not touch!!)
+ ca::FieldValue fTargB _fvecalignment {}; // field in the target point (modifiable, do not touch!!)
L1XYMeasurementInfo TargetXYInfo _fvecalignment {}; // target constraint [cm]
int fGhostSuppression {1}; // NOTE: Should be equal to 0 in TRACKS_FROM_TRIPLETS mode!
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.h b/reco/L1/L1Algo/L1BaseStationInfo.h
index 5c0966f39ccc1db5b2d220e441c462ef3a949072..2d79889236311d962d739b27da3aba9f54fc3f3e 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.h
+++ b/reco/L1/L1Algo/L1BaseStationInfo.h
@@ -54,7 +54,7 @@ public:
kZmin, ///< min z of the station
kZmax, ///< max z of the station
kThicknessMap, ///< thickness map of the station (optional?)
- kFieldSlice, ///< L1Station.L1FieldSlice object initialization
+ kFieldSlice, ///< L1Station.ca::FieldSlice object initialization
// The last item is equal to the number of bits in fInitFlags
kEnd
};
diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/reco/L1/L1Algo/L1BranchExtender.cxx
index 69dd750fcb86b29714848a7424d107298c3a77cb..d95994ce780653de1b02e17c58c381795e1a1232 100644
--- a/reco/L1/L1Algo/L1BranchExtender.cxx
+++ b/reco/L1/L1Algo/L1BranchExtender.cxx
@@ -86,16 +86,15 @@ void L1Algo::BranchFitterFast(const L1Branch& t, TrackParamV& Tout, const bool u
T.C10() = hit0.dxy;
T.C11() = hit0.dy2;
- L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
- L1FieldRegion fld _fvecalignment;
+ ca::FieldRegion fld _fvecalignment;
fvec fldZ0 = sta1.fZ; // suppose field is smoth
fvec fldZ1 = sta2.fZ;
fvec fldZ2 = sta0.fZ;
- sta1.fieldSlice.GetFieldValue(x1, y1, fldB0);
- sta2.fieldSlice.GetFieldValue(x2, y2, fldB1);
- sta0.fieldSlice.GetFieldValue(x0, y0, fldB2);
+ ca::FieldValue fldB0 = sta1.fieldSlice.GetFieldValue(x1, y1);
+ ca::FieldValue fldB1 = sta2.fieldSlice.GetFieldValue(x2, y2);
+ ca::FieldValue fldB2 = sta0.fieldSlice.GetFieldValue(x0, y0);
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
@@ -114,7 +113,7 @@ void L1Algo::BranchFitterFast(const L1Branch& t, TrackParamV& Tout, const bool u
fldB1 = fldB2;
fldZ0 = fldZ1;
fldZ1 = fldZ2;
- sta.fieldSlice.GetFieldValue(hit.x, hit.y, fldB2);
+ fldB2 = sta.fieldSlice.GetFieldValue(hit.x, hit.y);
fldZ2 = sta.fZ;
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
} // i
@@ -173,15 +172,14 @@ void L1Algo::FindMoreHits(L1Branch& t, TrackParamV& Tout, const bool upstream, c
fvec x2 = hit2.x;
fvec y2 = hit2.y;
- L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
- L1FieldRegion fld _fvecalignment;
+ ca::FieldRegion fld _fvecalignment;
fvec fldZ0 = sta1.fZ;
fvec fldZ1 = sta2.fZ;
fvec fldZ2 = sta0.fZ;
- sta1.fieldSlice.GetFieldValue(x1, y1, fldB0);
- sta2.fieldSlice.GetFieldValue(x2, y2, fldB1);
- sta0.fieldSlice.GetFieldValue(x0, y0, fldB2);
+ ca::FieldValue fldB0 = sta1.fieldSlice.GetFieldValue(x1, y1);
+ ca::FieldValue fldB1 = sta2.fieldSlice.GetFieldValue(x2, y2);
+ ca::FieldValue fldB2 = sta0.fieldSlice.GetFieldValue(x0, y0);
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
@@ -275,7 +273,7 @@ void L1Algo::FindMoreHits(L1Branch& t, TrackParamV& Tout, const bool upstream, c
fldB1 = fldB2;
fldZ0 = fldZ1;
fldZ1 = fldZ2;
- sta.fieldSlice.GetFieldValue(hit.x, hit.y, fldB2);
+ fldB2 = sta.fieldSlice.GetFieldValue(hit.x, hit.y);
fldZ2 = sta.fZ;
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
}
diff --git a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
index db87816857c3ea3d261bf23cedbd797a15c8ffb1..9468c04a9d2b6f0e21162046dfe67575c6d63931 100644
--- a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
+++ b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
@@ -326,7 +326,7 @@ void L1Algo::CaTrackFinderSlice()
// Select magnetic field. For primary tracks - fVtxFieldValue, for secondary tracks - st.fieldSlice
if (caIteration.GetPrimaryFlag()) { fTargB = fParameters.GetVertexFieldValue(); }
else {
- fParameters.GetStation(0).fieldSlice.GetFieldValue(0, 0, fTargB);
+ fTargB = fParameters.GetStation(0).fieldSlice.GetFieldValue(0, 0);
} // NOTE: calculates field fTargB in the center of 0th station
fIsTargetField = (fabs(fTargB.y[0]) > 0.001);
diff --git a/reco/L1/L1Algo/L1CloneMerger.cxx b/reco/L1/L1Algo/L1CloneMerger.cxx
index 9d5fefaea44f533c662d221b5b86addf468d23a5..1d5b9f44e4bcd943633aa9f19571e101da76cccb 100644
--- a/reco/L1/L1Algo/L1CloneMerger.cxx
+++ b/reco/L1/L1Algo/L1CloneMerger.cxx
@@ -92,8 +92,8 @@ void L1CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extRe
TrackParamV& Tb = fitB.Tr();
TrackParamV& Tf = fitF.Tr();
- L1FieldValue fBm, fBb, fBf _fvecalignment;
- L1FieldRegion fld _fvecalignment;
+ ca::FieldValue fBm, fBb, fBf _fvecalignment;
+ ca::FieldRegion fld _fvecalignment;
// Max length for merging
unsigned char maxLengthForMerge = static_cast<unsigned char>(frAlgo.GetParameters()->GetNstationsActive() - 3);
@@ -125,8 +125,8 @@ void L1CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extRe
unsigned short stam;
- frAlgo.GetParameters()->GetStation(staf).fieldSlice.GetFieldValue(Tf.X(), Tf.Y(), fBf);
- frAlgo.GetParameters()->GetStation(stab).fieldSlice.GetFieldValue(Tb.X(), Tb.Y(), fBb);
+ fBf = frAlgo.GetParameters()->GetStation(staf).fieldSlice.GetFieldValue(Tf.X(), Tf.Y());
+ fBb = frAlgo.GetParameters()->GetStation(stab).fieldSlice.GetFieldValue(Tb.X(), Tb.Y());
unsigned short dist = firstStation[iTr] - lastStation[jTr];
@@ -137,7 +137,7 @@ void L1CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extRe
fvec zm = frAlgo.GetParameters()->GetStation(stam).fZ;
fvec xm = fvec(0.5) * (Tf.GetX() + Tf.Tx() * (zm - Tf.Z()) + Tb.GetX() + Tb.Tx() * (zm - Tb.Z()));
fvec ym = fvec(0.5) * (Tf.Y() + Tf.Ty() * (zm - Tf.Z()) + Tb.Y() + Tb.Ty() * (zm - Tb.Z()));
- frAlgo.GetParameters()->GetStation(stam).fieldSlice.GetFieldValue(xm, ym, fBm);
+ fBm = frAlgo.GetParameters()->GetStation(stam).fieldSlice.GetFieldValue(xm, ym);
fld.Set(fBb, Tb.Z(), fBm, zm, fBf, Tf.Z());
fvec zMiddle = fvec(0.5) * (Tb.Z() + Tf.Z());
diff --git a/reco/L1/L1Algo/L1Field.h b/reco/L1/L1Algo/L1Field.h
deleted file mode 100644
index dbbb4e601a3287ce65d7eb6f220c2a676fa51102..0000000000000000000000000000000000000000
--- a/reco/L1/L1Algo/L1Field.h
+++ /dev/null
@@ -1,230 +0,0 @@
-/* Copyright (C) 2007-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer], Igor Kulakov, Maksym Zyzak, Sergei Zharko */
-
-#ifndef L1Field_h
-#define L1Field_h 1
-
-#include <boost/serialization/access.hpp>
-
-#include <string>
-
-#include "CaConstants.h"
-#include "CaSimd.h"
-#include "CaTrackParam.h"
-
-using namespace cbm::algo::ca; // TODO: remove
-
-class L1FieldValue {
-
-public:
- fvec x {0.f}; //< x-component of the field
- fvec y {0.f}; //< y-component of the field
- fvec z {0.f}; //< z-component of the field
-
- /// Combines the magnetic field with another field value object using weight
- /// \param B other field value to combine with
- /// \param w weight from 0 to 1
- void Combine(const L1FieldValue& B, const fvec& w);
-
- void Combine(const L1FieldValue& B, const fmask& w);
-
- /// Consistency checker
- void CheckConsistency() const;
-
- /// Operator << overloading
- friend std::ostream& operator<<(std::ostream& out, const L1FieldValue& B);
-
- /// String representation of class contents
- /// \param indentLevel number of indent characters in the output
- std::string ToString(int indentLevel) const;
-
- /// Serialization function
- friend class boost::serialization::access;
- template<class Archive>
- void serialize(Archive& ar, const unsigned int)
- {
- ar& x;
- ar& y;
- ar& z;
- }
-} _fvecalignment;
-
-[[gnu::always_inline]] inline void L1FieldValue::Combine(const L1FieldValue& B, const fvec& w)
-{
- x += w * (B.x - x);
- y += w * (B.y - y);
- z += w * (B.z - z);
-}
-
-[[gnu::always_inline]] inline void L1FieldValue::Combine(const L1FieldValue& B, const fmask& w)
-{
- x(w) = B.x;
- y(w) = B.y;
- z(w) = B.z;
-}
-
-/// Class represents a set of magnetic field approximation coefficients
-///
-// TODO: Crosscheck the default content (S.Zharko)
-class L1FieldSlice {
-public:
- /// Default constructor
- L1FieldSlice();
-
- /// Consistency checker
- void CheckConsistency() const;
-
- /// Gets field value from (x, y) fvec point
- /// \param x x-coordinate of input
- /// \param y y-coordinate of input
- /// \param B the L1FieldValue output
- void GetFieldValue(const fvec& x, const fvec& y, L1FieldValue& B) const;
-
- /// Gets field value for the intersection with a straight track
- /// \param t straight track
- /// \param B the L1FieldValue output
- void GetFieldValueForLine(const cbm::algo::ca::TrackParamV& t, L1FieldValue& B) const;
-
- /// String representation of class contents
- /// \param indentLevel number of indent characters in the output
- std::string ToString(int indentLevel = 0) const;
-
-public:
- // NOTE: We don't use an initialization of arrays here because we cannot be sure
- // if the underlying type (fvec) has a default constructor, but
- // we are sure, that it can be initialized with a float. (S.Zharko)
- static constexpr auto kMaxNFieldApproxCoefficients = cbm::algo::ca::constants::size::MaxNFieldApproxCoefficients;
- fvec cx[kMaxNFieldApproxCoefficients]; ///< Polynomial coefficients for x-component of the field value
- fvec cy[kMaxNFieldApproxCoefficients]; ///< Polynomial coefficients for y-component of the field value
- fvec cz[kMaxNFieldApproxCoefficients]; ///< Polynomial coefficients for z-component of the field value
-
- fvec z {constants::Undef<fscal>}; ///< z coordinate of the slice
-
- /// Serialization function
- friend class boost::serialization::access;
- template<class Archive>
- void serialize(Archive& ar, const unsigned int)
- {
- ar& cx;
- ar& cy;
- ar& cz;
- ar& z;
- }
-} _fvecalignment;
-
-
-class L1FieldRegion {
-public:
- // NOTE: When a custom constructor is defined, default constructor also should be provided (S.Zharko)
- L1FieldRegion() = default;
-
- L1FieldRegion(float reg[10]) noexcept;
-
- /// Consistency checker
- void CheckConsistency() const;
-
- /// Gets the field vector and writes it into B pointer
- /// \param x x-coordinate of the point to calculate the field
- /// \param y y-coordinate of the point to calculate the field
- /// \param z z-coordinate of the point to calculate the field
- /// \param B pointer to the output fvec array of the magnetic field
- void Get(const fvec& x, const fvec& y, const fvec& z, fvec* B) const;
-
- /// Interpolates the magnetic field by three nodal points and sets the result to this L1FieldRegion object
- /// The result is a quadratic interpolation of the field as a function of z
- /// \param b0 field value in the first nodal point
- /// \param b0z z-coordinate of the first nodal point
- /// \param b1 field value in the second nodal point
- /// \param b1z z-coordinate of the second nodal point
- /// \param b2 field value in the third nodal point
- /// \param b2z z-coordinate of the third nodal point
- /// TODO: does the sequence of b0z, b1z and b2z matter? If yes, probalby we need an assertion (S.Zharko)
- void Set(const L1FieldValue& b0, const fvec b0z, const L1FieldValue& b1, const fvec b1z, const L1FieldValue& b2,
- const fvec b2z);
-
- /// Interpolates the magnetic field by thwo nodal points and sets the result to this L1FieldRegion object.
- /// The result is a linear interpolation of the field as a function of z
- /// \param b0 field value in the first nodal point
- /// \param b0z z-coordinate of the first nodal point
- /// \param b1 field value in the second nodal point
- /// \param b1z z-coordinate of the second nodal point
- /// TODO: does the sequence of b0z and b1z matter? If yes, probalby we need an assertion (S.Zharko)
- void Set(const L1FieldValue& b0, const fvec b0z, const L1FieldValue& b1, const fvec b1z);
-
- /// Shifts the coefficients to new central point
- /// \param z z-coordinate of the new central point
- void Shift(fvec z);
-
- /// Replaces the selected layer of the coefficients with one from another
- /// L1FieldRegion object
- /// \param i0 the index of the fvect layer in this L1FieldRegion object
- /// \param fl the other L1FieldRegion object, which layer is going to be replaced
- /// \param i1 the index of the source fvect layer to copy
- void SetOneEntry(const int i0, const L1FieldRegion& f1, const int i1);
-
- /// Replaces all the layers of the coefficients with one selected layer from another
- /// L1FieldRegion object
- /// \param fl the other L1FieldRegion object, which layer is going to be replaced
- /// \param i1 the index of the source fvect layer to copy
- void SetOneEntry(const L1FieldRegion& f1, const int i1);
-
- /// Saves the contents of the particular layer of the coefficients into an array of floats
- /// \param reg output array of 10 floats
- /// \param iVec index of the input
- // TODO: Probably, it would be better to rename "reg" into "output" and make it the second
- // parameter of this function (S.Zharko)
- // TODO: Probably we need a const specifier here, because the method does not change the fields
- void GetOneEntry(float reg[10], const int iVec);
-
- /// String representation of class contents
- /// \param indentLevel number of indent characters in the output
- std::string ToString(int indentLevel = 0) const;
-
- /// Use original field instead of the field approximation
- void SetUseOriginalField(bool v = true) { fUseOriginalField = v; }
-
- /// Force using the original field
- static void ForceUseOfOriginalField(bool v = true) { fgkForceUseOfOriginalField = v; }
-
-public:
- static bool fgkForceUseOfOriginalField;
-
- // TODO: Probably it's better to have arrays instead of separate fvec values? (S.Zharko)
- // Bx(z) = cx0 + cx1*(z-z0) + cx2*(z-z0)^2
- fvec cx0 {0.f};
- fvec cx1 {0.f};
- fvec cx2 {0.f};
-
- // By(z) = cy0 + cy1*(z-z0) + cy2*(z-z0)^2
- fvec cy0 {0.f};
- fvec cy1 {0.f};
- fvec cy2 {0.f};
-
- // Bz(z) = cz0 + cz1*(z-z0) + cz2*(z-z0)^2
- fvec cz0 {0.f};
- fvec cz1 {0.f};
- fvec cz2 {0.f};
-
- fvec z0 {0.f}; ///< z-coordinate of the field region central point
-
- bool fUseOriginalField {false};
-
- /// Serialization function
- friend class boost::serialization::access;
- template<class Archive>
- void serialize(Archive& ar, const unsigned int)
- {
- ar& cx0;
- ar& cx1;
- ar& cx2;
- ar& cy0;
- ar& cy1;
- ar& cy2;
- ar& cz0;
- ar& cz1;
- ar& cz2;
- }
-} _fvecalignment;
-
-#endif
diff --git a/reco/L1/L1Algo/L1Fit.cxx b/reco/L1/L1Algo/L1Fit.cxx
index 993ceb2774cb8e60fc95bdc7c961650e75328e03..2b9b5b770b176d02f200fa4634f450d000d0d0b2 100644
--- a/reco/L1/L1Algo/L1Fit.cxx
+++ b/reco/L1/L1Algo/L1Fit.cxx
@@ -580,7 +580,7 @@ void L1Fit::FilterVi(fvec vi)
void L1Fit::Extrapolate // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
(cnst& z_out, // extrapolate to this z position
- const L1FieldRegion& F)
+ const ca::FieldRegion& F)
{
// use Q/p linearisation at fQp0
@@ -594,7 +594,7 @@ void L1Fit::Extrapolate // extrapolates track parameters and returns jacobian f
void L1Fit::ExtrapolateStepFull // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
(cnst& zOut, // extrapolate to this z position
- const L1FieldRegion& Field)
+ const ca::FieldRegion& Field)
{
// use Q/p linearisation at fQp0
// implementation of the Runge-Kutta method without optimization
@@ -641,7 +641,7 @@ void L1Fit::ExtrapolateStepFull // extrapolates track parameters and returns ja
cnst stepDz[5] = {0., 0., h * fvec(0.5), h * fvec(0.5), h};
- fvec f[5][7] = {{0.}}; // ( d*/dz ) [step]
+ fvec f[5][7] = {{0.}}; // ( d*/dz ) [step]
fvec F[5][7][7] = {{{0.}}}; // ( d *new [step] / d *old )
// Runge-Kutta steps
@@ -661,11 +661,10 @@ void L1Fit::ExtrapolateStepFull // extrapolates track parameters and returns ja
}
fvec z = fTr.GetZ() + stepDz[step];
- fvec B[3];
- cnst& Bx = B[0];
- cnst& By = B[1];
- cnst& Bz = B[2];
- Field.Get(rstep[0], rstep[1], z, B);
+ ca::FieldValue B = Field.Get(rstep[0], rstep[1], z);
+ cnst& Bx = B.x;
+ cnst& By = B.y;
+ cnst& Bz = B.z;
fvec tx = rstep[2];
fvec ty = rstep[3];
@@ -793,7 +792,7 @@ void L1Fit::ExtrapolateStepFull // extrapolates track parameters and returns ja
void L1Fit::ExtrapolateStep // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
(cnst& zOut, // extrapolate to this z position
- const L1FieldRegion& F)
+ const ca::FieldRegion& F)
{
// use Q/p linearisation at fQp0
@@ -866,12 +865,11 @@ void L1Fit::ExtrapolateStep // extrapolates track parameters and returns jacobi
}
}
- fvec B[3];
- cnst& Bx = B[0];
- cnst& By = B[1];
- cnst& Bz = B[2];
- F.Get(r[0], r[1], z, B);
+ ca::FieldValue B = F.Get(r[0], r[1], z);
+ cnst& Bx = B.x;
+ cnst& By = B.y;
+ cnst& Bz = B.z;
fvec tx = r[2];
fvec ty = r[3];
@@ -1133,7 +1131,7 @@ void L1Fit::ExtrapolateStep // extrapolates track parameters and returns jacobi
void
L1Fit::ExtrapolateStepAnalytic // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
(cnst& z_out, // extrapolate to this z position
- const L1FieldRegion& F)
+ const ca::FieldRegion& F)
{
//
// Part of the analytic extrapolation formula with error (c_light*B*dz)^4/4!
@@ -1397,7 +1395,7 @@ void
//cout<<"Extrapolation ok"<<endl;
}
-void L1Fit::ExtrapolateLine(cnst& z_out, const L1FieldRegion& F)
+void L1Fit::ExtrapolateLine(cnst& z_out, const ca::FieldRegion& F)
{
// extrapolate the track assuming fQp0 == 0
// TODO: write special simplified procedure
@@ -1654,7 +1652,7 @@ void L1Fit::EnergyLossCorrectionAl(cnst& radThick, cnst& direction)
EnergyLossCorrection<atomicZ>(atomicA, rho, radLen, radThick, direction);
}
-void L1Fit::GetExtrapolatedXYline(cnst& z, const L1FieldRegion& F, fvec& extrX, fvec& extrY, fvec Jx[6],
+void L1Fit::GetExtrapolatedXYline(cnst& z, const ca::FieldRegion& F, fvec& extrX, fvec& extrY, fvec Jx[6],
fvec Jy[6]) const
{
// extrapolate track assuming it is straight (qp==0)
@@ -1703,7 +1701,7 @@ void L1Fit::GetExtrapolatedXYline(cnst& z, const L1FieldRegion& F, fvec& extrX,
void L1Fit::AddTargetToLine(cnst& targX, cnst& targY, cnst& targZ, const L1XYMeasurementInfo& targXYInfo,
- const L1FieldRegion& F)
+ const ca::FieldRegion& F)
{
// Add the target constraint to a straight line track
diff --git a/reco/L1/L1Algo/L1Fit.h b/reco/L1/L1Algo/L1Fit.h
index 31ad3feef49c2e6001019a641b87759042480e0e..a5d77ebb223ffcae7f2fff79871413ead5fab9b0 100644
--- a/reco/L1/L1Algo/L1Fit.h
+++ b/reco/L1/L1Algo/L1Fit.h
@@ -11,9 +11,9 @@
#ifndef L1Fit_h
#define L1Fit_h
+#include "CaField.h"
#include "CaTrackParam.h"
#include "L1Def.h"
-#include "L1Field.h"
#include "L1UMeasurementInfo.h"
#include "L1XYMeasurementInfo.h"
@@ -108,13 +108,13 @@ public:
void MeasureVelocityWithQp();
- void Extrapolate(cnst& z_out, const L1FieldRegion& F);
+ void Extrapolate(cnst& z_out, const ca::FieldRegion& F);
- void ExtrapolateStep(cnst& z_out, const L1FieldRegion& F);
- void ExtrapolateStepFull(cnst& z_out, const L1FieldRegion& F);
- void ExtrapolateStepAnalytic(cnst& z_out, const L1FieldRegion& F);
+ void ExtrapolateStep(cnst& z_out, const ca::FieldRegion& F);
+ void ExtrapolateStepFull(cnst& z_out, const ca::FieldRegion& F);
+ void ExtrapolateStepAnalytic(cnst& z_out, const ca::FieldRegion& F);
void ExtrapolateLineNoField(cnst& z_out);
- void ExtrapolateLine(cnst& z_out, const L1FieldRegion& F);
+ void ExtrapolateLine(cnst& z_out, const ca::FieldRegion& F);
void EnergyLossCorrection(cnst& radThick, cnst& upstreamDirection);
@@ -129,7 +129,7 @@ public:
void AddMsInThickMaterial(cnst& radThick, cnst& thickness, bool fDownstream);
- void GetExtrapolatedXYline(cnst& z, const L1FieldRegion& F, fvec& extrX, fvec& extrY, fvec Jx[6], fvec Jy[6]) const;
+ void GetExtrapolatedXYline(cnst& z, const ca::FieldRegion& F, fvec& extrX, fvec& extrY, fvec Jx[6], fvec Jy[6]) const;
void ExtrapolateXC00Line(cnst& z_out, fvec& x, fvec& C00) const;
void ExtrapolateYC11Line(cnst& z_out, fvec& y, fvec& C11) const;
@@ -137,7 +137,7 @@ public:
void AddTargetToLine(cnst& targX, cnst& targY, cnst& targZ, const L1XYMeasurementInfo& targXYInfo,
- const L1FieldRegion& F);
+ const ca::FieldRegion& F);
static fvec ApproximateBetheBloch(cnst& bg2);
diff --git a/reco/L1/L1Algo/L1InitManager.cxx b/reco/L1/L1Algo/L1InitManager.cxx
index 73ff58facb345113fb8271cb92606409e5fdaead..0b931ce4a749fa7c9fa17d029036d7bb3192f069 100644
--- a/reco/L1/L1Algo/L1InitManager.cxx
+++ b/reco/L1/L1Algo/L1InitManager.cxx
@@ -61,8 +61,8 @@ void L1InitManager::ClearSetupInfo()
fInitController.SetFlag(EInitKey::kActiveDetectorIDs, false);
// Clear field info
- fParameters.fVertexFieldRegion = L1FieldRegion();
- fParameters.fVertexFieldValue = L1FieldValue();
+ fParameters.fVertexFieldRegion = ca::FieldRegion();
+ fParameters.fVertexFieldValue = ca::FieldValue();
fInitController.SetFlag(EInitKey::kPrimaryVertexField, false);
// Clear target position
@@ -258,7 +258,7 @@ void L1InitManager::InitTargetField(double zStep)
constexpr int nPointsNodal {3};
std::array<double, nPointsNodal> inputNodalZ {fTargetZ, fTargetZ + zStep, fTargetZ + 2. * zStep};
- std::array<L1FieldValue, nPointsNodal> B {};
+ std::array<ca::FieldValue, nPointsNodal> B {};
std::array<fvec, nPointsNodal> z {};
// loop over nodal points
for (int idx = 0; idx < nPointsNodal; ++idx) {
diff --git a/reco/L1/L1Algo/L1InitManager.h b/reco/L1/L1Algo/L1InitManager.h
index 7fee948313cfe55e7397cb80c4e6130e3a9d7786..2f3f47a55c4d76556cf7fe6b0b3cf5f044a486f5 100644
--- a/reco/L1/L1Algo/L1InitManager.h
+++ b/reco/L1/L1Algo/L1InitManager.h
@@ -18,11 +18,11 @@
#include <unordered_map>
#include "CaConstants.h"
+#include "CaField.h"
#include "CaSimd.h"
#include "L1BaseStationInfo.h"
#include "L1CAIteration.h"
#include "L1EArray.h"
-#include "L1Field.h"
#include "L1ObjectInitController.h"
#include "L1Parameters.h"
#include "L1Utils.h"
@@ -183,7 +183,7 @@ public:
/// available for modifications.
void InitStationLayout();
- /// @brief Calculates L1FieldValue and L1FieldReference values for a selected step in z-axis from the target position
+ /// @brief Calculates ca::FieldValue and L1FieldReference values for a selected step in z-axis from the target position
/// \param zStep step between nodal points
// TODO: Consider possibility for linear approximation (S.Zh.)
void InitTargetField(double zStep);
diff --git a/reco/L1/L1Algo/L1Parameters.h b/reco/L1/L1Algo/L1Parameters.h
index 719b4271d9a2ec675c8a3c448a47010072367403..2e450c97b205dace586c63bcfa93270f62dfc906 100644
--- a/reco/L1/L1Algo/L1Parameters.h
+++ b/reco/L1/L1Algo/L1Parameters.h
@@ -19,6 +19,7 @@
#include <utility>
#include "CaConstants.h"
+#include "CaField.h"
#include "CaVector.h"
#include "L1CAIteration.h"
#include "L1Material.h"
@@ -26,6 +27,8 @@
#include "L1Station.h"
#include "L1Utils.h" // IS DEPRECATED?
+using namespace cbm::algo; // TODO: remove
+
class L1InitManager;
enum class L1DetectorID;
@@ -177,11 +180,11 @@ public:
/// Gets Z component of target position
fvec GetTargetPositionZ() const { return fTargetPos[2]; }
- /// Gets L1FieldRegion object at primary vertex
- const L1FieldRegion& GetVertexFieldRegion() const { return fVertexFieldRegion; }
+ /// Gets ca::FieldRegion object at primary vertex
+ const ca::FieldRegion& GetVertexFieldRegion() const { return fVertexFieldRegion; }
- /// Gets L1FieldValue object at primary vertex
- const L1FieldValue& GetVertexFieldValue() const { return fVertexFieldValue; }
+ /// Gets ca::FieldValue object at primary vertex
+ const ca::FieldValue& GetVertexFieldValue() const { return fVertexFieldValue; }
/// @brief Gets random seed
///
@@ -231,10 +234,10 @@ private:
L1Utils::kNaN};
/// Field value object at primary vertex (between target and the first station)
- alignas(constants::misc::Alignment) L1FieldValue fVertexFieldValue {};
+ alignas(constants::misc::Alignment) ca::FieldValue fVertexFieldValue {};
/// Field region object at primary vertex (between target and the first station)
- alignas(constants::misc::Alignment) L1FieldRegion fVertexFieldRegion {};
+ alignas(constants::misc::Alignment) ca::FieldRegion fVertexFieldRegion {};
/// Array of stations
alignas(constants::misc::Alignment) L1StationsContainer_t fStations {};
diff --git a/reco/L1/L1Algo/L1Station.h b/reco/L1/L1Algo/L1Station.h
index 58c00f5e8d7363f0d39554381dcb6daffb9cdf38..22ceb1ed726c0878155a76eea96704e8a91b4145 100644
--- a/reco/L1/L1Algo/L1Station.h
+++ b/reco/L1/L1Algo/L1Station.h
@@ -9,10 +9,13 @@
#include <type_traits>
#include "CaConstants.h"
+#include "CaField.h"
#include "CaSimd.h"
-#include "L1Field.h"
#include "L1Utils.h"
+using namespace cbm::algo; // TODO: Remove this line
+using namespace cbm::algo::ca; // TODO: Remove this line
+
/// Structure L1Station
/// Contains a set of geometry parameters for a particular station
///
@@ -21,15 +24,15 @@ public:
// TODO: SZh 12.05.2022: Rewrite type into L1DetectorID, change detector indexing scheme
// TODO: SZh 12.05.2022: Provide getters to stations
- int type = cbm::algo::ca::constants::Undef<int>; // ? Detector type?
- int timeInfo = cbm::algo::ca::constants::Undef<int>; ///< flag: if time information can be used
+ int type = ca::constants::Undef<int>; // ? Detector type?
+ int timeInfo = ca::constants::Undef<int>; ///< flag: if time information can be used
int fieldStatus =
- cbm::algo::ca::constants::Undef<int>; ///< flag: 1 - station is INSIDE the field, 0 - station is OUTSIDE the field
- fvec fZ = cbm::algo::ca::constants::Undef<fvec>; ///< z position of station [cm]
- fvec Xmax = cbm::algo::ca::constants::Undef<fvec>; ///< min radius of the station [cm]
- fvec Ymax = cbm::algo::ca::constants::Undef<fvec>; ///< max radius of the station [cm]
+ ca::constants::Undef<int>; ///< flag: 1 - station is INSIDE the field, 0 - station is OUTSIDE the field
+ fvec fZ = ca::constants::Undef<fvec>; ///< z position of station [cm]
+ fvec Xmax = ca::constants::Undef<fvec>; ///< min radius of the station [cm]
+ fvec Ymax = ca::constants::Undef<fvec>; ///< max radius of the station [cm]
- L1FieldSlice fieldSlice {};
+ ca::FieldSlice fieldSlice {};
// Serialization block
friend class boost::serialization::access;
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index aea32e357e7ab8ea435ee206e2fe36892ab593b3..2e3f2d582e8b66e819256f8ada31b00196a98f88 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -26,14 +26,14 @@ void L1Algo::L1KFTrackFitter()
// cout << " Start L1 Track Fitter " << endl;
int start_hit = 0; // for interation in fSliceRecoHits[]
- // static L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
- // static L1FieldRegion fld _fvecalignment;
- L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
- L1FieldRegion fld _fvecalignment;
+ // static ca::FieldValue fldB0, fldB1, fldB2 _fvecalignment;
+ // static ca::FieldRegion fld _fvecalignment;
+ ca::FieldValue fldB0, fldB1, fldB2 _fvecalignment;
+ ca::FieldRegion fld _fvecalignment;
- L1FieldValue fldB01, fldB11, fldB21 _fvecalignment;
- L1FieldRegion fld1 _fvecalignment;
+ ca::FieldValue fldB01, fldB11, fldB21 _fvecalignment;
+ ca::FieldRegion fld1 _fvecalignment;
const int nStations = fParameters.GetNstationsActive();
int nTracks_SIMD = fvec::size();
@@ -86,7 +86,7 @@ void L1Algo::L1KFTrackFitter()
fvec z_start;
fvec z_end;
- L1FieldValue fB[constants::size::MaxNstations], fB_temp _fvecalignment;
+ ca::FieldValue fB[constants::size::MaxNstations], fB_temp _fvecalignment;
fvec ZSta[constants::size::MaxNstations];
@@ -143,7 +143,7 @@ void L1Algo::L1KFTrackFitter()
dt2[ista][iVec] = hit.dt2;
if (!sta[ista].timeInfo) { dt2[ista][iVec] = 1.e4; }
z[ista][iVec] = hit.z;
- sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista], fB_temp);
+ 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;
@@ -211,13 +211,13 @@ void L1Algo::L1KFTrackFitter()
fldZ1 = z[ista];
- sta[ista].fieldSlice.GetFieldValue(tr.X(), tr.Y(), fldB1);
+ fldB1 = sta[ista].fieldSlice.GetFieldValue(tr.X(), tr.Y());
fldB1.Combine(fB[ista], w[ista]);
fldZ2 = z[ista - 2];
fvec dz = fldZ2 - fldZ1;
- sta[ista].fieldSlice.GetFieldValue(tr.X() + tr.Tx() * dz, tr.Y() + tr.Ty() * dz, fldB2);
+ fldB2 = sta[ista].fieldSlice.GetFieldValue(tr.X() + tr.Tx() * dz, tr.Y() + tr.Ty() * dz);
fldB2.Combine(fB[ista - 2], w[ista - 2]);
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
@@ -225,7 +225,7 @@ void L1Algo::L1KFTrackFitter()
fldZ0 = z[ista];
dz = (fldZ1 - fldZ0);
- sta[ista].fieldSlice.GetFieldValue(tr.X() - tr.Tx() * dz, tr.Y() - tr.Ty() * dz, fldB0);
+ fldB0 = sta[ista].fieldSlice.GetFieldValue(tr.X() - tr.Tx() * dz, tr.Y() - tr.Ty() * dz);
fldB0.Combine(fB[ista], w[ista]);
fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2);
@@ -312,19 +312,19 @@ void L1Algo::L1KFTrackFitter()
fit.SetQp0(tr.Qp());
fldZ1 = z[ista];
- sta[ista].fieldSlice.GetFieldValue(tr.X(), tr.Y(), fldB1);
+ fldB1 = sta[ista].fieldSlice.GetFieldValue(tr.X(), tr.Y());
fldB1.Combine(fB[ista], w[ista]);
fldZ2 = z[ista + 2];
dz = fldZ2 - fldZ1;
- sta[ista].fieldSlice.GetFieldValue(tr.X() + tr.Tx() * dz, tr.Y() + tr.Ty() * dz, fldB2);
+ fldB2 = sta[ista].fieldSlice.GetFieldValue(tr.X() + tr.Tx() * dz, tr.Y() + tr.Ty() * dz);
fldB2.Combine(fB[ista + 2], w[ista + 2]);
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
for (++ista; ista < nStations; ista++) {
fldZ0 = z[ista];
dz = (fldZ1 - fldZ0);
- sta[ista].fieldSlice.GetFieldValue(tr.X() - tr.Tx() * dz, tr.Y() - tr.Ty() * dz, fldB0);
+ fldB0 = sta[ista].fieldSlice.GetFieldValue(tr.X() - tr.Tx() * dz, tr.Y() - tr.Ty() * dz);
fldB0.Combine(fB[ista], w[ista]);
fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2);
diff --git a/reco/L1/L1Algo/L1TripletConstructor.cxx b/reco/L1/L1Algo/L1TripletConstructor.cxx
index 729355fcf28fd9c8bd27356aaaa3a4ad5b79eb01..636f4b63e3eb82d7fe26b928dffabd889522c362 100644
--- a/reco/L1/L1Algo/L1TripletConstructor.cxx
+++ b/reco/L1/L1Algo/L1TripletConstructor.cxx
@@ -98,8 +98,8 @@ void L1TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar,
fFit.SetQp0(fAlgo->fMaxInvMom);
- L1FieldValue lB, mB, cB, rB _fvecalignment;
- L1FieldValue l_B_global, centB_global _fvecalignment;
+ ca::FieldValue lB, mB, cB, rB _fvecalignment;
+ ca::FieldValue l_B_global, centB_global _fvecalignment;
// get the magnetic field along the track.
// (suppose track is straight line with origin in the target)
@@ -142,20 +142,18 @@ void L1TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar,
// field made by the left hit, the target and the station istac in-between.
// is used for extrapolation to the target and to the middle hit
- L1FieldRegion fld0;
+ ca::FieldRegion fld0;
{
- L1FieldValue B0, B1;
- fFld0Sta[0]->fieldSlice.GetFieldValueForLine(T, B0);
- fFld0Sta[1]->fieldSlice.GetFieldValueForLine(T, B1);
+ ca::FieldValue B0 = fFld0Sta[0]->fieldSlice.GetFieldValueForLine(T);
+ ca::FieldValue B1 = fFld0Sta[1]->fieldSlice.GetFieldValueForLine(T);
fld0.Set(fAlgo->fTargB, fAlgo->fTargZ, B0, fFld0Sta[0]->fZ, B1, fFld0Sta[1]->fZ);
}
{ // field, made by the left hit, the middle station and the right station
// Will be used for extrapolation to the right hit
- L1FieldValue B0, B1, B2;
- fFld1Sta[0]->fieldSlice.GetFieldValueForLine(T, B0);
- fFld1Sta[1]->fieldSlice.GetFieldValueForLine(T, B1);
- fFld1Sta[2]->fieldSlice.GetFieldValueForLine(T, B2);
+ ca::FieldValue B0 = fFld1Sta[0]->fieldSlice.GetFieldValueForLine(T);
+ ca::FieldValue B1 = fFld1Sta[1]->fieldSlice.GetFieldValueForLine(T);
+ ca::FieldValue B2 = fFld1Sta[2]->fieldSlice.GetFieldValueForLine(T);
fFldL.Set(B0, fFld1Sta[0]->fZ, B1, fFld1Sta[1]->fZ, B2, fFld1Sta[2]->fZ);
}
@@ -435,15 +433,15 @@ void L1TripletConstructor::FitTriplets()
// find the field along the track
- L1FieldValue B[3] _fvecalignment;
- L1FieldRegion fld _fvecalignment;
- L1FieldRegion fldTarget _fvecalignment;
+ ca::FieldValue B[3] _fvecalignment;
+ ca::FieldRegion fld _fvecalignment;
+ ca::FieldRegion fldTarget _fvecalignment;
fvec tx[3] = {(x[1] - x[0]) / (z[1] - z[0]), (x[2] - x[0]) / (z[2] - z[0]), (x[2] - x[1]) / (z[2] - z[1])};
fvec ty[3] = {(y[1] - y[0]) / (z[1] - z[0]), (y[2] - y[0]) / (z[2] - z[0]), (y[2] - y[1]) / (z[2] - z[1])};
for (int ih = 0; ih < NHits; ++ih) {
fvec dz = (sta[ih].fZ - z[ih]);
- sta[ih].fieldSlice.GetFieldValue(x[ih] + tx[ih] * dz, y[ih] + ty[ih] * dz, B[ih]);
+ B[ih] = sta[ih].fieldSlice.GetFieldValue(x[ih] + tx[ih] * dz, y[ih] + ty[ih] * dz);
};
fld.Set(B[0], sta[0].fZ, B[1], sta[1].fZ, B[2], sta[2].fZ);
diff --git a/reco/L1/L1Algo/L1TripletConstructor.h b/reco/L1/L1Algo/L1TripletConstructor.h
index bfd88d1304e7874154e4aebc3ce34e5f6edbc600..c6d795717f4279fd3ccf369eda06154f6ad1b00b 100644
--- a/reco/L1/L1Algo/L1TripletConstructor.h
+++ b/reco/L1/L1Algo/L1TripletConstructor.h
@@ -5,10 +5,10 @@
#ifndef L1TripletConstructor_h
#define L1TripletConstructor_h
+#include "CaField.h"
#include "CaTrackParam.h"
#include "CaVector.h"
#include "L1Algo.h"
-#include "L1Field.h"
#include "L1Fit.h"
#include "L1HitPoint.h"
#include "L1Station.h"
@@ -103,7 +103,7 @@ private:
ca::HitIndex_t fIhitL;
TrackParamV fTrL;
- L1FieldRegion fFldL;
+ ca::FieldRegion fFldL;
Vector<ca::HitIndex_t> fHitsM_2 {"L1TripletConstructor::fHitsM_2"};
Vector<TrackParamV> fTracks_2 {"L1TripletConstructor::fTracks_2"};
diff --git a/reco/L1/L1Algo/L1Utils.h b/reco/L1/L1Algo/L1Utils.h
index 41ad13f8f5ee56ad8ac54a610224beb72dc94f3f..9764c4133abce8b6c328beaca27860cf20ee2978 100644
--- a/reco/L1/L1Algo/L1Utils.h
+++ b/reco/L1/L1Algo/L1Utils.h
@@ -24,6 +24,7 @@
#include <cmath>
#include "CaSimd.h"
+#include "CaUtils.h"
#include "L1Def.h"
// TODO: SZh 16.05.2023: Split this class into several headers and place them into L1Algo::utils
@@ -47,18 +48,7 @@ namespace L1Utils
|| fabs(lhs - rhs) < std::numeric_limits<T>::min();
}
- [[gnu::always_inline]] inline void CheckSimdVectorEquality(cbm::algo::ca::fvec v, const char* name)
- {
- bool ok = true;
- for (size_t i = 1; i < cbm::algo::ca::fvec::size(); i++) {
- ok = ok && (v[i] == v[0]);
- }
- if (!ok) {
- std::stringstream msg;
- msg << name << " SIMD vector is inconsistent, not all of the words are equal each other: " << v;
- throw std::logic_error(msg.str());
- }
- }
+ using cbm::algo::ca::utils::CheckSimdVectorEquality;
/// Hash for unordered_map with enum class keys
struct EnumClassHash {
diff --git a/reco/L1/L1Algo/utils/L1CADebug.h b/reco/L1/L1Algo/utils/L1CADebug.h
index e6144f242eb495244513b8fbd0728dde7ab79962..40b610167e7a79d0f81879b3ca2f3d9f6f943bfe 100644
--- a/reco/L1/L1Algo/utils/L1CADebug.h
+++ b/reco/L1/L1Algo/utils/L1CADebug.h
@@ -15,7 +15,7 @@ TH1F *h_pick_res_x[20][3][3], *h_pick_res_y[20][3][3], *h_pick_pull_x[20][3][3],
*h_pull_qp[20][3][3];
//TH2F *h_dyvsy[20], *h_dxvsx[20];
-void Pulls(int i, int j, int k, double* mc, TrackParamV& T, fvec qp0, L1FieldRegion& fld)
+void Pulls(int i, int j, int k, double* mc, TrackParamV& T, fvec qp0, ca::FieldRegion& fld)
{
TrackParamV tmp_T = T;
fvec z = mc[5];
diff --git a/reco/L1/qa/CbmCaTrackFitQa.cxx b/reco/L1/qa/CbmCaTrackFitQa.cxx
index 2f4125cf3347e95e7ef61fddadb24798e399f3f2..980dbf9e8d3dc93a3c6b18b27418356be655c92e 100644
--- a/reco/L1/qa/CbmCaTrackFitQa.cxx
+++ b/reco/L1/qa/CbmCaTrackFitQa.cxx
@@ -20,8 +20,8 @@
#include <algorithm>
+#include "CaField.h"
#include "CaToolsMCData.h"
-#include "L1Field.h"
#include "L1Fit.h"
#include "L1Utils.h"
@@ -33,6 +33,8 @@
using cbm::ca::TrackFitQa;
using cbm::ca::tools::MCPoint;
+using namespace cbm::algo; // TODO: SZh 16.05.2023: Remove this line
+
// ---------------------------------------------------------------------------------------------------------------------
//
TrackFitQa::TrackFitQa(const char* pointTag, const char* prefix, std::shared_ptr<ObjList_t> pObjList)
@@ -137,7 +139,7 @@ void TrackFitQa::Fill(const TrackParamV& trPar, const tools::MCPoint& mcPoint, b
fitter.SetMask(fmask::One());
fitter.SetDoFitVelocity(true);
fitter.SetTrack(trPar);
- L1FieldRegion fieldRegion;
+ cbm::algo::ca::FieldRegion fieldRegion;
fieldRegion.SetUseOriginalField(); // Precised magnetic field is used
fitter.Extrapolate(mcPoint.GetZOut(), fieldRegion);
diff --git a/reco/L1/qa/CbmCaTrackFitQa.h b/reco/L1/qa/CbmCaTrackFitQa.h
index 6d626c981369c722e3baf6cf6bc02095861e41fa..ba3e2b901bde6d7491075d5831c199f5d008d32c 100644
--- a/reco/L1/qa/CbmCaTrackFitQa.h
+++ b/reco/L1/qa/CbmCaTrackFitQa.h
@@ -15,8 +15,8 @@
#include <array>
#include "CaConstants.h"
+#include "CaField.h"
#include "L1EArray.h"
-#include "L1Field.h"
#include "L1Fit.h"
// Forward declarations