diff --git a/algo/ca/core/pars/CaField.cxx b/algo/ca/core/pars/CaField.cxx
index a4b24f596cb9f1c6064016a17ba08477457fe360..b2666bd454bd0fbdb2692cef8949267553037697 100644
--- a/algo/ca/core/pars/CaField.cxx
+++ b/algo/ca/core/pars/CaField.cxx
@@ -28,35 +28,6 @@ std::function<std::tuple<double, double, double>(double x, double y, double z)>
 template<typename DataT>
 bool FieldRegion<DataT>::fgIsZeroOriginalField = false;
 
-//----------------------------------------------------------------------------------------------------------------------
-//
-template<typename DataT>
-void FieldValue<DataT>::CheckConsistency() const
-{
-  /* Check SIMD data vectors for consistent initialization */
-  kfutils::CheckSimdVectorEquality(x, "FieldValue::x");
-  kfutils::CheckSimdVectorEquality(y, "FieldValue::y");
-  kfutils::CheckSimdVectorEquality(z, "FieldValue::z");
-
-  // TODO: Any other checks? (S.Zharko)
-}
-
-//----------------------------------------------------------------------------------------------------------------------
-// TODO:
-template<typename DataT>
-std::string FieldValue<DataT>::ToString(int) const
-{
-  std::stringstream aStream{};
-  aStream << "B = (" << x << ", " << y << ", " << z << ") [kG]";
-  return aStream.str();
-}
-
-namespace cbm::algo::ca
-{
-  template class FieldValue<fvec>;
-  template class FieldValue<float>;
-  template class FieldValue<double>;
-}  // namespace cbm::algo::ca
 
 //
 // FieldSlice methods
@@ -89,6 +60,7 @@ void FieldSlice<DataT>::CheckConsistency() const
   kfutils::CheckSimdVectorEquality(z, "FieldSlice: z");
 }
 
+
 //----------------------------------------------------------------------------------------------------------------------
 // TODO: Should it be inline? (S.Zharko)
 template<typename DataT>
@@ -116,19 +88,20 @@ FieldValue<DataT> FieldSlice<DataT>::GetFieldValue(const DataT& x, const DataT&
   DataT x3y2 = x3 * y2;
   DataT x4y  = x4 * y;
 
-  FieldValue<DataT> B;
-
-  B.x = cx[0] + cx[1] * x + cx[2] * y + cx[3] * x2 + cx[4] * xy + cx[5] * y2 + cx[6] * x3 + cx[7] * x2y + cx[8] * xy2
-        + cx[9] * y3 + cx[10] * x4 + cx[11] * x3y + cx[12] * x2y2 + cx[13] * xy3 + cx[14] * y4 + cx[15] * x5
-        + cx[16] * x4y + cx[17] * x3y2 + cx[18] * x2y3 + cx[19] * xy4 + cx[20] * y5;
-
-  B.y = cy[0] + cy[1] * x + cy[2] * y + cy[3] * x2 + cy[4] * xy + cy[5] * y2 + cy[6] * x3 + cy[7] * x2y + cy[8] * xy2
-        + cy[9] * y3 + cy[10] * x4 + cy[11] * x3y + cy[12] * x2y2 + cy[13] * xy3 + cy[14] * y4 + cy[15] * x5
-        + cy[16] * x4y + cy[17] * x3y2 + cy[18] * x2y3 + cy[19] * xy4 + cy[20] * y5;
+  FieldValue<DataT> B(
+    // Bx
+    cx[0] + cx[1] * x + cx[2] * y + cx[3] * x2 + cx[4] * xy + cx[5] * y2 + cx[6] * x3 + cx[7] * x2y + cx[8] * xy2
+      + cx[9] * y3 + cx[10] * x4 + cx[11] * x3y + cx[12] * x2y2 + cx[13] * xy3 + cx[14] * y4 + cx[15] * x5
+      + cx[16] * x4y + cx[17] * x3y2 + cx[18] * x2y3 + cx[19] * xy4 + cx[20] * y5,
+    // By
+    cy[0] + cy[1] * x + cy[2] * y + cy[3] * x2 + cy[4] * xy + cy[5] * y2 + cy[6] * x3 + cy[7] * x2y + cy[8] * xy2
+      + cy[9] * y3 + cy[10] * x4 + cy[11] * x3y + cy[12] * x2y2 + cy[13] * xy3 + cy[14] * y4 + cy[15] * x5
+      + cy[16] * x4y + cy[17] * x3y2 + cy[18] * x2y3 + cy[19] * xy4 + cy[20] * y5,
+    // Bz
+    cz[0] + cz[1] * x + cz[2] * y + cz[3] * x2 + cz[4] * xy + cz[5] * y2 + cz[6] * x3 + cz[7] * x2y + cz[8] * xy2
+      + cz[9] * y3 + cz[10] * x4 + cz[11] * x3y + cz[12] * x2y2 + cz[13] * xy3 + cz[14] * y4 + cz[15] * x5
+      + cz[16] * x4y + cz[17] * x3y2 + cz[18] * x2y3 + cz[19] * xy4 + cz[20] * y5);
 
-  B.z = cz[0] + cz[1] * x + cz[2] * y + cz[3] * x2 + cz[4] * xy + cz[5] * y2 + cz[6] * x3 + cz[7] * x2y + cz[8] * xy2
-        + cz[9] * y3 + cz[10] * x4 + cz[11] * x3y + cz[12] * x2y2 + cz[13] * xy3 + cz[14] * y4 + cz[15] * x5
-        + cz[16] * x4y + cz[17] * x3y2 + cz[18] * x2y3 + cz[19] * xy4 + cz[20] * y5;
   return B;
 }
 
