diff --git a/algo/kf/core/geo/KfField.cxx b/algo/kf/core/geo/KfField.cxx
index 888971655eea6e73566e0a7f148e390000cdd79c..29ecaebb24bdef96f7d74d818b829b4238b0a4c6 100644
--- a/algo/kf/core/geo/KfField.cxx
+++ b/algo/kf/core/geo/KfField.cxx
@@ -20,34 +20,61 @@ using cbm::algo::kf::detail::FieldBase;
// ---------------------------------------------------------------------------------------------------------------------
//
-template<typename T, EFieldMode FldMode>
-std::string Field<T, FldMode>::ToString(int indentLevel, int verbose) const
+template<typename T>
+std::string FieldBase<T, EFieldMode::Intrpl>::ToString(int indentLevel, int verbose) const
+{
+ constexpr char IndentChar = '\t';
+ std::stringstream msg;
+ std::string indent(indentLevel, IndentChar);
+ msg << indent << "Field slices:";
+ for (int iSlice = 0; iSlice < GetNofFieldSlices(); ++iSlice) {
+ const auto& fldSlice = fvFieldSlices[iSlice];
+ msg << indent << "\n - " << fldSlice.ToString(indentLevel + 1, verbose);
+ }
+ return msg.str();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+std::string FieldBase<T, EFieldMode::Orig>::ToString(int indentLevel, int verbose) const
+{
+ constexpr char IndentChar = '\t';
+ std::stringstream msg;
+ std::string indent(indentLevel, IndentChar);
+ msg << indent << "Original field function";
+ msg << indent << "Field slice z-positions: ";
+ for (int iSlice = 0; iSlice < GetNofFieldSlices(); ++iSlice) {
+ msg << indent << "\n " << iSlice << ") " << fvFieldSliceZ[iSlice];
+ }
+ return msg.str();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+std::string Field<T>::ToString(int indentLevel, int verbose) const
{
constexpr char IndentChar = '\t';
std::stringstream msg;
std::string indent(indentLevel, IndentChar);
msg << indent << "Field near primary vertex:\n" << this->fPrimVertexField.ToString(indentLevel + 1, verbose) << '\n';
msg << indent << "Field type: " << static_cast<int>(this->fFieldType) << '\n';
- if constexpr (FldMode == EFieldMode::Orig) {
- //msg << indent << "Original field function";
+ msg << indent << "Field mode: " << static_cast<int>(this->fFieldMode) << '\n';
+ if (fFieldMode == EFieldMode::Orig) {
+ msg << foFldOrig->ToString(indentLevel, verbose);
}
- else if constexpr (FldMode == EFieldMode::Intrpl) {
- msg << indent << "Field slices:";
- for (const auto& fldSlice : this->fvFieldSlices) {
- msg << indent << "\n - " << fldSlice.ToString(indentLevel + 1, verbose);
- }
+ else if (fFieldMode == EFieldMode::Intrpl) {
+ msg << foFldIntrpl->ToString(indentLevel, verbose);
}
return 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>;
+ template class Field<float>;
+ template class Field<double>;
+ template class Field<fvec>;
} // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfField.h b/algo/kf/core/geo/KfField.h
index c0bb840be05b85ad8df200e2a59f0130d3bdd63a..1e07c1afb20a1ca60fb63d3533c13a205b4f11e7 100644
--- a/algo/kf/core/geo/KfField.h
+++ b/algo/kf/core/geo/KfField.h
@@ -13,7 +13,11 @@
#include "KfFieldRegion.h"
#include "KfFieldSlice.h"
+#include <boost/serialization/access.hpp>
+#include <boost/serialization/split_free.hpp>
+
#include <array>
+#include <optional>
#include <set>
namespace cbm::algo::kf
@@ -37,17 +41,6 @@ namespace cbm::algo::kf
friend class FieldBase;
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(fFieldType, fFieldFn); }
-
- protected:
/// \brief Default constructor
FieldBase() = default;
@@ -56,22 +49,56 @@ namespace cbm::algo::kf
template<typename I>
FieldBase(const FieldBase<I, EFieldMode::Orig>& other) : fFieldFn(other.fFieldFn)
{
+ fvFieldSliceZ.reserve(other.fvFieldSliceZ.size());
+ for (const auto& slice : other.fvFieldSliceZ) {
+ fvFieldSliceZ.emplace_back(utils::simd::Cast<I, T>(slice));
+ }
}
/// \brief Copy assignment operator
- FieldBase& operator=(const FieldBase& other) = default;
+ FieldBase& operator=(const FieldBase& other)
+ {
+ if (this != &other) {
+ fvFieldSliceZ = other.fvFieldSliceZ;
+ 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
+ /// \brief Sets a field slice (z-ref)
+ /// \param fieldSliceZref Reference z-coordinate of a field slice
+ void AddFieldSlice(const T& fieldSliceZref) { fvFieldSliceZ.push_back(fieldSliceZref); }
- private:
- /// Serialization function
- friend class boost::serialization::access;
- template<class Archive>
- void serialize(Archive& ar, const unsigned int)
+ /// \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
{
- ar& fFieldType;
+ return GlobalField::GetFieldValue(fFieldFn, x, y, fvFieldSliceZ[sliceID]);
}
+
+ /// \brief Gets field function
+ const FieldFn_t& GetFieldFunction() const { return fFieldFn; }
+
+ /// \brief Gets number of field slices in the instance
+ int GetNofFieldSlices() const { return fvFieldSliceZ.size(); }
+
+ /// \brief Makes field region
+ FieldRegion<T> MakeFieldRegion(EFieldType fldType) const { return FieldRegion(fldType, fFieldFn); }
+
+ /// \brief Sets magnetic field function
+ /// \param fieldFn Magnetic field function (KF-format)
+ void SetFieldFunction(const FieldFn_t& fieldFn) { fFieldFn = fieldFn; }
+
+ /// \brief String representation of the class
+ /// \param indentLevel Indent level of the string output
+ /// \param verbose Verbosity level
+ std::string ToString(int indentLevel, int verbose) const;
+
+ private:
+ std::vector<T> fvFieldSliceZ{}; ///< z-positions of field slices to emulate functionality
+ FieldFn_t fFieldFn{defs::ZeroFieldFn}; ///< Field function: (x,y,z) [cm] -> (Bx,By,Bz) [kG]
};
@@ -86,8 +113,28 @@ namespace cbm::algo::kf
using SlicesContainer_t = std::vector<FieldSlice<T>>;
public:
- using FieldRegion_t = FieldRegion<T, EFieldMode::Intrpl>;
- using FieldValue_t = FieldValue<T>;
+ /// \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)
+ {
+ fvFieldSlices.reserve(other.fvFieldSlices.size());
+ for (const auto& slice : other.fvFieldSlices) {
+ fvFieldSlices.emplace_back(FieldSlice<T>(slice));
+ }
+ }
+
+ /// \brief Copy assignment operator
+ FieldBase& operator=(const FieldBase& other)
+ {
+ if (this != &other) {
+ fvFieldSlices = other.fvFieldSlices;
+ }
+ return *this;
+ }
/// \brief Sets a field slice
/// \param fieldSlice Field slice object
@@ -97,7 +144,7 @@ namespace cbm::algo::kf
/// \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
+ FieldValue<T> GetFieldValue(int sliceID, const T& x, const T& y) const
{
return fvFieldSlices[sliceID].GetFieldValue(x, y);
}
@@ -106,39 +153,26 @@ namespace cbm::algo::kf
int GetNofFieldSlices() const { return fvFieldSlices.size(); }
/// \brief Creates field region object
+ /// \brief fldType FieldType
/// \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
+ FieldRegion<T> MakeFieldRegion(EFieldType fldType, 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(fFieldType, b0, z0, b1, z1, b2, z2);
+ return FieldRegion<T>(fldType, b0, z0, b1, z1, b2, z2);
}
- protected:
- /// \brief Default constructor
- FieldBase() = default;
+ /// \brief String representation of the class
+ /// \param indentLevel Indent level of the string output
+ /// \param verbose Verbosity level
+ std::string ToString(int indentLevel, int verbose) const;
- /// \brief Copy constructor
- /// \tparam I Underlying floating type of the source
- template<typename I>
- FieldBase(const FieldBase<I, EFieldMode::Intrpl>& other)
- {
- fvFieldSlices.clear();
- fvFieldSlices.reserve(other.fvFieldSlices.size());
- for (const auto& slice : other.fvFieldSlices) {
- fvFieldSlices.emplace_back(FieldSlice<T>(slice));
- }
- }
-
- /// \brief Copy assignment operator
- FieldBase& operator=(const FieldBase& other) = default;
-
- SlicesContainer_t fvFieldSlices; ///< Array of field slices
- EFieldType fFieldType{EFieldType::Null}; ///< Field type
+ protected:
+ SlicesContainer_t fvFieldSlices{}; ///< Array of field slices
private:
/// \brief Serialization function
@@ -147,7 +181,6 @@ namespace cbm::algo::kf
void serialize(Archive& ar, const unsigned int)
{
ar& fvFieldSlices;
- ar& fFieldType;
}
};
} // namespace detail
@@ -156,25 +189,39 @@ namespace cbm::algo::kf
/// \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> {
- template<typename, EFieldMode>
+ template<typename T>
+ class alignas(VcMemAlign) Field {
+ template<typename>
friend class Field;
+ friend class FieldFactory;
+ friend class boost::serialization::access;
public:
- static constexpr EFieldMode FieldType = FldMode;
-
- using FieldRegion_t = FieldRegion<T, FldMode>;
- using Real_t = T;
-
- /// \brief Default constructor
- Field() = default;
+ /// \brief Constructor
+ /// \param fldMode Field mode
+ Field(EFieldMode fldMode, EFieldType fldType)
+ : foFldIntrpl(fldMode == EFieldMode::Intrpl ? std::make_optional(detail::FieldBase<T, EFieldMode::Intrpl>())
+ : std::nullopt)
+ , foFldOrig(fldMode == EFieldMode::Orig ? std::make_optional(detail::FieldBase<T, EFieldMode::Orig>())
+ : std::nullopt)
+ , fPrimVertexField(FieldRegion<T>(fldMode, fldType))
+ , fFieldType(EFieldType::Normal)
+ , fFieldMode(fldMode)
+ {
+ }
/// \brief Copy constructor
template<typename I>
- Field(const Field<I, FldMode>& other)
- : detail::FieldBase<T, FldMode>(other)
- , fPrimVertexField(FieldRegion<I, FldMode>(other.fPrimVertexField))
+ Field(const Field<I>& other)
+ : foFldIntrpl(other.foFldIntrpl.has_value()
+ ? std::make_optional(detail::FieldBase<T, EFieldMode::Intrpl>(*other.foFldIntrpl))
+ : std::nullopt)
+ , foFldOrig(other.foFldOrig.has_value()
+ ? std::make_optional(detail::FieldBase<T, EFieldMode::Orig>(*other.foFldOrig))
+ : std::nullopt)
+ , fPrimVertexField(FieldRegion<I>(other.fPrimVertexField))
+ , fFieldType(other.fFieldType)
+ , fFieldMode(other.fFieldMode)
{
}
@@ -182,24 +229,48 @@ namespace cbm::algo::kf
~Field() = default;
/// \brief Copy assignment operator
- Field& operator=(const Field& other) = default;
+ Field& operator=(const Field& other)
+ {
+ if (this != &other) {
+ fPrimVertexField = other.fPrimVertexField;
+ foFldIntrpl = other.foFldIntrpl;
+ foFldOrig = other.foFldOrig;
+ fFieldType = other.fFieldType;
+ fFieldMode = other.fFieldMode;
+ }
+ return *this;
+ }
+
+ /// \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 fFieldMode == EFieldMode::Intrpl ? foFldIntrpl->GetFieldValue(sliceID, x, y)
+ : foFldOrig->GetFieldValue(sliceID, x, y);
+ }
/// \brief Gets field type
EFieldType GetFieldType() const { return this->fFieldType; }
/// \brief Gets field region near primary vertex
- const FieldRegion_t& GetPrimVertexField() const { return fPrimVertexField; }
-
- /// \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; }
+ const FieldRegion<T>& GetPrimVertexField() const { return fPrimVertexField; }
+
+ /// \brief Makes 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]
+ /// \note Parameters b0-b2, z0-z2 are ignored, if fFieldMode == EFieldMode::Orig
+ FieldRegion<T> GetFieldRegion(const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1, const T& z1,
+ const FieldValue<T>& b2, const T& z2) const
+ {
+ return (fFieldMode == EFieldMode::Intrpl ? FieldRegion<T>(fFieldType, b0, z0, b1, z1, b2, z2)
+ : FieldRegion<T>(fFieldType, foFldOrig->GetFieldFunction()));
+ }
/// \brief String representation of the class
/// \param indentLevel Indent level of the string output
@@ -207,16 +278,49 @@ namespace cbm::algo::kf
std::string ToString(int indentLevel, int verbose) const;
private:
- /// \brief Serialization function
- friend class boost::serialization::access;
+ /// \brief Default constructor
+ Field() = default;
+
+ /// \brief Serialization load method
+ template<class Archive>
+ void load(Archive& ar, const unsigned int /*version*/)
+ {
+ auto field = detail::FieldBase<T, EFieldMode::Intrpl>{};
+ ar >> field;
+ foFldIntrpl = std::make_optional(field);
+ ar >> fPrimVertexField;
+ ar >> fFieldType;
+ ar >> fFieldMode;
+ foFldOrig = std::nullopt; // Note: original field cannot be serialized
+ }
+
+ /// \brief Serialization save method
template<class Archive>
- void serialize(Archive& ar, const unsigned int)
+ void save(Archive& ar, const unsigned int /*version*/) const
{
- ar& boost::serialization::base_object<detail::FieldBase<T, FldMode>>(*this);
- ar& fPrimVertexField;
+ if (fFieldMode == EFieldMode::Intrpl) {
+ ar << (*foFldIntrpl);
+ ar << fPrimVertexField;
+ ar << fFieldType;
+ ar << fFieldMode;
+ }
+ else if (fFieldMode == EFieldMode::Orig) {
+ throw std::logic_error("Attempt to serialize a kf::FieldRegion object with fFieldMode == EFieldMode::Orig");
+ }
}
- FieldRegion_t fPrimVertexField; ///< Field region near primary vertex (constant field is assumed)
+ /// \brief Serialization method
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int version)
+ {
+ boost::serialization::split_member(ar, *this, version);
+ }
+
+ std::optional<detail::FieldBase<T, EFieldMode::Intrpl>> foFldIntrpl{std::nullopt}; ///< Interpolated field
+ std::optional<detail::FieldBase<T, EFieldMode::Orig>> foFldOrig{std::nullopt}; ///< Original field
+ FieldRegion<T> fPrimVertexField; ///< Field region near primary vertex (constant field is assumed)
+ EFieldType fFieldType{EFieldType::Normal}; ///< Field type
+ EFieldMode fFieldMode; ///< Field mode
};
@@ -259,25 +363,34 @@ namespace cbm::algo::kf
/// \param refZ Reference z-position of the slice [cm]
void AddSliceReference(double halfSizeX, double halfSizeY, double refZ);
+ /// \brief Gets field mode
+ EFieldMode GetFieldMode() const { return fFieldMode; }
+
/// \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 Create field
+ /// \tparam T Underlying floating-point data type
+ template<typename T>
+ Field<T> 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; }
+ void SetFieldFunction(const FieldFn_t& fieldFn, EFieldType fldType)
+ {
+ fFieldFn = fieldFn;
+ fFieldType = fldType;
+ }
- /// \brief Sets field type
- /// \param fieldType Field type
- void SetFieldType(EFieldType fieldType) { fFieldType = fieldType; }
+ /// \brief Sets field mode
+ void SetFieldMode(EFieldMode fldMode) { fFieldMode = fldMode; }
+
+ /// \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; }
/// \brief Sets target
/// \param x x-coordinate of the target position [cm]
@@ -285,15 +398,12 @@ namespace cbm::algo::kf
/// \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
+ EFieldMode fFieldMode{EFieldMode::Intrpl}; ///< FieldMode
/// \brief Target position
std::array<double, 3> fTarget{{defs::Undef<double>, defs::Undef<double>, defs::Undef<double>}};
@@ -301,10 +411,10 @@ namespace cbm::algo::kf
// -------------------------------------------------------------------------------------------------------------------
//
- template<typename T, EFieldMode FldMode>
- Field<T, FldMode> FieldFactory::MakeField() const
+ template<typename T>
+ Field<T> FieldFactory::MakeField() const
{
- auto field = Field<T, FldMode>{};
+ auto field = Field<T>(fFieldMode, fFieldType);
// Check initialization
if (std::any_of(fTarget.begin(), fTarget.end(), [](double x) { return !utils::IsFinite(x); })) {
@@ -318,14 +428,17 @@ namespace cbm::algo::kf
}
// Initialize the Field object
- field.SetFieldType(fFieldType);
- if constexpr (FldMode == EFieldMode::Orig) {
- field.SetFieldFunction(fFieldFn);
- field.SetPrimVertexField(field.MakeFieldRegion());
+ if (fFieldMode == EFieldMode::Orig) {
+ field.foFldOrig->SetFieldFunction(fFieldFn);
+ field.fPrimVertexField = FieldRegion<T>(fFieldType, fFieldFn);
+ for (const auto& sliceRef : fSliceReferences) {
+ field.foFldOrig->AddFieldSlice(sliceRef.fRefZ);
+ }
}
- else if constexpr (FldMode == EFieldMode::Intrpl) { // A version with the interpolated field
+ else if (fFieldMode == EFieldMode::Intrpl) { // A version with the interpolated field
for (const auto& sliceRef : fSliceReferences) {
- field.AddFieldSlice(FieldSlice<T>(fFieldFn, sliceRef.fHalfSizeX, sliceRef.fHalfSizeY, sliceRef.fRefZ));
+ field.foFldIntrpl->AddFieldSlice(
+ FieldSlice<T>(fFieldFn, sliceRef.fHalfSizeX, sliceRef.fHalfSizeY, sliceRef.fRefZ));
}
{
// PV field initialization
@@ -338,7 +451,8 @@ namespace cbm::algo::kf
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]));
+ field.fPrimVertexField = FieldRegion<T>(fFieldType, fld[0], pos[0], fld[1], pos[1], fld[2], pos[2]);
+ // TODO: Provide alternative method for Orig
}
}
return field;
diff --git a/algo/kf/core/geo/KfFieldRegion.cxx b/algo/kf/core/geo/KfFieldRegion.cxx
index fc892352463994563cd665a5adca6256732a1b8d..1f824eef447606d404b5c0d695d398520b222c81 100644
--- a/algo/kf/core/geo/KfFieldRegion.cxx
+++ b/algo/kf/core/geo/KfFieldRegion.cxx
@@ -55,6 +55,20 @@ 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) {
+ 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])};
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T>
@@ -70,66 +84,94 @@ void FieldRegionBase<T, EFieldMode::Intrpl>::Shift(const T& z)
fZfirst = z;
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+std::string FieldRegionBase<T, EFieldMode::Intrpl>::ToString(int indentLevel, int) const
+{
+ constexpr char IndentChar = '\t';
+ std::stringstream msg;
+ std::string indent(indentLevel, IndentChar);
+ using fmt::format;
+ auto Cnv = [&](const auto& val) { return utils::simd::Cast<T, Literal_t<T>>(val); }; // alias for the conversion fn
+ const auto& coeff = this->fCoeff;
+ msg << indent << format("Field Region: dz = z - ({})", Cnv(this->fZfirst));
+ msg << '\n'
+ << indent << IndentChar
+ << format("Bx(dz) = {:>12} + {:>12} dz + {:>12} dz2", Cnv(coeff[0][0]), Cnv(coeff[0][1]), Cnv(coeff[0][2]));
+ msg << '\n'
+ << indent << IndentChar
+ << format("Bx(dz) = {:>12} + {:>12} dz + {:>12} dz2", Cnv(coeff[1][0]), Cnv(coeff[1][1]), Cnv(coeff[1][2]));
+ 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";
+ }
+ return msg.str();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+std::string FieldRegionBase<T, EFieldMode::Orig>::ToString(int indentLevel, int) const
+{
+ constexpr char IndentChar = '\t';
+ std::stringstream msg;
+ std::string indent(indentLevel, IndentChar);
+ msg << indent << "Field region: created from the original field function";
+ return msg.str();
+}
+
+
namespace cbm::algo::kf::detail
{
template class FieldRegionBase<float, EFieldMode::Intrpl>;
template class FieldRegionBase<double, EFieldMode::Intrpl>;
template class FieldRegionBase<fvec, EFieldMode::Intrpl>;
+ template class FieldRegionBase<float, EFieldMode::Orig>;
+ template class FieldRegionBase<double, EFieldMode::Orig>;
+ template class FieldRegionBase<fvec, EFieldMode::Orig>;
} // namespace cbm::algo::kf::detail
+
// ---------------------------------------------------------------------------------------------------------------------
//
-template<typename T, EFieldMode FldMode>
-FieldValue<T> FieldRegion<T, FldMode>::Get(const T& x, const T& y, const T& z) const
+template<typename T>
+void FieldRegion<T>::Set(const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1, const T& z1,
+ const FieldValue<T>& b2, const T& z2)
{
- if 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])};
+ if (fFieldMode == EFieldMode::Intrpl) {
+ foFldIntrpl->Set(b0, z0, b1, z1, b2, z2);
}
- else if constexpr (FldMode == EFieldMode::Orig) {
- return GlobalField::GetFieldValue(this->fFieldFn, x, y, z);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename T>
+void FieldRegion<T>::Shift(const T& z)
+{
+ if (fFieldMode == EFieldMode::Intrpl) {
+ foFldIntrpl->Shift(z);
}
}
// ---------------------------------------------------------------------------------------------------------------------
//
-template<typename T, EFieldMode FldMode>
-std::string FieldRegion<T, FldMode>::ToString(int indentLevel, int) const
+template<typename T>
+std::string FieldRegion<T>::ToString(int indentLevel, int) const
{
constexpr char IndentChar = '\t';
std::stringstream msg;
std::string indent(indentLevel, IndentChar);
- if constexpr (FldMode == EFieldMode::Intrpl) {
- using fmt::format;
- auto Cnv = [&](const auto& val) { return utils::simd::Cast<T, Literal_t<T>>(val); }; // alias for the conversion fn
- const auto& coeff = this->fCoeff;
- msg << indent << format("Field Region: dz = z - ({})", Cnv(this->fZfirst));
- msg << '\n'
- << indent << IndentChar
- << format("Bx(dz) = {:>12} + {:>12} dz + {:>12} dz2", Cnv(coeff[0][0]), Cnv(coeff[0][1]), Cnv(coeff[0][2]));
- msg << '\n'
- << indent << IndentChar
- << format("Bx(dz) = {:>12} + {:>12} dz + {:>12} dz2", Cnv(coeff[1][0]), Cnv(coeff[1][1]), Cnv(coeff[1][2]));
- 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";
- }
- }
- else if constexpr (FldMode == EFieldMode::Orig) {
- msg << indent << "Field region: created from the original field function";
- }
+ msg << indent << "FieldType: " << static_cast<int>(fFieldType) << '\n';
+ msg << indent << "FieldMode: " << static_cast<int>(fFieldMode) << '\n';
+ msg << indent
+ << (fFieldMode == EFieldMode::Intrpl ? foFldIntrpl->ToString(indentLevel) : foFldOrig->ToString(indentLevel));
return msg.str();
}
@@ -147,10 +189,7 @@ std::tuple<double, double, double> GlobalField::UndefField(double, double, doubl
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>;
+ template class FieldRegion<float>;
+ template class FieldRegion<double>;
+ template class FieldRegion<fvec>;
} // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfFieldRegion.h b/algo/kf/core/geo/KfFieldRegion.h
index 9e8ceec5e44f2744a66fc03b480992f8c8a52ceb..e31f1ce5af3655dba7f2c0bd1611bfbc3bd07b97 100644
--- a/algo/kf/core/geo/KfFieldRegion.h
+++ b/algo/kf/core/geo/KfFieldRegion.h
@@ -17,9 +17,11 @@
#include <boost/serialization/access.hpp>
#include <boost/serialization/array.hpp>
-#include <boost/serialization/base_object.hpp>
+#include <boost/serialization/optional.hpp>
+#include <boost/serialization/split_free.hpp>
#include <array>
+#include <optional>
#include <string>
#include <type_traits>
@@ -63,7 +65,6 @@ namespace cbm::algo::kf
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)
@@ -73,8 +74,6 @@ namespace cbm::algo::kf
friend class FieldRegionBase;
public:
- using FieldValue_t = FieldValue<T>;
-
/// \brief Default constructor
FieldRegionBase() = default;
@@ -94,16 +93,23 @@ namespace cbm::algo::kf
/// \brief Copy assignment operator
FieldRegionBase& operator=(const FieldRegionBase& other) = default;
- 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)
+ /// \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
{
+ return GlobalField::GetFieldValue(fFieldFn, x, y, z);
}
+
+ /// \brief String representation of the class content
+ /// \param intdentLevel Indent level
+ /// \param verbose Verbosity level
+ std::string ToString(int indentLevel = 0, int verbose = 1) const;
+
+ protected:
+ FieldFn_t fFieldFn{defs::ZeroFieldFn}; ///< Field function: (x,y,z) [cm] -> (Bx,By,Bz) [kG]
};
@@ -154,6 +160,12 @@ namespace cbm::algo::kf
/// \brief Copy assignment operator
FieldRegionBase& operator=(const FieldRegionBase& other) = default;
+ /// \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]
+ FieldValue<T> Get(const T& x, const T& y, const T& z) const;
+
/// \brief Gets the first z-coordinate of the interpolation
const T& GetZfirst() const { return fZfirst; }
@@ -172,6 +184,11 @@ namespace cbm::algo::kf
/// \param z new z-coordinate [cm]
void Shift(const T& z);
+ /// \brief String representation of the class content
+ /// \param intdentLevel Indent level
+ /// \param verbose Verbosity level
+ std::string ToString(int indentLevel = 0, int verbose = 1) const;
+
protected:
CoeffArray_t fCoeff{{T(0.)}}; ///< Approximation coefficients
T fZfirst{Literal_t<T>()}; ///< Reference z-coordinate of the field region
@@ -192,35 +209,56 @@ namespace cbm::algo::kf
/// \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> {
- template<typename, EFieldMode>
+ template<typename T>
+ class alignas(VcMemAlign) FieldRegion {
+ template<typename>
friend class FieldRegion;
public:
- static constexpr EFieldMode FieldType = FldMode;
-
- using FieldValue_t = FieldValue<T>;
- using Real_t = T;
-
- /// \brief Default constructor
- FieldRegion() = default;
+ /// \brief Constructor of the generic field
+ FieldRegion(EFieldMode fldMode, EFieldType fldType)
+ : foFldIntrpl(fldMode == EFieldMode::Intrpl ? std::make_optional<detail::FieldRegionBase<T, EFieldMode::Intrpl>>()
+ : std::nullopt)
+ , foFldOrig(fldMode == EFieldMode::Orig ? std::make_optional<detail::FieldRegionBase<T, EFieldMode::Orig>>()
+ : std::nullopt)
+ , fFieldType(fldType)
+ , fFieldMode(fldMode)
+ {
+ }
/// \brief Copy constructor
/// \tparam I Underlying floating point type of the source
template<typename I>
- FieldRegion(const FieldRegion<I, FldMode>& other) : detail::FieldRegionBase<T, FldMode>(other)
+ FieldRegion(const FieldRegion<I>& other)
+ : foFldIntrpl(other.foFldIntrpl.has_value()
+ ? std::make_optional(detail::FieldRegionBase<T, EFieldMode::Intrpl>(*other.foFldIntrpl))
+ : std::nullopt)
+ , foFldOrig(other.foFldOrig.has_value()
+ ? std::make_optional(detail::FieldRegionBase<T, EFieldMode::Orig>(*other.foFldOrig))
+ : std::nullopt)
+ , fFieldType(other.fFieldType)
+ , fFieldMode(other.fFieldMode)
{
}
- /// \brief Constructor from the arguments
+ /// \brief Constructor for the interpolated field region
/// \param fieldType Type of the field
- /// \param args Arguments, which have to be passed to the base class constructors
- template<typename... Args>
- FieldRegion(EFieldType fieldType, Args... args)
- : detail::FieldRegionBase<T, FldMode>(args...)
+ FieldRegion(EFieldType fieldType, const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1, const T& z1,
+ const FieldValue<T>& b2, const T& z2)
+ : foFldIntrpl(std::make_optional<detail::FieldRegionBase<T, EFieldMode::Intrpl>>(b0, z0, b1, z1, b2, z2))
+ , foFldOrig(std::nullopt)
, fFieldType(fieldType)
+ , fFieldMode(EFieldMode::Intrpl)
+ {
+ }
+
+ /// \brief Constructor for the field region with the original field function
+ /// \param fieldType Type of the field
+ FieldRegion(EFieldType fieldType, const FieldFn_t& fieldFn)
+ : foFldIntrpl(std::nullopt)
+ , foFldOrig(std::make_optional<detail::FieldRegionBase<T, EFieldMode::Orig>>(fieldFn))
+ , fFieldType(fieldType)
+ , fFieldMode(EFieldMode::Orig)
{
}
@@ -235,28 +273,83 @@ namespace cbm::algo::kf
/// \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;
+ FieldValue<T> Get(const T& x, const T& y, const T& z) const
+ {
+ return fFieldMode == EFieldMode::Intrpl ? foFldIntrpl->Get(x, y, z) : foFldOrig->Get(x, y, z);
+ }
+
+ /// \brief Gets field mode
+ EFieldMode GetFieldMode() const { return fFieldMode; }
/// \brief Gets the field type
EFieldType GetFieldType() const { return fFieldType; }
+ /// \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 Sets the field type
+ void SetFieldType(EFieldType fldType) { fFieldType = fldType; }
+
+ /// \brief Shifts the coefficients to another first z-coordinate
+ /// \param z new z-coordinate [cm]
+ void Shift(const T& z);
+
/// \brief String representation of the class content
/// \param intdentLevel Indent level
/// \param verbose Verbosity level
std::string ToString(int indentLevel = 0, int verbose = 1) const;
private:
- /// \brief Serialization method
friend class boost::serialization::access;
+
+ /// \brief Default constructor
+ FieldRegion() = default;
+
+ /// \brief Serialization load method
template<class Archive>
- void serialize(Archive& ar, const unsigned int)
+ void load(Archive& ar, const unsigned int /*version*/)
{
- using boost::serialization::base_object;
- ar& base_object<detail::FieldRegionBase<T, FldMode>>(*this);
- ar& fFieldType;
+ auto field = detail::FieldRegionBase<T, EFieldMode::Intrpl>{};
+ ar >> field;
+ foFldIntrpl = std::make_optional(field);
+ ar >> fFieldType;
+ ar >> fFieldMode;
+ foFldOrig = std::nullopt; // Note: original field cannot be serialized
}
+ /// \brief Serialization save method
+ template<class Archive>
+ void save(Archive& ar, const unsigned int /*version*/) const
+ {
+ if (fFieldMode == EFieldMode::Intrpl) {
+ ar << (*foFldIntrpl);
+ ar << fFieldType;
+ ar << fFieldMode;
+ }
+ else if (fFieldMode == EFieldMode::Orig) {
+ throw std::logic_error("Attempt to serialize a kf::FieldRegion object with fFieldMode == EFieldMode::Orig");
+ }
+ }
+
+ /// \brief Serialization method
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int version)
+ {
+ boost::serialization::split_member(ar, *this, version);
+ }
+
+ std::optional<detail::FieldRegionBase<T, EFieldMode::Intrpl>> foFldIntrpl{std::nullopt};
+ std::optional<detail::FieldRegionBase<T, EFieldMode::Orig>> foFldOrig{std::nullopt};
EFieldType fFieldType{EFieldType::Null}; ///< Field type in a given region
+ EFieldMode fFieldMode; ///< Field mode
};
@@ -265,8 +358,8 @@ namespace cbm::algo::kf
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;
+ using utils::simd::Cast;
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);
diff --git a/algo/kf/core/geo/KfFieldSlice.h b/algo/kf/core/geo/KfFieldSlice.h
index 295051c2423da08d86d1fe574f0b18f876ea65bc..6be46bad7ce19241551da0f348dea12fecc890e7 100644
--- a/algo/kf/core/geo/KfFieldSlice.h
+++ b/algo/kf/core/geo/KfFieldSlice.h
@@ -90,7 +90,7 @@ namespace cbm::algo::kf
}
/// \brief Gets reference z-coordinate of the slice [cm]
- const T& GetZref() { return fZref; }
+ const T& GetZref() const { return fZref; }
/// \brief String representation of the class content
/// \param intdentLevel Indent level
diff --git a/algo/kf/core/geo/KfSetup.cxx b/algo/kf/core/geo/KfSetup.cxx
index ec8862775685506164e4325da359e7d930ba9096..dd4fa54a174614b424df709588116575d3458434 100644
--- a/algo/kf/core/geo/KfSetup.cxx
+++ b/algo/kf/core/geo/KfSetup.cxx
@@ -31,13 +31,7 @@ std::string Setup<T>::ToString(int verbosity, int indentLevel) const
for (const auto& layer : fvMaterialLayers) {
msg << layer.ToString(indentLevel + 1, verbosity) << '\n';
}
- msg << indent << "\nINTERPOLATED MAGNETIC FIELD:\n" << fIntrplField.ToString(indentLevel + 1, verbosity);
- if (foOrigField) {
- msg << indent << "\nORIGINAL MAGNETIC FIELD:\n" << foOrigField->ToString(indentLevel + 1, verbosity);
- }
- else {
- msg << indent << "\nORIGINAL MANGETIC FIELD: not provided";
- }
+ msg << indent << "\nMAGNETIC FIELD:\n" << fField.ToString(indentLevel + 1, verbosity);
}
return msg.str();
}
diff --git a/algo/kf/core/geo/KfSetup.h b/algo/kf/core/geo/KfSetup.h
index 75829aa2e64b4f53d9b0a2938a60d1a9c175e7c4..2f1a8a3d846024352d9ca5b04fe495ba5554c17f 100644
--- a/algo/kf/core/geo/KfSetup.h
+++ b/algo/kf/core/geo/KfSetup.h
@@ -37,8 +37,14 @@ namespace cbm::algo::kf
public:
- /// \brief Default constructor
- Setup() = default;
+ /// \brief Constructor
+ /// \param fldMode Field mode
+ Setup(EFieldMode fldMode)
+ : fvMaterialLayers({})
+ , fField(Field<T>(fldMode, EFieldType::Normal))
+ , fTarget(Target<T>())
+ {
+ }
/// \brief Destructor
~Setup() = default;
@@ -48,8 +54,7 @@ namespace cbm::algo::kf
template<typename I>
Setup(const Setup<I>& other)
: fvMaterialLayers(other.fvMaterialLayers)
- , fIntrplField(other.fIntrplField)
- , foOrigField(other.foOrigField.has_value() ? std::make_optional(*other.foOrigField) : std::nullopt)
+ , fField(other.fField)
, fTarget(other.fTarget)
{
}
@@ -68,10 +73,8 @@ namespace cbm::algo::kf
void AddMaterial(const MaterialMap& map) { fvMaterialLayers.emplace_back(map); }
/// \brief Makes an instance of the field depending on the template parameter
- /// \tparam FldMode Field mode of the returned instance
/// \throw std::logic_error If the particular field member is undefined (nullopt)
- template<EFieldMode FldMode>
- const Field<T, FldMode>& GetField() const;
+ const Field<T>& GetField() const { return fField; }
/// \brief Gets material layer
/// \param iLayer Index of layer
@@ -89,8 +92,7 @@ namespace cbm::algo::kf
std::string ToString(int verbosity = 1, int indentLevel = 0) const;
/// \brief Sets field
- template<EFieldMode FldMode>
- void SetField(const Field<double, FldMode>& field);
+ void SetField(const Field<double>& field) { fField = Field<T>(field); }
/// \brief Sets target (moves from the source)
/// \param target Target instance
@@ -102,22 +104,16 @@ namespace cbm::algo::kf
void load(Archive& ar, const unsigned int /*version*/)
{
ar >> fvMaterialLayers;
- ar >> fvMaterialLayerIds;
ar >> fTarget;
- ar >> fIntrplField;
- foOrigField = std::nullopt;
-
- // NOTE: The original field cannot be serialized, because it contains a pointer to the external magnetic field
- // function.
+ ar >> fField;
}
template<class Archive>
void save(Archive& ar, const unsigned int /*version*/) const
{
ar << fvMaterialLayers;
- ar << fvMaterialLayerIds;
ar << fTarget;
- ar << fIntrplField;
+ ar << fField;
}
template<class Archive>
@@ -127,45 +123,8 @@ namespace cbm::algo::kf
}
std::vector<MaterialMap> fvMaterialLayers = {}; ///< Container of the material maps
- std::vector<int> fvMaterialLayerIds = {}; ///< Map of the layer ids: original -> actual
- Field<T, EFieldMode::Intrpl> fIntrplField; ///< Interpolated field (NOTE: maybe make optional)
- std::optional<Field<T, EFieldMode::Orig>> foOrigField{std::nullopt}; ///< Optional original field
- Target<T> fTarget; ///< Target layer
+ Field<T> fField; ///< Interpolated field (NOTE: maybe make optional)
+ Target<T> fTarget; ///< Target layer
};
- // -------------------------------------------------------------------------------------------------------------------
- //
- template<typename T>
- template<EFieldMode FldMode>
- const Field<T, FldMode>& Setup<T>::GetField() const
- {
- if constexpr (FldMode == EFieldMode::Orig) {
- if (foOrigField) {
- return *foOrigField;
- }
- else {
- throw std::logic_error("kf::Setup::MakeField(): attempt to create a Field instance in the original field mode "
- "in the unsuitable scenario (for example, in the online reconstruction code). Please, "
- "use the interpolated magnetic field calling Setup::MakeField<T, EFieldMode::Intrpl>() "
- "instead.");
- }
- }
- else if constexpr (FldMode == EFieldMode::Intrpl) {
- return fIntrplField;
- }
- }
-
- // -------------------------------------------------------------------------------------------------------------------
- //
- template<typename T>
- template<EFieldMode FldMode>
- void Setup<T>::SetField(const Field<double, FldMode>& field)
- {
- if constexpr (FldMode == EFieldMode::Intrpl) {
- fIntrplField = field;
- }
- else if constexpr (FldMode == EFieldMode::Orig) {
- foOrigField = std::make_optional(field);
- }
- }
} // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfSetupInitializer.h b/algo/kf/core/geo/KfSetupInitializer.h
index 7fb84e9f6bd413d0da7eafe3c2393eabf79dff46..1f2f51e51c2f01ed03a7fb3a90da3a9778c94a61 100644
--- a/algo/kf/core/geo/KfSetupInitializer.h
+++ b/algo/kf/core/geo/KfSetupInitializer.h
@@ -65,20 +65,20 @@ namespace cbm::algo::kf
virtual bool InitProperties() { return false; }
/// \brief Creates a setup instance
- /// \param bProvideOrigField Flag: true - defines optional variable of the original field
template<typename T>
- Setup<T> MakeSetup(bool bProvideOrigField = true) const;
+ Setup<T> MakeSetup() const;
/// \brief Resets the instance
void Reset();
/// \brief Sets magnetic field function
+ /// \param fieldMode Magnetic field mode:
/// \param fieldFn Magnetic field function
/// \param fieldType Magnetic field type
- void SetFieldFunction(const FieldFn_t& fieldFn, EFieldType fieldType)
+ void SetFieldProperty(EFieldMode fieldMode, const FieldFn_t& fieldFn, EFieldType fieldType)
{
- fFieldFactory.SetFieldFunction(fieldFn);
- fFieldFactory.SetFieldType(fieldType);
+ fFieldFactory.SetFieldFunction(fieldFn, fieldType);
+ fFieldFactory.SetFieldMode(fieldMode);
}
/// \brief Sets target initialization properties
@@ -120,17 +120,14 @@ namespace cbm::algo::kf
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T>
- Setup<T> SetupInitializer::MakeSetup(bool bProvideOrigField) const
+ Setup<T> SetupInitializer::MakeSetup() const
{
- Setup<T> setup;
+ Setup<T> setup(fFieldFactory.GetFieldMode());
// Target initialization
setup.SetTarget(fTarget);
// Field initialization
- setup.SetField(fFieldFactory.MakeField<double, EFieldMode::Intrpl>());
- if (bProvideOrigField) {
- setup.SetField(fFieldFactory.MakeField<double, EFieldMode::Orig>());
- }
+ setup.SetField(fFieldFactory.MakeField<double>());
// Material layers initialization
for (const auto& material : fMaterialLayers) {
@@ -145,7 +142,7 @@ namespace cbm::algo::kf
template<typename T>
Setup<T> SetupInitializer::Load(const std::string& fileName)
{
- Setup<double> setup;
+ Setup<double> setup(EFieldMode::Intrpl);
std::ifstream ifs(fileName, std::ios::binary);
if (!ifs) {
std::stringstream msg;
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index d8e95834592e437f43167b4676b04ee25200e423..09568aa6a63e73ceb9afcd13a15e11dc4c1dc795 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -28,6 +28,7 @@
#include "CbmStsTrackingInterface.h"
#include "CbmTofTrackingInterface.h"
#include "CbmTrdTrackingInterface.h"
+#include "KfDefs.h"
#include <boost/filesystem.hpp>
// TODO: include of CbmSetup.h creates problems on Mac
@@ -544,7 +545,7 @@ try {
auto pSetupInitializer = std::make_unique<TrackingSetupInitializer>();
pSetupInitializer->Use(fUseMVD, fUseSTS, fUseMUCH, fUseTRD, fUseTOF);
pSetupInitializer->InitProperties();
- auto trackerSetup = pSetupInitializer->MakeSetup<double>(/*provide original field*/ true);
+ auto trackerSetup = pSetupInitializer->MakeSetup<double>();
LOG(info) << "!!!!\n!!!!\n!!!!\n!!!! KF SETUP:\n" << trackerSetup.ToString(0) << "\n!!!!\n!!!!\n!!!!";
// Storing setup
@@ -558,12 +559,13 @@ try {
{
// TEST: magnetic field
- const auto& fld = trackerSetupVec.GetField<cbm::algo::kf::EFieldMode::Intrpl>();
+ const auto& fld = trackerSetupVec.GetField();
struct Node {
int iSt;
fvec x;
fvec y;
};
+
{
Node n{.iSt = 1, .x = 2.3, .y = 1.2};
auto kfFldValue = fld.GetFieldValue(n.iSt, n.x, n.y);
@@ -583,6 +585,27 @@ try {
LOG(info) << "B: " << kfFldValue.ToString() << " --- " << caFldValue.ToString();
}
+ // TEST: field region building
+ {
+ using cbm::algo::kf::EFieldType;
+ fvec x = 1.2;
+ fvec y = 1.3;
+ fvec z0 = -32.;
+ fvec z1 = -28.;
+ fvec z2 = -24.;
+ auto kfFldVal_0 = fld.GetFieldValue(0, x, y);
+ auto kfFldVal_1 = fld.GetFieldValue(1, x, y);
+ auto kfFldVal_2 = fld.GetFieldValue(2, x, y);
+ auto kfFldReg = fld.GetFieldRegion(kfFldVal_0, z0, kfFldVal_1, z1, kfFldVal_2, z2);
+ auto caFldVal_0 = fpParameters->GetStation(0).fieldSlice.GetFieldValue(x, y);
+ auto caFldVal_1 = fpParameters->GetStation(1).fieldSlice.GetFieldValue(x, y);
+ auto caFldVal_2 = fpParameters->GetStation(2).fieldSlice.GetFieldValue(x, y);
+ auto caFldReg = cbm::algo::ca::FieldRegion<fvec>{};
+ caFldReg.Set(caFldVal_0, z0, caFldVal_1, z1, caFldVal_2, z2);
+ LOG(info) << "CA ----> " << caFldReg.ToString();
+ LOG(info) << "KF ----> " << kfFldReg.ToString();
+ }
+
// TEST: material maps
{
Node n{.iSt = 1, .x = 2.3, .y = 1.2};
diff --git a/reco/kfnew/CbmKfTrackingSetupInitializer.cxx b/reco/kfnew/CbmKfTrackingSetupInitializer.cxx
index 9650ac15884b69c2976eda0f5f5553d516acaa75..b93a1306112883c814397a94ca6e3ba81342942e 100644
--- a/reco/kfnew/CbmKfTrackingSetupInitializer.cxx
+++ b/reco/kfnew/CbmKfTrackingSetupInitializer.cxx
@@ -32,6 +32,7 @@ using kf::tools::MaterialMapCreator;
//
bool TrackingSetupInitializer::InitProperties()
try {
+ using cbm::algo::kf::EFieldMode;
using cbm::algo::kf::EFieldType;
// Check, if a subsystem exists in the setup
@@ -43,11 +44,10 @@ try {
// Magnetic field initialization
if (FairRunAna::Instance()->GetField()) {
- this->SetFieldFunction(OriginalField(), EFieldType::Normal);
+ this->SetFieldProperty(EFieldMode::Intrpl, OriginalField(), EFieldType::Normal);
}
else {
- //this->SetFieldFunction(cbm::algo::kf::defs::ZeroFieldFn, EFieldType::Null);
- this->SetFieldFunction(ZeroField(), EFieldType::Null);
+ this->SetFieldProperty(EFieldMode::Intrpl, ZeroField(), EFieldType::Null);
}
// Tracking station property initialization