diff --git a/algo/ca/core/pars/CaField.cxx b/algo/ca/core/pars/CaField.cxx
index 426555698ddb24bf41da76b468a50b8b47dc5316..c5bddcc6c1542aec6c0df77939ac54574d8f6f2f 100644
--- a/algo/ca/core/pars/CaField.cxx
+++ b/algo/ca/core/pars/CaField.cxx
@@ -4,7 +4,7 @@
 
 #include "CaField.h"
 
-#include "KfSimdUtils.h"
+#include "KfUtils.h"
 
 #include <iomanip>
 #include <iostream>
@@ -217,6 +217,7 @@ FieldValue<DataT> FieldRegion<DataT>::Get(const DataT& x, const DataT& y, 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] =
         fgOdiginalField(kfutils::simd::Cast<DataT, double>(x, i), kfutils::simd::Cast<DataT, double>(y, i),
@@ -274,23 +275,23 @@ void FieldRegion<DataT>::Set(const FieldValue<DataT>& b0, const DataT b0z, const
 
 //----------------------------------------------------------------------------------------------------------------------
 // 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;
-}
+//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)
diff --git a/algo/ca/core/pars/CaField.h b/algo/ca/core/pars/CaField.h
index e511b6d2635c4b6069da5ed49c76d3216006c1fe..dedb95fe0a6d3fe064bf2da67ac3df14308a6082 100644
--- a/algo/ca/core/pars/CaField.h
+++ b/algo/ca/core/pars/CaField.h
@@ -241,7 +241,7 @@ namespace cbm::algo::ca
     /// \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);
+    //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
@@ -262,10 +262,10 @@ namespace cbm::algo::ca
     /// \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);
-    }
+    //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
diff --git a/algo/ca/core/tracking/CaTrackFit.cxx b/algo/ca/core/tracking/CaTrackFit.cxx
index af0745e0f3bb7501822ee6d0aeef0d70189b532c..50007a34f77685e58522235abce58a5469ab53fa 100644
--- a/algo/ca/core/tracking/CaTrackFit.cxx
+++ b/algo/ca/core/tracking/CaTrackFit.cxx
@@ -693,6 +693,7 @@ namespace cbm::algo::ca
       fvec Bx          = B.x;
       fvec By          = B.y;
       fvec Bz          = B.z;
+      // NOTE: SZh 13.08.2024: The rstep[0] and rstep[1] make no effect on the B value, if Field is an approximation
 
       fvec tx    = rstep[2];
       fvec ty    = rstep[3];
diff --git a/algo/ca/core/utils/CaSimd.h b/algo/ca/core/utils/CaSimd.h
index 2cca35b687ba36d03aee4c5390588e7fb45c6d3b..60a77a3231117685312577eac4129cd7cb36625b 100644
--- a/algo/ca/core/utils/CaSimd.h
+++ b/algo/ca/core/utils/CaSimd.h
@@ -5,7 +5,7 @@
 #pragma once  // include this header only once per compilation unit
 
 #include "KfSimd.h"
-#include "KfSimdUtils.h"
+#include "KfUtils.h"
 
 namespace cbm::algo::ca
 {
diff --git a/algo/kf/core/CMakeLists.txt b/algo/kf/core/CMakeLists.txt
index e34fb10bac91bf9a45b5d149a070a88a2d513514..55646bf8c5e58f8105264a5aa4678c65fccac269 100644
--- a/algo/kf/core/CMakeLists.txt
+++ b/algo/kf/core/CMakeLists.txt
@@ -1,17 +1,25 @@
 set(INCLUDE_DIRECTORIES
   ${CMAKE_CURRENT_SOURCE_DIR}
+  ${CMAKE_CURRENT_SOURCE_DIR}/algo
   ${CMAKE_CURRENT_SOURCE_DIR}/data
-  ${CMAKE_CURRENT_SOURCE_DIR}/pars 
+  ${CMAKE_CURRENT_SOURCE_DIR}/geo
+  ${CMAKE_CURRENT_SOURCE_DIR}/pars
   ${CMAKE_CURRENT_SOURCE_DIR}/utils
 )
 
 set(SRCS 
   ${CMAKE_CURRENT_SOURCE_DIR}/KfFramework.cxx
   ${CMAKE_CURRENT_SOURCE_DIR}/data/KfTrackParam.cxx
-  ${CMAKE_CURRENT_SOURCE_DIR}/pars/KfMaterialMap.cxx
+  ${CMAKE_CURRENT_SOURCE_DIR}/geo/KfMaterialMap.cxx
+  ${CMAKE_CURRENT_SOURCE_DIR}/geo/KfGeoLayer.cxx
+  ${CMAKE_CURRENT_SOURCE_DIR}/geo/KfTarget.cxx
+  ${CMAKE_CURRENT_SOURCE_DIR}/geo/KfField.cxx
+  ${CMAKE_CURRENT_SOURCE_DIR}/geo/KfFieldValue.cxx
+  ${CMAKE_CURRENT_SOURCE_DIR}/geo/KfFieldSlice.cxx
+  ${CMAKE_CURRENT_SOURCE_DIR}/geo/KfFieldRegion.cxx
+  ${CMAKE_CURRENT_SOURCE_DIR}/geo/KfSetup.cxx
   ${CMAKE_CURRENT_SOURCE_DIR}/pars/KfParameters.cxx
-  ${CMAKE_CURRENT_SOURCE_DIR}/pars/KfSetup.cxx
-  ${CMAKE_CURRENT_SOURCE_DIR}/utils/KfSimdUtils.cxx
+  ${CMAKE_CURRENT_SOURCE_DIR}/utils/KfUtils.cxx
 )
 
 SET_SOURCE_FILES_PROPERTIES(${SRCS} PROPERTIES COMPILE_FLAGS "-O3")
@@ -28,7 +36,9 @@ add_library(KfCore SHARED ${SRCS})
 list(APPEND HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/utils/KfSimd.h)
 
 target_include_directories(KfCore
-  PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/data
+  PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/algo
+         ${CMAKE_CURRENT_SOURCE_DIR}/data
+         ${CMAKE_CURRENT_SOURCE_DIR}/geo
          ${CMAKE_CURRENT_SOURCE_DIR}/pars
          ${CMAKE_CURRENT_SOURCE_DIR}/utils
          ${CMAKE_CURRENT_SOURCE_DIR}
@@ -49,6 +59,8 @@ install(TARGETS KfCore DESTINATION lib)
 install(DIRECTORY kf TYPE INCLUDE FILES_MATCHING PATTERN "*.h")
 install(DIRECTORY kf/utils TYPE INCLUDE FILES_MATCHING PATTERN "*.h")
 install(DIRECTORY kf/data TYPE INCLUDE FILES_MATCHING PATTERN "*.h")
+install(DIRECTORY kf/geo TYPE INCLUDE FILES_MATCHING PATTERN "*.h")
+install(DIRECTORY kf/algo TYPE INCLUDE FILES_MATCHING PATTERN "*.h")
 install(DIRECTORY kf/pars TYPE INCLUDE FILES_MATCHING PATTERN "*.h")
 
 install(
@@ -56,12 +68,16 @@ install(
     KfFramework.h
     KfDefs.h
     data/KfTrackParam.h
-    pars/KfMaterialMap.h 
+    geo/KfMaterialMap.h 
+    geo/KfGeoLayer.h
+    geo/KfSetup.h
+    geo/KfTarget.h
     pars/KfParameters.h 
-    pars/KfSetup.h
     utils/KfVector.h
     utils/KfSimd.h
-    utils/KfSimdUtils.h
+    utils/KfSimdVc.h 
+    utils/KfSimdPseudo.h
+    utils/KfUtils.h
   DESTINATION
     include/
 )
diff --git a/algo/kf/core/KfDefs.h b/algo/kf/core/KfDefs.h
index 2282b08afe145e763fff37015b5480606e7aa54a..96588834c608aef76760318b618fef4def0f610f 100644
--- a/algo/kf/core/KfDefs.h
+++ b/algo/kf/core/KfDefs.h
@@ -7,14 +7,76 @@
 /// @since  28.03.2024
 /// @author Sergei Zharko <s.zharko@gsi.de>
 
-#pragma once
+#ifndef KF_CORE_KfDefs_h
+#define KF_CORE_KfDefs_h 1
 
 #include "KfSimd.h"
 
+#include <array>
+#include <functional>
 #include <limits>
+#include <tuple>
+
+namespace cbm::algo::kf
+{
+  // ---- Enumerations -------------------------------------------------------------------------------------------------
+  //
+  /// \enum  EFieldMode
+  /// \brief Enumiration for the magnetic field representation variants in the track fitting algorithm
+  enum class EFieldMode
+  {
+    Intrpl,  ///< Interpolated magnetic field
+    Orig     ///< Original magnetic field function
+  };
+
+  /// \enum  EFieldType
+  /// \brief Magnetic field type in different setup regions
+  // TODO: SZh 21.08.2024: Develope the concept of field types at different setup regions
+  enum class EFieldType
+  {
+    Target,   ///< Field near the target
+    Tracker,  ///< Field near the tracker subsystem
+    Low,      ///< Magnetic field in the intermediate region (from tracker to null)
+    Null      ///< No magnetic field
+  };
+
+  /// \struct Literal
+  /// \brief  Replaces the type T with the Literal::type to handle the constant expressions for different constants
+  template<typename T>
+  struct Literal {
+    using type = T;
+  };
+
+  template<>
+  struct Literal<fvec> {
+    using type = fscal;
+  };
+
+  template<typename T>
+  using Literal_t = typename Literal<T>::type;
+
+
+  // ---- Common data types --------------------------------------------------------------------------------------------
+  //
+  /// \brief Geometry (spatial) vector
+  template<class T>
+  using GeoVector_t = std::array<T, 3>;
+
+  /// \brief Magnetic field function type
+  /// Signature: tuple<Bx, By, Bz>(x, y, z);
+  using FieldFn_t = std::function<std::tuple<double, double, double>(double, double, double)>;
+
+  /// \brief Magnetic field function type
+  /// \note  The field function follows the FairField::GetFieldFunction signature
+  using FieldFnFair_t = std::function<void(const double (&r)[3], double*)>;
+}  // namespace cbm::algo::kf
 
 namespace cbm::algo::kf::defs
 {
+  // ----- Array sizes -------------------------------------------------------------------------------------------------
+  constexpr int MaxNofFieldSlices = 64;  ///< Max number of field slices
+
+
   // ----- Control -----------------------------------------------------------------------------------------------------
   constexpr int DebugLvl     = 0;     ///< Level of debug output
   constexpr bool GetterCheck = true;  ///< Bound check in getters
@@ -22,38 +84,32 @@ namespace cbm::algo::kf::defs
 
   // ----- Math --------------------------------------------------------------------------------------------------------
   template<class T = double>
-  constexpr T Pi = T(3.14159265358979323846L);  ///< Value of PI, used in ROOT TMath
-
+  constexpr auto Pi = Literal_t<T>(3.14159265358979323846L);  ///< Value of PI, used in ROOT TMath
 
   // ----- Physics -----------------------------------------------------------------------------------------------------
   template<class T = double>
-  constexpr T MuonMass = T(0.105658375523L);  ///< Muon mass [GeV/c2]
+  constexpr auto MuonMass = Literal_t<T>(0.105658375523L);  ///< Muon mass [GeV/c2]
 
   template<class T = double>
-  constexpr T PionMass = T(0.1395703918L);  ///< Pion mass [GeV/c2]
+  constexpr auto PionMass = Literal_t<T>(0.1395703918L);  ///< Pion mass [GeV/c2]
 
   template<class T = double>
-  constexpr T KaonMass = T(0.493677L);  ///< Kaon mass [GeV/c2] (PDG 22.08.2023)
+  constexpr auto KaonMass = Literal_t<T>(0.493677L);  ///< Kaon mass [GeV/c2] (PDG 22.08.2023)
 
   template<class T = double>
-  constexpr T ElectronMass = T(0.0005109989500015L);  ///< Electron mass [GeV/c2]
+  constexpr auto ElectronMass = Literal_t<T>(0.0005109989500015L);  ///< Electron mass [GeV/c2]
 
   template<class T = double>
-  constexpr T ProtonMass = T(0.93827208816L);  ///< Proton mass [GeV/c2] (PDG 11.08.2022)
+  constexpr auto ProtonMass = Literal_t<T>(0.93827208816L);  ///< Proton mass [GeV/c2] (PDG 11.08.2022)
 
   template<class T = double>
-  constexpr T SpeedOfLight = T(29.9792458L);  ///< Speed of light [cm/ns]
+  constexpr auto SpeedOfLight = Literal_t<T>(29.9792458L);  ///< Speed of light [cm/ns]
 
   template<class T = double>
-  constexpr T SpeedOfLightInv = T(1.L) / SpeedOfLight<T>;  ///< Inverse speed of light [ns/cm]
+  constexpr auto SpeedOfLightInv = Literal_t<T>(1.L) / SpeedOfLight<T>;  ///< Inverse speed of light [ns/cm]
 
   template<class T = double>
-  constexpr T MinFeild = T(1.e-4L);  ///< Minimal (negligible) magnetic field value [kG]
-
-
-  // ----- Misc --------------------------------------------------------------------------------------------------------
-  constexpr int MemAlign = 16;  ///< Default alignment of memory (bytes)
-
+  constexpr auto MinField = Literal_t<T>(1.e-4L);  ///< Minimal (negligible) magnetic field value [kG]
 
   // ----- Undefined values --------------------------------------------------------------------------------------------
   /// \brief Undefined values
@@ -75,4 +131,17 @@ namespace cbm::algo::kf::defs
   template<>
   inline constexpr fscal Undef<fvec> = std::numeric_limits<fscal>::signaling_NaN();
 
+  // ----- Other -------------------------------------------------------------------------------------------------------
+  /// \brief Zero magnetic field function
+  constexpr auto ZeroFieldFn = [](double, double, double) constexpr { return std::make_tuple(0., 0., 0.); };
+
+  // ----- Debug flags (NOTE: to be modified only by the experts) ------------------------------------------------------
+  //
+  namespace dbg
+  {
+    /// \brief Replaces the magnetic field extrapolation with the true field function in the ALL FieldRegion objects
+    constexpr bool ForceOriginalField = false;
+  }  // namespace dbg
 }  // namespace cbm::algo::kf::defs
+
+#endif  // KF_CORE_KfDefs_h
diff --git a/algo/kf/core/data/KfTrackParam.h b/algo/kf/core/data/KfTrackParam.h
index 949a49da9d5b30f496abecd2a8b43076f6ea8f00..fb5fd53a77f115a403f84bed6e28eafa26e7f3db 100644
--- a/algo/kf/core/data/KfTrackParam.h
+++ b/algo/kf/core/data/KfTrackParam.h
@@ -13,7 +13,7 @@
 
 #include "KfDefs.h"
 #include "KfSimd.h"
-#include "KfSimdUtils.h"
+#include "KfUtils.h"
 
 #include <boost/serialization/access.hpp>
 
@@ -760,7 +760,6 @@ namespace cbm::algo::kf
     SetC66(dvi * dvi);
   }
 
-}  // namespace cbm::algo::ca
+}  // namespace cbm::algo::kf
 
 #endif  // CA_CORE_CaTrackParam_h