@@ -139,6 +112,7 @@ FieldValue<DataT> FieldSlice<DataT>::GetFieldValueForLine(const TrackParamBase<D
   return GetFieldValue(t.GetX() + t.GetTx() * dz, t.GetY() + t.GetTy() * dz);
 }
 
+
 //----------------------------------------------------------------------------------------------------------------------
 //
 template<typename DataT>
@@ -225,9 +199,9 @@ FieldValue<DataT> FieldRegion<DataT>::Get(const DataT& x, const DataT& y, const
   else {
     DataT dz  = z - z0;
     DataT dz2 = dz * dz;
-    B.x       = cx0 + cx1 * dz + cx2 * dz2;
-    B.y       = cy0 + cy1 * dz + cy2 * dz2;
-    B.z       = cz0 + cz1 * dz + cz2 * dz2;
+    B.Set(cx0 + cx1 * dz + cx2 * dz2,   // Bx
+          cy0 + cy1 * dz + cy2 * dz2,   // By
+          cz0 + cz1 * dz + cz2 * dz2);  // Bz
   }
   return B;
 }
@@ -239,7 +213,7 @@ template<typename DataT>
 void FieldRegion<DataT>::Set(const FieldValue<DataT>& b0, const DataT b0z, const FieldValue<DataT>& b1, const DataT b1z,
                              const FieldValue<DataT>& b2, const DataT b2z)
 {
-  fIsZeroField = (b0.IsZeroField() && b1.IsZeroField() && b2.IsZeroField());
+  fIsZeroField = (b0.IsZero() && b1.IsZero() && b2.IsZero());
 
   z0        = b0z;
   DataT dz1 = b1z - b0z;
@@ -251,21 +225,21 @@ void FieldRegion<DataT>::Set(const FieldValue<DataT>& b0, const DataT b0z, const
   DataT w11 = -dz2 * w21;
   DataT w12 = -dz1 * w22;
 
-  DataT db1 = b1.x - b0.x;
-  DataT db2 = b2.x - b0.x;
-  cx0       = b0.x;
+  DataT db1 = b1.GetBx() - b0.GetBx();
+  DataT db2 = b2.GetBx() - b0.GetBx();
+  cx0       = b0.GetBx();
   cx1       = db1 * w11 + db2 * w12;
   cx2       = db1 * w21 + db2 * w22;
 
-  db1 = b1.y - b0.y;
-  db2 = b2.y - b0.y;
-  cy0 = b0.y;
+  db1 = b1.GetBy() - b0.GetBy();
+  db2 = b2.GetBy() - b0.GetBy();
+  cy0 = b0.GetBy();
   cy1 = db1 * w11 + db2 * w12;
   cy2 = db1 * w21 + db2 * w22;
 
-  db1 = b1.z - b0.z;
-  db2 = b2.z - b0.z;
-  cz0 = b0.z;
+  db1 = b1.GetBz() - b0.GetBz();
+  db2 = b2.GetBz() - b0.GetBz();
+  cz0 = b0.GetBz();
   cz1 = db1 * w11 + db2 * w12;
   cz2 = db1 * w21 + db2 * w22;
 }
diff --git a/algo/ca/core/pars/CaField.h b/algo/ca/core/pars/CaField.h
index 5aaeadbe9fe113961d39aa956bdd49a4ffd6067f..a12ae18e5b5010ce506ff44e6cd993ab5d8aa055 100644
--- a/algo/ca/core/pars/CaField.h
+++ b/algo/ca/core/pars/CaField.h
@@ -16,110 +16,7 @@
 
 namespace cbm::algo::ca
 {
-  template<typename DataT>
-  class FieldValue {
-
-   public:
-    DataT x{0.f};  //< x-component of the field
-    DataT y{0.f};  //< y-component of the field
-    DataT z{0.f};  //< z-component of the field
-
-    FieldValue() = default;
-
-    /// \brief Copy constructor with type conversion
-    template<typename DataIn>
-    FieldValue(const FieldValue<DataIn>& other)
-      : x(kfutils::simd::Cast<DataIn, DataT>(other.x))
-      , y(kfutils::simd::Cast<DataIn, DataT>(other.y))
-      , z(kfutils::simd::Cast<DataIn, DataT>(other.z))
-    {
-    }
-
-    /// Combines the magnetic field with another field value object using weight
-    /// \param B  other field value to combine with
-    /// \param w  weight from 0 to 1
-    void Combine(const FieldValue<DataT>& B, const DataT w);
-
-    void Combine(const FieldValue<DataT>& B, const bool w);
-
-    void Combine(const FieldValue<fvec>& B, const fmask w);
-
-    /// Is the field value zero?
-    bool IsZeroField() const
-    {
-      constexpr auto e = constants::misc::NegligibleFieldkG;
-      auto isZero      = (x * x + y * y + z * z <= e * e);
-      if constexpr (std::is_same<DataT, fvec>::value) {
-        return isZero.isFull();
-      }
-      else {
-        return isZero;
-      }
-    }
-
-    /// Consistency checker
-    void CheckConsistency() const;
-
-    /// Operator << overloading
-    friend std::ostream& operator<<(std::ostream& out, const FieldValue<DataT>& B)
-    {
-      return out << B.x << " | " << B.y << " | " << B.z;
-    };
-
-    /// String representation of class contents
-    /// \param indentLevel      number of indent characters in the output
-    std::string ToString(int indentLevel = 0) const;
-
-    /// \brief Sets the entries of the field components into the i-th elements of the SIMD vectors
-    /// \details If DataT is not a SIMD type, the function simply overwrites values ignoring the index i
-    /// \param bx The value to set for the x-component of the field
-    /// \param by The value to set for the y-component of the field
-    /// \param bz The value to set for the z-component of the field
-    /// \param i The index at which the SIMD entries are set (ignored if DataT is not a SIMD type)
-    [[gnu::always_inline]] inline void SetSimdEntry(double bx, double by, double bz, size_t i)
-    {
-      kfutils::simd::SetEntry(x, bx, i);  //B.x[i]            = bx;
-      kfutils::simd::SetEntry(y, by, i);  //B.y[i]            = by;
-      kfutils::simd::SetEntry(z, bz, i);  //B.z[i]            = bz;
-    }
-
-    /// Serialization function
-    friend class boost::serialization::access;
-    template<class Archive>
-    void serialize(Archive& ar, const unsigned int)
-    {
-      ar& x;
-      ar& y;
-      ar& z;
-    }
-  } _fvecalignment;
-
-  template<typename DataT>
-  [[gnu::always_inline]] inline void FieldValue<DataT>::Combine(const FieldValue<DataT>& B, const DataT w)
-  {
-    x += w * (B.x - x);
-    y += w * (B.y - y);
-    z += w * (B.z - z);
-  }
-
-  template<typename DataT>
-  [[gnu::always_inline]] inline void FieldValue<DataT>::Combine(const FieldValue<DataT>& B, const bool w)
-  {
-    if (w) {
-      x = B.x;
-      y = B.y;
-      z = B.z;
-    }
-  }
-
-  // TODO: Add specific masking implementation for vir (G.Kozlov)
-  template<>
-  [[gnu::always_inline]] inline void FieldValue<fvec>::Combine(const FieldValue<fvec>& B, const fmask w)
-  {
-    x(w) = B.x;
-    y(w) = B.y;
-    z(w) = B.z;
-  }
+  using cbm::algo::kf::FieldValue;
 
   /// Class represents a set of magnetic field approximation coefficients
   ///
diff --git a/algo/ca/core/pars/CaInitManager.cxx b/algo/ca/core/pars/CaInitManager.cxx
index 067c91d0c8e826f505d91f2fbabb42b7db6cea9e..443226f132504c799258d986d33808a71119dff9 100644
--- a/algo/ca/core/pars/CaInitManager.cxx
+++ b/algo/ca/core/pars/CaInitManager.cxx
@@ -301,10 +301,8 @@ void InitManager::InitTargetField(double zStep)
     double point[nDimensions]{0., 0., inputNodalZ[idx]};
     double field[nDimensions]{};
     fFieldFunction(point, field);
-    z[idx]   = inputNodalZ[idx];
-    B[idx].x = field[0];
-    B[idx].y = field[1];
-    B[idx].z = field[2];
+    z[idx] = inputNodalZ[idx];
+    B[idx].Set(field[0], field[1], field[2]);
   }  // loop over nodal points: end
   fParameters.fVertexFieldRegion.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
   fParameters.fVertexFieldValue = B[0];
diff --git a/algo/ca/core/pars/CaParameters.cxx b/algo/ca/core/pars/CaParameters.cxx
index 7b09d411bedde6ecf2112b430df5306f26c12342..af9522bd4fd2edcfef48acd44e1456bb81fa3bd9 100644
--- a/algo/ca/core/pars/CaParameters.cxx
+++ b/algo/ca/core/pars/CaParameters.cxx
@@ -266,11 +266,11 @@ std::string Parameters<DataT>::ToString(int verbosity, int indentLevel) const
   msg << indent << indentCh << clrs::CLb << "FIELD:\n" << clrs::CL;
   msg << indent << indentCh << indentCh << "Target field:\n";
   msg << indent << indentCh << indentCh << indentCh
-      << "Bx = " << kfutils::simd::Cast<DataT, float>(fVertexFieldValue.x, 0) << " Kg\n";
+      << "Bx = " << kfutils::simd::Cast<DataT, float>(fVertexFieldValue.GetBx(), 0) << " Kg\n";
   msg << indent << indentCh << indentCh << indentCh
-      << "By = " << kfutils::simd::Cast<DataT, float>(fVertexFieldValue.y, 0) << " Kg\n";
+      << "By = " << kfutils::simd::Cast<DataT, float>(fVertexFieldValue.GetBy(), 0) << " Kg\n";
   msg << indent << indentCh << indentCh << indentCh
-      << "Bz = " << kfutils::simd::Cast<DataT, float>(fVertexFieldValue.z, 0) << " Kg\n";
+      << "Bz = " << kfutils::simd::Cast<DataT, float>(fVertexFieldValue.GetBz(), 0) << " Kg\n";
 
 
   msg << indent << indentCh << clrs::CLb << "NUMBER OF STATIONS:\n" << clrs::CL;
diff --git a/algo/ca/core/tracking/CaTrackFit.cxx b/algo/ca/core/tracking/CaTrackFit.cxx
index 8d6899a8ef196ae9b6d7a825a8f5f1abecd300c5..f90285b3191ac957dbf62f1c781fecba11fe4b0c 100644
--- a/algo/ca/core/tracking/CaTrackFit.cxx
+++ b/algo/ca/core/tracking/CaTrackFit.cxx
@@ -705,9 +705,9 @@ namespace cbm::algo::ca
       DataT z = fTr.GetZ() + stepDz[step];
 
       ca::FieldValue B = Field.Get(rstep[0], rstep[1], z);
-      DataT Bx         = B.x;
-      DataT By         = B.y;
-      DataT Bz         = B.z;
+      DataT Bx         = B.GetBx();
+      DataT By         = B.GetBy();
+      DataT Bz         = B.GetBz();
       // NOTE: SZh 13.08.2024: The rstep[0] and rstep[1] make no effect on the B value, if Field is an approximation
 
       DataT tx    = rstep[2];
diff --git a/algo/ca/core/tracking/CaTrackFitter.cxx b/algo/ca/core/tracking/CaTrackFitter.cxx
index 41afdee5b5b9ab844f44ea5488129eda37563437..6b3b8407459195a32185744940b6aaf9d6eb1083 100644
--- a/algo/ca/core/tracking/CaTrackFitter.cxx
+++ b/algo/ca/core/tracking/CaTrackFitter.cxx
@@ -166,9 +166,8 @@ void TrackFitter::FitCaTracks(const ca::InputData& input, WindowData& wData)
         mxy[ista].NdfX()[iVec] = 1.;
         mxy[ista].NdfY()[iVec] = 1.;
 
-        fB[ista].x[iVec] = fB_temp.x[iVec];
-        fB[ista].y[iVec] = fB_temp.y[iVec];
-        fB[ista].z[iVec] = fB_temp.z[iVec];
+        fB[ista].SetSimdEntry(fB_temp.GetBx()[iVec], fB_temp.GetBy()[iVec], fB_temp.GetBz()[iVec], iVec);
+
         if (ih == 0) {
           z_start[iVec]          = z[ista][iVec];
           x_first[iVec]          = x[ista][iVec];
@@ -241,12 +240,12 @@ void TrackFitter::FitCaTracks(const ca::InputData& input, WindowData& wData)
 
       fldB1 = sta[ista].fieldSlice.GetFieldValue(tr.X(), tr.Y());
 
-      fldB1.Combine(fB[ista], w[ista]);
+      fldB1.SetSimdEntries(fB[ista], w[ista]);
 
       fldZ2   = z[ista - 2];
       fvec dz = fldZ2 - fldZ1;
       fldB2   = sta[ista].fieldSlice.GetFieldValue(tr.X() + tr.Tx() * dz, tr.Y() + tr.Ty() * dz);
-      fldB2.Combine(fB[ista - 2], w[ista - 2]);
+      fldB2.SetSimdEntries(fB[ista - 2], w[ista - 2]);
       fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
 
       for (--ista; ista >= 0; ista--) {
@@ -254,7 +253,7 @@ void TrackFitter::FitCaTracks(const ca::InputData& input, WindowData& wData)
         fldZ0 = z[ista];
         dz    = (fldZ1 - fldZ0);
         fldB0 = sta[ista].fieldSlice.GetFieldValue(tr.X() - tr.Tx() * dz, tr.Y() - tr.Ty() * dz);
-        fldB0.Combine(fB[ista], w[ista]);
+        fldB0.SetSimdEntries(fB[ista], w[ista]);
         fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2);
 
         fmask initialised = (z[ista] < z_end) & (z_start <= z[ista]);
@@ -348,19 +347,19 @@ void TrackFitter::FitCaTracks(const ca::InputData& input, WindowData& wData)
 
       fldZ1 = z[ista];
       fldB1 = sta[ista].fieldSlice.GetFieldValue(tr.X(), tr.Y());
-      fldB1.Combine(fB[ista], w[ista]);
+      fldB1.SetSimdEntries(fB[ista], w[ista]);
 
       fldZ2 = z[ista + 2];
       dz    = fldZ2 - fldZ1;
       fldB2 = sta[ista].fieldSlice.GetFieldValue(tr.X() + tr.Tx() * dz, tr.Y() + tr.Ty() * dz);
-      fldB2.Combine(fB[ista + 2], w[ista + 2]);
+      fldB2.SetSimdEntries(fB[ista + 2], w[ista + 2]);
       fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
 
       for (++ista; ista < nStations; ista++) {
         fldZ0 = z[ista];
         dz    = (fldZ1 - fldZ0);
         fldB0 = sta[ista].fieldSlice.GetFieldValue(tr.X() - tr.Tx() * dz, tr.Y() - tr.Ty() * dz);
-        fldB0.Combine(fB[ista], w[ista]);
+        fldB0.SetSimdEntries(fB[ista], w[ista]);
         fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2);
 
         fmask initialised = (z[ista] <= z_end) & (z_start < z[ista]);
diff --git a/algo/ca/core/tracking/CaTripletConstructor.cxx b/algo/ca/core/tracking/CaTripletConstructor.cxx
index 04069c0e77131f486cf1a96a39b28c3017ec838f..53e843a733e354faa92a3904be9405c1076b0c88 100644
--- a/algo/ca/core/tracking/CaTripletConstructor.cxx
+++ b/algo/ca/core/tracking/CaTripletConstructor.cxx
@@ -35,7 +35,7 @@ TripletConstructor::TripletConstructor(const ca::Parameters<fvec>& pars, WindowD
                                      [](const ca::Station<fvec>& s, int edge) { return bool(s.fieldStatus) > edge; })
                     - fParameters.GetStations().cbegin();
 
-  fIsTargetField = (fabs(frWData.TargB().y[0]) > 0.001);
+  fIsTargetField = !frWData.TargB().IsZero();
 }
 
 
diff --git a/algo/kf/core/geo/KfFieldValue.cxx b/algo/kf/core/geo/KfFieldValue.cxx
index ddbcb98bf02ba204b13367f5b1ad1225f22d408f..3054b8579e7a54c2951d0656a81683454b633cad 100644
--- a/algo/kf/core/geo/KfFieldValue.cxx
+++ b/algo/kf/core/geo/KfFieldValue.cxx
@@ -11,10 +11,21 @@
 
 #include <sstream>
 
-using cbm::algo::kf::FieldValue;
+using namespace cbm::algo::kf;
 
 // ---------------------------------------------------------------------------------------------------------------------
 //
+template<typename T>
+void FieldValue<T>::CheckConsistency() const
+{
+  // Check SIMD data vectors for consistent initialization
+  utils::CheckSimdVectorEquality(fB[0], "FieldValue::x");
+  utils::CheckSimdVectorEquality(fB[1], "FieldValue::y");
+  utils::CheckSimdVectorEquality(fB[2], "FieldValue::z");
+
+  // TODO: Any other checks? (S.Zharko)
+}
+
 template<typename T>
 std::string FieldValue<T>::ToString(int indentLevel) const
 {
diff --git a/algo/kf/core/geo/KfFieldValue.h b/algo/kf/core/geo/KfFieldValue.h
index b4b1a8fc22c7832f11c2ae89f2fdd66faef6cd42..d9898579adaa832efda09797298a4751cb288123 100644
--- a/algo/kf/core/geo/KfFieldValue.h
+++ b/algo/kf/core/geo/KfFieldValue.h
@@ -58,21 +58,28 @@ namespace cbm::algo::kf
     /// \brief Copy assignment operator
     FieldValue& operator=(const FieldValue& other) = default;
 
-    /// \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 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>
+    void Set(const I& bx, const I& by, const I& bz)
+    {
+      fB[0] = utils::simd::Cast<I, T>(bx);
+      fB[1] = utils::simd::Cast<I, T>(by);
+      fB[2] = utils::simd::Cast<I, T>(bz);
+    }
 
     /// \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);
+    void SetSimdEntries(const FieldValue& other, const kf::utils::masktype<T>& mask)
+    {
+      for (size_t iD = 0; iD < fB.size(); ++iD) {
+        this->fB[iD] = kf::utils::iif(mask, other.fB[iD], this->fB[iD]);
+      }
+    }
 
     /// \brief Gets magnetic flux density x, y, z-components [kG]
     /// \return  tuple (Bx, By, Bz)
@@ -89,13 +96,13 @@ namespace cbm::algo::kf
     [[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]; }
+    [[gnu::always_inline]] T GetBx() const { return fB[0]; }
 
     /// \brief Gets magnetic flux density y-component [kG]
-    [[gnu::always_inline]] T GetY() const { return fB[1]; }
+    [[gnu::always_inline]] T GetBy() const { return fB[1]; }
 
     /// \brief Gets magnetic flux density z-component [kG]
-    [[gnu::always_inline]] T GetZ() const { return fB[2]; }
+    [[gnu::always_inline]] T GetBz() const { return fB[2]; }
 
     /// \brief Checks, if the field value is zero (negligible)
     [[gnu::always_inline]] bool IsZero() const
@@ -109,6 +116,9 @@ namespace cbm::algo::kf
       }
     }
 
