diff --git a/algo/ca/core/pars/CaField.cxx b/algo/ca/core/pars/CaField.cxx
index 99b8ca5f21d171921d287aa6f25d4825692bab69..dcfcbdb625f536d4739d001e2cf8284501c311ee 100644
--- a/algo/ca/core/pars/CaField.cxx
+++ b/algo/ca/core/pars/CaField.cxx
@@ -3,225 +3,3 @@
Authors: Sergey Gorbunov, Sergei Zharko [committer] */
#include "CaField.h"
-
-#include "KfUtils.h"
-
-#include <iomanip>
-#include <iostream>
-#include <sstream>
-
-using namespace cbm::algo::ca;
-
-
-template<typename DataT>
-bool FieldRegion<DataT>::fgForceUseOfOriginalField = false;
-
-template<typename DataT>
-std::function<std::tuple<double, double, double>(double x, double y, double z)>
- FieldRegion<DataT>::fgOriginalField([](double, double, double) {
- assert(false && "ca::FieldRegion: The original field is not set, use cbm::ca::tools::SetCbmField()");
- return std::tuple(kfdefs::Undef<double>, kfdefs::Undef<double>, kfdefs::Undef<double>);
- });
-
-template<typename DataT>
-bool FieldRegion<DataT>::fgIsZeroOriginalField = false;
-
-
-//
-// FieldRegion methdos
-//
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-template<typename DataT>
-FieldRegion<DataT>::FieldRegion(float reg[10]) noexcept
- : cx0(reg[0])
- , cx1(reg[1])
- , cx2(reg[2])
- , cy0(reg[3])
- , cy1(reg[4])
- , cy2(reg[5])
- , cz0(reg[6])
- , cz1(reg[7])
- , cz2(reg[8])
- , z0(reg[9])
-{
-}
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-template<typename DataT>
-void FieldRegion<DataT>::CheckConsistency() const
-{
- // Check SIMD data vectors for consistent initialization
- kfutils::CheckSimdVectorEquality(cx0, "FieldRegion::cx0");
- kfutils::CheckSimdVectorEquality(cx1, "FieldRegion::cx1");
- kfutils::CheckSimdVectorEquality(cx2, "FieldRegion::cx2");
- kfutils::CheckSimdVectorEquality(cy0, "FieldRegion::cy0");
- kfutils::CheckSimdVectorEquality(cy1, "FieldRegion::cy1");
- kfutils::CheckSimdVectorEquality(cy2, "FieldRegion::cy2");
- kfutils::CheckSimdVectorEquality(cz0, "FieldRegion::cz0");
- kfutils::CheckSimdVectorEquality(cz1, "FieldRegion::cz1");
- kfutils::CheckSimdVectorEquality(cz2, "FieldRegion::cz2");
- kfutils::CheckSimdVectorEquality(z0, "FieldRegion::z0");
- // TODO: Any other checks? (S.Zharko)
-}
-
-
-//----------------------------------------------------------------------------------------------------------------------
-// TODO: Should it be inline? (S.Zharko)
-template<typename DataT>
-FieldValue<DataT> FieldRegion<DataT>::Get(const DataT& x, const DataT& y, const DataT& z) const
-{
- FieldValue<DataT> B;
- if (fgForceUseOfOriginalField || fUseOriginalField) {
- // TODO: ? Why the simd words are set indvidually?
- for (size_t i = 0; i < kfutils::simd::Size<DataT>(); i++) {
- auto [bx, by, bz] =
- fgOriginalField(kfutils::simd::Cast<DataT, double>(x, i), kfutils::simd::Cast<DataT, double>(y, i),
- kfutils::simd::Cast<DataT, double>(z, i));
- B.SetSimdEntry(bx, by, bz, i);
- }
- }
- else {
- DataT dz = z - z0;
- DataT dz2 = dz * dz;
- B.Set(cx0 + cx1 * dz + cx2 * dz2, // Bx
- cy0 + cy1 * dz + cy2 * dz2, // By
- cz0 + cz1 * dz + cz2 * dz2); // Bz
- }
- return B;
-}
-
-
-//----------------------------------------------------------------------------------------------------------------------
-template<typename T>
-std::tuple<T, T, T> FieldRegion<T>::GetDoubleIntegrals(const T& /*x1*/, const T& /*y1*/, const T& z1, //
- const T& /*x2*/, const T& /*y2*/, const T& z2) const
-{
- // double integral of the field along z
-
- if (fUseOriginalField) {
- // TODO: implement the double integral for the original field
- assert(false && "FieldRegion::GetDoubleIntegrals() is not implemented for the original field");
- }
-
- auto fld = *this;
- fld.Shift(z1);
- T dz = z2 - z1;
- T dz2 = dz * dz;
- T c0 = dz2 * T(1. / 2.);
- T c1 = dz2 * dz * T(1. / 6.);
- T c2 = dz2 * dz2 * T(1. / 12.);
- return {c0 * fld.cx0 + c1 * fld.cx1 + c2 * fld.cx2, //
- c0 * fld.cy0 + c1 * fld.cy1 + c2 * fld.cy2, //
- c0 * fld.cz0 + c1 * fld.cz1 + c2 * fld.cz2};
-}
-
-
-//----------------------------------------------------------------------------------------------------------------------
-// TODO: Should it be inline? (S.Zharko)
-template<typename DataT>
-void FieldRegion<DataT>::Set(const FieldValue<DataT>& b0, const DataT b0z, const FieldValue<DataT>& b1, const DataT b1z,
- const FieldValue<DataT>& b2, const DataT b2z)
-{
- fIsZeroField = (b0.IsZero() && b1.IsZero() && b2.IsZero());
-
- z0 = b0z;
- DataT dz1 = b1z - b0z;
- DataT dz2 = b2z - b0z;
- DataT det = kfutils::simd::One<DataT>() / (dz1 * dz2 * (b2z - b1z));
-
- DataT w21 = -dz2 * det;
- DataT w22 = dz1 * det;
- DataT w11 = -dz2 * w21;
- DataT w12 = -dz1 * w22;
-
- DataT db1 = b1.GetBx() - b0.GetBx();
- DataT db2 = b2.GetBx() - b0.GetBx();
- cx0 = b0.GetBx();
- cx1 = db1 * w11 + db2 * w12;
- cx2 = db1 * w21 + db2 * w22;
-
- db1 = b1.GetBy() - b0.GetBy();
- db2 = b2.GetBy() - b0.GetBy();
- cy0 = b0.GetBy();
- cy1 = db1 * w11 + db2 * w12;
- cy2 = db1 * w21 + db2 * w22;
-
- db1 = b1.GetBz() - b0.GetBz();
- db2 = b2.GetBz() - b0.GetBz();
- cz0 = b0.GetBz();
- cz1 = db1 * w11 + db2 * w12;
- cz2 = db1 * w21 + db2 * w22;
-}
-
-//----------------------------------------------------------------------------------------------------------------------
-// TODO: Should it be inline? (S.Zharko)
-//template<typename DataT>
-//void FieldRegion<DataT>::Set(const FieldValue<DataT>& b0, const DataT b0z, const FieldValue<DataT>& b1, const DataT b1z)
-//{
-// fIsZeroField = (b0.IsZeroField() && b1.IsZeroField());
-//
-// z0 = b0z;
-// DataT dzi = kfutils::simd::One<DataT>() / (b1z - b0z);
-// cx0 = b0.x;
-// cy0 = b0.y;
-// cz0 = b0.z;
-// cx1 = (b1.x - b0.x) * dzi;
-// cy1 = (b1.y - b0.y) * dzi;
-// cz1 = (b1.z - b0.z) * dzi;
-// cx2 = 0.f;
-// cy2 = 0.f;
-// cz2 = 0.f;
-//}
-
-//----------------------------------------------------------------------------------------------------------------------
-// TODO: Should it be inline? (S.Zharko)
-template<typename DataT>
-void FieldRegion<DataT>::Shift(DataT z)
-{
- DataT dz = z - z0;
- DataT cx2dz = cx2 * dz;
- DataT cy2dz = cy2 * dz;
- DataT cz2dz = cz2 * dz;
- z0 = z;
- cx0 += (cx1 + cx2dz) * dz;
- cy0 += (cy1 + cy2dz) * dz;
- cz0 += (cz1 + cz2dz) * dz;
- cx1 += cx2dz + cx2dz;
- cy1 += cy2dz + cy2dz;
- cz1 += cz2dz + cz2dz;
-}
-
-//----------------------------------------------------------------------------------------------------------------------
-//
-template<typename DataT>
-std::string FieldRegion<DataT>::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 << "Bx(z) [kG] components: " << '\n';
- aStream << indent << indentChar << "cx0: " << cx0 << '\n';
- aStream << indent << indentChar << "cx1: " << cx1 << '\n';
- aStream << indent << indentChar << "cx2: " << cx2 << '\n';
- aStream << indent << "By(z) [kG] components: " << '\n';
- aStream << indent << indentChar << "cy0: " << cy0 << '\n';
- aStream << indent << indentChar << "cy1: " << cy1 << '\n';
- aStream << indent << indentChar << "cy2: " << cy2 << '\n';
- aStream << indent << "Bz(z) [kG] components: " << '\n';
- aStream << indent << indentChar << "cz0: " << cz0 << '\n';
- aStream << indent << indentChar << "cz1: " << cz1 << '\n';
- aStream << indent << indentChar << "cz2: " << cz2 << '\n';
- aStream << indent << "z0 [cm]: " << z0;
- return aStream.str();
-}
-
-namespace cbm::algo::ca
-{
- template class FieldRegion<fvec>;
- template class FieldRegion<float>;
- template class FieldRegion<double>;
-} // namespace cbm::algo::ca
diff --git a/algo/ca/core/pars/CaField.h b/algo/ca/core/pars/CaField.h
index 07d3cee5c1bb41e808203b2af7ee4cc9fa7e91b9..18c73b1ad6946e40cab4ccfcac3b8331be8d661f 100644
--- a/algo/ca/core/pars/CaField.h
+++ b/algo/ca/core/pars/CaField.h
@@ -7,221 +7,13 @@
#include "CaDefs.h"
#include "CaSimd.h"
+#include "KfFieldRegion.h"
#include "KfFieldValue.h"
#include "KfTrackParam.h"
-#include <boost/serialization/access.hpp>
-
-#include <functional>
-#include <string>
-
namespace cbm::algo::ca
{
+ using cbm::algo::kf::FieldRegion;
using cbm::algo::kf::FieldSlice;
using cbm::algo::kf::FieldValue;
-
- //template<typename T>
- //using FieldRegion = cbm::algo::kf::FieldRegion<T, cbm::algo::kf::EFieldMode::Intrpl>;
-
- template<typename DataT>
- 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;
-
- /// \brief Copy constructor with type conversion
- template<typename DataIn>
- FieldRegion(const FieldRegion<DataIn>& other)
- : cx0(kfutils::simd::Cast<DataIn, DataT>(other.cx0))
- , cx1(kfutils::simd::Cast<DataIn, DataT>(other.cx1))
- , cx2(kfutils::simd::Cast<DataIn, DataT>(other.cx2))
- , cy0(kfutils::simd::Cast<DataIn, DataT>(other.cy0))
- , cy1(kfutils::simd::Cast<DataIn, DataT>(other.cy1))
- , cy2(kfutils::simd::Cast<DataIn, DataT>(other.cy2))
- , cz0(kfutils::simd::Cast<DataIn, DataT>(other.cz0))
- , cz1(kfutils::simd::Cast<DataIn, DataT>(other.cz1))
- , cz2(kfutils::simd::Cast<DataIn, DataT>(other.cz2))
- , z0(kfutils::simd::Cast<DataIn, DataT>(other.z0))
- , fUseOriginalField(other.fUseOriginalField)
- {
- }
-
- /// 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<DataT> Get(const DataT& x, const DataT& y, const DataT& z) const;
-
- std::tuple<DataT, DataT, DataT> GetDoubleIntegrals(const DataT& x1, const DataT& y1, const DataT& z1, //
- const DataT& x2, const DataT& y2, const DataT& z2) const;
-
- /// Is the field region empty?
- bool IsZeroField() const
- {
- if (fgForceUseOfOriginalField || fUseOriginalField) {
- return fgIsZeroOriginalField;
- }
- return fIsZeroField;
- }
-
- /// 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<DataT>& b0, const DataT b0z, const FieldValue<DataT>& b1, const DataT b1z,
- const FieldValue<DataT>& b2, const DataT 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<DataT>& b0, const DataT b0z, const FieldValue<DataT>& b1, const DataT b1z);
-
- /// Shifts the coefficients to new central point
- /// \param z z-coordinate of the new central point
- void Shift(DataT z);
-
- /// Set one SIMD entry from one SIMD entry of the other class
- /// It only works when DataT is fvec, TdataB is fvec
- /// \param fb the other FieldRegion object, which layer is going to be replaced
- /// \param ib the index of the source fvec layer to copy
- template<typename TdataB>
- void SetSimdEntry(const FieldRegion<TdataB>& fb, const int ib)
- {
- CopyBase<TdataB, true, false>(0, fb, ib);
- }
-
- /// Set one SIMD entry from one SIMD entry of the other class
- /// It only works when DataT is fvec, TdataB is fvec
- /// \param ia the index of the fvec layer in this FieldRegion object
- /// \param fb the other FieldRegion object, which layer is going to be replaced
- /// \param ib the index of the source fvec layer to copy
- //void SetSimdEntry(const int ia, const FieldRegion<fvec>& fb, const int ib)
- //{
- // CopyBase<fvec, false, false>(ia, fb, ib);
- //}
-
- /// 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,
- bool isZeroField)
- {
- fgOriginalField = f;
- fgIsZeroOriginalField = isZeroField;
- }
-
- /// \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 FieldRegion<TdataB>& Tb, const int ib);
-
- public:
- static bool fgForceUseOfOriginalField;
- static std::function<std::tuple<double, double, double>(double x, double y, double z)> fgOriginalField;
- static bool fgIsZeroOriginalField;
-
- // 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
- DataT cx0{0.f};
- DataT cx1{0.f};
- DataT cx2{0.f};
-
- // By(z) = cy0 + cy1*(z-z0) + cy2*(z-z0)^2
- DataT cy0{0.f};
- DataT cy1{0.f};
- DataT cy2{0.f};
-
- // Bz(z) = cz0 + cz1*(z-z0) + cz2*(z-z0)^2
- DataT cz0{0.f};
- DataT cz1{0.f};
- DataT cz2{0.f};
-
- DataT z0{0.f}; ///< z-coordinate of the field region central point
-
- bool fUseOriginalField{false}; ///< Use original field instead of the field approximation
-
- bool fIsZeroField{false}; ///< Is the field region empty?
-
- /// 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;
- ar& z0;
- ar& fUseOriginalField;
- ar& fIsZeroField;
- }
- } _fvecalignment;
-
- // ---------------------------------------------------------------------------------------------------------------------
- //
- template<typename TdataA>
- template<typename TdataB, bool TDoAllA, bool TDoAllB>
- inline void FieldRegion<TdataA>::CopyBase(const int ia, const FieldRegion<TdataB>& Tb, const int ib)
- {
- auto copy = [&](TdataA& a, const TdataB& b) {
- kfutils::VecCopy<TdataA, TdataB, TDoAllA, TDoAllB>::CopyEntries(a, ia, b, ib);
- };
-
- copy(cx0, Tb.cx0);
- copy(cx1, Tb.cx1);
- copy(cx2, Tb.cx2);
-
- copy(cy0, Tb.cy0);
- copy(cy1, Tb.cy1);
- copy(cy2, Tb.cy2);
-
- copy(cz0, Tb.cz0);
- copy(cz1, Tb.cz1);
- copy(cz2, Tb.cz2);
-
- copy(z0, Tb.z0);
-
- fUseOriginalField = Tb.fUseOriginalField;
- fIsZeroField = Tb.fIsZeroField;
-
- } // CopyBase
-
- /// \brief Explicit instantiation declarations for FieldRegion class with specific template types to ensure proper fgOriginalField definition
- extern template class FieldRegion<fvec>;
- extern template class FieldRegion<float>;
- extern template class FieldRegion<double>;
-
} // namespace cbm::algo::ca
diff --git a/algo/ca/core/tracking/CaFramework.cxx b/algo/ca/core/tracking/CaFramework.cxx
index eae683f55b19c6a8ee7a38788714dfc0ba198f37..422f282f51704d1985bf7002af116d82322dc115 100644
--- a/algo/ca/core/tracking/CaFramework.cxx
+++ b/algo/ca/core/tracking/CaFramework.cxx
@@ -61,7 +61,6 @@ void Framework::ReceiveParameters(Parameters<fvec>&& parameters)
fParameters = std::move(parameters);
fNstationsBeforePipe = fParameters.GetNstationsActive(static_cast<EDetectorID>(0));
- ca::FieldRegion<fvec>::ForceUseOfOriginalField(fParameters.DevIsUseOfOriginalField());
kf::GlobalField::ForceUseOfOriginalField(fParameters.DevIsUseOfOriginalField());
}
diff --git a/algo/ca/core/tracking/CaTrackFit.cxx b/algo/ca/core/tracking/CaTrackFit.cxx
index d31d4aa2c1ea1f5f5280d6bcd0fc0af3c527d39d..1ab9c9809312e14e4ffc5cccf7a10af5b5da2296 100644
--- a/algo/ca/core/tracking/CaTrackFit.cxx
+++ b/algo/ca/core/tracking/CaTrackFit.cxx
@@ -619,7 +619,7 @@ namespace cbm::algo::ca
{
// use Q/p linearisation at fQp0
- if (F.IsZeroField()) {
+ if (F.GetFieldType() == kf::EFieldType::Null) {
ExtrapolateLineNoField(z_out);
}
else {
diff --git a/algo/ca/core/tracking/CaTrackFitter.cxx b/algo/ca/core/tracking/CaTrackFitter.cxx
index a2f3214b9b1a92a3f15b238c5f4627491c693e47..7f49cc9927ad4554054729709e84328354e198aa 100644
--- a/algo/ca/core/tracking/CaTrackFitter.cxx
+++ b/algo/ca/core/tracking/CaTrackFitter.cxx
@@ -288,8 +288,7 @@ void TrackFitter::FitCaTracks(const ca::InputData& input, WindowData& wData)
vtxInfo.SetDy2(1.e-8);
if (ca::TrackingMode::kGlobal == fTrackingMode) {
- ca::FieldRegion<fvec> fldFull _fvecalignment;
- fldFull.SetUseOriginalField(true);
+ ca::FieldRegion<fvec> fldFull(kf::GlobalField::fgOriginalFieldType, kf::GlobalField::fgOriginalField);
fitpv.SetMaxExtrapolationStep(1.);
for (int vtxIter = 0; vtxIter < 2; vtxIter++) {
fitpv.SetQp0(fitpv.Tr().Qp());
diff --git a/algo/kf/core/geo/KfField.h b/algo/kf/core/geo/KfField.h
index 1e07c1afb20a1ca60fb63d3533c13a205b4f11e7..e1227faa5abe9246ad616484052f4578bcf4e8d6 100644
--- a/algo/kf/core/geo/KfField.h
+++ b/algo/kf/core/geo/KfField.h
@@ -268,7 +268,7 @@ namespace cbm::algo::kf
FieldRegion<T> GetFieldRegion(const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1, const T& z1,
const FieldValue<T>& b2, const T& z2) const
{
- return (fFieldMode == EFieldMode::Intrpl ? FieldRegion<T>(fFieldType, b0, z0, b1, z1, b2, z2)
+ return (fFieldMode == EFieldMode::Intrpl ? FieldRegion<T>(b0, z0, b1, z1, b2, z2)
: FieldRegion<T>(fFieldType, foFldOrig->GetFieldFunction()));
}
diff --git a/algo/kf/core/geo/KfFieldRegion.cxx b/algo/kf/core/geo/KfFieldRegion.cxx
index f87839f61eea8e777a7fff03ef57dab0b5b98ba5..340e22a18d6b5378857c434fd3bf61473ef04bde 100644
--- a/algo/kf/core/geo/KfFieldRegion.cxx
+++ b/algo/kf/core/geo/KfFieldRegion.cxx
@@ -24,9 +24,12 @@ using cbm::algo::kf::FieldValue;
using cbm::algo::kf::GlobalField;
using cbm::algo::kf::detail::FieldRegionBase;
+using namespace cbm::algo;
+FieldFn_t GlobalField::fgOriginalField = &GlobalField::UndefField;
+
+kf::EFieldType GlobalField::fgOriginalFieldType{kf::EFieldType::Null};
-FieldFn_t GlobalField::gOriginalField = &GlobalField::UndefField;
bool GlobalField::fgForceUseOfOriginalField = false;
@@ -63,7 +66,7 @@ template<typename T>
FieldValue<T> FieldRegionBase<T, EFieldMode::Intrpl>::Get(const T& x, const T& y, const T& z) const
{
if (GlobalField::IsUsingOriginalFieldForced()) {
- return GlobalField::GetFieldValue(GlobalField::gOriginalField, x, y, z);
+ return GlobalField::GetFieldValue(GlobalField::fgOriginalField, x, y, z);
}
auto dz = z - this->fZfirst;
return FieldValue<T>{this->fCoeff[0][0] + dz * (this->fCoeff[0][1] + dz * this->fCoeff[0][2]),
@@ -71,6 +74,30 @@ FieldValue<T> FieldRegionBase<T, EFieldMode::Intrpl>::Get(const T& x, const T& y
this->fCoeff[2][0] + dz * (this->fCoeff[2][1] + dz * this->fCoeff[2][2])};
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+std::tuple<T, T, T>
+FieldRegionBase<T, EFieldMode::Intrpl>::GetDoubleIntegrals(const T& /*x1*/, const T& /*y1*/, const T& z1, //
+ const T& /*x2*/, const T& /*y2*/, const T& z2) const
+{
+ // double integral of the field along z
+
+ if (GlobalField::IsUsingOriginalFieldForced()) {
+ // TODO: implement the double integral for the original field
+ }
+ auto fld = *this;
+ fld.Shift(z1);
+ T dz = z2 - z1;
+ T dz2 = dz * dz;
+ T c0 = dz2 * T(1. / 2.);
+ T c1 = dz2 * dz * T(1. / 6.);
+ T c2 = dz2 * dz2 * T(1. / 12.);
+ return {c0 * fld.fCoeff[0][0] + c1 * fld.fCoeff[0][1] + c2 * fld.fCoeff[0][2],
+ c0 * fld.fCoeff[1][0] + c1 * fld.fCoeff[1][1] + c2 * fld.fCoeff[1][2],
+ c0 * fld.fCoeff[2][0] + c1 * fld.fCoeff[2][1] + c2 * fld.fCoeff[2][2]};
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T>
@@ -166,17 +193,6 @@ namespace cbm::algo::kf::detail
} // namespace cbm::algo::kf::detail
-// ---------------------------------------------------------------------------------------------------------------------
-//
-template<typename T>
-void FieldRegion<T>::Set(const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1, const T& z1,
- const FieldValue<T>& b2, const T& z2)
-{
- if (fFieldMode == EFieldMode::Intrpl) {
- foFldIntrpl->Set(b0, z0, b1, z1, b2, z2);
- }
-}
-
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T>
diff --git a/algo/kf/core/geo/KfFieldRegion.h b/algo/kf/core/geo/KfFieldRegion.h
index 0ca1f66f9ca1fa46e8c552725ffa1d2a911ec73b..299aa5fb853fbe5f941135daa1e9897098792248 100644
--- a/algo/kf/core/geo/KfFieldRegion.h
+++ b/algo/kf/core/geo/KfFieldRegion.h
@@ -33,14 +33,21 @@ namespace cbm::algo::kf
class GlobalField {
public:
/// \brief Global variable to store the fielf funciton (x,y,z)->(Bx,By,Bz)
- static FieldFn_t gOriginalField;
+ static FieldFn_t fgOriginalField;
+
+ /// \brief Global field type
+ static EFieldType fgOriginalFieldType; ///< global field type
/// \brief Undefined magnetic field function
static std::tuple<double, double, double> UndefField(double, double, double);
/// \brief Sets global field function
/// \param fn Field function of the FairRoot format
- static void SetFieldFunction(const FieldFn_t& fn) { gOriginalField = fn; }
+ static void SetFieldFunction(EFieldType fldType, const FieldFn_t& fn)
+ {
+ fgOriginalFieldType = fldType;
+ fgOriginalField = fn;
+ }
/// \brief Creates a field value object in a spactial point of a generic floating-point type
/// \tparam T Floating point data-type of the involved variables
@@ -112,6 +119,14 @@ namespace cbm::algo::kf
return GlobalField::GetFieldValue(fFieldFn, x, y, z);
}
+ /// \brief Gets the double integrals of the field along the track
+ std::tuple<T, T, T> GetDoubleIntegrals(const T& /*x1*/, const T& /*y1*/, const T& /*z1*/, //
+ const T& /*x2*/, const T& /*y2*/, const T& /*z2*/) const
+ {
+ assert("cbm::algo::kf::FieldRegionBase<T, EFieldMode::Orig>::GetDoubleIntegrals() is not implemented");
+ return {defs::Undef<T>, defs::Undef<T>, defs::Undef<T>};
+ }
+
/// Consistency checker
void CheckConsistency() const;
@@ -178,9 +193,16 @@ namespace cbm::algo::kf
/// \param z z-coordinate of the point [cm]
FieldValue<T> Get(const T& x, const T& y, const T& z) const;
+ /// \brief Gets the double integrals of the field along the track
+ std::tuple<T, T, T> GetDoubleIntegrals(const T& x1, const T& y1, const T& z1, //
+ const T& x2, const T& y2, const T& z2) const;
+
/// \brief Gets the first z-coordinate of the interpolation
const T& GetZfirst() const { return fZfirst; }
+ /// \brief Gets the coefficients of the polynomial approximation
+ const auto& GetCoefficients() const { return fCoeff; }
+
/// \brief Interpolates the field by three nodal points with the second order polynomial
/// \note The interpolation is carried out vs. z-coordinates of the space points of the field values
/// \param b0 Field value at z0 [kG]
@@ -192,6 +214,13 @@ namespace cbm::algo::kf
void Set(const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1, const T& z1, const FieldValue<T>& b2,
const T& z2);
+ /// \brief Set the coefficients of the polynomial approximation
+ void Set(const CoeffArray_t& coeff, const T& z0)
+ {
+ fCoeff = coeff;
+ fZfirst = z0;
+ }
+
/// \brief Shifts the coefficients to another first z-coordinate
/// \param z new z-coordinate [cm]
void Shift(const T& z);
@@ -231,7 +260,7 @@ namespace cbm::algo::kf
public:
/// \brief Constructor of the generic field
- FieldRegion(EFieldMode fldMode, EFieldType fldType)
+ FieldRegion(EFieldMode fldMode = EFieldMode::Intrpl, EFieldType fldType = EFieldType::Normal)
: foFldIntrpl(fldMode == EFieldMode::Intrpl ? std::make_optional<detail::FieldRegionBase<T, EFieldMode::Intrpl>>()
: std::nullopt)
, foFldOrig(fldMode == EFieldMode::Orig ? std::make_optional<detail::FieldRegionBase<T, EFieldMode::Orig>>()
@@ -256,11 +285,22 @@ namespace cbm::algo::kf
{
}
+ /// \brief Constructor for the interpolated field region
+ FieldRegion(const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1, const T& z1, const FieldValue<T>& b2,
+ const T& z2)
+ : foFldIntrpl(std::make_optional<detail::FieldRegionBase<T, EFieldMode::Intrpl>>())
+ , foFldOrig(std::nullopt)
+ , fFieldType(EFieldType::Normal)
+ , fFieldMode(EFieldMode::Intrpl)
+ {
+ // Set the field type automatically
+ Set(b0, z0, b1, z1, b2, z2);
+ }
+
/// \brief Constructor for the interpolated field region
/// \param fieldType Type of the field
- FieldRegion(EFieldType fieldType, const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1, const T& z1,
- const FieldValue<T>& b2, const T& z2)
- : foFldIntrpl(std::make_optional<detail::FieldRegionBase<T, EFieldMode::Intrpl>>(b0, z0, b1, z1, b2, z2))
+ FieldRegion(EFieldType fieldType, const detail::FieldRegionBase<T, EFieldMode::Intrpl>& fld)
+ : foFldIntrpl(std::make_optional(fld))
, foFldOrig(std::nullopt)
, fFieldType(fieldType)
, fFieldMode(EFieldMode::Intrpl)
@@ -293,6 +333,14 @@ namespace cbm::algo::kf
return fFieldMode == EFieldMode::Intrpl ? foFldIntrpl->Get(x, y, z) : foFldOrig->Get(x, y, z);
}
+ /// \brief Gets the double integrals of the field along the track
+ std::tuple<T, T, T> GetDoubleIntegrals(const T& x1, const T& y1, const T& z1, const T& x2, const T& y2,
+ const T& z2) const
+ {
+ return fFieldMode == EFieldMode::Intrpl ? foFldIntrpl->GetDoubleIntegrals(x1, y1, z1, x2, y2, z2)
+ : foFldOrig->GetDoubleIntegrals(x1, y1, z1, x2, y2, z2);
+ }
+
/// \brief Gets field mode
EFieldMode GetFieldMode() const { return fFieldMode; }
@@ -308,7 +356,16 @@ namespace cbm::algo::kf
/// \param b2 Field value at z2 [kG]
/// \param z2 Third nodal point in z-axis direction [cm]
void Set(const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1, const T& z1, const FieldValue<T>& b2,
- const T& z2);
+ const T& z2)
+ {
+ fFieldMode = EFieldMode::Intrpl;
+ fFieldType = EFieldType::Normal;
+ // Set the field type automatically
+ if (b0.IsZero() && b1.IsZero() && b2.IsZero()) {
+ fFieldType = EFieldType::Null;
+ }
+ foFldIntrpl->Set(b0, z0, b1, z1, b2, z2);
+ }
/// \brief Sets the field type
void SetFieldType(EFieldType fldType) { fFieldType = fldType; }
@@ -325,12 +382,11 @@ namespace cbm::algo::kf
/// \param verbose Verbosity level
std::string ToString(int indentLevel = 0, int verbose = 1) const;
+ const auto& GetIntrplField() const { return foFldIntrpl; }
+
private:
friend class boost::serialization::access;
- /// \brief Default constructor
- FieldRegion() = default;
-
/// \brief Serialization load method
template<class Archive>
void load(Archive& ar, const unsigned int /*version*/)
diff --git a/reco/KF/CbmKfTrackFitter.cxx b/reco/KF/CbmKfTrackFitter.cxx
index 1a5afc8f9385d1bd443c71823065f456a4ea0724..7b3f9652b3d6755a01783dfd62f52609d516a73c 100644
--- a/reco/KF/CbmKfTrackFitter.cxx
+++ b/reco/KF/CbmKfTrackFitter.cxx
@@ -39,10 +39,7 @@ using namespace std;
using ca::fmask;
using ca::fvec;
-namespace
-{
- using namespace cbm::algo;
-}
+using namespace cbm::algo;
void CbmKfTrackFitter::Track::OrderNodesInZ()
{
@@ -529,8 +526,7 @@ bool CbmKfTrackFitter::FitTrack(CbmKfTrackFitter::Track& t)
fFit.SetParticleMass(fMass);
- ca::FieldRegion<double> field;
- field.SetUseOriginalField();
+ ca::FieldRegion<double> field(kf::GlobalField::fgOriginalFieldType, kf::GlobalField::fgOriginalField);
{ // fit downstream
diff --git a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
index 0b4613a5cbcfae9fd844f8e990ce3f7ce88cd58a..b0464a6c9cc994ae9a31c700009dbc5fcac7f763 100644
--- a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
+++ b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
@@ -55,38 +55,40 @@ namespace NS_L1TrackFitter
using namespace NS_L1TrackFitter;
-inline void CbmL1PFFitter::PFFieldRegion::setFromL1FieldRegion(const ca::FieldRegion<fvec>& fld, int i)
+inline void CbmL1PFFitter::PFFieldRegion::setFromL1FieldRegion(const ca::FieldRegion<fvec>& fld, int ind)
{
- fP[0] = fld.cx0[i];
- fP[1] = fld.cx1[i];
- fP[2] = fld.cx2[i];
-
- fP[3] = fld.cy0[i];
- fP[4] = fld.cy1[i];
- fP[5] = fld.cy2[i];
-
- fP[6] = fld.cz0[i];
- fP[7] = fld.cz1[i];
- fP[8] = fld.cz2[i];
-
- fP[9] = fld.z0[i];
+ assert(fld.GetFieldMode() == kf::EFieldMode::Intrpl);
+ const auto& coeff = fld.GetIntrplField()->GetCoefficients();
+
+ int i = 0;
+ for (int j = 0; j < 3; j++) {
+ for (int k = 0; k < 3; k++) {
+ fP[i] = coeff[j][k][ind];
+ i++;
+ }
+ }
+ fP[9] = fld.GetIntrplField()->GetZfirst()[ind];
}
-inline void CbmL1PFFitter::PFFieldRegion::getL1FieldRegion(ca::FieldRegion<fvec>& fld, int i)
+inline void CbmL1PFFitter::PFFieldRegion::getL1FieldRegion(ca::FieldRegion<fvec>& fld, int ind)
{
- fld.cx0[i] = fP[0];
- fld.cx1[i] = fP[1];
- fld.cx2[i] = fP[2];
+ assert(fld.GetFieldMode() == kf::EFieldMode::Intrpl);
- fld.cy0[i] = fP[3];
- fld.cy1[i] = fP[4];
- fld.cy2[i] = fP[5];
+ auto coeff = fld.GetIntrplField()->GetCoefficients();
+ auto z = fld.GetIntrplField()->GetZfirst();
- fld.cz0[i] = fP[6];
- fld.cz1[i] = fP[7];
- fld.cz2[i] = fP[8];
+ int i = 0;
+ for (int j = 0; j < 3; j++) {
+ for (int k = 0; k < 3; k++) {
+ coeff[j][k][ind] = fP[i];
+ i++;
+ }
+ }
+ z[ind] = fP[9];
- fld.z0[i] = fP[9];
+ kf::detail::FieldRegionBase<fvec, kf::EFieldMode::Intrpl> fldInterpolated;
+ fldInterpolated.Set(coeff, fP[9]);
+ fld = kf::FieldRegion<fvec>(kf::EFieldType::Normal, fldInterpolated);
}
inline CbmL1PFFitter::PFFieldRegion::PFFieldRegion(const ca::FieldRegion<fvec>& fld, int i)
@@ -170,7 +172,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
Initialize();
ca::FieldValue<fvec> b0, b1, b2 _fvecalignment;
- ca::FieldRegion<fvec> fld _fvecalignment;
+ ca::FieldRegion<fvec> fld;
// fld.SetUseOriginalField();
static int nHits = CbmL1::Instance()->fpAlgo->GetParameters().GetNstationsActive();
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 92ac7e05a2b809388037707af92b4a064cd3ec10..9cdf511f0b72ccf25f9f9bba71642e9b2691ba5a 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -59,6 +59,7 @@
#include <vector>
using namespace cbm::algo::ca;
+using namespace cbm::algo;
using cbm::algo::kf::TrackParamV;
using cbm::ca::tools::Debugger;
@@ -1310,8 +1311,7 @@ void CbmL1::TrackFitPerformance()
if (!mc.IsPrimary()) { // secondary
{ // extrapolate to vertex
- ca::FieldRegion<fvec> fld _fvecalignment;
- fld.SetUseOriginalField();
+ ca::FieldRegion<fvec> fld(kf::GlobalField::fgOriginalFieldType, kf::GlobalField::fgOriginalField);
fit.Extrapolate(mc.GetStartZ(), fld);
// add material
const int fSta = fvHitDebugInfo[it->Hits[0]].iStation;
@@ -1378,8 +1378,7 @@ void CbmL1::TrackFitPerformance()
else { // primary
{ // extrapolate to vertex
- ca::FieldRegion<fvec> fld _fvecalignment;
- fld.SetUseOriginalField();
+ ca::FieldRegion<fvec> fld(kf::GlobalField::fgOriginalFieldType, kf::GlobalField::fgOriginalField);
// add material
const int fSta = fvHitDebugInfo[it->Hits[0]].iStation;
@@ -1527,8 +1526,7 @@ void CbmL1::FillFitHistos(TrackParamV& track, const cbm::ca::tools::MCPoint& mcP
//fit.SetMaxExtrapolationStep(10.);
fit.SetDoFitVelocity(true);
fit.SetTrack(track);
- ca::FieldRegion<fvec> fld _fvecalignment;
- fld.SetUseOriginalField();
+ ca::FieldRegion<fvec> fld(kf::EFieldType::Normal, kf::GlobalField::fgOriginalField);
fit.Extrapolate(mcP.GetZOut(), fld);
track = fit.Tr();
diff --git a/reco/L1/catools/CaToolsField.h b/reco/L1/catools/CaToolsField.h
index b77604318cf976018dd38925896569fe4bb2e406..9757b5034558fecfcce0b0944b8bb2ee630e8e5f 100644
--- a/reco/L1/catools/CaToolsField.h
+++ b/reco/L1/catools/CaToolsField.h
@@ -37,20 +37,10 @@ namespace cbm::ca::tools
return std::tuple(0., 0., 0.);
}
};
- if (FairRunAna::Instance()->GetField()) {
- // NOTE: SZh 28.08.2024:
- // A template class instantiates an independent static variable for each class specification, so each specification must
- // be addressed.
- ca::FieldRegion<ca::fvec>::SetOriginalField(fld, false);
- ca::FieldRegion<float>::SetOriginalField(fld, false);
- ca::FieldRegion<double>::SetOriginalField(fld, false);
- }
- else {
- ca::FieldRegion<ca::fvec>::SetOriginalField(fld, true);
- ca::FieldRegion<float>::SetOriginalField(fld, true);
- ca::FieldRegion<double>::SetOriginalField(fld, true);
- }
- cbm::algo::kf::GlobalField::SetFieldFunction(fld);
+
+ cbm::algo::kf::GlobalField::SetFieldFunction( //
+ (FairRunAna::Instance()->GetField() ? cbm::algo::kf::EFieldType::Normal : cbm::algo::kf::EFieldType::Null), //
+ fld);
}
} // namespace cbm::ca::tools
diff --git a/reco/L1/qa/CbmCaTrackFitQa.cxx b/reco/L1/qa/CbmCaTrackFitQa.cxx
index 357bae8080ea908c543c9b2b2026455871cc09c2..0de6cb97e8d09ad488b3b0a643750e517b85500e 100644
--- a/reco/L1/qa/CbmCaTrackFitQa.cxx
+++ b/reco/L1/qa/CbmCaTrackFitQa.cxx
@@ -165,8 +165,8 @@ void TrackFitQa::Fill(const TrackParamV& trPar, const tools::MCPoint& mcPoint, b
fitter.SetMask(fmask::One());
fitter.SetDoFitVelocity(true);
fitter.SetTrack(trPar);
- cbm::algo::ca::FieldRegion<ca::fvec> fieldRegion;
- fieldRegion.SetUseOriginalField(); // Precised magnetic field is used
+ cbm::algo::ca::FieldRegion<ca::fvec> fieldRegion(kf::GlobalField::fgOriginalFieldType,
+ kf::GlobalField::fgOriginalField);
fitter.Extrapolate(mcPoint.GetZ(), fieldRegion);
const TrackParamV& trParExtr = fitter.Tr(); // Track parameters extrapolated to given MC point
diff --git a/reco/L1/qa/CbmCaTrackTypeQa.cxx b/reco/L1/qa/CbmCaTrackTypeQa.cxx
index 8d152adb885c4d5019e877dac9a0650d8a798c5e..bb07084d91763bced95b79da1f0ac32f97d4b4f9 100644
--- a/reco/L1/qa/CbmCaTrackTypeQa.cxx
+++ b/reco/L1/qa/CbmCaTrackTypeQa.cxx
@@ -19,6 +19,7 @@ using cbm::ca::TrackTypeQa;
using cbm::ca::tools::MCPoint;
using cbm::ca::tools::MCTrack;
+using namespace cbm::algo;
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -240,7 +241,7 @@ void TrackTypeQa::Init()
// Track fitter
fTrackFit.SetDoFitVelocity(true); // TODO: set flag to the configuration
fTrackFit.SetMask(fmask::One());
- fFieldRegion.SetUseOriginalField();
+ fFieldRegion = kf::FieldRegion<fvec>(kf::GlobalField::fgOriginalFieldType, kf::GlobalField::fgOriginalField);
}
// ---------------------------------------------------------------------------------------------------------------------