diff --git a/algo/ca/core/pars/CaField.cxx b/algo/ca/core/pars/CaField.cxx
index 676a53135bedacfde382cdc1c72422105e0aef49..99b8ca5f21d171921d287aa6f25d4825692bab69 100644
--- a/algo/ca/core/pars/CaField.cxx
+++ b/algo/ca/core/pars/CaField.cxx
@@ -53,7 +53,7 @@ FieldRegion<DataT>::FieldRegion(float reg[10]) noexcept
template<typename DataT>
void FieldRegion<DataT>::CheckConsistency() const
{
- /* Check SIMD data vectors for consistent initialization */
+ // Check SIMD data vectors for consistent initialization
kfutils::CheckSimdVectorEquality(cx0, "FieldRegion::cx0");
kfutils::CheckSimdVectorEquality(cx1, "FieldRegion::cx1");
kfutils::CheckSimdVectorEquality(cx2, "FieldRegion::cx2");
@@ -94,6 +94,31 @@ FieldValue<DataT> FieldRegion<DataT>::Get(const DataT& x, const DataT& y, const
}
+//----------------------------------------------------------------------------------------------------------------------
+template<typename T>
+std::tuple<T, T, T> FieldRegion<T>::GetDoubleIntegrals(const T& /*x1*/, const T& /*y1*/, const T& z1, //
+ const T& /*x2*/, const T& /*y2*/, const T& z2) const
+{
+ // double integral of the field along z
+
+ if (fUseOriginalField) {
+ // TODO: implement the double integral for the original field
+ assert(false && "FieldRegion::GetDoubleIntegrals() is not implemented for the original field");
+ }
+
+ auto fld = *this;
+ fld.Shift(z1);
+ T dz = z2 - z1;
+ T dz2 = dz * dz;
+ T c0 = dz2 * T(1. / 2.);
+ T c1 = dz2 * dz * T(1. / 6.);
+ T c2 = dz2 * dz2 * T(1. / 12.);
+ return {c0 * fld.cx0 + c1 * fld.cx1 + c2 * fld.cx2, //
+ c0 * fld.cy0 + c1 * fld.cy1 + c2 * fld.cy2, //
+ c0 * fld.cz0 + c1 * fld.cz1 + c2 * fld.cz2};
+}
+
+
//----------------------------------------------------------------------------------------------------------------------
// TODO: Should it be inline? (S.Zharko)
template<typename DataT>
diff --git a/algo/ca/core/pars/CaField.h b/algo/ca/core/pars/CaField.h
index 023eed73858fec5ceadba4c81e17b6738865ed36..07d3cee5c1bb41e808203b2af7ee4cc9fa7e91b9 100644
--- a/algo/ca/core/pars/CaField.h
+++ b/algo/ca/core/pars/CaField.h
@@ -20,6 +20,8 @@ namespace cbm::algo::ca
using cbm::algo::kf::FieldSlice;
using cbm::algo::kf::FieldValue;
+ //template<typename T>
+ //using FieldRegion = cbm::algo::kf::FieldRegion<T, cbm::algo::kf::EFieldMode::Intrpl>;
template<typename DataT>
class FieldRegion {
@@ -56,6 +58,9 @@ namespace cbm::algo::ca
/// \return the magnetic field value
FieldValue<DataT> Get(const DataT& x, const DataT& y, const DataT& z) const;
+ std::tuple<DataT, DataT, DataT> GetDoubleIntegrals(const DataT& x1, const DataT& y1, const DataT& z1, //
+ const DataT& x2, const DataT& y2, const DataT& z2) const;
+
/// Is the field region empty?
bool IsZeroField() const
{
diff --git a/algo/ca/core/tracking/CaFramework.cxx b/algo/ca/core/tracking/CaFramework.cxx
index 1366844f5ec0da1ee78d3434cac338716b12bc95..eae683f55b19c6a8ee7a38788714dfc0ba198f37 100644
--- a/algo/ca/core/tracking/CaFramework.cxx
+++ b/algo/ca/core/tracking/CaFramework.cxx
@@ -62,6 +62,7 @@ void Framework::ReceiveParameters(Parameters<fvec>&& parameters)
fNstationsBeforePipe = fParameters.GetNstationsActive(static_cast<EDetectorID>(0));
ca::FieldRegion<fvec>::ForceUseOfOriginalField(fParameters.DevIsUseOfOriginalField());
+ kf::GlobalField::ForceUseOfOriginalField(fParameters.DevIsUseOfOriginalField());
}
// ---------------------------------------------------------------------------------------------------------------------
diff --git a/algo/ca/core/tracking/CaTrackFit.cxx b/algo/ca/core/tracking/CaTrackFit.cxx
index f90285b3191ac957dbf62f1c781fecba11fe4b0c..d31d4aa2c1ea1f5f5280d6bcd0fc0af3c527d39d 100644
--- a/algo/ca/core/tracking/CaTrackFit.cxx
+++ b/algo/ca/core/tracking/CaTrackFit.cxx
@@ -1169,28 +1169,26 @@ namespace cbm::algo::ca
// extrapolate track assuming it is straight (qp==0)
// return the extrapolated X, Y and the derivatives of the extrapolated X and Y
- cnst c_light(0.000299792458), c1(1.), c2i(0.5), c6i(1. / 6.), c12i(1. / 12.);
+ cnst c_light(0.000299792458);
cnst tx = fTr.GetTx();
cnst ty = fTr.GetTy();
DataT dz = z - fTr.GetZ();
- DataT dz2 = dz * dz;
- DataT dzc6i = dz * c6i;
- DataT dz2c12i = dz2 * c12i;
-
DataT xx = tx * tx;
DataT yy = ty * ty;
DataT xy = tx * ty;
- DataT Ay = -xx - c1;
- DataT Bx = yy + c1;
+ DataT Ay = -xx - DataT(1.);
+ DataT Bx = yy + DataT(1.);
+
+ DataT ct = c_light * sqrt(DataT(1.) + xx + yy);
+
- DataT ctdz2 = c_light * sqrt(c1 + xx + yy) * dz2;
+ DataT Sx, Sy, Sz;
+ std::tie(Sx, Sy, Sz) = F.GetDoubleIntegrals(fTr.GetX(), fTr.GetY(), fTr.GetZ(), //
+ fTr.GetX() + dz * tx, fTr.GetY() + dz * ty, z);
- DataT Sx = F.cx0 * c2i + F.cx1 * dzc6i + F.cx2 * dz2c12i;
- DataT Sy = F.cy0 * c2i + F.cy1 * dzc6i + F.cy2 * dz2c12i;
- DataT Sz = F.cz0 * c2i + F.cz1 * dzc6i + F.cz2 * dz2c12i;
extrX = fTr.GetX() + tx * dz;
extrY = fTr.GetY() + ty * dz;
@@ -1202,7 +1200,7 @@ namespace cbm::algo::ca
Jx[1] = DataT(0.);
Jx[2] = dz;
Jx[3] = DataT(0.);
- Jx[4] = ctdz2 * (Sx * xy + Sy * Ay + Sz * ty);
+ Jx[4] = ct * (Sx * xy + Sy * Ay + Sz * ty);
Jx[5] = DataT(0.);
Jx[6] = DataT(0.);
@@ -1210,7 +1208,7 @@ namespace cbm::algo::ca
Jy[1] = DataT(1.);
Jy[2] = DataT(0.);
Jy[3] = dz;
- Jy[4] = ctdz2 * (Sx * Bx - Sy * xy - Sz * tx);
+ Jy[4] = ct * (Sx * Bx - Sy * xy - Sz * tx);
Jy[5] = DataT(0.);
Jx[6] = DataT(0.);
}
diff --git a/algo/kf/core/KfDefs.h b/algo/kf/core/KfDefs.h
index 7a8a795b29da3eb98cac677b7d538b838f99c7f9..a71a871d19b55ffaefcc99eba821d763477a45a0 100644
--- a/algo/kf/core/KfDefs.h
+++ b/algo/kf/core/KfDefs.h
@@ -134,8 +134,6 @@ namespace cbm::algo::kf::defs
//
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
diff --git a/algo/kf/core/geo/KfField.cxx b/algo/kf/core/geo/KfField.cxx
index 29ecaebb24bdef96f7d74d818b829b4238b0a4c6..74172f541c1864ed0f3c855074a27ff8df8ce21a 100644
--- a/algo/kf/core/geo/KfField.cxx
+++ b/algo/kf/core/geo/KfField.cxx
@@ -37,7 +37,7 @@ std::string FieldBase<T, EFieldMode::Intrpl>::ToString(int indentLevel, int verb
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T>
-std::string FieldBase<T, EFieldMode::Orig>::ToString(int indentLevel, int verbose) const
+std::string FieldBase<T, EFieldMode::Orig>::ToString(int indentLevel, int /*verbose*/) const
{
constexpr char IndentChar = '\t';
std::stringstream msg;
diff --git a/algo/kf/core/geo/KfFieldRegion.cxx b/algo/kf/core/geo/KfFieldRegion.cxx
index 1f824eef447606d404b5c0d695d398520b222c81..f87839f61eea8e777a7fff03ef57dab0b5b98ba5 100644
--- a/algo/kf/core/geo/KfFieldRegion.cxx
+++ b/algo/kf/core/geo/KfFieldRegion.cxx
@@ -28,6 +28,8 @@ using cbm::algo::kf::detail::FieldRegionBase;
FieldFn_t GlobalField::gOriginalField = &GlobalField::UndefField;
+bool GlobalField::fgForceUseOfOriginalField = false;
+
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -60,7 +62,7 @@ void FieldRegionBase<T, EFieldMode::Intrpl>::Set(const FieldValue<T>& b0, const
template<typename T>
FieldValue<T> FieldRegionBase<T, EFieldMode::Intrpl>::Get(const T& x, const T& y, const T& z) const
{
- if (defs::dbg::ForceOriginalField) {
+ if (GlobalField::IsUsingOriginalFieldForced()) {
return GlobalField::GetFieldValue(GlobalField::gOriginalField, x, y, z);
}
auto dz = z - this->fZfirst;
@@ -84,6 +86,20 @@ void FieldRegionBase<T, EFieldMode::Intrpl>::Shift(const T& z)
fZfirst = z;
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+void FieldRegionBase<T, EFieldMode::Intrpl>::CheckConsistency() const
+{
+ // Check SIMD data vectors for consistent initialization
+ for (int iD = 0; iD < 3; ++iD) {
+ for (int iC = 0; iC < 3; ++iC) {
+ utils::CheckSimdVectorEquality(fCoeff[iD][iC], fmt::format("FieldRegion::fCoeff[{}][{}]", iD, iC).c_str());
+ }
+ }
+ utils::CheckSimdVectorEquality(fZfirst, "FieldRegion::fZfirst");
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T>
@@ -105,16 +121,27 @@ std::string FieldRegionBase<T, EFieldMode::Intrpl>::ToString(int indentLevel, in
msg << '\n'
<< indent << IndentChar
<< format("Bx(dz) = {:>12} + {:>12} dz + {:>12} dz2", Cnv(coeff[2][0]), Cnv(coeff[2][1]), Cnv(coeff[2][2]));
- if constexpr (defs::dbg::ForceOriginalField) {
- msg << indent << "\nWARNING: the defs::dbg::ForceOriginalField is enabled, so the magnetic field interpolation "
- << indent
- << "\nis replaced with the global magnetic function for debugging purposes. If you see this message and are "
- << indent
- << "\nnot sure, what is going on, please set the defs::dbg::ForceOriginalField to false and re-run the routine";
+ if (GlobalField::IsUsingOriginalFieldForced()) {
+ msg
+ << indent
+ << "\nWARNING: the GlobalField::ForceUseOfOriginalField() is enabled, so the magnetic field interpolation "
+ << indent
+ << "\nis replaced with the global magnetic function for debugging purposes. If you see this message and are "
+ << indent
+ << "\nnot sure, what is going on, please call GlobalField::ForceUseOfOriginalField(false) and re-run the routine";
}
return msg.str();
}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+void FieldRegionBase<T, EFieldMode::Orig>::CheckConsistency() const
+{
+ // Check SIMD data vectors for consistent initialization
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T>
@@ -175,11 +202,26 @@ std::string FieldRegion<T>::ToString(int indentLevel, int) const
return msg.str();
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+void FieldRegion<T>::CheckConsistency() const
+{
+ // Check SIMD data vectors for consistent initialization
+
+ if (fFieldMode == EFieldMode::Intrpl) {
+ foFldIntrpl->CheckConsistency();
+ }
+ else {
+ foFldOrig->CheckConsistency();
+ }
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
std::tuple<double, double, double> GlobalField::UndefField(double, double, double)
{
- assert(!defs::dbg::ForceOriginalField
+ assert(!GlobalField::IsUsingOriginalFieldForced()
&& "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.");
diff --git a/algo/kf/core/geo/KfFieldRegion.h b/algo/kf/core/geo/KfFieldRegion.h
index e31f1ce5af3655dba7f2c0bd1611bfbc3bd07b97..0ca1f66f9ca1fa46e8c552725ffa1d2a911ec73b 100644
--- a/algo/kf/core/geo/KfFieldRegion.h
+++ b/algo/kf/core/geo/KfFieldRegion.h
@@ -52,6 +52,15 @@ namespace cbm::algo::kf
template<typename T>
[[gnu::always_inline]] static FieldValue<T> GetFieldValue(const FieldFn_t& fieldFn, const T& x, const T& y,
const T& z);
+
+ /// \brief Forces the use of the original field function
+ static void ForceUseOfOriginalField(bool v = true) { fgForceUseOfOriginalField = v; }
+
+ /// \brief Checks if the original field function is used
+ static bool IsUsingOriginalFieldForced() { return fgForceUseOfOriginalField; }
+
+ private:
+ static bool fgForceUseOfOriginalField; ///< Flag to force the use of the original field in all field classes
};
namespace detail
@@ -103,6 +112,9 @@ namespace cbm::algo::kf
return GlobalField::GetFieldValue(fFieldFn, x, y, z);
}
+ /// Consistency checker
+ void CheckConsistency() const;
+
/// \brief String representation of the class content
/// \param intdentLevel Indent level
/// \param verbose Verbosity level
@@ -184,6 +196,9 @@ namespace cbm::algo::kf
/// \param z new z-coordinate [cm]
void Shift(const T& z);
+ /// Consistency checker
+ void CheckConsistency() const;
+
/// \brief String representation of the class content
/// \param intdentLevel Indent level
/// \param verbose Verbosity level
@@ -302,6 +317,9 @@ namespace cbm::algo::kf
/// \param z new z-coordinate [cm]
void Shift(const T& z);
+ /// Consistency checker
+ void CheckConsistency() const;
+
/// \brief String representation of the class content
/// \param intdentLevel Indent level
/// \param verbose Verbosity level
diff --git a/reco/KF/ParticleFitter/CbmL1PFFitter.h b/reco/KF/ParticleFitter/CbmL1PFFitter.h
index 06e81d6d130b392c179b576c836eedbf6f7cbcfb..759fde6768f0aaa1e6ce6fb1b3b48373140f7e85 100644
--- a/reco/KF/ParticleFitter/CbmL1PFFitter.h
+++ b/reco/KF/ParticleFitter/CbmL1PFFitter.h
@@ -20,6 +20,7 @@
#ifndef _CbmL1PFFitter_h_
#define _CbmL1PFFitter_h_
+#include "CaField.h"
#include "CaSimd.h"
#include <vector>
@@ -27,11 +28,6 @@
class CbmMvdHit;
class CbmStsHit;
class CbmStsTrack;
-namespace cbm::algo::ca
-{
- template<typename DataT>
- class FieldRegion;
-}
using namespace cbm::algo;
diff --git a/reco/L1/catools/CaToolsField.h b/reco/L1/catools/CaToolsField.h
index 88de2e69699dade09c09a985811da62d1f585343..b77604318cf976018dd38925896569fe4bb2e406 100644
--- a/reco/L1/catools/CaToolsField.h
+++ b/reco/L1/catools/CaToolsField.h
@@ -50,6 +50,7 @@ namespace cbm::ca::tools
ca::FieldRegion<float>::SetOriginalField(fld, true);
ca::FieldRegion<double>::SetOriginalField(fld, true);
}
+ cbm::algo::kf::GlobalField::SetFieldFunction(fld);
}
} // namespace cbm::ca::tools