+    /// Consistency checker
+    void CheckConsistency() const;
+
     // TODO (?): Provide setters/accessors
 
     /// \brief Sets magnetic flux density components to the field function
@@ -140,35 +150,4 @@ namespace cbm::algo::kf
       {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/reco/KF/ParticleFitter/CbmL1PFFitter.cxx b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
index 288c725b0db898f055de1362860f41144dfd3a5e..0b4613a5cbcfae9fd844f8e990ce3f7ce88cd58a 100644
--- a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
+++ b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
@@ -262,10 +262,8 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
       mxy[i].SetDxy(0.);
       mxy[i].SetNdfX(1.);
       mxy[i].SetNdfY(1.);
-      dt2[i]  = 1.;
-      fB[i].x = 0.;
-      fB[i].y = 0.;
-      fB[i].z = 0.;
+      dt2[i] = 1.;
+      fB[i].Set(0., 0., 0.);
     }
 
     for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
@@ -311,10 +309,8 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
         w[ista][iVec] = true;
 
 
-        fB_temp          = sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista]);
-        fB[ista].x[iVec] = fB_temp.x[iVec];
-        fB[ista].y[iVec] = fB_temp.y[iVec];
-        fB[ista].z[iVec] = fB_temp.z[iVec];
+        fB_temp = sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista]);
+        fB[ista].SetSimdEntry(fB_temp.GetBx()[iVec], fB_temp.GetBy()[iVec], fB_temp.GetBz()[iVec], iVec);
         if (i == 0) {
           z_start[iVec]          = z[ista][iVec];
           x_first[iVec]          = x[ista][iVec];
@@ -355,17 +351,17 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
 
     z1 = z[i];
     b1 = sta[i].fieldSlice.GetFieldValue(T.X(), T.Y());
-    b1.Combine(fB[i], w[i]);
+    b1.SetSimdEntries(fB[i], w[i]);
     z2 = z[i + 2];
     dz = z2 - z1;
     b2 = sta[i].fieldSlice.GetFieldValue(T.X() + T.Tx() * dz, T.Y() + T.Ty() * dz);
-    b2.Combine(fB[i + 2], w[i + 2]);
+    b2.SetSimdEntries(fB[i + 2], w[i + 2]);
     fld.Set(b2, z2, b1, z1, b0, z0);
     for (++i; i < nHits; i++) {
       z0 = z[i];
       dz = (z1 - z0);
       b0 = sta[i].fieldSlice.GetFieldValue(T.X() - T.Tx() * dz, T.Y() - T.Ty() * dz);
-      b0.Combine(fB[i], w[i]);
+      b0.SetSimdEntries(fB[i], w[i]);
       fld.Set(b0, z0, b1, z1, b2, z2);
 
       fmask initialised = (z[i] <= z_end) & (z_start < z[i]);
@@ -415,18 +411,18 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
 
     z1 = z[i];
     b1 = sta[i].fieldSlice.GetFieldValue(T.X(), T.Y());
-    b1.Combine(fB[i], w[i]);
+    b1.SetSimdEntries(fB[i], w[i]);
 
     z2 = z[i - 2];
     dz = z2 - z1;
     b2 = sta[i].fieldSlice.GetFieldValue(T.X() + T.Tx() * dz, T.Y() + T.Ty() * dz);
-    b2.Combine(fB[i - 2], w[i - 2]);
+    b2.SetSimdEntries(fB[i - 2], w[i - 2]);
     fld.Set(b2, z2, b1, z1, b0, z0);
     for (--i; i >= 0; i--) {
       z0 = z[i];
       dz = (z1 - z0);
       b0 = sta[i].fieldSlice.GetFieldValue(T.X() - T.Tx() * dz, T.Y() - T.Ty() * dz);
-      b0.Combine(fB[i], w[i]);
+      b0.SetSimdEntries(fB[i], w[i]);
       fld.Set(b0, z0, b1, z1, b2, z2);
 
       fmask initialised = (z[i] < z_end) & (z_start <= z[i]);
@@ -584,10 +580,8 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
           }
         }
 
-        fB_temp              = sta[ista].fieldSlice.GetFieldValue(posx, posy);
-        fB[iH + 1].x[iVec]   = fB_temp.x[iVec];
-        fB[iH + 1].y[iVec]   = fB_temp.y[iVec];
-        fB[iH + 1].z[iVec]   = fB_temp.z[iVec];
+        fB_temp = sta[ista].fieldSlice.GetFieldValue(posx, posy);
+        fB[iH + 1].SetSimdEntry(fB_temp.GetBx()[iVec], fB_temp.GetBy()[iVec], fB_temp.GetBz()[iVec], iVec);
         zField[iH + 1][iVec] = posz;
       }
     }