-
diff --git a/algo/kf/core/geo/KfField.cxx b/algo/kf/core/geo/KfField.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c9ab45a6b6b70f69f18d845a2b6fcdb2c917c1ec
--- /dev/null
+++ b/algo/kf/core/geo/KfField.cxx
@@ -0,0 +1,73 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   KfField.cxx
+/// \brief  Magnetic field representation in KF (source)
+/// \author Sergei Zharko <s.zharko@gsi.de>
+/// \since  22.07.2024
+
+#include "KfField.h"
+
+#include <mutex>
+#include <sstream>
+#include <tuple>
+
+using cbm::algo::kf::EFieldMode;
+using cbm::algo::kf::FieldFactory;
+using cbm::algo::kf::FieldFnFair_t;
+using cbm::algo::kf::detail::FieldBase;
+
+template<typename T>
+using FieldBaseIntrpl_t = FieldBase<T, EFieldMode::Intrpl>;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+template<typename I>
+FieldBaseIntrpl_t<T>::FieldBase(const FieldBaseIntrpl_t<I>& other)
+{
+  for (size_t i = 0; i < defs::MaxNofFieldSlices; ++i) {
+    fvFieldSlices[i] = utils::simd::Cast<I, T>(other.fvFieldSlices[i]);
+  }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+template<typename I>
+FieldBaseIntrpl_t<T>& FieldBaseIntrpl_t<T>::operator=(const FieldBaseIntrpl_t<I>& other)
+{
+  if (this != &other) {
+    for (size_t i = 0; i < defs::MaxNofFieldSlices; ++i) {
+      fvFieldSlices[i] = utils::simd::Cast<I, T>(other.fvFieldSlices[i]);
+    }
+  }
+  return *this;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void FieldFactory::AddSliceReference(double halfSizeX, double halfSizeY, double zRef)
+{
+  if (!fSliceReferences.emplace().second) {
+    std::stringstream msg;
+    msg << "FieldContainerFactory::AddReference: attempt of adding another slice reference with zRef = " << zRef
+        << "(halfSizeX = " << halfSizeX << ", halfSizeY = " << halfSizeY << ").\nThe next slice references were "
+        << "added:";
+    for (const auto& el : fSliceReferences) {
+      msg << "\n\t- halfSizeX = " << el.fHalfSizeX << ", halfSizeY = " << el.fHalfSizeY << ", zRef = " << el.fRefZ;
+    }
+    throw std::runtime_error(msg.str());
+  }
+}
+
+namespace cbm::algo::kf
+{
+  template class Field<float, EFieldMode::Orig>;
+  template class Field<float, EFieldMode::Intrpl>;
+  template class Field<double, EFieldMode::Orig>;
+  template class Field<double, EFieldMode::Intrpl>;
+  template class Field<fvec, EFieldMode::Orig>;
+  template class Field<fvec, EFieldMode::Intrpl>;
+}  // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfField.h b/algo/kf/core/geo/KfField.h
new file mode 100644
index 0000000000000000000000000000000000000000..5a18714ebaec9ce08cb955f87ad486db63289270
--- /dev/null
+++ b/algo/kf/core/geo/KfField.h
@@ -0,0 +1,366 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   KfField.h
+/// \brief  Magnetic field representation in KF (header)
+/// \author Sergei Zharko <s.zharko@gsi.de>
+/// \since  22.07.2024
+
+#pragma once
+
+#include "KfDefs.h"
+#include "KfFieldRegion.h"
+#include "KfFieldSlice.h"
+
+#include <array>
+#include <set>
+
+namespace cbm::algo::kf
+{
+  namespace detail
+  {
+    /// \class  FieldBase
+    /// \brief  Data members of the Field class (primary template)
+    /// \tparam T        Underlying floating point type (float/double/fvec)
+    /// \tparam FldMode  Field type, defined by EFieldMode
+    template<typename T, EFieldMode FldMode>
+    class alignas(VcMemAlign) FieldBase {
+    };
+
+    /// \class  FieldBase<T, Orig>
+    /// \brief  FieldBase implementation for the original field
+    /// \tparam T  Underlying floating point type (float/double/fvec)
+    template<typename T>
+    class alignas(VcMemAlign) FieldBase<T, EFieldMode::Orig> {
+     public:
+      using FieldRegion_t = FieldRegion<T, EFieldMode::Orig>;
+      using FieldValue_t  = FieldValue<T>;
+
+      /// \brief Sets magnetic field function
+      /// \param fieldFn  Magnetic field function (KF-format)
+      void SetFieldFunction(const FieldFn_t& fieldFn) { fFieldFn = fieldFn; }
+
+      /// \brief Creates field region object
+      FieldRegion_t MakeFieldRegion() const { return FieldRegion_t(fFieldFn, fFieldType); }
+
+     protected:
+      /// \brief Default constructor
+      FieldBase() = default;
+
+      /// \brief  Copy constructor
+      /// \tparam I  Underlying floating type of the source
+      template<typename I>
+      FieldBase(const FieldBase<I, EFieldMode::Orig>& other) : fFieldFn(other.fFieldFn)
+      {
+      }
+
+      /// \brief  Copy assignment operator
+      /// \tparam I  Underlying floating type of the source
+      template<typename I>
+      FieldBase& operator=(const FieldBase<I, EFieldMode::Orig>& other)
+      {
+        if (this != &other) {
+          fFieldFn = other.fFieldFn;
+        }
+        return *this;
+      }
+
+      FieldFn_t fFieldFn{defs::ZeroFieldFn};    ///< Field function: (x,y,z) [cm] -> (Bx,By,Bz) [kG]
+      EFieldType fFieldType{EFieldType::Null};  ///< Field type
+
+     private:
+      /// Serialization function
+      friend class boost::serialization::access;
+      template<class Archive>
+      void serialize(Archive& ar, const unsigned int)
+      {
+        ar& fFieldType;
+      }
+    };
+
+
+    /// \class  FieldBase<T, EFieldMode::Intrpl>
+    /// \brief  Data members of the Field class with the interpolation of the magnetic field
+    /// \tparam T  Underlying floating point type (float/double/fvec)
+    template<typename T>
+    class alignas(VcMemAlign) FieldBase<T, EFieldMode::Intrpl> {
+      using SlicesContainer_t = std::array<FieldSlice<T>, defs::MaxNofFieldSlices>;
+
+     public:
+      using FieldRegion_t = FieldRegion<T, EFieldMode::Intrpl>;
+      using FieldValue_t  = FieldValue<T>;
+
+      /// \brief Creates field value object
+      /// \param sliceID  Index of slice
+      /// \param x        x-coordinate of the point [cm]
+      /// \param y        y-coordinate of the point [cm]
+      FieldValue_t GetFieldValue(int sliceID, const T& x, const T& y) const
+      {
+        return fvFieldSlices[sliceID].GetFieldValue(x, y);
+      }
+
+      /// \brief Creates field region object
+      /// \param b0  Field value in the first node [kG]
+      /// \param z0  First node z-coordinate [cm]
+      /// \param b1  Field value in the first node [kG]
+      /// \param z1  Second node z-coordinate [cm]
+      /// \param b2  Field value in the first node [kG]
+      /// \param z2  Third node z-coordinate [cm]
+      FieldRegion_t MakeFieldRegion(const FieldValue_t& b0, const T& z0, const FieldValue_t& b1, const T& z1,
+                                    const FieldValue_t& b2, const T& z2) const
+      {
+        return FieldRegion_t(b0, z0, b1, z1, b2, z2, fFieldType);
+      }
+
+      /// \brief Sets a field slice
+      /// \param id          Index of the layer
+      /// \param fieldSlice  Field slice object
+      void SetFieldSlice(int id, const FieldSlice<T>& fieldSlice) { fvFieldSlices[id] = fieldSlice; }
+
+     protected:
+      /// \brief Default constructor
+      FieldBase() = default;
+
+      /// \brief  Copy constructor
+      /// \tparam I  Underlying floating type of the source
+      template<typename I>
+      FieldBase(const FieldBase<I, EFieldMode::Intrpl>& other);
+
+      /// \brief  Copy assignment operator
+      /// \tparam I  Underlying floating type of the source
+      template<typename I>
+      FieldBase& operator=(const FieldBase<I, EFieldMode::Intrpl>& other);
+
+      SlicesContainer_t fvFieldSlices;          ///< Array of field slices
+      EFieldType fFieldType{EFieldType::Null};  ///< Field type
+
+     private:
+      /// \brief Serialization function
+      friend class boost::serialization::access;
+      template<class Archive>
+      void serialize(Archive& ar, const unsigned int)
+      {
+        ar& fvFieldSlices;
+        ar& fFieldType;
+      }
+    };
+  }  // namespace detail
+
+
+  /// \class  Field
+  /// \brief  Magnetic field manager class
+  /// \tparam T  Underlying floating point type (float/double/fvec)
+  template<typename T, EFieldMode FldMode>
+  class alignas(VcMemAlign) Field : public detail::FieldBase<T, FldMode> {
+    using FieldBase_t = detail::FieldBase<T, FldMode>;
+
+   public:
+    static constexpr EFieldMode FieldType = FldMode;
+
+    using FieldRegion_t = typename FieldBase_t::FieldRegion_t;
+
+    /// \brief Default constructor
+    Field() = default;
+
+    /// \brief Copy constructor
+    template<typename I>
+    Field(const Field<I, FldMode>& other)
+      : FieldBase_t(other)
+      , fPrimVertexField(FieldRegion<I, FldMode>(other.fPrimVertexField))
+    {
+    }
+
+    /// \brief Destructor
+    ~Field() = default;
+
+    /// \brief Copy assignment operator
+    template<typename I>
+    Field& operator=(const Field<I, FldMode>& other)
+    {
+      if (this == &other) {
+        return *this;
+      }
+      fPrimVertexField            = other.fPrimVertexField;
+      return FieldBase_t::operator=(other);
+    }
+
+    /// \brief Gets field type
+    EFieldType GetFieldType() const { return this->fFieldType; }
+
+    /// \brief Sets field type
+    /// \param fieldType  Field type
+    // TODO: SZh 21.08.2024:
+    //   The concept of the field type depending on the z-regions should be developed. At the moment the full setup can
+    //   have either EFieldType::Null (mCBM), or something else (CBM).
+    void SetFieldType(EFieldType fieldType) { this->fFieldType = fieldType; }
+
+    /// \brief Sets field region near primary vertex
+    /// \param fld  Field region near vertex
+    void SetPrimVertexField(const FieldRegion_t& primVertexField) { fPrimVertexField = primVertexField; }
+
+    /// \brief Gets field region near primary vertex
+    const FieldRegion_t& GetPrimVertexField() const { return fPrimVertexField; }
+
+   private:
+    /// \brief Serialization function
+    friend class boost::serialization::access;
+    template<class Archive>
+    void serialize(Archive& ar, const unsigned int)
+    {
+      ar& boost::serialization::base_object<FieldBase_t>(*this);
+      ar& fPrimVertexField;
+    }
+
+    FieldRegion_t fPrimVertexField;  ///< Field region near primary vertex (constant field is assumed)
+  };
+
+
+  /// \class FieldFactory
+  /// \brief A helper class to instantiate a Field object
+  class FieldFactory {
+    /// \struct SliceRef
+    /// \brief  A helper structure for field slices initialization
+    struct SliceRef {
+      double fHalfSizeX{defs::Undef<double>};  ///< Half-size of the slice in x-direction [cm]
+      double fHalfSizeY{defs::Undef<double>};  ///< Half-size of the slice in y-direction [cm]
+      double fRefZ{defs::Undef<double>};       ///< Reference z-position of the slice [cm]
+
+      /// \brief Default constructor
+      SliceRef() = default;
+
+      /// \brief Destructor
+      //~SliceRef() = default;
+
+      /// \brief Constructor
+      /// \param halfX  Half-size of the slice in x-direction [cm]
+      /// \param halfY  Half-size of the slice in y-direction [cm]
+      /// \param refZ   Reference z-position of the slice [cm]
+      SliceRef(double halfX, double halfY, double refZ) : fHalfSizeX(halfX), fHalfSizeY(halfY), fRefZ(refZ) {}
+
+      /// \brief Comparision operator
+      friend constexpr bool operator<(const SliceRef& l, const SliceRef& r) { return l.fRefZ < r.fRefZ; }
+    };
+
+   public:
+    /// \brief Default constructor
+    FieldFactory() = default;
+
+    /// \brief Copy constructor
+    FieldFactory(const FieldFactory&) = default;
+
+    /// \brief Move constructor
+    FieldFactory(FieldFactory&&) = default;
+
+    /// \brief Destructor
+    ~FieldFactory() = default;
+
+    /// \brief Copy assignment operator
+    FieldFactory& operator=(const FieldFactory&) = default;
+
+    /// \brief Move assignment operator
+    FieldFactory& operator=(FieldFactory&&) = default;
+
+    /// \brief Adds a slice reference
+    /// \param halfSizeX  Half-size of the slice in x-direction [cm]
+    /// \param halfSizeY  Half-size of the slice in y-direction [cm]
+    /// \param refZ       Reference z-position of the slice [cm]
+    void AddSliceReference(double halfSizeX, double halfSizeY, double refZ);
+
+    /// \brief Gets field type
+    EFieldType GetFieldType() const { return fFieldType; }
+
+    /// \brief  Create field manager
+    /// \tparam T          Underlying floating-point data type
+    /// \tparam OrigField  If orig field is used
+    template<typename T, EFieldMode FldMode>
+    Field<T, FldMode> MakeField() const;
+
+    /// \brief Resets the instance
+    void Reset() { *this = FieldFactory(); }
+
+    /// \brief Sets magnetic field function
+    /// \param fieldFn  Magnetic field function (KF-format)
+    void SetFieldFunction(const FieldFn_t& fieldFn) { fFieldFn = fieldFn; }
+
+    /// \brief Sets magnetic field function
+    /// \param fieldFn  Magnetic field function (FairRoot format)
+    void SetFieldFunction(const FieldFnFair_t& fieldFn) { fFieldFn = GlobalField::ConvertFairFieldFunction(fieldFn); }
+
+    /// \brief Sets field type
+    /// \param fieldType  Field type
+    void SetFieldType(EFieldType fieldType) { fFieldType = fieldType; }
+
+    /// \brief Sets target
+    /// \param x  x-coordinate of the target position [cm]
+    /// \param y  y-coordinate of the target position [cm]
+    /// \param z  z-coordinate of the target position [cm]
+    void SetTarget(double x, double y, double z) { fTarget = {x, y, z}; }
+
+    /// \brief Sets a step for the primary vertex field region estimation
+    /// \param step  A step between nodal points in z-axis direction [cm]
+    void SetStep(double step = 2.5) { fTargetStep = step; }
+
+   private:
+    std::set<SliceRef> fSliceReferences;      ///< Set of slice references
+    FieldFn_t fFieldFn{defs::ZeroFieldFn};    ///< Field function (x, y, z) [cm] -> (Bx, By, Bz) [kG]
+    double fTargetStep{2.5};                  ///< Step between nodal points for the primary vertex field estimation
+    EFieldType fFieldType{EFieldType::Null};  ///< Field type
+
+    /// \brief Target position
+    std::array<double, 3> fTarget{{defs::Undef<double>, defs::Undef<double>, defs::Undef<double>}};
+  };
+
+  // -------------------------------------------------------------------------------------------------------------------
+  //
+  template<typename T, EFieldMode FldMode>
+  Field<T, FldMode> FieldFactory::MakeField() const
+  {
+    auto field = Field<T, FldMode>{};
+
+    // Check initialization
+    if (std::any_of(fTarget.begin(), fTarget.end(), [](double x) { return !utils::IsFinite(x); })) {
+      throw std::runtime_error("FieldFactory::MakeField: target is undefined");
+    }
+    if (!fSliceReferences.size()) {  // TODO: Remove requirement of slice references
+      throw std::runtime_error("FieldFactory::MakeField: no slice references were provided");
+    }
+    else if (fSliceReferences.size() > defs::MaxNofFieldSlices) {
+      std::stringstream msg;
+      msg << "FieldFactory::MakeField: to many slice references are provided (" << fSliceReferences.size() << "), "
+          << "the maximum allowed number of field slices is " << defs::MaxNofFieldSlices;
+      throw std::runtime_error(msg.str());
+    }
+    if (!fFieldFn) {
+      throw std::runtime_error("FieldFactory::CreateField: no field function is provided");
+    }
+
+    // Initialize the Field object
+    field.SetFieldType(fFieldType);
+    if constexpr (FldMode == EFieldMode::Orig) {
+      field.SetFieldFunction(fFieldFn);
+      field.SetPrimVertexField(field.MakeFieldRegion());
+    }
+    else if constexpr (FldMode == EFieldMode::Intrpl) {  // A version with the interpolated field
+      int layerId = 0;
+      for (const auto& sliceRef : fSliceReferences) {
+        field.SetFieldSlice(layerId++,
+                            FieldSlice<T>(fFieldFn, sliceRef.fHalfSizeX, sliceRef.fHalfSizeY, sliceRef.fRefZ));
+      }
+      {
+        // PV field initialization
+        constexpr size_t nNodes = 3;
+        double z                = fTarget[2];
+        std::array<FieldValue<T>, nNodes> fld;
+        std::array<T, nNodes> pos;
+        for (int iNode = 0; iNode < nNodes; ++iNode) {
+          pos[iNode] = utils::simd::Cast<double, T>(z);
+          fld[iNode] = GlobalField::GetFieldValue<double>(fFieldFn, fTarget[0], fTarget[1], z);
+          z += fTargetStep;
+        }
+        field.SetPrimVertexField(field.MakeFieldRegion(fld[0], pos[0], fld[1], pos[1], fld[2], pos[2]));
+      }
+    }
+    return field;
+  }
+}  // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfFieldRegion.cxx b/algo/kf/core/geo/KfFieldRegion.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..54d3cc53e6ed1696662c1967aea8d43b36a707db
--- /dev/null
+++ b/algo/kf/core/geo/KfFieldRegion.cxx
@@ -0,0 +1,157 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   KfFieldRegion.cxx
+/// \brief  Magnetic flux density interpolation along the track vs. z-coordinate (source)
+/// \author Sergei Zharko <s.zharko@gsi.de>
+/// \since  22.07.2024
+///
+
+#include "KfFieldRegion.h"
+
+#include "KfUtils.h"
+
+#include <mutex>
+#include <sstream>
+
+using cbm::algo::kf::EFieldMode;
+using cbm::algo::kf::FieldFn_t;
+using cbm::algo::kf::FieldFnFair_t;
+using cbm::algo::kf::FieldRegion;
+using cbm::algo::kf::FieldValue;
+using cbm::algo::kf::GlobalField;
+
+template<typename T>
+using FieldRegionIntrpl_t = cbm::algo::kf::detail::FieldRegionBase<T, EFieldMode::Intrpl>;
+
+
+FieldFn_t GlobalField::gOriginalField = &GlobalField::UndefField;
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+void FieldRegionIntrpl_t<T>::Set(const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1, const T& z1,
+                                 const FieldValue<T>& b2, const T& z2)
+{
+  fZfirst = z0;
+
+  auto dz1 = z1 - z0;
+  auto dz2 = z2 - z0;
+  auto det = utils::simd::One<T>() / (dz1 * dz2 * (z2 - z1));
+  auto w21 = -dz2 * det;
+  auto w22 = dz1 * det;
+  auto w11 = -dz2 * w21;
+  auto w12 = -dz1 * w22;
+
+  for (int iD = 0; iD < 3; ++iD) {
+    auto db1    = b1.GetComponent(iD) - b0.GetComponent(iD);
+    auto db2    = b2.GetComponent(iD) - b0.GetComponent(iD);
+    auto& coeff = fCoeff[iD];
+    coeff[0]    = b0.GetComponent(iD);
+    coeff[1]    = db1 * w11 + db2 * w12;
+    coeff[2]    = db1 * w21 + db2 * w22;
+  }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+void FieldRegionIntrpl_t<T>::Shift(const T& z)
+{
+  auto dz = z - fZfirst;
+  for (int iD = 0; iD < 3; ++iD) {
+    auto& coeff = fCoeff[iD];
+    auto c2dz   = coeff[2] * dz;
+    coeff[0] += (coeff[1] + c2dz) * dz;
+    coeff[1] += (2 * c2dz);
+  }
+  fZfirst = z;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T, EFieldMode FldMode>
+FieldValue<T> FieldRegion<T, FldMode>::Get(const T& x, const T& y, const T& z) const
+{
+  if constexpr (FldMode == EFieldMode::Intrpl) {
+    if constexpr (defs::dbg::ForceOriginalField) {
+      return GlobalField::GetFieldValue(GlobalField::gOriginalField, 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]),
+                         this->fCoeff[1][0] + dz * (this->fCoeff[1][1] + dz * this->fCoeff[1][2]),
+                         this->fCoeff[2][0] + dz * (this->fCoeff[2][1] + dz * this->fCoeff[2][2])};
+  }
+  else if constexpr (FldMode == EFieldMode::Orig) {
+    return GlobalField::GetFieldValue(this->fFieldFn, x, y, z);
+  }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T, EFieldMode FldMode>
+std::string FieldRegion<T, FldMode>::ToString(int indentLevel) const
+{
+  constexpr char IndentChar = '\t';
+  std::stringstream msg;
+  std::string indent(indentLevel, IndentChar);
+  if constexpr (FldMode == EFieldMode::Intrpl) {
+    const auto& coeff = this->fCoeff;
+    msg << indent << "Field region: dz = z - " << this->fZfirst << " [cm]\n";
+    msg << indent << "  Bx(dz) = " << coeff[0][0] << " + " << coeff[0][1] << "dz + " << coeff[0][2] << "dz2\n";
+    msg << indent << "  By(dz) = " << coeff[1][0] << " + " << coeff[1][1] << "dz + " << coeff[1][2] << "dz2\n";
+    msg << indent << "  Bz(dz) = " << coeff[2][0] << " + " << coeff[2][1] << "dz + " << coeff[2][2] << "dz2";
+    if constexpr (defs::dbg::ForceOriginalField) {
+      msg << indent << "\nWARNING: the defs::dbg::ForceOriginalField is enabled, so the magnetic field interpolation "
+          << "is replaced with the global magnetic function for debugging purposes. If you see this message and are "
+          << "not sure, what is going on, please set the defs::dbg::ForceOriginalField to false and re-run the routine";
+    }
+  }
+  else if constexpr (FldMode == EFieldMode::Orig) {
+    msg << indent << "Field region: created from the original field function";
+  }
+  return msg.str();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+std::tuple<double, double, double> GlobalField::UndefField(double, double, double)
+{
+  assert(!defs::dbg::ForceOriginalField
+         && "cbm::algo::kf::GlobalField: The original globa magnetic field is required by "
+            "kf::defs::dbg::ForceOriginalField = true. Please provide it with the "
+            "kf::GlobalField::SetGlobalField() function.");
+  return std::make_tuple(defs::Undef<double>, defs::Undef<double>, defs::Undef<double>);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+FieldFn_t GlobalField::ConvertFairFieldFunction(const FieldFnFair_t& fn)
+{
+  if (!fn) {
+    throw std::runtime_error("GlobalField::ConvertFairFieldFunction(): attempt to set an undefined field function");
+  }
+
+  return [&](double x, double y, double z) -> std::tuple<double, double, double> {
+    double pos[3] = {x, y, z};
+    double fld[3] = {0.};
+    static std::mutex mtx;
+    mtx.lock();
+    fn(pos, fld);
+    mtx.unlock();
+    return std::make_tuple(fld[0], fld[1], fld[2]);
+  };
+}
+
+
+namespace cbm::algo::kf
+{
+  template class FieldRegion<float, EFieldMode::Orig>;
+  template class FieldRegion<float, EFieldMode::Intrpl>;
+  template class FieldRegion<double, EFieldMode::Orig>;
+  template class FieldRegion<double, EFieldMode::Intrpl>;
+  template class FieldRegion<fvec, EFieldMode::Orig>;
+  template class FieldRegion<fvec, EFieldMode::Intrpl>;
+}  // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfFieldRegion.h b/algo/kf/core/geo/KfFieldRegion.h
new file mode 100644
index 0000000000000000000000000000000000000000..6906a1f9bb832ad76200124ca43a0ff4f89fc540
--- /dev/null
+++ b/algo/kf/core/geo/KfFieldRegion.h
@@ -0,0 +1,293 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   KfFieldRegion.h
+/// \brief  Magnetic flux density interpolation along the track vs. z-coordinate (header)
+/// \author Sergei Zharko <s.zharko@gsi.de>
+/// \since  19.07.2024
+///
+
+#pragma once
+
+#include "KfDefs.h"
+#include "KfFieldValue.h"
+#include "KfMath.h"
+#include "KfUtils.h"
+
+#include <boost/serialization/access.hpp>
+#include <boost/serialization/array.hpp>
+#include <boost/serialization/base_object.hpp>
+
+#include <array>
+#include <string>
+#include <type_traits>
+
+namespace cbm::algo::kf
+{
+
+  /// \class GlobalField
+  /// \brief Handler for the global magnetic field related functions
+  class GlobalField {
+   public:
+    /// \brief Global variable to store the fielf funciton (x,y,z)->(Bx,By,Bz)
+    static FieldFn_t gOriginalField;
+
+    /// \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 FieldFnFair_t& fn) { gOriginalField = ConvertFairFieldFunction(fn); }
+
+    /// \brief Converts a FairRoot field function to the KF format
+    /// \param fn  Field function of the FairRoot format
+    static FieldFn_t ConvertFairFieldFunction(const FieldFnFair_t& 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
+    /// \param fieldFn  Field function (x,y,z)->(Bx,By,Bz)
+    /// \param x        The point x-coordinate [cm]
+    /// \param y        The point y-coordinate [cm]
+    /// \param z        The point z-coordinate [cm]
+    /// \return  FieldValue<T> object
+    template<typename T>
+    [[gnu::always_inline]] static FieldValue<T> GetFieldValue(const FieldFn_t& fieldFn, const T& x, const T& y,
+                                                              const T& z);
+  };
+
+
+  namespace detail
+  {
+    /// \class  FieldRegionBase
+    /// \brief  Base class of the FieldRegion class (primary template)
+    /// \tparam T          Underlying floating point type (float/double/fvec)
+    /// \tparam OrigField  When true, the original CBM field function is used. When false, the polynomial
+    ///                    interpolation is used instead.
+    template<typename T, EFieldMode FldMode>
+    class alignas(VcMemAlign) FieldRegionBase {
+    };
+
+
+    /// \class  FieldRegionBase<T, EFieldMode::Orig>
+    /// \brief  Data members and specific methods of the FieldRegion with the original magnetic field
+    /// \tparam T  Underlying floating point type (float/double/fvec)
+    template<typename T>
+    class alignas(VcMemAlign) FieldRegionBase<T, EFieldMode::Orig> {
+     public:
+      using FieldValue_t = FieldValue<T>;
+
+      /// \brief Default constructor
+      FieldRegionBase() = default;
+
+      /// \brief Constructor
+      FieldRegionBase(const FieldFn_t& fieldFn) : fFieldFn(fieldFn) {}
+
+      /// \brief  Copy constructor
+      /// \tparam I  Underlying floating point type of the source
+      template<typename I>
+      FieldRegionBase(const FieldRegionBase<I, EFieldMode::Orig>& other) : fFieldFn(other.fFieldFn)
+      {
+      }
+
+      /// \brief Destructor
+      ~FieldRegionBase() = default;
+
+      /// \brief Copy assignment operator
+      /// \tparam I  Underlying floating point type of the source
+      template<typename I>
+      FieldRegionBase& operator=(const FieldRegionBase<I, EFieldMode::Orig>& other)
+      {
+        if (this != &other) {
+          fFieldFn = other.fFieldFn;
+        }
+        return *this;
+      }
+
+     protected:
+      FieldFn_t fFieldFn{defs::ZeroFieldFn};  ///< Field function: (x,y,z) [cm] -> (Bx,By,Bz) [kG]
+
+     private:
+      /// \brief Serialization function
+      friend class boost::serialization::access;
+      template<class Archive>
+      void serialize(Archive&, const unsigned int)
+      {
+      }
+    };
+
+
+    /// \class  FieldRegionBase<T, Intrpl>
+    /// \brief  Data members and specific mehtods of the FieldRegion class with the interpolated magnetic field
+    /// \tparam T  Underlying floating  point type (float/double/fvec)
+    template<typename T>
+    class alignas(VcMemAlign) FieldRegionBase<T, EFieldMode::Intrpl> {
+      using CoeffArray_t = std::array<std::array<T, 3>, 3>;  ///< Approximation coefficients array
+
+     public:
+      using FieldValue_t = FieldValue<T>;
+
+      /// \brief Default constructor
+      FieldRegionBase() = default;
+
+      /// \brief Constructor from parameters
+      /// \param b0  Field value in the first node [kG]
+      /// \param z0  First node z-coordinate [cm]
+      /// \param b1  Field value in the first node [kG]
+      /// \param z1  Second node z-coordinate [cm]
+      /// \param b2  Field value in the first node [kG]
+      /// \param z2  Third node z-coordinate [cm]
+      FieldRegionBase(const FieldValue_t& b0, const T& z0, const FieldValue_t& b1, const T& z1, const FieldValue_t& b2,
+                      const T& z2)
+      {
+        this->Set(b0, z0, b1, z1, b2, z2);
+      }
+
+      /// \brief Copy constructor
+      template<typename I>
+      FieldRegionBase(const FieldRegionBase<I, EFieldMode::Intrpl>& other)
+        : fCoeff(utils::simd::Cast<I, T>(other.fCoeff))
+        , fZfirst(utils::simd::Cast<I, T>(other.fZfirst))
+      {
+      }
+
+      /// \brief Destructor
+      ~FieldRegionBase() = default;
+
+      /// \brief Copy assignment operator
+      /// \tparam I  Underlying floating point type of the source
+      template<typename I>
+      FieldRegionBase& operator=(const FieldRegionBase<I, EFieldMode::Orig>& other)
+      {
+        if (this != &other) {
+          this->fCoeff  = utils::simd::Cast<I, T>(other.fCoeff);
+          this->fZfirst = utils::simd::Cast<I, T>(other.fZfirst);
+        }
+        return *this;
+      }
+
+      /// \brief Gets the first z-coordinate of the interpolation
+      const T& GetZfirst() const { return fZfirst; }
+
+      /// \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]
+      /// \param z0  First nodal point in z-axis direction [cm]
+      /// \param b1  Field value at z1 [kG]
+      /// \param z1  Second nodal point in z-axis direction [cm]
+      /// \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);
+
+      /// \brief Shifts the coefficients to another first z-coordinate
+      /// \param z  new z-coordinate [cm]
+      void Shift(const T& z);
+
+     protected:
+      CoeffArray_t fCoeff{{T(0.)}};  ///< Approximation coefficients
+      T fZfirst{Literal_t<T>()};     ///< Reference z-coordinate of the field region
+
+     private:
+      /// Serialization function
+      friend class boost::serialization::access;
+      template<class Archive>
+      void serialize(Archive& ar, const unsigned int)
+      {
+        ar& fCoeff;
+        ar& fZfirst;
+      }
+    };
+  }  // namespace detail
+
+
+  /// \class  FieldRegion
+  /// \brief  Magnetic field region, corresponding to a hit triplet
+  /// \tparam T        Underlying data-type
+  /// \tparam FldMode  Field type
+  template<typename T, EFieldMode FldMode>
+  class alignas(VcMemAlign) FieldRegion : public detail::FieldRegionBase<T, FldMode> {
+    using FieldRegionBase_t = detail::FieldRegionBase<T, FldMode>;
+
+   public:
+    static constexpr EFieldMode FieldType = FldMode;
+
+    using FieldValue_t = FieldValue<T>;
+
+    /// \brief  Copy constructor
+    /// \tparam I  Underlying floating point type of the source
+    template<typename I>
+    FieldRegion(const FieldRegion<I, FldMode>& other) : FieldRegionBase_t(other)
+    {
+    }
+
+    /// \brief  Constructor from the arguments
+    /// \param  args       Arguments, which have to be passed to the base class constructors
+    /// \param  fieldType  Type of the field
+    template<typename... Args>
+    FieldRegion(Args... args, EFieldType fieldType) : FieldRegionBase_t(args...)
+                                                    , fFieldType(fieldType)
+    {
+    }
+
+    /// \brief Destructor
+    ~FieldRegion() = default;
+
+    /// \brief Copy assignment operator
+    /// \tparam I  Underlying floating point type of the source
+    template<typename I>
+    FieldRegion& operator=(const FieldRegion<I, FldMode>& other)
+    {
+      if (this == &other) {
+        return *this;
+      }
+      fFieldType                        = other.fFieldType;
+      return FieldRegionBase_t::operator=(other);
+    }
+
+    /// \brief Gets the field value in a spatial point on the track
+    /// \param x  x-coordinate of the point [cm]
+    /// \param y  y-coordinate of the point [cm]
+    /// \param z  z-coordinate of the point [cm]
+    /// \note  The x and y coordinates are ignored, if the interpolated field is used
+    FieldValue_t Get(const T& x, const T& y, const T& z) const;
+
+    /// \brief Gets the field type
+    EFieldType GetFieldType() const { return fFieldType; }
+
+    /// \brief String representation of the class content
+    /// \param intdentLevel  Indent level
+    std::string ToString(int indentLevel = 0) const;
+
+   private:
+    /// \brief Default constructor
+    FieldRegion() = default;
+
+    /// \brief Serialization method
+    friend class boost::serialization::access;
+    template<class Archive>
+    void serialize(Archive& ar, const unsigned int)
+    {
+      ar& boost::serialization::base_object<FieldRegionBase_t>;
+      ar& fFieldType;
+    }
+
+    EFieldType fFieldType{EFieldType::Null};  ///< Field type in a given region
+  };
+
+
+  // -------------------------------------------------------------------------------------------------------------------
+  //
+  template<typename T>
+  inline FieldValue<T> GlobalField::GetFieldValue(const FieldFn_t& fieldFn, const T& x, const T& y, const T& z)
+  {
+    using utils::simd::Cast;
+    FieldValue<T> B;
+    for (size_t i = 0; i < utils::simd::Size<T>(); ++i) {
+      auto [bx, by, bz] = fieldFn(Cast<T, double>(x, i), Cast<T, double>(y, i), Cast<T, double>(z, i));
+      B.SetSimdEntry(bx, by, bz, i);
+    }
+    return B;
+  }
+
+}  // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfFieldSlice.cxx b/algo/kf/core/geo/KfFieldSlice.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b5f94447294ab8898fe5a37a97b43ddade9d4a01
--- /dev/null
+++ b/algo/kf/core/geo/KfFieldSlice.cxx
@@ -0,0 +1,146 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   KfFieldSlice.h
+/// \brief  A class for a magnetic field approximation on a transverse plane (source)
+/// \author Sergei Zharko <s.zharko@gsi.de>
+/// \since  13.08.2024
+
+#include "KfFieldSlice.h"
+
+#include <iomanip>
+#include <sstream>
+
+using cbm::algo::kf::FieldFn_t;
+using cbm::algo::kf::FieldSlice;
+using cbm::algo::kf::FieldValue;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+FieldSlice<T>::FieldSlice(const FieldFn_t& fieldFn, double xMax, double yMax, double zRef) : fZref(zRef)
+{
+  // SLE initialization
+  // Augmented matrix N x (N + 3), where N - number of coefficients, defining the polynomial
+  //
+  std::array<std::array<double, kNofCoeff + 3>, kNofCoeff> A = {{}};
+
+  const double dx = (xMax < 2. * kNofCoeff) ? (xMax / kNofCoeff / 4.) : 1.;
+  const double dy = (yMax < 2. * kNofCoeff) ? (yMax / kNofCoeff / 4.) : 1.;
+
+  // Point and field arrays to access the field
+  std::array<double, 3> field;
+
+  // Monomial values for each coefficient
+  std::array<double, kNofCoeff> m = {{}};
+
+  // Fill the augmented matrix
+  for (double x = -xMax; x <= xMax; x += dx) {
+    for (double y = -yMax; y <= yMax; y += dy) {
+      const double r2 = x * x / (xMax * xMax) + y * y / (yMax * yMax);
+      if (r2 > 1.) {
+        continue;
+      }
+
+      std::tie(field[0], field[1], field[2]) = fieldFn(x, y, zRef);
+      // Fill the monomial array
+      {
+        m[kNofCoeff - 1] = 1.;
+        for (int i = kNofCoeff - 2; i >= kNofCoeff - 1 - kPolDegree; --i) {
+          m[i] = y * m[i + 1];
+        }
+        double xFactor = x;
+        int k          = ((kPolDegree - 1) * (kPolDegree + 2)) / 2;  // index of the monomial vector element
+        for (int i = kPolDegree; i > 0; --i) {                       // loop over y powers
+          for (int j = kNofCoeff; j >= kNofCoeff - i; --j) {         // loop over x powers
+            m[k] = xFactor * m[j];
+          }
+          xFactor *= x;
+        }
+      }
+
+      {
+        double w = 1. / (r2 + 1.);
+        for (int i = 0; i < kNofCoeff; ++i) {
+          // Fill the left part
+          for (int j = 0; j < kNofCoeff; ++j) {
+            A[i][j] += w * m[i] * m[j];
+          }
+          // Fill the right part
+          for (int j = 0; j < 3; ++j) {
+            A[i][kNofCoeff + j] += w * field[j] * m[i];
+          }
+        }
+      }
+    }
+  }
+
+  // SLE solution using Gaussian elimination
+  {
+    for (int k = 0; k < kNofCoeff - 1; ++k) {
+      for (int j = k + 1; j < kNofCoeff; ++j) {
+        double factor = A[j][k] / A[k][k];
+        for (int i = 0; i < kNofCoeff + 3; ++i) {
+          A[j][i] -= factor * A[k][i];
+        }
+      }
+    }
+    for (int k = kNofCoeff - 1; k > 0; --k) {
+      for (int j = k - 1; j >= 0; --j) {
+        double factor = A[j][k] / A[k][k];
+        for (int i = k; i < kNofCoeff + 3; ++i) {
+          A[j][i] -= factor * A[k][i];
+        }
+      }
+    }
+  }
+
+  for (int iD = 0; iD < 3; ++iD) {
+    auto& coeff = fCoeff[iD];
+    for (int j = 0; j < kNofCoeff; ++j) {
+      coeff[j] = utils::simd::Cast<double, T>(A[j][kNofCoeff + iD] / A[j][j]);
+    }
+  }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+std::string FieldSlice<T>::ToString(int indentLevel) const
+{
+  using std::setw;
+  constexpr char indentChar = '\t';
+  std::stringstream msg;
+  std::string indent(indentLevel, indentChar);
+
+  std::array<std::string, kNofCoeff> m;
+  {
+    m[kNofCoeff - 1] = "";
+    for (int i = kNofCoeff - 2; i >= kNofCoeff - 1 - kPolDegree; --i) {
+      m[i] = std::string("y") + std::to_string(kNofCoeff - i - 1);
+    }
+    int k = ((kPolDegree - 1) * (kPolDegree + 2)) / 2;  // index of the monomial vector element
+    for (int i = kPolDegree; i > 0; --i) {              // loop over y powers
+      std::string x = std::string("x") + std::to_string(kPolDegree - i + 1);
+      for (int j = kNofCoeff; j >= kNofCoeff - i; --j) {  // loop over x powers
+        m[k] = x + m[j];
+      }
+    }
+    m[kNofCoeff - 1] = "1";
+  }
+  msg << indent << setw(15) << "Monomial" << ' ' << setw(15) << "Bx" << ' ' << setw(15) << "By" << ' ' << setw(15)
+      << "Bz\n";
+  for (int i = 0; i < kNofCoeff; ++i) {
+    msg << indent << setw(15) << m[i] << ' ' << setw(15) << fCoeff[0][i] << ' ' << setw(15) << fCoeff[1][i] << ' '
+        << setw(15) << fCoeff[2][i] << '\n';
+  }
+  return msg.str();
+}
+
+namespace cbm::algo::kf
+{
+  template class FieldSlice<float>;
+  template class FieldSlice<double>;
+  template class FieldSlice<fvec>;
+}  // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfFieldSlice.h b/algo/kf/core/geo/KfFieldSlice.h
new file mode 100644
index 0000000000000000000000000000000000000000..11fc02ad6b5df4a5378b4e99e9541573180e31eb
--- /dev/null
+++ b/algo/kf/core/geo/KfFieldSlice.h
@@ -0,0 +1,116 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   KfFieldSlice.h
+/// \brief  A class for a magnetic field approximation on a transverse plane (header)
+/// \author Sergei Zharko <s.zharko@gsi.de>
+/// \since  13.08.2024
+
+#pragma once
+
+#include "KfDefs.h"
+#include "KfFieldValue.h"
+#include "KfMath.h"
+#include "KfTrackParam.h"
+#include "KfUtils.h"
+
+#include <boost/serialization/access.hpp>
+#include <boost/serialization/array.hpp>
+
+#include <array>
+#include <functional>
+#include <string>
+
+namespace cbm::algo::kf
+{
+  /// \class FieldSlice
+  /// \brief A magnetic field approximation on the two-dimensional plane
+  /// \tparam T Underlying data-type
+  template<typename T>
+  class alignas(VcMemAlign) FieldSlice {
+    static constexpr int kPolDegree{5};                                         ///< Approximation polynomial degree
+    static constexpr int kNofCoeff{((kPolDegree + 1) * (kPolDegree + 2)) / 2};  ///< Number of coefficients
+
+    /// \brief Array of the approximation coefficients [<dimension>][<index>]
+    using CoeffArray_t = std::array<std::array<T, kNofCoeff>, 3>;
+
+   public:
+    /// \brief Default constructor
+    FieldSlice() = default;
+
+    /// \brief Constructor from parameters
+    /// \param fieldFn  Magnetic field function (point xyz, field xyz)
+    /// \param xMax     Half-size of the slice in x-direction [cm]
+    /// \param yMax     Half-size of the slice in y-direction [cm]
+    /// \param zRef     Reference z-coordinate of the slice [cm]
+    FieldSlice(const FieldFn_t& fieldFn, double xMax, double yMax, double zRef);
+
+    /// \brief Copy constructor
+    /// \tparam I The underlying data type of the source object
+    template<typename I>
+    FieldSlice(const FieldSlice<I>& other)
+      : fCoeff(utils::simd::Cast(other.fCoeff))
+      , fZref(utils::simd::Cast(other.fZref))
+    {
+    }
+
+    /// \brief Destructor
+    ~FieldSlice() = default;
+
+    /// \brief Copy assignment operator
+    /// \tparam I The underlying data type of the source object
+    template<typename I>
+    FieldSlice& operator=(const FieldSlice<I>& other)
+    {
+      if (this != &other) {
+        this->fCoeff = utils::simd::Cast<I, T>(other.fCoeff);
+        this->fZref  = utils::simd::Cast<I, T>(other.fZref);
+      }
+      return *this;
+    }
+
+    /// \brief Gets field value at a point on the transverse plane
+    /// \param x  x-coordinate of the point [cm]
+    /// \param y  y-coordinate of the point [cm]
+    /// \return  Field value [kG]
+    constexpr FieldValue<T> GetFieldValue(const T& x, const T& y) const
+    {
+      using math::Horner;
+      /* clang-format off */
+      return FieldValue<T>(Horner<kPolDegree>(fCoeff[0].cbegin(), x, y), 
+                           Horner<kPolDegree>(fCoeff[1].cbegin(), x, y), 
+                           Horner<kPolDegree>(fCoeff[2].cbegin(), x, y));
+      /* clang-format on */
+    }
+
+    /// \brief Gets field value for the intersection with a straight track
+    /// \param trackParam  Track parameter set
+    /// \return  Field value [kG]
+    constexpr FieldValue<T> GetFieldValueForLine(const TrackParamBase<T>& trkPar) const
+    {
+      T dz = fZref - trkPar.GetZ();
+      return GetFieldValue(trkPar.GetX() + trkPar.GetTx() * dz, trkPar.GetY() + trkPar.GetTy() * dz);
+    }
+
+    /// \brief Gets reference z-coordinate of the slice [cm]
+    const T& GetZref() { return fZref; }
+
+    /// \brief String representation of the class content
+    /// \param intdentLevel  Indent level
+    std::string ToString(int indentLevel = 0) const;
+
+   private:
+    /// \brief Serialization function
+    friend class boost::serialization::access;
+    template<class Archive>
+    void serialize(Archive& ar, const unsigned int)
+    {
+      ar& fCoeff;
+      ar& fZref;
+    }
+
+    CoeffArray_t fCoeff{{T(0.)}};  ///< Approximation coefficients for the field slice
+    T fZref{defs::Undef<T>};       ///< Reference z of the coefficients [cm]
+  };
+}  // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfFieldValue.cxx b/algo/kf/core/geo/KfFieldValue.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..ddbcb98bf02ba204b13367f5b1ad1225f22d408f
--- /dev/null
+++ b/algo/kf/core/geo/KfFieldValue.cxx
@@ -0,0 +1,33 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   KfFieldValue.cxx
+/// \brief  Magnetic flux density vector representation (source)
+/// \author Sergei Zharko <s.zharko@gsi.de>
+/// \since  13.08.2024
+
+#include "KfFieldValue.h"
+
+#include <sstream>
+
+using cbm::algo::kf::FieldValue;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+std::string FieldValue<T>::ToString(int indentLevel) const
+{
+  constexpr char IndentChar = '\t';
+  std::stringstream msg;
+  std::string indent(indentLevel, IndentChar);
+  msg << indent << "B = {" << fB[0] << ", " << fB[1] << ", " << fB[2] << "} [kG]";
+  return msg.str();
+}
+
+namespace cbm::algo::kf
+{
+  template class FieldValue<float>;
+  template class FieldValue<double>;
+  template class FieldValue<fvec>;
+}  // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfFieldValue.h b/algo/kf/core/geo/KfFieldValue.h
new file mode 100644
index 0000000000000000000000000000000000000000..b26a3791b5c4ff3ea8d701857a9c64db899d3c0d
--- /dev/null
+++ b/algo/kf/core/geo/KfFieldValue.h
@@ -0,0 +1,182 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file   KfFieldValue.h
+/// \brief  Magnetic flux density vector representation
+/// \author Sergei Zharko <s.zharko@gsi.de>
+/// \since  13.08.2024
+
+#pragma once
+
+#include "KfDefs.h"
+#include "KfMath.h"
+#include "KfUtils.h"
+
+#include <boost/serialization/access.hpp>
+#include <boost/serialization/array.hpp>
+
+#include <array>
+#include <functional>
+#include <string>
+#include <tuple>
+
+namespace cbm::algo::kf
+{
+
+  /// \class  FieldValue
+  /// \brief  Magnetic flux density vector
+  /// \tparam T  Underlying data-type (float, double, fvec)
+  template<typename T>
+  class alignas(VcMemAlign) FieldValue {
+   public:
+    /// \brief Default constructor
+    /// \note  By default field is zero
+    FieldValue() = default;
+
+    /// \brief Constructor from components
+    /// \tparam I  Data type of the input parameters
+    /// \param bx  x-component of magnetic field vector [kG]
+    /// \param by  y-component of magnetic field vector [kG]
+    /// \param bz  z-component of magnetic field vector [kG]
+    template<typename I>
+    FieldValue(const I& bx, const I& by, const I& bz)
+      : fB({utils::simd::Cast<I, T>(bx), utils::simd::Cast<I, T>(by), utils::simd::Cast<I, T>(bz)})
+    {
+    }
+
+    /// \brief Copy constructor
+    /// \tparam I  Underlying data-type of the other object
+    template<typename I>
+    FieldValue(const FieldValue<I>& other) : fB(utils::simd::Cast<I, T, 3>(other.fB))
+    {
+    }
+
+    /// \brief Destructor
+    ~FieldValue() = default;
+
+    /// \brief Copy assignment operator
+    /// \tparam I  Underlying data-type of the other object
+    template<typename I>
+    FieldValue& operator=(const FieldValue<I>& other)
+    {
+      if (this != &other) {
+        this->fB = utils::simd::Cast<I, T, 3>(other.fB);
+      }
+      return *this;
+    }
+
+    /// \brief Combines the current magnetic field value with another one using a mask
+    /// \param other Other flux value
+    /// \param mask  Mask to perform a combination
+    void Combine(const FieldValue& other, const T& mask);
+
+    /// \brief Combines the current magnetic field value with another one using a mask
+    /// \param other Other flux value
+    /// \param mask  Mask to perform a combination
+    void Combine(const FieldValue& other, const bool mask);
+
+    /// \brief Combines the current magnetic field value with another one using a mask
+    /// \note  Only for T == fvec
+    /// \param other Other flux value
+    /// \param mask  Mask to perform a combination
+    void Combine(const FieldValue<fvec>& other, const fmask mask);
+
+    /// \brief Gets magnetic flux density x, y, z-components [kG]
+    /// \return  tuple (Bx, By, Bz)
+    [[gnu::always_inline]] std::tuple<T, T, T> Get() const { return std::make_tuple(fB[0], fB[1], fB[2]); }
+
+    /// \brief Gets squared absolute magnetic flux [kG2]
+    [[gnu::always_inline]] T GetAbsSq() const { return fB[0] * fB[0] + fB[1] * fB[1] + fB[2] * fB[2]; }
+
+    /// \brief Gets absolute magnetic flux [kG]
+    [[gnu::always_inline]] T GetAbs() const { return std::sqrt(this->GetAbsSq()); }
+
+    /// \brief Gets component by index
+    /// \param iD  index of the component (i = 0: x, 1: y, 2: z)
+    [[gnu::always_inline]] T GetComponent(int iD) const { return fB[iD]; }
+
+    /// \brief Gets magnetic flux x-component [kG]
+    [[gnu::always_inline]] T GetX() const { return fB[0]; }
+
+    /// \brief Gets magnetic flux density y-component [kG]
+    [[gnu::always_inline]] T GetY() const { return fB[1]; }
+
+    /// \brief Gets magnetic flux density z-component [kG]
+    [[gnu::always_inline]] T GetZ() const { return fB[2]; }
+
+    /// \brief Checks, if the field value is zero (negligible)
+    [[gnu::always_inline]] bool IsZero() const
+    {
+      auto bZero = this->GetAbsSq() <= defs::MinField<T> * defs::MinField<T>;
+      if constexpr (std::is_same_v<T, fvec>) {
+        return bZero.isFull();
+      }
+      else {
+        return bZero;
+      }
+    }
+
+    // TODO (?): Provide setters/accessors
+
+    /// \brief Sets magnetic flux density components to the field function
+    /// \param bx x-component of the magnetic flux density [kG]
+    /// \param by y-component of the magnetic flux density [kG]
+    /// \param bz z-component of the magnetic flux density [kG]
+    /// \param i  index of the SIMD entry (ignored, if T is not a SIMD vector)
+    [[gnu::always_inline]] void SetSimdEntry(double bx, double by, double bz, size_t i)
+    {
+      utils::simd::SetEntry(fB[0], bx, i);  //B.x[i] = bx;
+      utils::simd::SetEntry(fB[1], by, i);  //B.y[i] = by;
+      utils::simd::SetEntry(fB[2], bz, i);  //B.z[i] = bz;
+    }
+
+    /// \brief String representation of the field
+    std::string ToString(int indentLevel = 0) const;
+
+   private:
+    /// Serialization function
+    friend class boost::serialization::access;
+    template<class Archive>
+    void serialize(Archive& ar, const unsigned int)
+    {
+      ar& fB;
+    }
+
+    /// \brief Magnetic flux vector [kG]
+    GeoVector_t<T> fB{
+      {utils::simd::Cast<float, T>(0.F), utils::simd::Cast<float, T>(0.F), utils::simd::Cast<float, T>(0.F)}};
+  };
+
+  // -------------------------------------------------------------------------------------------------------------------
+  //
+  template<typename T>
+  [[gnu::always_inline]] inline void FieldValue<T>::Combine(const FieldValue& other, const T& mask)
+  {
+    for (size_t iD = 0; iD < fB.size(); ++iD) {
+      this->fB[iD] += mask * (other.fB[iD] - this->fB[iD]);
+    }
+  }
+
+  // -------------------------------------------------------------------------------------------------------------------
+  //
+  template<typename T>
+  [[gnu::always_inline]] inline void FieldValue<T>::Combine(const FieldValue& other, const bool mask)
+  {
+    if (mask) {
+      for (size_t iD = 0; iD < fB.size(); ++iD) {
+        this->fB[iD] += mask * (other.fB[iD] - this->fB[iD]);
+      }
+    }
+  }
+
+  // -------------------------------------------------------------------------------------------------------------------
+  //
+  template<>
+  [[gnu::always_inline]] inline void FieldValue<fvec>::Combine(const FieldValue<fvec>& other, const fmask mask)
+  {
+    for (size_t iD = 0; iD < fB.size(); ++iD) {
+      this->fB[iD](mask) = other.fB[iD];
+    }
+  }
+}  // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfGeoLayer.cxx b/algo/kf/core/geo/KfGeoLayer.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..eaa8e9214b6ca46b71fdbaa166bd7044598fc86e
--- /dev/null
+++ b/algo/kf/core/geo/KfGeoLayer.cxx
@@ -0,0 +1,30 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// @file   KfGeoLayer.cxx
+/// @brief  A logical layer of the KF-setup (implementation)
+/// @since  18.07.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#include "KfGeoLayer.h"
+
+using cbm::algo::kf::GeoLayer;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+GeoLayer::GeoLayer(kf::MaterialMap&& materialMap) : fMaterialMap(std::move(materialMap)) {}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+std::string GeoLayer::ToString(int verbose)
+{
+  if (verbose < 1) {
+    return std::string{};
+  }
+
+  // TODO: SZh. 19.07.2024: Implement
+  std::stringstream msg;
+  msg << "";
+  return msg.str();
+}
diff --git a/algo/kf/core/geo/KfGeoLayer.h b/algo/kf/core/geo/KfGeoLayer.h
new file mode 100644
index 0000000000000000000000000000000000000000..60bd954410c0eb33585fc6450a5155e1ac4ccd12
--- /dev/null
+++ b/algo/kf/core/geo/KfGeoLayer.h
@@ -0,0 +1,69 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// @file   KfGeoLayer.h
+/// @brief  A logical layer of the KF-setup (header)
+/// @since  18.07.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include "KfMaterialMap.h"
+
+namespace cbm::algo::kf
+{
+  /// \class GeoLayer
+  /// \brief A logical layer of the KF-setup for a given [zMin, zMax] interval
+  ///
+  /// A layer includes a material budget map and a magnetic field, estimated for a reference z-position.
+  class GeoLayer {
+   public:
+    /// \brief Default constructor
+    GeoLayer() = default;
+
+    /// \brief Constructor from the parameters
+    /// \param materialMap  Instance of the material budget map
+    ///
+    /// \note  The constructor moves the material budget map from the source
+    GeoLayer(kf::MaterialMap&& materialMap);
+
+    /// \brief Copy constructor
+    GeoLayer(const GeoLayer&) = default;
+
+    /// \brief Move constructor
+    GeoLayer(GeoLayer&&) = default;
+
+    /// \brief Copy assignment operator
+    GeoLayer& operator=(const GeoLayer&) = default;
+
+    /// \brief Move assignment operator
+    GeoLayer& operator=(GeoLayer&&) = default;
+
+    /// \brief Destructor
+    ~GeoLayer() = default;
+
+    /// \brief Material map accessor (mutable)
+    kf::MaterialMap& MaterialMap() { return fMaterialMap; }
+
+    /// \brief Material map accessor (constant)
+    const kf::MaterialMap& MaterialMap() const { return fMaterialMap; }
+
+    /// \brief String representation of the class
+    /// \param verbose  Verbosity level
+    std::string ToString(int verbose = 1);
+
+   private:
+    /// \brief Serialization method
+    friend class boost::serialization::access;
+    template<class Archive>
+    void serialize(Archive& ar, const unsigned int /*version*/)
+    {
+      ar& fMaterialMap;
+    }
+
+    kf::MaterialMap fMaterialMap;  ///< Material map instance for the geometry layer
+  };
+
+
+}  // namespace cbm::algo::kf
diff --git a/algo/kf/core/pars/KfMaterialMap.cxx b/algo/kf/core/geo/KfMaterialMap.cxx
similarity index 70%
rename from algo/kf/core/pars/KfMaterialMap.cxx
rename to algo/kf/core/geo/KfMaterialMap.cxx
index 3490bd3d387162d93a4ce8ad29c0a0320e55301c..9128811ff379a415d1772a2b84f0401a13acbf3d 100644
--- a/algo/kf/core/pars/KfMaterialMap.cxx
+++ b/algo/kf/core/geo/KfMaterialMap.cxx
@@ -3,17 +3,17 @@
    Authors: Sergey Gorbunov, Sergei Zharko [committer] */
 
 #include "KfMaterialMap.h"
-#include "KfSimd.h"
 
 #include "AlgoFairloggerCompat.h"
+#include "KfSimd.h"
 
+#include <cmath>
 #include <iomanip>
 #include <sstream>
 #include <vector>
-#include <cmath>
 
-using cbm::algo::kf::MaterialMap;
 using cbm::algo::kf::fvec;
+using cbm::algo::kf::MaterialMap;
 
 // ---------------------------------------------------------------------------------------------------------------------
 //
@@ -37,15 +37,15 @@ void MaterialMap::Add(const MaterialMap& other, float zTarg)
   // The function allows to add a material layer either from the left or from the right to the station
   // NOTE: A symmetry of x-y is assumed
   constexpr int nRays    = 3;                   // Number of rays in a bin in a dimension
-  const bool bRadialRays = !std::isnan(zTarg);  // Are rays radial (true, coming from target) or parallel (false)?  
-  const auto scaleFactor = bRadialRays ? ((other.fZref - zTarg) / (this->fZref - zTarg)) : 1.;  
-  const auto binSize     = 2. * scaleFactor * this->fXYmax / this->fNbins; // Size of each bin dimension [cm]
-  const auto stepSize    = binSize / static_cast<float>(nRays);            // Step between two neighboring rays [cm]
+  const bool bRadialRays = !std::isnan(zTarg);  // Are rays radial (true, coming from target) or parallel (false)?
+  const auto scaleFactor = bRadialRays ? ((other.fZref - zTarg) / (this->fZref - zTarg)) : 1.;
+  const auto binSize     = 2. * scaleFactor * this->fXYmax / this->fNbins;  // Size of each bin dimension [cm]
+  const auto stepSize    = binSize / static_cast<float>(nRays);             // Step between two neighboring rays [cm]
 
   // The coordinates of the first ray intersection with the other material layer [cm]
   float xBinOther = -this->fXYmax * scaleFactor + stepSize * 0.5;
-  float yBinOther = xBinOther; 
-  
+  float yBinOther = xBinOther;
+
   // Loop over bins of the active (this)
   for (int iBinY = 0; iBinY < this->fNbins; ++iBinY) {
     for (int iBinX = 0; iBinX < this->fNbins; ++iBinX) {
@@ -53,7 +53,7 @@ void MaterialMap::Add(const MaterialMap& other, float zTarg)
       float avgThickness = 0;  // Collected average thickness
       for (int iRayY = 0; iRayY < nRays; ++iRayY) {
         for (int iRayX = 0; iRayX < nRays; ++iRayX) {
-          avgThickness += other.GetRadThickScal(xBinOther + iBinX * stepSize, yBinOther + iBinY * stepSize);
+          avgThickness += other.GetRadThick(xBinOther + iBinX * stepSize, yBinOther + iBinY * stepSize);
         }
       }
       this->fTable[iBinX + this->fNbins * iBinY] += avgThickness / (nRays * nRays);
@@ -77,53 +77,6 @@ int MaterialMap::GetBin(float x, float y) const
   return i + j * fNbins;
 }
 
-// ---------------------------------------------------------------------------------------------------------------------
-//
-float MaterialMap::GetRadThickScal(float x, float y) const
-{
-  x     = (x < fXYmax && x >= -fXYmax) ? x : 0;
-  y     = (y < fXYmax && y >= -fXYmax) ? y : 0;
-  int i = static_cast<int>((x + fXYmax) * fFactor);
-  int j = static_cast<int>((y + fXYmax) * fFactor);
-  i     = (i < fNbins && i >= 0) ? i : fNbins / 2;
-  j     = (j < fNbins && j >= 0) ? j : fNbins / 2;
-  return fTable[i + j * fNbins];
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-fvec MaterialMap::GetRadThickVec(fvec x, fvec y) const
-{
-  fvec r;
-  for (size_t i = 0; i < fvec::size(); i++)
-    r[i] = GetRadThickScal(x[i], y[i]);
-  return r;
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void MaterialMap::CheckConsistency() const
-{
-  /* (i) Check, if the object was initialized */
-  // TODO: provide undef checking utility
-  //if (IsNaN()) {
-  //  throw std::logic_error("MaterialMap: class was not initialized");
-  //}
-
-  /* (ii) Check if the thickness values correct (non-negative) */
-  bool isThicknessOk = true;
-  for (int i = 0; i < int(fTable.size()); ++i) {
-    if (fTable[i] < 0) {
-      isThicknessOk = false;
-      LOG(error) << "MaterialMap: found illegal negative station thickness value " << fTable[i]
-                 << " at iBinX = " << (i % fNbins) << ", iBin = " << (i / fNbins);
-    }
-  }
-  if (!isThicknessOk) {
-    throw std::logic_error("MaterialMap: incorrect station thickness values found in the thickness map");
-  }
-}
-
 // ---------------------------------------------------------------------------------------------------------------------
 //
 void MaterialMap::Initialize(int nBins, float xyMax, float zRef, float zMin, float zMax)
diff --git a/algo/kf/core/pars/KfMaterialMap.h b/algo/kf/core/geo/KfMaterialMap.h
similarity index 76%
rename from algo/kf/core/pars/KfMaterialMap.h
rename to algo/kf/core/geo/KfMaterialMap.h
index e8ea76f508adc8c312930a091d13baac60adbb20..d336d84eadbbe602a86b3f95f489d93858990bea 100644
--- a/algo/kf/core/pars/KfMaterialMap.h
+++ b/algo/kf/core/geo/KfMaterialMap.h
@@ -9,19 +9,22 @@
 //#include "CaUtils.h"
 #include "AlgoFairloggerCompat.h"
 #include "KfDefs.h"
+#include "KfUtils.h"
 
 #include <boost/serialization/vector.hpp>
 
 #include <iomanip>
 #include <sstream>
 #include <string>
+#include <type_traits>
 #include <vector>
 
+//TODO: rewrite
 namespace cbm::algo::kf
 {
   /// \class MaterialMap
   /// \brief A map of station thickness in units of radiation length (X0) to the specific point in XY plane
-  class MaterialMap {
+  class alignas(VcMemAlign) MaterialMap {
    public:
     /// \brief Default constructor
     MaterialMap() = default;
@@ -49,16 +52,16 @@ namespace cbm::algo::kf
     /// \brief Gets number of bins (rows or columns) of the material table
     int GetNbins() const { return fNbins; }
 
-    /// \brief Gets radius in cm of the material table
+    /// \brief Gets radius of the material table [cm]
     float GetXYmax() const { return fXYmax; }
 
-    /// \brief Gets reference Z of the material in cm
+    /// \brief Gets reference Z of the material [cm]
     float GetZref() const { return fZref; }
 
-    /// \brief Gets minimal Z of the collected material in cm
+    /// \brief Gets minimal Z of the collected material [cm]
     float GetZmin() const { return fZmin; }
 
-    /// \brief Gets maximal Z of the collected material in cm
+    /// \brief Gets maximal Z of the collected material [cm]
     float GetZmax() const { return fZmax; }
 
     /// \brief Gets value of X/X0 in a given cell of the material table by the indeces of the row and column
@@ -70,28 +73,37 @@ namespace cbm::algo::kf
     /// \param iBin  Index of the bin in 2d table
     float GetRadThickBin(int iBin) const { return fTable[iBin]; }
 
-    /// \brief Gets material thickness in units of X0 in (x,y) point of the station
-    /// \param x  X coordinate of the point [cm]
-    /// \param y  Y coordinate of the point [cm]
-    float GetRadThickScal(float x, float y) const;
-
-    /// \brief Gets material thickness in units of X0 in (x,y) point of the station
-    /// cbm::algo::ca::fvec type can be float, that is why "Vec" and "Scal" specifications
-    /// \param x  X coordinate of the point [cm] (SIMDized vector)
-    /// \param y  Y coordinate of the point [cm] (SIMDized veotor)
-    // TODO: provide SIMD
-    fvec GetRadThickVec(fvec x, fvec y) const;
+    /// \brief   Gets material thickness in units of X0 in (x,y) point of the station
+    /// \tparam  I  Type of the x and y (floating point)
+    /// \param   x  X coordinate of the point [cm]
+    /// \param   y  Y coordinate of the point [cm]
+    template<typename I>
+    I GetRadThick(const I& x, const I& y) const
+    {
+      if constexpr (std::is_same_v<I, fvec>) {
+        fvec res;
+        for (size_t i = 0; i < utils::simd::Size<I>(); ++i) {
+          res[i] = GetRadThick(x[i], y[i]);
+        }
+        return res;
+      }
+      else {
+        I xNew = (x < fXYmax && x >= -fXYmax) ? x : 0;
+        I yNew = (y < fXYmax && y >= -fXYmax) ? y : 0;
+        int i  = static_cast<int>((xNew + fXYmax) * fFactor);
+        int j  = static_cast<int>((yNew + fXYmax) * fFactor);
+        i      = (i < fNbins && i >= 0) ? i : fNbins / 2;
+        j      = (j < fNbins && j >= 0) ? j : fNbins / 2;
+        return fTable[i + j * fNbins];
+      }
+    }
 
     /// \brief Checks, if the fields are NaN
     bool IsNaN() const
     {
-      return defs::IsUndefined(fNbins) || defs::IsUndefined(fXYmax) || defs::IsUndefined(fFactor)
-             || defs::IsUndefined(fZref) || defs::IsUndefined(fZmin) || defs::IsUndefined(fZmax);
+      return utils::IsUndefined(fNbins) || utils::IsUndefined(fXYmax * fFactor * fZref * fZmin * fZmax);
     }
 
-    /// \brief Verifies class invariant consistency
-    void CheckConsistency() const;
-
     /// \brief Sets value of material thickness in units of X0 for a given cell of the material table
     /// \param iBinX      Index of table column
     /// \param iBinY      Index of table row
diff --git a/algo/kf/core/pars/KfSetup.cxx b/algo/kf/core/geo/KfSetup.cxx
similarity index 100%
rename from algo/kf/core/pars/KfSetup.cxx
rename to algo/kf/core/geo/KfSetup.cxx
diff --git a/algo/kf/core/pars/KfSetup.h b/algo/kf/core/geo/KfSetup.h
similarity index 66%
rename from algo/kf/core/pars/KfSetup.h
rename to algo/kf/core/geo/KfSetup.h
index 1a5e57f59894e16010087c6e56f847fa2df0bb52..755b02ab5fd97142466fdee938346d6c27a4dfb6 100644
--- a/algo/kf/core/pars/KfSetup.h
+++ b/algo/kf/core/geo/KfSetup.h
@@ -10,7 +10,7 @@
 #pragma once  // include this header only once per compilation unit
 
 #include "KfDefs.h"
-#include "KfMaterialMap.h"
+#include "KfGeoLayer.h"
 #include "KfVector.h"
 
 #include <boost/serialization/access.hpp>
@@ -22,37 +22,12 @@ namespace cbm::algo::kf
 {
   using MaterialContainer_t = Vector<MaterialMap>;
 
-  /// \class  CollisionPoint
-  /// \brief  KF-framework representation of the nominal collision point region
-  /// \tparam DataT Underlying data type
-  template<typename DataT>
-  class alignas(defs::MemAlign) CollisionPoint {
-    /// \brief Default constructor
-    CollisionPoint() = default;
-
-    /// \brief Copy constructor
-    CollisionPoint(const CollisionPoint&) = default;
-
-    /// \brief Move constructor
-    CollisionPoint(CollisionPoint&&) = default;
-
-    /// \brief Destructor
-    ~CollisionPoint() = default;
-
-    /// \brief Copy assignment operator
-    CollisionPoint& operator=(const CollisionPoint&) = default;
-    
-    /// \brief Move assignment operator
-    CollisionPoint& operator=(CollisionPoint&&) = default;
-
-
-  };
-
   /// \class  Setup
   /// \brief  KF-framework representation of the detector setup
-  /// \tparam DataT Underlying data type
-  template<typename DataT>
-  class alignas(defs::MemAlign) Setup {
+  /// \tparam T  Underlying floating point data-type (float, double or fvec)
+  // TODO: Should the class have a template parameter?
+  template<typename T>
+  class alignas(VcMemAlign) Setup {
    public:
     /// \brief Default constructor
     Setup() = default;
@@ -60,15 +35,9 @@ namespace cbm::algo::kf
     /// \brief Destructor
     ~Setup() = default;
 
-    /// \brief Copy constructor
-    Setup(const Setup&) = default;
-
     /// \brief Move constructor
     Setup(Setup&&) = default;
 
-    /// \brief Copy assignment operator
-    Setup& operator=(const Setup&) = default;
-
     /// \brief Move assignment operator
     Setup& operator=(Setup&&) = default;
 
@@ -83,7 +52,7 @@ namespace cbm::algo::kf
     void CheckConsistency() const;
 
     /// \brief Gets material layer
-    [[gnu::always_inline]] const MaterialMap& GetMaterialLayer(int iLayer) const
+    [[gnu::always_inline]] const MaterialMap& GetLayer(int iLayer) const
     {
       return fvMatLayers.at<defs::GetterCheck>(iLayer);
     }
@@ -103,6 +72,16 @@ namespace cbm::algo::kf
     std::string ToString(int verbosity = 0, int indentLevel = 0) const;
 
    private:
+    /// \brief  Copy constructor
+    /// \tparam I  Underlying floating-point data type of the source
+    //template<typename I>
+    Setup(const Setup& other) = default;
+
+    /// \brief  Copy assignment operator
+    /// \tparam I  Underlying floating-point data type of the source
+    //template<typename I>
+    Setup& operator=(const Setup& other) = default;
+
     /// \brief Serialization method
     friend class boost::serialization::access;
     template<class Archive>
@@ -114,6 +93,7 @@ namespace cbm::algo::kf
 
     MaterialContainer_t fvMatLayers = {};  ///< Material layers map
 
-    bool fbInitialized = false;  ///< Is instance initialized
+    bool fbInitialized{false};         ///< Is instance initialized
+    bool fbOrigFieldAvailable{false};  ///< Is original field available
   };
 }  // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfTarget.cxx b/algo/kf/core/geo/KfTarget.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..cb564fbd11f854b8bac9c1afab76f790ab676252
--- /dev/null
+++ b/algo/kf/core/geo/KfTarget.cxx
@@ -0,0 +1,24 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// @file   KfTarget.cxx
+/// @brief  A target layer in the KF-setup (implementation)
+/// @since  19.07.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#include "KfTarget.h"
+
+using cbm::algo::kf::GeoLayer;
+using cbm::algo::kf::MaterialMap;
+using cbm::algo::kf::Target;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+Target::Target(kf::MaterialMap&& materialMap, float x, float y, float z)
+  : GeoLayer(std::move(materialMap))
+  , fX(x)
+  , fY(y)
+  , fZ(z)
+{
+}
diff --git a/algo/kf/core/geo/KfTarget.h b/algo/kf/core/geo/KfTarget.h
new file mode 100644
index 0000000000000000000000000000000000000000..f37bca6eac969c53db01ba6fb5b2c4c3b722030d
--- /dev/null
+++ b/algo/kf/core/geo/KfTarget.h
@@ -0,0 +1,86 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// @file   KfTarget.h
+/// @brief  A target layer in the KF-setup (header)
+/// @since  18.07.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include "KfDefs.h"
+#include "KfGeoLayer.h"
+
+#include <boost/serialization/access.hpp>
+#include <boost/serialization/base_object.hpp>
+
+namespace cbm::algo::kf
+{
+  /// \class Target
+  /// \brief A geometry layer in the target region
+  class Target : public GeoLayer {
+   public:
+    /// \brief Default constructor
+    Target() = default;
+
+    /// \brief Constructor from parameters
+    /// \param materialMap  Instance of the material budget map
+    /// \param x            x-coordinate of the nominal target center [cm]
+    /// \param y            y-coordinate of the nominal target center [cm]
+    /// \param z            z-coordinate of the nominal target center [cm]
+    Target(kf::MaterialMap&& materialMap, float x, float y, float z);
+
+    /// \brief Copy constructor
+    Target(const Target&) = default;
+
+    /// \brief Move constructor
+    Target(Target&&) = default;
+
+    /// \brief Copy assignment operator
+    Target& operator=(const Target&) = default;
+
+    /// \brief Move assignment operator
+    Target& operator=(Target&&) = default;
+
+    /// \brief Destructor
+    ~Target() = default;
+
+    /// \brief Gets x-coordinate of the nominal target center
+    float GetX() const { return fX; }
+
+    /// \brief Gets x-coordinate of the nominal target center
+    float GetY() const { return fY; }
+
+    /// \brief Gets x-coordinate of the nominal target center
+    float GetZ() const { return fZ; }
+
+    /// \brief Sets x-coordinate of the nominal target center
+    /// \param x  x-coordinate [cm]
+    void SetX(float x) { fX = x; }
+
+    /// \brief Sets y-coordinate of the nominal target center
+    /// \param y  y-coordinate [cm]
+    void SetY(float y) { fY = y; }
+
+    /// \brief Sets x-coordinate of the nominal target center
+    /// \param z  x-coordinate [cm]
+    void SetZ(float z) { fZ = z; }
+
+   private:
+    /// \brief Serialization method
+    friend class boost::serialization::access;
+    template<class Archive>
+    void serialize(Archive& ar, const unsigned int /*version*/)
+    {
+      ar& boost::serialization::base_object<GeoLayer>(*this);
+      ar& fX;
+      ar& fY;
+      ar& fZ;
+    }
+
+    float fX = defs::Undef<float>;  ///< x-coordinate of the nominal target center [cm]
+    float fY = defs::Undef<float>;  ///< y-coordinate of the nominal target center [cm]
+    float fZ = defs::Undef<float>;  ///< z-coordinate of the nominal target center [cm]
+  };
+}  // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/README.md b/algo/kf/core/geo/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..c787296acbe7d29f15b8e0978c510813e2f4fe69
--- /dev/null
+++ b/algo/kf/core/geo/README.md
@@ -0,0 +1 @@
+Directory kf/core/geo: contains classes representing the setup geometry in the KF-fitter package
diff --git a/algo/kf/core/pars/README.md b/algo/kf/core/pars/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..cfbfddc6188ec314d847365714f86942ac7dbd8d
--- /dev/null
+++ b/algo/kf/core/pars/README.md
@@ -0,0 +1,2 @@
+Directory kf/core/par: parameter and regimes definition
+
diff --git a/algo/kf/core/utils/KfMath.h b/algo/kf/core/utils/KfMath.h
new file mode 100644
index 0000000000000000000000000000000000000000..00c158916ea9033106859726f8019f6b95c144c7
--- /dev/null
+++ b/algo/kf/core/utils/KfMath.h
@@ -0,0 +1,86 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// @file   KfMath.h
+/// @brief  Collection of generic mathematical methods
+/// @since  30.07.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include <cstddef>
+
+namespace cbm::algo::kf::math
+{
+  /// \brief Number of coefficients in a polynomial
+  /// \param N  Degree of the polynomial
+  /// \param M  Number of dimensions
+  constexpr size_t NofPolCoefficients(size_t N, size_t M)
+  {
+    return M == 1 ? N + 1 : (NofPolCoefficients(N, M - 1) * (M + N)) / M;
+  }
+
+  /// \brief Horner's scheme for a 1D-polynomial estimation
+  /// \tparam  T  Underlying data type
+  /// \tparam  N  Degree of the polynomial
+  /// \param   c  Pointer to the first element of the polynomial coefficients
+  /// \param   x  Variable, for which the polynomial is to be estimated
+  ///
+  /// The polynomial coefficients are indexed as follows:
+  ///
+  /// P<N>(x) = x * (... * (x * (x * c[0] + c[1]) + c[2]) + ...) + c[N]
+  ///         = c[0] * x^N + c[1] * x^(N-1) + ... + c[N - 1] * x + c[N]
+  template<size_t N, typename T>
+  constexpr T Horner(const T* c, const T& x)
+  {
+    if constexpr (N == 0) {
+      return c[0];
+    }
+    else {
+      return Horner<N - 1, T>(c, x) * x + c[N];
+    }
+  }
+
+  /// \brief Horner's scheme for a multidvariable polynomial estimation
+  /// \tparam  T  Underlying data type
+  /// \tparam  N  Degree of the polynomial
+  /// \param   c  Pointer to the first element of the polynomial coefficients
+  /// \param  x1  First variable of the polynomial
+  /// \param  xI  Other variables: x2, ..., xM
+  ///
+  /// Example of the polynomial coefficient indexing for N = 4, M = 3:
+  ///
+  /// P(x,y,z) =
+  ///            c[0]  * x4 +
+  ///
+  ///            c[1]  * x3y +
+  ///            c[2]  * x3z  + c[3]  * x3 +
+  ///
+  ///            c[4]  * x2y2 +
+  ///            c[5]  * x2yz + c[6]  * x2y +
+  ///            c[7]  * x2z2 + c[8]  * x2z + c[9]  * x2 +
+  ///
+  ///            c[10] * xy3  +
+  ///            c[11] * xy2z + c[12] * xy2 +
+  ///            c[13] * xyz2 + c[14] * xyz + c[15] * xy +
+  ///            c[16] * xz3  + c[17] * xz2 + c[18] * xz + c[19] * x +
+  ///
+  ///            c[20] * y4   +
+  ///            c[21] * y3z  + c[22] * y3  +
+  ///            c[23] * y2z2 + c[24] * y2z + c[25] * y2 +
+  ///            c[26] * yz3  + c[27] * yz2 + c[28] * yz + c[29] * y +
+  ///            c[30] * z4   + c[31] * z3  + c[32] * z2 + c[33] * z + c[34];
+  template<size_t N, typename T, typename... Args>
+  constexpr T Horner(const T* c, const T& x1, Args... xI)
+  {
+    if constexpr (N == 0) {
+      return c[0];
+    }
+    else {
+      constexpr size_t M = sizeof...(Args) + 1;
+      constexpr size_t k = NofPolCoefficients(N - 1, M);
+      return Horner<N - 1, T>(c, x1, xI...) * x1 + Horner<N, T>(c + k, xI...);
+    }
+  }
+}  // namespace cbm::algo::kf::math
diff --git a/algo/kf/core/utils/KfSimdPseudo.h b/algo/kf/core/utils/KfSimdPseudo.h
index e6783ab61dc8339457edd5975954a05fd4edb602..bfdb741408860539cf22cde6499571ca63db4e7d 100644
--- a/algo/kf/core/utils/KfSimdPseudo.h
+++ b/algo/kf/core/utils/KfSimdPseudo.h
@@ -27,6 +27,8 @@ namespace cbm::algo::kf
   fscal sgn(fscal x);
 
 
+  constexpr auto VcMemAlign = 16;
+
   class fmask {
 
    public:
diff --git a/algo/kf/core/utils/KfSimdVc.h b/algo/kf/core/utils/KfSimdVc.h
index 20394e3fa4fd532b9533a930e41a3a29157c2ca4..65a664c933b9dcd92f2004e4f1cac893676f5aee 100644
--- a/algo/kf/core/utils/KfSimdVc.h
+++ b/algo/kf/core/utils/KfSimdVc.h
@@ -16,10 +16,13 @@ namespace cbm::algo::kf
   //using fvec = Vc::double_v;
   //using fvec = Vc::Vector<float, Vc::VectorAbi::Scalar>;
   //using fvec = Vc::SimdArray<float, 4>;
-  
+
   using fscal = fvec::EntryType;
   using fmask = fvec::MaskType;
 
+  constexpr auto VcMemAlign = Vc::VectorAlignment;
+
+// TODO: redefine via C++11 alignas
 #define _fvecalignment __attribute__((aligned(Vc::VectorAlignment)))
 }  // namespace cbm::algo::kf
 
diff --git a/algo/kf/core/utils/KfSimdUtils.cxx b/algo/kf/core/utils/KfUtils.cxx
similarity index 99%
rename from algo/kf/core/utils/KfSimdUtils.cxx
rename to algo/kf/core/utils/KfUtils.cxx
index 3a94f5ffc539b25fedac1bbf857813df91a89095..11e2b9a34b34520f09dffe084217e059adfad7a5 100644
--- a/algo/kf/core/utils/KfSimdUtils.cxx
+++ b/algo/kf/core/utils/KfUtils.cxx
@@ -2,7 +2,7 @@
    SPDX-License-Identifier: GPL-3.0-only
    Authors: Sergey Gorbunov, Sergei Zharko [committer] */
 
-#include "KfSimdUtils.h"
+#include "KfUtils.h"
 
 /// Namespace contains compile-time constants definition for the CA tracking algorithm
 ///
@@ -47,6 +47,7 @@ namespace cbm::algo::kf::utils::math
 
   void CholeskyFactorization(const double a[], const int n, int nn, double u[], int* nullty, int* ifault)
   {
+    // TODO: move to KfMath.h
     //
     //  Purpose:
     //
@@ -210,6 +211,7 @@ namespace cbm::algo::kf::utils::math
 
   void SymInv(const double a[], const int n, double c[], double w[], int* nullty, int* ifault)
   {
+    // TODO: move to KfMath.h
     //
     //  Purpose:
     //
diff --git a/algo/kf/core/utils/KfSimdUtils.h b/algo/kf/core/utils/KfUtils.h
similarity index 92%
rename from algo/kf/core/utils/KfSimdUtils.h
rename to algo/kf/core/utils/KfUtils.h
index 45de4a8af2fc20ea7efaaa4b21068c9243209601..560bad4cc62129e26329ee16a30124ec2d38ea79 100644
--- a/algo/kf/core/utils/KfSimdUtils.h
+++ b/algo/kf/core/utils/KfUtils.h
@@ -12,6 +12,7 @@
 #include "KfDefs.h"
 #include "KfSimd.h"
 
+#include <array>
 #include <sstream>
 
 /// Namespace contains compile-time constants definition for the CA tracking algorithm
@@ -46,7 +47,7 @@ namespace cbm::algo::kf::utils
   /// \brief Checks whether a variable of a particular type defined
   /// \param val Value to be checked
   template<typename T>
-  inline bool IsUndefined(T val)
+  inline bool IsUndefined(const T& val)
   {
     if constexpr (std::is_same_v<T, float>) {
       return std::isnan(val);
@@ -68,7 +69,7 @@ namespace cbm::algo::kf::utils
   /// \brief Checks whether a variable of a particular type is finite
   /// \param val  Value to check
   template<typename T>
-  inline bool IsFinite(T val)
+  inline bool IsFinite(const T& val)
   {
     if constexpr (std::is_same_v<T, float>) {
       return std::isfinite(val);
@@ -98,6 +99,7 @@ namespace cbm::algo::kf::utils
   }
 
   /// \brief Checks, if a SIMD vector horizontally equal
+  /// \note  Throws std::logic_error, if check fails
   /// TODO: Find this method in the VC!
   template<>
   void CheckSimdVectorEquality(fvec v, const char* name);
@@ -239,6 +241,24 @@ namespace cbm::algo::kf::utils::simd
     return static_cast<double>(val[i]);
   }
 
+  /// \brief  Converts a value of type DataT at a specific index to type DataOut
+  /// \details This function is a specialization for the array conversions
+  /// \tparam DataT Input value data type
+  /// \tparam DataOut Converted value data type
+  /// \tparam N Size of the array
+  /// \param val Input value of type DataT to be converted
+  /// \param i The index indicating which element to convert, ignored for non-SIMD DataT
+  /// \return Return val as DataOut
+  template<typename DataT, typename DataOut, size_t N>
+  inline std::array<DataOut, N> Cast(const std::array<DataT, N>& in, size_t i)
+  {
+    std::array<DataOut, N> ret;
+    for (size_t j = 0; j < N; ++j) {
+      ret[j] = Cast<DataT, DataOut>(in[j], i);
+    }
+    return ret;
+  }
+
   /// \brief Returns the number of elements available for a given data type
   /// \details This function is a template that provides a default implementation returning 1
   /// \tparam DataT The data type for which the size is determined
diff --git a/algo/test/CMakeLists.txt b/algo/test/CMakeLists.txt
index ee412e2f2507529fe976d43233d4dac40d43ada7..a0c7cde46b4f6342e1c3e940c7ceae02897c34d6 100644
--- a/algo/test/CMakeLists.txt
+++ b/algo/test/CMakeLists.txt
@@ -14,6 +14,7 @@ Function(AddBasicTest name)
   )
 EndFunction()
 
+#AddBasicTest(_GTestKf)
 AddBasicTest(_GTestTimeClusterTrigger)
 AddBasicTest(_GTestEventBuilder)
 AddBasicTest(_GTestDigiEventSelector)
diff --git a/reco/L1/qa/CbmCaOutputQa.cxx b/reco/L1/qa/CbmCaOutputQa.cxx
index 0e6ff2e369880228bb97169794f648023ada8313..e4a646149649aeef6d93fe04f24a7913c88635dd 100644
--- a/reco/L1/qa/CbmCaOutputQa.cxx
+++ b/reco/L1/qa/CbmCaOutputQa.cxx
@@ -445,7 +445,7 @@ void OutputQa::ExecQa()
   // ** Fill distributions for reconstructed tracks (including ghost ones) **
   // ************************************************************************
 
-  for (size_t iTrkReco = 0; iTrkReco < nRecoTracks; ++iTrkReco) {
+  for (size_t iTrkReco = 0; iTrkReco < static_cast<size_t>(nRecoTracks); ++iTrkReco) {
     const auto& recoTrk = fvRecoTracks[iTrkReco];
 
     // Reject tracks, which do not contain hits
diff --git a/reco/detectors/rich/finder/CbmL1RichENNRingFinderParallel.cxx b/reco/detectors/rich/finder/CbmL1RichENNRingFinderParallel.cxx
index 83446dbf3246fa980d3c5054e4e5e04adc254758..3d2dfd24a9b539e0742e28a33f71435f54df5745 100644
--- a/reco/detectors/rich/finder/CbmL1RichENNRingFinderParallel.cxx
+++ b/reco/detectors/rich/finder/CbmL1RichENNRingFinderParallel.cxx
@@ -26,7 +26,7 @@
 #include "CbmEvent.h"
 #include "CbmRichHit.h"
 #include "CbmRichRing.h"
-#include "KfSimdUtils.h"
+#include "KfUtils.h"
 #include "TClonesArray.h"
 #include "TStopwatch.h"
 #include "assert.h"