@@ -723,10 +717,8 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<PFF
           }
         }
 
-        fB_temp              = sta[ista].fieldSlice.GetFieldValue(posx, posy);
-        fB[iH + 1].x[iVec]   = fB_temp.x[iVec];
-        fB[iH + 1].y[iVec]   = fB_temp.y[iVec];
-        fB[iH + 1].z[iVec]   = fB_temp.z[iVec];
+        fB_temp = sta[ista].fieldSlice.GetFieldValue(posx, posy);
+        fB[iH + 1].SetSimdEntry(fB_temp.GetBx()[iVec], fB_temp.GetBy()[iVec], fB_temp.GetBz()[iVec], iVec);
         zField[iH + 1][iVec] = posz;
       }
     }
@@ -809,9 +801,7 @@ void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks,
 
         fB_temp = sta[ista].fieldSlice.GetFieldValue(posx, posy);
 
-        fB[iH].x[iVec]   = fB_temp.x[iVec];
-        fB[iH].y[iVec]   = fB_temp.y[iVec];
-        fB[iH].z[iVec]   = fB_temp.z[iVec];
+        fB[iH].SetSimdEntry(fB_temp.GetBx()[iVec], fB_temp.GetBy()[iVec], fB_temp.GetBz()[iVec], iVec);
         zField[iH][iVec] = posz;
       }
     }
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 53df76f1dcc4e763a1445600db9b76bf4465739d..3449bf4656639f9d40bb3dd92f0e465b044f014d 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -1641,12 +1641,12 @@ void CbmL1::FieldApproxCheck()
         bbb = sqrt(B[0] * B[0] + B[1] * B[1] + B[2] * B[2]);
 
         B_L1   = FSl.GetFieldValue(x, y);
-        bbb_L1 = sqrt(B_L1.x[0] * B_L1.x[0] + B_L1.y[0] * B_L1.y[0] + B_L1.z[0] * B_L1.z[0]);
+        bbb_L1 = B_L1.GetAbs()[0];
 
         stB->SetBinContent(ii, jj, (bbb - bbb_L1));
-        stBx->SetBinContent(ii, jj, (B[0] - B_L1.x[0]));
-        stBy->SetBinContent(ii, jj, (B[1] - B_L1.y[0]));
-        stBz->SetBinContent(ii, jj, (B[2] - B_L1.z[0]));
+        stBx->SetBinContent(ii, jj, (B[0] - B_L1.GetBx()[0]));
+        stBy->SetBinContent(ii, jj, (B[1] - B_L1.GetBy()[0]));
+        stBz->SetBinContent(ii, jj, (B[2] - B_L1.GetBz()[0]));
         j++;
       }
       i++;