diff --git a/algo/ca/core/pars/CaField.cxx b/algo/ca/core/pars/CaField.cxx
index fe9b2d71cf7daeef5933870abdc40a89a0d6233d..a4b24f596cb9f1c6064016a17ba08477457fe360 100644
--- a/algo/ca/core/pars/CaField.cxx
+++ b/algo/ca/core/pars/CaField.cxx
@@ -44,16 +44,10 @@ void FieldValue<DataT>::CheckConsistency() const
//----------------------------------------------------------------------------------------------------------------------
// TODO:
template<typename DataT>
-std::string FieldValue<DataT>::ToString(int indentLevel) const
+std::string FieldValue<DataT>::ToString(int) const
{
std::stringstream aStream{};
- constexpr char indentChar = '\t';
- std::string indent(indentLevel, indentChar);
- aStream << indent << "Bx [kG]: " << std::setw(12) << std::setfill(' ') << kfutils::simd::Cast<DataT, float>(x)
- << '\n';
- aStream << indent << "By [kG]: " << std::setw(12) << std::setfill(' ') << kfutils::simd::Cast<DataT, float>(y)
- << '\n';
- aStream << indent << "Bz [kG]: " << std::setw(12) << std::setfill(' ') << kfutils::simd::Cast<DataT, float>(z);
+ aStream << "B = (" << x << ", " << y << ", " << z << ") [kG]";
return aStream.str();
}
diff --git a/algo/ca/core/pars/CaField.h b/algo/ca/core/pars/CaField.h
index 3921f298b8162f4f8c7c3b017265111c647e1d72..5aaeadbe9fe113961d39aa956bdd49a4ffd6067f 100644
--- a/algo/ca/core/pars/CaField.h
+++ b/algo/ca/core/pars/CaField.h
@@ -68,7 +68,7 @@ namespace cbm::algo::ca
/// String representation of class contents
/// \param indentLevel number of indent characters in the output
- std::string ToString(int indentLevel) const;
+ 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
diff --git a/algo/kf/core/CMakeLists.txt b/algo/kf/core/CMakeLists.txt
index 9247f6cf8c976e8da792e8ba2f333e5cb081e242..f412291aba2ff71cfc6c854e0ae73dc39a6dd4ba 100644
--- a/algo/kf/core/CMakeLists.txt
+++ b/algo/kf/core/CMakeLists.txt
@@ -17,6 +17,7 @@ set(SRCS
${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}/geo/KfSetupInitializer.cxx
${CMAKE_CURRENT_SOURCE_DIR}/utils/KfUtils.cxx
)
@@ -45,11 +46,12 @@ target_include_directories(KfCore
target_compile_definitions(KfCore PUBLIC NO_ROOT)
target_link_libraries(KfCore
- PUBLIC Vc::Vc
- Boost::serialization
+ PUBLIC Vc::Vc
OnlineDataLog # needed for the logger
external::fles_logging # needed for the logger
- external::fles_ipc # needed for the logger
+ external::fles_ipc # needed for the logger
+ PRIVATE Boost::serialization
+ fmt::fmt
external::yaml-cpp
)
@@ -72,6 +74,7 @@ install(
geo/KfFieldValue.h
geo/KfMaterialMap.h
geo/KfSetup.h
+ geo/KfSetupInitializer.h
geo/KfTarget.h
geo/KfField.h
pars/KfParticlePDG.h
diff --git a/algo/kf/core/KfDefs.h b/algo/kf/core/KfDefs.h
index abed6a7468806ba882536440240e55162b3d3c4b..7a8a795b29da3eb98cac677b7d538b838f99c7f9 100644
--- a/algo/kf/core/KfDefs.h
+++ b/algo/kf/core/KfDefs.h
@@ -64,18 +64,13 @@ namespace cbm::algo::kf
/// \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
- constexpr int MaxNofMaterialLayers = 32; ///< Max number of the material layers
+ //constexpr int MaxNofFieldSlices = 64; ///< Max number of field slices
+ //constexpr int MaxNofMaterialLayers = 32; ///< Max number of the material layers
// ----- Control -----------------------------------------------------------------------------------------------------
constexpr int DebugLvl = 0; ///< Level of debug output
diff --git a/algo/kf/core/geo/KfField.cxx b/algo/kf/core/geo/KfField.cxx
index fc5ae8cd980673c5d8bc84e6329ce85d290172dc..888971655eea6e73566e0a7f148e390000cdd79c 100644
--- a/algo/kf/core/geo/KfField.cxx
+++ b/algo/kf/core/geo/KfField.cxx
@@ -16,45 +16,41 @@
using cbm::algo::kf::EFieldMode;
using cbm::algo::kf::Field;
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, EFieldMode FldMode>
-std::string Field<T, FldMode>::ToString(int indentLevel) const
+std::string Field<T, FldMode>::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) << '\n';
+ 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 << "Original field function";
}
else if constexpr (FldMode == EFieldMode::Intrpl) {
msg << indent << "Field slices:";
for (const auto& fldSlice : this->fvFieldSlices) {
- msg << indent << "\n - " << fldSlice.ToString(indentLevel + 1);
+ msg << indent << "\n - " << fldSlice.ToString(indentLevel + 1, 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>;
+} // namespace cbm::algo::kf
+
+
// ---------------------------------------------------------------------------------------------------------------------
//
void FieldFactory::AddSliceReference(double halfSizeX, double halfSizeY, double zRef)
@@ -71,12 +67,3 @@ void FieldFactory::AddSliceReference(double halfSizeX, double halfSizeY, double
}
}
-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
index 98882d7c2286c4415860c7e1d292c3034aaf7896..c0bb840be05b85ad8df200e2a59f0130d3bdd63a 100644
--- a/algo/kf/core/geo/KfField.h
+++ b/algo/kf/core/geo/KfField.h
@@ -33,6 +33,9 @@ namespace cbm::algo::kf
/// \tparam T Underlying floating point type (float/double/fvec)
template<typename T>
class alignas(VcMemAlign) FieldBase<T, EFieldMode::Orig> {
+ template<typename, EFieldMode>
+ friend class FieldBase;
+
public:
using FieldRegion_t = FieldRegion<T, EFieldMode::Orig>;
using FieldValue_t = FieldValue<T>;
@@ -42,7 +45,7 @@ namespace cbm::algo::kf
void SetFieldFunction(const FieldFn_t& fieldFn) { fFieldFn = fieldFn; }
/// \brief Creates field region object
- FieldRegion_t MakeFieldRegion() const { return FieldRegion_t(fFieldFn, fFieldType); }
+ FieldRegion_t MakeFieldRegion() const { return FieldRegion_t(fFieldType, fFieldFn); }
protected:
/// \brief Default constructor
@@ -56,18 +59,7 @@ namespace cbm::algo::kf
}
/// \brief Copy assignment operator
- /// \tparam I Underlying floating type of the source
- template<typename I>
- FieldBase& operator=(const FieldBase<I, EFieldMode::Orig>& other)
- {
- if constexpr (std::is_same_v<I, T>) {
- if (this == &other) {
- return *this;
- }
- }
- fFieldFn = other.fFieldFn;
- return *this;
- }
+ FieldBase& operator=(const FieldBase& other) = default;
FieldFn_t fFieldFn{defs::ZeroFieldFn}; ///< Field function: (x,y,z) [cm] -> (Bx,By,Bz) [kG]
EFieldType fFieldType{EFieldType::Null}; ///< Field type
@@ -88,12 +80,19 @@ namespace cbm::algo::kf
/// \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>;
+ template<typename, EFieldMode>
+ friend class FieldBase;
+
+ using SlicesContainer_t = std::vector<FieldSlice<T>>;
public:
using FieldRegion_t = FieldRegion<T, EFieldMode::Intrpl>;
using FieldValue_t = FieldValue<T>;
+ /// \brief Sets a field slice
+ /// \param fieldSlice Field slice object
+ void AddFieldSlice(const FieldSlice<T>& fieldSlice) { fvFieldSlices.push_back(fieldSlice); }
+
/// \brief Creates field value object
/// \param sliceID Index of slice
/// \param x x-coordinate of the point [cm]
@@ -103,6 +102,9 @@ namespace cbm::algo::kf
return fvFieldSlices[sliceID].GetFieldValue(x, y);
}
+ /// \brief Gets number of field slices in the instance
+ int GetNofFieldSlices() const { return fvFieldSlices.size(); }
+
/// \brief Creates field region object
/// \param b0 Field value in the first node [kG]
/// \param z0 First node z-coordinate [cm]
@@ -113,14 +115,9 @@ namespace cbm::algo::kf
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);
+ return FieldRegion_t(fFieldType, b0, z0, b1, z1, b2, z2);
}
- /// \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;
@@ -128,24 +125,18 @@ namespace cbm::algo::kf
/// \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)
+ FieldBase(const FieldBase<I, EFieldMode::Intrpl>& other)
{
- if constexpr (std::is_same_v<I, T>) {
- if (this == &other) {
- return *this;
- }
+ fvFieldSlices.clear();
+ fvFieldSlices.reserve(other.fvFieldSlices.size());
+ for (const auto& slice : other.fvFieldSlices) {
+ fvFieldSlices.emplace_back(FieldSlice<T>(slice));
}
- for (size_t i = 0; i < defs::MaxNofFieldSlices; ++i) {
- fvFieldSlices[i] = utils::simd::Cast<I, T>(other.fvFieldSlices[i]);
- }
- return *this;
}
+ /// \brief Copy assignment operator
+ FieldBase& operator=(const FieldBase& other) = default;
+
SlicesContainer_t fvFieldSlices; ///< Array of field slices
EFieldType fFieldType{EFieldType::Null}; ///< Field type
@@ -167,12 +158,13 @@ namespace cbm::algo::kf
/// \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>;
+ template<typename, EFieldMode>
+ friend class Field;
public:
static constexpr EFieldMode FieldType = FldMode;
- using FieldRegion_t = typename FieldBase_t::FieldRegion_t;
+ using FieldRegion_t = FieldRegion<T, FldMode>;
using Real_t = T;
/// \brief Default constructor
@@ -181,7 +173,7 @@ namespace cbm::algo::kf
/// \brief Copy constructor
template<typename I>
Field(const Field<I, FldMode>& other)
- : FieldBase_t(other)
+ : detail::FieldBase<T, FldMode>(other)
, fPrimVertexField(FieldRegion<I, FldMode>(other.fPrimVertexField))
{
}
@@ -190,21 +182,14 @@ namespace cbm::algo::kf
~Field() = default;
/// \brief Copy assignment operator
- template<typename I>
- Field& operator=(const Field<I, FldMode>& other)
- {
- if constexpr (std::is_same_v<I, T>) {
- if (this == &other) {
- return *this;
- }
- }
- fPrimVertexField = other.fPrimVertexField;
- return FieldBase_t::operator=(other);
- }
+ Field& operator=(const Field& other) = default;
/// \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:
@@ -216,12 +201,10 @@ namespace cbm::algo::kf
/// \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; }
-
/// \brief String representation of the class
/// \param indentLevel Indent level of the string output
- std::string ToString(int indentLevel) const;
+ /// \param verbose Verbosity level
+ std::string ToString(int indentLevel, int verbose) const;
private:
/// \brief Serialization function
@@ -229,7 +212,7 @@ namespace cbm::algo::kf
template<class Archive>
void serialize(Archive& ar, const unsigned int)
{
- ar& boost::serialization::base_object<FieldBase_t>(*this);
+ ar& boost::serialization::base_object<detail::FieldBase<T, FldMode>>(*this);
ar& fPrimVertexField;
}
@@ -292,10 +275,6 @@ namespace cbm::algo::kf
/// \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; }
@@ -334,12 +313,6 @@ namespace cbm::algo::kf
if (!fSliceReferences.size()) { // TODO: Remove requirement of slice references
throw std::logic_error("FieldFactory::MakeField: no slice references were provided");
}
- else if (fSliceReferences.size() > defs::MaxNofFieldSlices) {
- std::stringstream msg;
- msg << "FieldFactory::MakeField: too many slice references are provided (" << fSliceReferences.size() << "), "
- << "the maximum allowed number of field slices is kf::defs::MaxNofFieldSlices = " << defs::MaxNofFieldSlices;
- throw std::logic_error(msg.str());
- }
if (!fFieldFn) {
throw std::logic_error("FieldFactory::CreateField: no field function is provided");
}
@@ -351,10 +324,8 @@ namespace cbm::algo::kf
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));
+ field.AddFieldSlice(FieldSlice<T>(fFieldFn, sliceRef.fHalfSizeX, sliceRef.fHalfSizeY, sliceRef.fRefZ));
}
{
// PV field initialization
@@ -362,7 +333,7 @@ namespace cbm::algo::kf
double z = fTarget[2];
std::array<FieldValue<T>, nNodes> fld;
std::array<T, nNodes> pos;
- for (int iNode = 0; iNode < nNodes; ++iNode) {
+ for (size_t 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;
diff --git a/algo/kf/core/geo/KfFieldRegion.cxx b/algo/kf/core/geo/KfFieldRegion.cxx
index 54d3cc53e6ed1696662c1967aea8d43b36a707db..fc892352463994563cd665a5adca6256732a1b8d 100644
--- a/algo/kf/core/geo/KfFieldRegion.cxx
+++ b/algo/kf/core/geo/KfFieldRegion.cxx
@@ -15,15 +15,15 @@
#include <mutex>
#include <sstream>
+#include <fmt/format.h>
+
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>;
+using cbm::algo::kf::detail::FieldRegionBase;
FieldFn_t GlobalField::gOriginalField = &GlobalField::UndefField;
@@ -32,8 +32,8 @@ 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)
+void FieldRegionBase<T, EFieldMode::Intrpl>::Set(const FieldValue<T>& b0, const T& z0, const FieldValue<T>& b1,
+ const T& z1, const FieldValue<T>& b2, const T& z2)
{
fZfirst = z0;
@@ -58,7 +58,7 @@ void FieldRegionIntrpl_t<T>::Set(const FieldValue<T>& b0, const T& z0, const Fie
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T>
-void FieldRegionIntrpl_t<T>::Shift(const T& z)
+void FieldRegionBase<T, EFieldMode::Intrpl>::Shift(const T& z)
{
auto dz = z - fZfirst;
for (int iD = 0; iD < 3; ++iD) {
@@ -70,6 +70,13 @@ void FieldRegionIntrpl_t<T>::Shift(const T& z)
fZfirst = z;
}
+namespace cbm::algo::kf::detail
+{
+ template class FieldRegionBase<float, EFieldMode::Intrpl>;
+ template class FieldRegionBase<double, EFieldMode::Intrpl>;
+ template class FieldRegionBase<fvec, EFieldMode::Intrpl>;
+} // namespace cbm::algo::kf::detail
+
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T, EFieldMode FldMode>
@@ -92,21 +99,32 @@ FieldValue<T> FieldRegion<T, FldMode>::Get(const T& x, const T& y, const T& z) c
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T, EFieldMode FldMode>
-std::string FieldRegion<T, FldMode>::ToString(int indentLevel) const
+std::string FieldRegion<T, FldMode>::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 << "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";
+ 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 "
- << "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";
+ 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) {
@@ -126,25 +144,6 @@ std::tuple<double, double, double> GlobalField::UndefField(double, double, doubl
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
{
diff --git a/algo/kf/core/geo/KfFieldRegion.h b/algo/kf/core/geo/KfFieldRegion.h
index d7f40d39b117084427671f1f039770572c7cde14..9e8ceec5e44f2744a66fc03b480992f8c8a52ceb 100644
--- a/algo/kf/core/geo/KfFieldRegion.h
+++ b/algo/kf/core/geo/KfFieldRegion.h
@@ -38,11 +38,7 @@ namespace cbm::algo::kf
/// \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);
+ static void SetFieldFunction(const FieldFn_t& fn) { gOriginalField = 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
@@ -56,7 +52,6 @@ namespace cbm::algo::kf
const T& z);
};
-
namespace detail
{
/// \class FieldRegionBase
@@ -74,6 +69,9 @@ namespace cbm::algo::kf
/// \tparam T Underlying floating point type (float/double/fvec)
template<typename T>
class alignas(VcMemAlign) FieldRegionBase<T, EFieldMode::Orig> {
+ template<typename, EFieldMode>
+ friend class FieldRegionBase;
+
public:
using FieldValue_t = FieldValue<T>;
@@ -94,18 +92,7 @@ namespace cbm::algo::kf
~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 constexpr (std::is_same_v<I, T>) {
- if (this == &other) {
- return *this;
- }
- }
- fFieldFn = other.fFieldFn;
- return *this;
- }
+ FieldRegionBase& operator=(const FieldRegionBase& other) = default;
protected:
FieldFn_t fFieldFn{defs::ZeroFieldFn}; ///< Field function: (x,y,z) [cm] -> (Bx,By,Bz) [kG]
@@ -125,6 +112,9 @@ namespace cbm::algo::kf
/// \tparam T Underlying floating point type (float/double/fvec)
template<typename T>
class alignas(VcMemAlign) FieldRegionBase<T, EFieldMode::Intrpl> {
+ template<typename, EFieldMode>
+ friend class FieldRegionBase;
+
using CoeffArray_t = std::array<std::array<T, 3>, 3>; ///< Approximation coefficients array
public:
@@ -149,28 +139,20 @@ namespace cbm::algo::kf
/// \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))
+ : fZfirst(utils::simd::Cast<I, T>(other.fZfirst))
{
+ for (size_t iDim = 0; iDim < fCoeff.size(); ++iDim) {
+ for (size_t iPow = 0; iPow < fCoeff[iDim].size(); ++iPow) {
+ fCoeff[iDim][iPow] = utils::simd::Cast<I, T>(other.fCoeff[iDim][iPow]);
+ }
+ }
}
/// \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 constexpr (std::is_same_v<I, T>) {
- if (this == &other) {
- return *this;
- }
- }
- this->fCoeff = utils::simd::Cast<I, T>(other.fCoeff);
- this->fZfirst = utils::simd::Cast<I, T>(other.fZfirst);
- return *this;
- }
+ FieldRegionBase& operator=(const FieldRegionBase& other) = default;
/// \brief Gets the first z-coordinate of the interpolation
const T& GetZfirst() const { return fZfirst; }
@@ -213,7 +195,8 @@ namespace cbm::algo::kf
/// \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>;
+ template<typename, EFieldMode>
+ friend class FieldRegion;
public:
static constexpr EFieldMode FieldType = FldMode;
@@ -227,16 +210,17 @@ namespace cbm::algo::kf
/// \brief Copy constructor
/// \tparam I Underlying floating point type of the source
template<typename I>
- FieldRegion(const FieldRegion<I, FldMode>& other) : FieldRegionBase_t(other)
+ FieldRegion(const FieldRegion<I, FldMode>& other) : detail::FieldRegionBase<T, FldMode>(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
+ /// \param args Arguments, which have to be passed to the base class constructors
template<typename... Args>
- FieldRegion(Args... args, EFieldType fieldType) : FieldRegionBase_t(args...)
- , fFieldType(fieldType)
+ FieldRegion(EFieldType fieldType, Args... args)
+ : detail::FieldRegionBase<T, FldMode>(args...)
+ , fFieldType(fieldType)
{
}
@@ -244,18 +228,7 @@ namespace cbm::algo::kf
~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 constexpr (std::is_same_v<I, T>) {
- if (this == &other) {
- return *this;
- }
- }
- fFieldType = other.fFieldType;
- return FieldRegionBase_t::operator=(other);
- }
+ FieldRegion& operator=(const FieldRegion& other) = default;
/// \brief Gets the field value in a spatial point on the track
/// \param x x-coordinate of the point [cm]
@@ -269,7 +242,8 @@ namespace cbm::algo::kf
/// \brief String representation of the class content
/// \param intdentLevel Indent level
- std::string ToString(int indentLevel = 0) const;
+ /// \param verbose Verbosity level
+ std::string ToString(int indentLevel = 0, int verbose = 1) const;
private:
/// \brief Serialization method
@@ -277,7 +251,8 @@ namespace cbm::algo::kf
template<class Archive>
void serialize(Archive& ar, const unsigned int)
{
- ar& boost::serialization::base_object<FieldRegionBase_t>(*this);
+ using boost::serialization::base_object;
+ ar& base_object<detail::FieldRegionBase<T, FldMode>>(*this);
ar& fFieldType;
}
diff --git a/algo/kf/core/geo/KfFieldSlice.cxx b/algo/kf/core/geo/KfFieldSlice.cxx
index b5f94447294ab8898fe5a37a97b43ddade9d4a01..8e8989ea1c8d0c212853abe841efafc92e7ae243 100644
--- a/algo/kf/core/geo/KfFieldSlice.cxx
+++ b/algo/kf/core/geo/KfFieldSlice.cxx
@@ -12,6 +12,8 @@
#include <iomanip>
#include <sstream>
+#include <fmt/format.h>
+
using cbm::algo::kf::FieldFn_t;
using cbm::algo::kf::FieldSlice;
using cbm::algo::kf::FieldValue;
@@ -53,8 +55,8 @@ FieldSlice<T>::FieldSlice(const FieldFn_t& fieldFn, double xMax, double yMax, do
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];
+ for (int j = kNofCoeff - 1; j >= kNofCoeff - i; --j) { // loop over x powers
+ m[k--] = xFactor * m[j];
}
xFactor *= x;
}
@@ -96,44 +98,47 @@ FieldSlice<T>::FieldSlice(const FieldFn_t& fieldFn, double xMax, double yMax, do
}
}
- 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]);
- }
+ for (int j = 0; j < kNofCoeff; ++j) {
+ fBx[j] = utils::simd::Cast<double, T>(A[j][kNofCoeff] / A[j][j]);
+ fBy[j] = utils::simd::Cast<double, T>(A[j][kNofCoeff + 1] / A[j][j]);
+ fBz[j] = utils::simd::Cast<double, T>(A[j][kNofCoeff + 2] / A[j][j]);
}
}
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T>
-std::string FieldSlice<T>::ToString(int indentLevel) const
+std::string FieldSlice<T>::ToString(int indentLevel, int verbose) const
{
- using std::setw;
+ using fmt::format;
constexpr char indentChar = '\t';
std::stringstream msg;
std::string indent(indentLevel, indentChar);
+ auto Cnv = [&](const auto& val) { return utils::simd::Cast<T, Literal_t<T>>(val); }; // alias for the conversion fn
- 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];
+ if (verbose > 0) {
+ msg << indent << format("FieldSlice: zRef = {} cm", Cnv(fZref));
+ if (verbose > 1) {
+ std::array<std::string, kNofCoeff> m;
+ {
+ m[kNofCoeff - 1] = "";
+ for (int i = kNofCoeff - 2; i >= kNofCoeff - 1 - kPolDegree; --i) {
+ m[i] = format("y{}", 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 = format("x{}", kPolDegree - i + 1);
+ for (int j = kNofCoeff - 1; j >= kNofCoeff - i; --j) { // loop over x powers
+ m[k--] = x + m[j];
+ }
+ }
+ m[kNofCoeff - 1] = "1";
+ }
+ msg << '\n' << indent << format("{:>15} {:>15} {:>15} {:>15}", "Monomial", "Bx", "By", "Bz");
+ for (int i = 0; i < kNofCoeff; ++i) {
+ msg << '\n' << indent << format("{:>15} {:>15} {:>15} {:>15}", m[i], Cnv(fBx[i]), Cnv(fBy[i]), Cnv(fBz[i]));
}
}
- 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();
}
diff --git a/algo/kf/core/geo/KfFieldSlice.h b/algo/kf/core/geo/KfFieldSlice.h
index 11fc02ad6b5df4a5378b4e99e9541573180e31eb..295051c2423da08d86d1fe574f0b18f876ea65bc 100644
--- a/algo/kf/core/geo/KfFieldSlice.h
+++ b/algo/kf/core/geo/KfFieldSlice.h
@@ -29,11 +29,14 @@ namespace cbm::algo::kf
/// \tparam T Underlying data-type
template<typename T>
class alignas(VcMemAlign) FieldSlice {
+ template<typename>
+ friend class 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>;
+ /// \brief Array of the approximation coefficients [<monomial>]
+ using CoeffArray_t = std::array<T, kNofCoeff>;
public:
/// \brief Default constructor
@@ -50,8 +53,10 @@ namespace cbm::algo::kf
/// \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))
+ : fBx(utils::simd::Cast<I, T, kNofCoeff>(other.fBx))
+ , fBy(utils::simd::Cast<I, T, kNofCoeff>(other.fBy))
+ , fBz(utils::simd::Cast<I, T, kNofCoeff>(other.fBz))
+ , fZref(utils::simd::Cast<I, T>(other.fZref))
{
}
@@ -59,16 +64,7 @@ namespace cbm::algo::kf
~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;
- }
+ FieldSlice& operator=(const FieldSlice& other) = default;
/// \brief Gets field value at a point on the transverse plane
/// \param x x-coordinate of the point [cm]
@@ -78,9 +74,9 @@ namespace cbm::algo::kf
{
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));
+ return FieldValue<T>(Horner<kPolDegree>(fBx.cbegin(), x, y),
+ Horner<kPolDegree>(fBy.cbegin(), x, y),
+ Horner<kPolDegree>(fBz.cbegin(), x, y));
/* clang-format on */
}
@@ -98,7 +94,8 @@ namespace cbm::algo::kf
/// \brief String representation of the class content
/// \param intdentLevel Indent level
- std::string ToString(int indentLevel = 0) const;
+ /// \param verbose Verbosity level
+ std::string ToString(int indentLevel = 0, int verbose = 1) const;
private:
/// \brief Serialization function
@@ -106,11 +103,15 @@ namespace cbm::algo::kf
template<class Archive>
void serialize(Archive& ar, const unsigned int)
{
- ar& fCoeff;
+ ar& fBx;
+ ar& fBy;
+ ar& fBz;
ar& fZref;
}
- CoeffArray_t fCoeff{{T(0.)}}; ///< Approximation coefficients for the field slice
- T fZref{defs::Undef<T>}; ///< Reference z of the coefficients [cm]
+ CoeffArray_t fBx{{T(0.)}}; ///< Approximation coefficients for the x-component of the field
+ CoeffArray_t fBy{{T(0.)}}; ///< Approximation coefficients for the y-component of the field
+ CoeffArray_t fBz{{T(0.)}}; ///< Approximation coefficients for the z-component of the field
+ T fZref{defs::Undef<T>}; ///< Reference z of the coefficients [cm]
};
} // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfFieldValue.h b/algo/kf/core/geo/KfFieldValue.h
index b26a3791b5c4ff3ea8d701857a9c64db899d3c0d..b4b1a8fc22c7832f11c2ae89f2fdd66faef6cd42 100644
--- a/algo/kf/core/geo/KfFieldValue.h
+++ b/algo/kf/core/geo/KfFieldValue.h
@@ -56,15 +56,7 @@ namespace cbm::algo::kf
~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;
- }
+ FieldValue& operator=(const FieldValue& other) = default;
/// \brief Combines the current magnetic field value with another one using a mask
/// \param other Other flux value
diff --git a/algo/kf/core/geo/KfMaterialMap.cxx b/algo/kf/core/geo/KfMaterialMap.cxx
index e6f69ce164128064de0ad604d36060374053e48d..5708155e66ef5d3d1ad01b096ccfa021b7bb73da 100644
--- a/algo/kf/core/geo/KfMaterialMap.cxx
+++ b/algo/kf/core/geo/KfMaterialMap.cxx
@@ -12,6 +12,8 @@
#include <sstream>
#include <vector>
+#include <fmt/format.h>
+
using cbm::algo::kf::fvec;
using cbm::algo::kf::MaterialMap;
@@ -50,24 +52,24 @@ 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]
+ constexpr int nRays{100}; // 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.F};
+ const auto binSize{2.F * 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 xBinOther{-this->fXYmax * scaleFactor + stepSize * 0.5F};
+ 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) {
+ for (int iBinY{0}; iBinY < this->fNbins; ++iBinY) {
+ for (int iBinX{0}; iBinX < this->fNbins; ++iBinX) {
// Collect material using ray shooting
- float avgThickness = 0; // Collected average thickness
- for (int iRayY = 0; iRayY < nRays; ++iRayY) {
- for (int iRayX = 0; iRayX < nRays; ++iRayX) {
- avgThickness += other.GetRadThick(xBinOther + iBinX * stepSize, yBinOther + iBinY * stepSize);
+ float avgThickness{0}; // Collected average thickness
+ for (int iRayY{0}; iRayY < nRays; ++iRayY) {
+ for (int iRayX{0}; iRayX < nRays; ++iRayX) {
+ avgThickness += other.GetThickness(xBinOther + iBinX * stepSize, yBinOther + iBinY * stepSize);
}
}
this->fTable[iBinX + this->fNbins * iBinY] += avgThickness / (nRays * nRays);
@@ -83,14 +85,49 @@ void MaterialMap::Add(const MaterialMap& other, float zTarg)
//
int MaterialMap::GetBin(float x, float y) const
{
- int i = static_cast<int>((x + fXYmax) * fFactor);
- int j = static_cast<int>((y + fXYmax) * fFactor);
+ int i{static_cast<int>((x + fXYmax) * fFactor)};
+ int j{static_cast<int>((y + fXYmax) * fFactor)};
if (i < 0 || j < 0 || i >= fNbins || j >= fNbins) {
return -1;
}
return i + j * fNbins;
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void MaterialMap::Rebin(int nGroups)
+{
+ if (nGroups < 1 && nGroups > fNbins) {
+ int nGroupsSet{nGroups};
+ nGroups = std::clamp(nGroups, 1, fNbins);
+ LOG(warn) << "kf::MaterialMap::Rebin(): incorrect input parameter nGroups = " << nGroupsSet << " is set to "
+ << nGroups;
+ }
+
+ int nBinsNew = static_cast<int>(static_cast<bool>(fNbins % nGroups)) + fNbins / nGroups;
+ std::vector<float> table(nBinsNew * nBinsNew);
+ for (int iX{0}; iX < nBinsNew; ++iX) {
+ for (int iY{0}; iY < nBinsNew; ++iY) {
+ float value{0.F};
+ int nEntries{0};
+ int iOldXmin{iX * nGroups};
+ int iOldXmax{std::max((iX + 1) * nGroups, fNbins)};
+ int iOldYmin{iY * nGroups};
+ int iOldYmax{std::max((iY + 1) * nGroups, fNbins)};
+ for (int iOldX{iOldXmin}; iOldX < iOldXmax; ++iOldX) {
+ for (int iOldY{iOldYmin}; iOldY < iOldYmax; ++iOldY) {
+ value += fTable[iOldX + fNbins * iOldY];
+ ++nEntries;
+ }
+ }
+ table[iX + nBinsNew * iY] = value / nEntries;
+ }
+ }
+ fTable = std::move(table);
+ fFactor = fFactor / fNbins * nBinsNew;
+ fNbins = nBinsNew;
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
void MaterialMap::Swap(MaterialMap& other) noexcept
@@ -106,12 +143,38 @@ void MaterialMap::Swap(MaterialMap& other) noexcept
// ---------------------------------------------------------------------------------------------------------------------
//
-std::string MaterialMap::ToString() const
+std::string MaterialMap::ToString(int indentLevel, int verbose) const
{
- std::stringstream msg;
+ using fmt::format;
using std::setw;
- msg << "[mat.map] nBins " << setw(6) << fNbins << ", XYmax " << setw(12) << fXYmax << " , z (ref, min, max) "
- << setw(12) << fZref << ' ' << setw(12) << fZmin << ' ' << setw(12) << fZmax;
+ std::stringstream msg;
+ if (verbose > 0) {
+ constexpr char indentCh = '\t';
+ std::string indent(indentLevel, indentCh);
+ msg << indent
+ << format("zRef = {:<12} cm, range [{:<12}, {:<12}] cm, nBins = {:<8}, XYmax = {:<12} cm", fZref, fZmin, fZmax,
+ fNbins, fXYmax);
+ if (verbose > 1) {
+ msg << indent << indentCh << "\nContent(rebinned to 10 bins):\n";
+ msg.precision(3);
+ auto mapTmp{*this};
+ mapTmp.Rebin(fNbins / 10);
+ msg << indent << indentCh << setw(15) << "y [cm] \\ x [cm]" << ' ';
+ for (int iX{0}; iX < mapTmp.fNbins; ++iX) {
+ float xRef{-mapTmp.fXYmax + (mapTmp.fXYmax * (2 * iX + 1)) / mapTmp.fNbins};
+ msg << indent << indentCh << setw(8) << xRef << ' ';
+ }
+ msg << '\n';
+ for (int iY{0}; iY < mapTmp.fNbins; ++iY) {
+ float yRef{-mapTmp.fXYmax + (mapTmp.fXYmax * (2 * iY + 1)) / mapTmp.fNbins};
+ msg << indent << indentCh << setw(15) << yRef << ' ';
+ for (int iX{0}; iX < mapTmp.fNbins; ++iX) {
+ msg << indent << indentCh << setw(8) << mapTmp.fTable[iX + mapTmp.fNbins * iY] << ' ';
+ }
+ msg << '\n';
+ }
+ }
+ }
return msg.str();
}
diff --git a/algo/kf/core/geo/KfMaterialMap.h b/algo/kf/core/geo/KfMaterialMap.h
index 30289c882bdc9628c2bd2d945f86a3f2f550a26f..46e40f00922f9b6ce8e409bbb51f895c8cc463a1 100644
--- a/algo/kf/core/geo/KfMaterialMap.h
+++ b/algo/kf/core/geo/KfMaterialMap.h
@@ -72,26 +72,17 @@ namespace cbm::algo::kf
/// \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
- /// \param iBinX Index of table column
- /// \param iBinY Index of table row
- float GetRadThickBin(int iBinX, int iBinY) const { return fTable[iBinX + fNbins * iBinY]; }
-
- /// \brief Gets value of X/X0 in a given cell of the material table by the index of the bin
- /// \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
+ /// \brief Gets material thickness in units of radiational length X0
/// \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
+ I GetThickness(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]);
+ res[i] = GetThickness(x[i], y[i]);
}
return res;
}
@@ -112,6 +103,10 @@ namespace cbm::algo::kf
return utils::IsUndefined(fNbins) || utils::IsUndefined(fXYmax * fFactor * fZref * fZmin * fZmax);
}
+ /// \brief Reduces number of bins by a given factor
+ /// \param factor Number of bins in a new bin
+ void Rebin(int factor);
+
/// \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
@@ -124,11 +119,10 @@ namespace cbm::algo::kf
/// \brief Swap method
void Swap(MaterialMap& other) noexcept;
- /// \brief Get bin index for (x,y). Returns -1 when outside of the map
- int GetBin(float x, float y) const;
-
/// \brief String representation of the object
- std::string ToString() const;
+ /// \param indentLevel Indent level of the string output
+ /// \param verbose Verbosity level
+ std::string ToString(int indentLevel = 0, int verbose = 1) const;
/// \brief Comparison operator (material map ordering by fZref)
friend bool operator<(const MaterialMap& lhs, const MaterialMap& rhs) { return lhs.fZref < rhs.fZref; }
@@ -138,6 +132,18 @@ namespace cbm::algo::kf
/// \throw std::logic_error If the object is in non-valid mode
void CheckConsistency() const;
+ /// \brief Get bin index for (x,y). Returns -1 when outside of the map
+ int GetBin(float x, float y) const;
+
+ /// \brief Gets material thickness in units of radiational length X0
+ /// \param iBin Index of the bin in 2d table
+ float GetThickness(int iBin) const { return fTable[iBin]; }
+
+ /// \brief Gets material thickness in units of radiational length X0
+ /// \param iBinX Index of table column
+ /// \param iBinY Index of table row
+ float GetThickness(int iBinX, int iBinY) const { return fTable[iBinX + fNbins * iBinY]; }
+
int fNbins = defs::Undef<int>; ///< Number of rows (== N columns) in the material budget table
float fXYmax = defs::Undef<float>; ///< Size of the station in x and y dimensions [cm]
float fFactor = defs::Undef<float>; ///< Util. var. for the conversion of point coordinates to row/column id
diff --git a/algo/kf/core/geo/KfSetup.cxx b/algo/kf/core/geo/KfSetup.cxx
index af7b923f2b59382013524fa515920584f0e4f35a..ec8862775685506164e4325da359e7d930ba9096 100644
--- a/algo/kf/core/geo/KfSetup.cxx
+++ b/algo/kf/core/geo/KfSetup.cxx
@@ -10,14 +10,9 @@
#include "KfSetup.h"
#include <sstream>
+#include <typeinfo>
using cbm::algo::kf::Setup;
-using cbm::algo::kf::SetupInitializer;
-
-//
-// Setup methods definition
-//
-
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -28,112 +23,28 @@ std::string Setup<T>::ToString(int verbosity, int indentLevel) const
if (verbosity > 0) {
constexpr char indentCh = '\t';
std::string indent(indentLevel, indentCh);
- msg << indent << "----- KF setup -----\n";
- msg << indent << "Target:\n" << fTarget.ToString(indentLevel + 1);
- msg << indent << "Material layers:\n";
+ msg << indent << "\n ----- KF Setup -----\n";
+ msg << indent << "\nFloating-point type: " << typeid(T).name();
+
+ msg << indent << "\nTARGET:\n" << fTarget.ToString(indentLevel + 1, verbosity);
+ msg << indent << "\nMATERIAL LAYERS:\n";
for (const auto& layer : fvMaterialLayers) {
- msg << indent << indentCh << layer.ToString() << '\n';
+ msg << layer.ToString(indentLevel + 1, verbosity) << '\n';
}
- msg << indent << "Interpolated magnetic field:\n" << fIntrplField.ToString(indentLevel + 1);
+ msg << indent << "\nINTERPOLATED MAGNETIC FIELD:\n" << fIntrplField.ToString(indentLevel + 1, verbosity);
if (foOrigField) {
- msg << indent << "Original magnetic field:\n" << foOrigField->ToString(indentLevel + 1);
+ msg << indent << "\nORIGINAL MAGNETIC FIELD:\n" << foOrigField->ToString(indentLevel + 1, verbosity);
}
else {
- msg << indent << "Original magnetic field: not provided";
+ msg << indent << "\nORIGINAL MANGETIC FIELD: not provided";
}
}
return msg.str();
}
-
namespace cbm::algo::kf
{
template class Setup<float>;
template class Setup<double>;
template class Setup<fvec>;
} // namespace cbm::algo::kf
-
-//
-// SetupInitializer methods definition
-//
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void SetupInitializer::AddMaterial(const MaterialMap& material)
-{
- if (!fMaterialLayers.emplace(material).second) {
- std::stringstream msg;
- msg << "SetupInitializer::AddMaterial: attempt of adding a duplicating material layer " << material.ToString()
- << ".\nThe next material layers were already added:";
- for (const auto& el : fMaterialLayers) {
- msg << "\n\t- " << el.ToString();
- }
- throw std::logic_error(msg.str());
- }
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-template<typename T>
-Setup<T> SetupInitializer::MakeSetup(bool bProvideOrigField) const
-{
- Setup<T> setup;
- // Target initialization
- setup.SetTarget(fTarget);
-
- // Field initialization
- setup.SetField(fFieldFactory.MakeField<T, EFieldMode::Intrpl>());
- if (bProvideOrigField) {
- setup.SetField(fFieldFactory.MakeField<T, EFieldMode::Orig>());
- }
-
- // Material layers initialization
- if (fMaterialLayers.size() > defs::MaxNofMaterialLayers) {
- std::stringstream msg;
- msg << "kf::SetupInitializer::MakeSetup(): too many material layers are provided (" << fMaterialLayers.size()
- << "), the maximum allowed number of the layers kf::defs::MaxNofMaterialLayers = "
- << defs::MaxNofMaterialLayers;
- throw std::logic_error(msg.str());
- }
- int layerId = 0;
- for (const auto& material : fMaterialLayers) {
- setup.SetMaterial(layerId, material);
- }
-
- return setup;
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void SetupInitializer::Reset()
-{
- fFieldFactory.Reset();
- fMaterialLayers.clear();
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void SetupInitializer::SetTargetProperties(double x, double y, double z, double fieldInitStep,
- const MaterialMap& material)
-{
- fTarget.SetX(x);
- fTarget.SetY(y);
- fTarget.SetZ(z);
- fTarget.SetMaterial(material);
- fFieldFactory.SetTarget(x, y, z);
- fFieldFactory.SetStep(fieldInitStep);
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void SetupInitializer::Store(const Setup<double>& setup, const std::string& fileName)
-{
- std::ofstream ofs(fileName, std::ios::binary);
- if (!ofs) {
- std::stringstream msg;
- msg << "kf::SetupInitializer::Store: failed openning file \"" << fileName << "\" to store the setup";
- throw std::runtime_error(msg.str());
- }
- boost::archive::binary_oarchive oa(ofs);
- oa << setup;
-}
diff --git a/algo/kf/core/geo/KfSetup.h b/algo/kf/core/geo/KfSetup.h
index 7cced174893dfce18170b7ba611abaecd15b0d78..75829aa2e64b4f53d9b0a2938a60d1a9c175e7c4 100644
--- a/algo/kf/core/geo/KfSetup.h
+++ b/algo/kf/core/geo/KfSetup.h
@@ -15,8 +15,6 @@
#include "KfTarget.h"
#include "KfVector.h"
-#include <boost/archive/binary_iarchive.hpp>
-#include <boost/archive/binary_oarchive.hpp>
#include <boost/serialization/access.hpp>
#include <boost/serialization/split_free.hpp>
@@ -34,10 +32,9 @@ namespace cbm::algo::kf
template<typename T>
class alignas(VcMemAlign) Setup {
friend class boost::serialization::access; // Boost serializer methods
- template<typename I>
+ template<typename>
friend class Setup;
- using MaterialContainer_t = std::array<MaterialMap, defs::MaxNofMaterialLayers>;
public:
/// \brief Default constructor
@@ -51,8 +48,8 @@ namespace cbm::algo::kf
template<typename I>
Setup(const Setup<I>& other)
: fvMaterialLayers(other.fvMaterialLayers)
- , foOrigField(other.foOrigField.has_value() ? std::make_optional(*other.foOrigField) : std::nullopt)
, fIntrplField(other.fIntrplField)
+ , foOrigField(other.foOrigField.has_value() ? std::make_optional(*other.foOrigField) : std::nullopt)
, fTarget(other.fTarget)
{
}
@@ -61,25 +58,15 @@ namespace cbm::algo::kf
//Setup(Setup&&) noexcept;
/// \brief Copy assignment operator
- /// \tparam I Underlying floating-point data type of the source
- template<typename I>
- Setup& operator=(const Setup<I>& other)
- {
- if constexpr (std::is_same_v<I, T>) {
- if (this == &other) {
- return *this;
- }
- }
- fvMaterialLayers = other.fvMaterialLayers;
- foOrigField = other.foOrigField;
- fIntrplField = other.fIntrplField;
- fTarget = other.fTarget;
- return *this;
- }
+ Setup& operator=(const Setup& other) = default;
/// \brief Move assignment operator
//Setup& operator=(Setup&&) noexcept;
+ /// \brief Sets material layer
+ /// \param map Material map
+ 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)
@@ -90,23 +77,21 @@ namespace cbm::algo::kf
/// \param iLayer Index of layer
const MaterialMap& GetMaterial(int iLayer) const { return fvMaterialLayers[iLayer]; }
+ /// \brief Gets number of defined material layers
+ int GetNofMaterialLayers() const { return static_cast<int>(fvMaterialLayers.size()); }
+
/// \brief Gets target
const Target<T>& GetTarget() const { return fTarget; }
/// \brief String representation of the class contents
/// \param verbosity A verbose level for output
/// \param indentLevel Indent level of the string output
- std::string ToString(int verbosity = 0, int indentLevel = 0) const;
+ std::string ToString(int verbosity = 1, int indentLevel = 0) const;
/// \brief Sets field
template<EFieldMode FldMode>
void SetField(const Field<double, FldMode>& field);
- /// \brief Sets material layer
- /// \param iLayer Index of layer
- /// \param map Material map
- void SetMaterial(int iLayer, const MaterialMap& map) { fvMaterialLayers[iLayer] = map; }
-
/// \brief Sets target (moves from the source)
/// \param target Target instance
void SetTarget(const Target<double>& target) { fTarget = target; }
@@ -117,6 +102,7 @@ namespace cbm::algo::kf
void load(Archive& ar, const unsigned int /*version*/)
{
ar >> fvMaterialLayers;
+ ar >> fvMaterialLayerIds;
ar >> fTarget;
ar >> fIntrplField;
foOrigField = std::nullopt;
@@ -129,6 +115,7 @@ namespace cbm::algo::kf
void save(Archive& ar, const unsigned int /*version*/) const
{
ar << fvMaterialLayers;
+ ar << fvMaterialLayerIds;
ar << fTarget;
ar << fIntrplField;
}
@@ -139,109 +126,13 @@ namespace cbm::algo::kf
boost::serialization::split_member(ar, *this, version);
}
- MaterialContainer_t fvMaterialLayers = {}; ///< Container of the material maps
- Field<T, EFieldMode::Intrpl> fIntrplField; ///< Interpolated field (NOTE: maybe make optional)
+ 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
};
-
- /// \class SetupInitializer
- /// \brief Creates a valid initialized Setup instance
- class SetupInitializer {
- public:
- /// \brief Default constructor
- SetupInitializer() = default;
-
- /// \brief Copy constructor
- SetupInitializer(const SetupInitializer&) = delete;
-
- /// \brief Move constructor
- SetupInitializer(SetupInitializer&&) = delete;
-
- /// \brief Destructor
- virtual ~SetupInitializer() = default;
-
- /// \brief Copy assignment operator
- SetupInitializer& operator=(const SetupInitializer&) = delete;
-
- /// \brief Move assignment operator
- SetupInitializer& operator=(SetupInitializer&&) = delete;
-
- /// \brief Adds magnetic field 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 AddFieldSliceReference(double halfSizeX, double halfSizeY, double refZ)
- {
- fFieldFactory.AddSliceReference(halfSizeX, halfSizeY, refZ);
- }
-
- /// \brief Adds material layer
- /// \param material Material layer
- void AddMaterial(const MaterialMap& material);
-
- /// \brief Virtual function, which initializes a KF-setup
- /// \return true Initialization was successful
- /// \return false Initialization failed
- ///
- /// This function can be re-implemented for a particular use-case...
- 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;
-
- /// \brief Resets the instance
- void Reset();
-
- /// \brief Sets magnetic field function
- /// \param fieldFn Magnetic field function (either KF, or FairRoot format)
- /// \param fieldType Magnetic field type
- template<typename FieldFn>
- void SetFieldFunction(const FieldFn& fieldFn, EFieldType fieldType)
- {
- fFieldFactory.SetFieldFunction(fieldFn);
- fFieldFactory.SetFieldType(fieldType);
- }
-
- /// \brief Sets target initialization properties
- /// \param x Target x-coordinate [cm]
- /// \param y Target y-coordinate [cm]
- /// \param z Target z-coordinate [cm]
- /// \param fieldInitStep Step between the nodal points in the target field region initialization [cm]
- /// \param material Target material layer
- void SetTargetProperties(double x, double y, double z, double fieldInitStep, const MaterialMap& material);
-
- /// \brief Stores a serialized setup to a file
- /// \param setup A reference to the Setup (NOTE: a double-Setup must be passed explicitly)
- /// \param fileName Output file name
- ///
- /// Only a Setup with T = double can be serialized, so if a particular setup instance has another floating-point
- /// type, it will be converted to the double-version. In the serialization only the interpolated magnetic field
- /// variant is stored, the original field version will be set to nullopt during the loading from the file.
- static void Store(const Setup<double>& setup, const std::string& fileName);
-
- /// \brief Loads a serialized setup from a file
- /// \tparam T Underlying floating-point type for the output setup object
- /// \param fileName Input file name
- /// \throw std::runtime_error If the fileName cannot be opened
- template<typename T>
- static Setup<T> Load(const std::string& fileName);
-
- private:
- std::set<MaterialMap> fMaterialLayers{}; ///< Set of material maps
- FieldFactory fFieldFactory; ///< Instance of field factory
- Target<double> fTarget; ///< Field to keep target properties
- bool fbPropertiesInitialized{false}; ///< Flag: if all the necessary fields are initialized
- };
-
-
- // ********************************
- // ** Template method definition **
- // ********************************
-
// -------------------------------------------------------------------------------------------------------------------
//
template<typename T>
@@ -277,31 +168,4 @@ namespace cbm::algo::kf
foOrigField = std::make_optional(field);
}
}
-
- // -------------------------------------------------------------------------------------------------------------------
- //
- template<typename T>
- Setup<T> SetupInitializer::Load(const std::string& fileName)
- {
- Setup<double> setup;
- std::ifstream ifs(fileName, std::ios::binary);
- if (!ifs) {
- std::stringstream msg;
- msg << "kf::SetupInitializer::Load: intput setup file \"" << fileName << "\" was not found";
- throw std::runtime_error(msg.str());
- }
- try {
- boost::archive::binary_iarchive ia(ifs);
- ia >> setup;
- }
- catch (const std::exception& err) {
- std::stringstream msg;
- msg << "kf::SetupInitializer::Load: input setup file \"" << fileName
- << "\" has inconsistent format or was "
- "corrupted. The exception message: "
- << err.what();
- throw std::runtime_error(msg.str());
- }
- return Setup<T>(setup);
- }
} // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfSetupInitializer.cxx b/algo/kf/core/geo/KfSetupInitializer.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..3148dc428d71061c8d647bc3a383da7f9c94d333
--- /dev/null
+++ b/algo/kf/core/geo/KfSetupInitializer.cxx
@@ -0,0 +1,63 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file KfSetupInitializer.h
+/// @brief A base KF-Setup initialization class (source)
+/// @since 28.03.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#include "KfSetupInitializer.h"
+
+using cbm::algo::kf::SetupInitializer;
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetupInitializer::AddMaterial(const MaterialMap& material)
+{
+ if (!fMaterialLayers.emplace(material).second) {
+ std::stringstream msg;
+ msg << "SetupInitializer::AddMaterial: attempt of adding a duplicating material layer " << material.ToString()
+ << ".\nThe next material layers were already added:";
+ for (const auto& el : fMaterialLayers) {
+ msg << "\n\t- " << el.ToString();
+ }
+ throw std::logic_error(msg.str());
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetupInitializer::Reset()
+{
+ fFieldFactory.Reset();
+ fMaterialLayers.clear();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetupInitializer::SetTargetProperties(double x, double y, double z, double fieldInitStep,
+ const MaterialMap& material)
+{
+ fTarget.SetX(x);
+ fTarget.SetY(y);
+ fTarget.SetZ(z);
+ fTarget.SetMaterial(material);
+ fFieldFactory.SetTarget(x, y, z);
+ fFieldFactory.SetStep(fieldInitStep);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetupInitializer::Store(const Setup<double>& setup, const std::string& fileName)
+{
+ std::ofstream ofs(fileName, std::ios::binary);
+ if (!ofs) {
+ std::stringstream msg;
+ msg << "kf::SetupInitializer::Store: failed openning file \"" << fileName << "\" to store the setup";
+ throw std::runtime_error(msg.str());
+ }
+ boost::archive::binary_oarchive oa(ofs);
+ oa << setup;
+}
diff --git a/algo/kf/core/geo/KfSetupInitializer.h b/algo/kf/core/geo/KfSetupInitializer.h
new file mode 100644
index 0000000000000000000000000000000000000000..7fb84e9f6bd413d0da7eafe3c2393eabf79dff46
--- /dev/null
+++ b/algo/kf/core/geo/KfSetupInitializer.h
@@ -0,0 +1,169 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file KfSetupInitializer.h
+/// @brief A base KF-Setup initialization class (header)
+/// @since 28.03.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once // include this header only once per compilation unit
+
+#include "KfSetup.h"
+
+#include <boost/archive/binary_iarchive.hpp>
+#include <boost/archive/binary_oarchive.hpp>
+
+#include <fstream>
+#include <optional>
+#include <sstream>
+#include <string>
+#include <vector>
+
+namespace cbm::algo::kf
+{
+ /// \class SetupInitializer
+ /// \brief Creates a valid initialized Setup instance
+ class SetupInitializer {
+ public:
+ /// \brief Default constructor
+ SetupInitializer() = default;
+
+ /// \brief Copy constructor
+ SetupInitializer(const SetupInitializer&) = delete;
+
+ /// \brief Move constructor
+ SetupInitializer(SetupInitializer&&) = delete;
+
+ /// \brief Destructor
+ virtual ~SetupInitializer() = default;
+
+ /// \brief Copy assignment operator
+ SetupInitializer& operator=(const SetupInitializer&) = delete;
+
+ /// \brief Move assignment operator
+ SetupInitializer& operator=(SetupInitializer&&) = delete;
+
+ /// \brief Adds magnetic field 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 AddFieldSliceReference(double halfSizeX, double halfSizeY, double refZ)
+ {
+ fFieldFactory.AddSliceReference(halfSizeX, halfSizeY, refZ);
+ }
+
+ /// \brief Adds material layer
+ /// \param material Material layer
+ void AddMaterial(const MaterialMap& material);
+
+ /// \brief Virtual function, which initializes a KF-setup
+ /// \return true Initialization was successful
+ /// \return false Initialization failed
+ ///
+ /// This function can be re-implemented for a particular use-case...
+ 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;
+
+ /// \brief Resets the instance
+ void Reset();
+
+ /// \brief Sets magnetic field function
+ /// \param fieldFn Magnetic field function
+ /// \param fieldType Magnetic field type
+ void SetFieldFunction(const FieldFn_t& fieldFn, EFieldType fieldType)
+ {
+ fFieldFactory.SetFieldFunction(fieldFn);
+ fFieldFactory.SetFieldType(fieldType);
+ }
+
+ /// \brief Sets target initialization properties
+ /// \param x Target x-coordinate [cm]
+ /// \param y Target y-coordinate [cm]
+ /// \param z Target z-coordinate [cm]
+ /// \param fieldInitStep Step between the nodal points in the target field region initialization [cm]
+ /// \param material Target material layer
+ void SetTargetProperties(double x, double y, double z, double fieldInitStep, const MaterialMap& material);
+
+ /// \brief Stores a serialized setup to a file
+ /// \param setup A reference to the Setup (NOTE: a double-Setup must be passed explicitly)
+ /// \param fileName Output file name
+ ///
+ /// Only a Setup with T = double can be serialized, so if a particular setup instance has another floating-point
+ /// type, it will be converted to the double-version. In the serialization only the interpolated magnetic field
+ /// variant is stored, the original field version will be set to nullopt during the loading from the file.
+ static void Store(const Setup<double>& setup, const std::string& fileName);
+
+ /// \brief Loads a serialized setup from a file
+ /// \tparam T Underlying floating-point type for the output setup object
+ /// \param fileName Input file name
+ /// \throw std::runtime_error If the fileName cannot be opened
+ template<typename T>
+ static Setup<T> Load(const std::string& fileName);
+
+ private:
+ std::set<MaterialMap> fMaterialLayers{}; ///< Set of material maps
+ FieldFactory fFieldFactory; ///< Instance of field factory
+ Target<double> fTarget; ///< Field to keep target properties
+ bool fbPropertiesInitialized{false}; ///< Flag: if all the necessary fields are initialized
+ };
+
+
+ // ********************************
+ // ** Template method definition **
+ // ********************************
+
+ // ---------------------------------------------------------------------------------------------------------------------
+ //
+ template<typename T>
+ Setup<T> SetupInitializer::MakeSetup(bool bProvideOrigField) const
+ {
+ Setup<T> setup;
+ // Target initialization
+ setup.SetTarget(fTarget);
+
+ // Field initialization
+ setup.SetField(fFieldFactory.MakeField<double, EFieldMode::Intrpl>());
+ if (bProvideOrigField) {
+ setup.SetField(fFieldFactory.MakeField<double, EFieldMode::Orig>());
+ }
+
+ // Material layers initialization
+ for (const auto& material : fMaterialLayers) {
+ setup.AddMaterial(material);
+ }
+
+ return setup;
+ }
+
+ // -------------------------------------------------------------------------------------------------------------------
+ //
+ template<typename T>
+ Setup<T> SetupInitializer::Load(const std::string& fileName)
+ {
+ Setup<double> setup;
+ std::ifstream ifs(fileName, std::ios::binary);
+ if (!ifs) {
+ std::stringstream msg;
+ msg << "kf::SetupInitializer::Load: intput setup file \"" << fileName << "\" was not found";
+ throw std::runtime_error(msg.str());
+ }
+ try {
+ boost::archive::binary_iarchive ia(ifs);
+ ia >> setup;
+ }
+ catch (const std::exception& err) {
+ std::stringstream msg;
+ msg << "kf::SetupInitializer::Load: input setup file \"" << fileName
+ << "\" has inconsistent format or was "
+ "corrupted. The exception message: "
+ << err.what();
+ throw std::runtime_error(msg.str());
+ }
+ return Setup<T>(setup);
+ }
+} // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfTarget.cxx b/algo/kf/core/geo/KfTarget.cxx
index 1df0856d6dee7072171f69d1cc9b84afad73e81e..6f4864d1253ea565f5845588e2a29941f4007a41 100644
--- a/algo/kf/core/geo/KfTarget.cxx
+++ b/algo/kf/core/geo/KfTarget.cxx
@@ -39,13 +39,14 @@ void Target<T>::SetMaterial(const MaterialMap& material)
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T>
-std::string Target<T>::ToString(int indentLevel) const
+std::string Target<T>::ToString(int indentLevel, int verbose) const
{
constexpr char IndentChar = '\t';
std::stringstream msg;
std::string indent(indentLevel, IndentChar);
- msg << indent << "position: {" << fX << ", " << fY << ", " << fZ << "} [cm]\n";
- msg << indent << "material: " << fMaterial.ToString();
+ auto Cnv = [&](const auto& val) { return utils::simd::Cast<T, Literal_t<T>>(val); }; // alias for the conversion fn
+ msg << indent << "position: {" << Cnv(fX) << ", " << Cnv(fY) << ", " << Cnv(fZ) << "} [cm]\n";
+ msg << indent << "material:\n" << indent << IndentChar << " " << fMaterial.ToString(verbose);
return msg.str();
}
diff --git a/algo/kf/core/geo/KfTarget.h b/algo/kf/core/geo/KfTarget.h
index 60734745845fab50562d28b8841e832a7657aa01..51ab42ed1c2db490cfd251d647e8f38399cf24b3 100644
--- a/algo/kf/core/geo/KfTarget.h
+++ b/algo/kf/core/geo/KfTarget.h
@@ -52,21 +52,7 @@ namespace cbm::algo::kf
~Target() = default;
/// \brief Copy assignment operator
- /// \tparam I Underlying floating-point type of the source
- template<typename I>
- Target& operator=(const Target<I>& other)
- {
- if constexpr (std::is_same_v<I, T>) {
- if (this == &other) {
- return *this;
- }
- }
- fMaterial = other.fMaterial;
- fX = utils::simd::Cast<I, T>(other.fX);
- fY = utils::simd::Cast<I, T>(other.fY);
- fZ = utils::simd::Cast<I, T>(other.fZ);
- return *this;
- }
+ Target& operator=(const Target& other) = default;
/// \brief Gets x-coordinate of the nominal target center
const T& GetX() const { return fX; }
@@ -98,7 +84,8 @@ namespace cbm::algo::kf
/// \brief String representation of the class
/// \param indentLevel Indent level of the string output
- std::string ToString(int indentLevel = 0) const;
+ /// \param verbose Verbosity level
+ std::string ToString(int indentLevel = 0, int vebose = 1) const;
private:
/// \brief Serialization method
diff --git a/algo/kf/core/utils/KfUtils.h b/algo/kf/core/utils/KfUtils.h
index 2f4c277f46804f90ed57b18b623ea0f16f17b41e..7d88186ec1d38e90669d218ebabb6927e7da6a0b 100644
--- a/algo/kf/core/utils/KfUtils.h
+++ b/algo/kf/core/utils/KfUtils.h
@@ -247,6 +247,26 @@ namespace cbm::algo::kf::utils::simd
return static_cast<DataOut>(val);
}
+ template<typename DataT, typename DataOut, size_t N>
+ inline std::array<DataOut, N> Cast(const std::array<DataT, N>& arr, size_t i)
+ {
+ std::array<DataOut, N> res;
+ for (size_t iEl = 0; iEl < arr.size(); ++iEl) {
+ res[iEl] = Cast<DataT, DataOut>(arr[iEl], i);
+ }
+ return res;
+ }
+
+ template<typename DataT, typename DataOut, size_t N>
+ inline std::array<DataOut, N> Cast(const std::array<DataT, N>& arr)
+ {
+ std::array<DataOut, N> res;
+ for (size_t iEl = 0; iEl < arr.size(); ++iEl) {
+ res[iEl] = Cast<DataT, DataOut>(arr[iEl]);
+ }
+ return res;
+ }
+
/// \brief Converts a value of type DataT at a specific index to type DataOut
/// \details This specialization extracts the element at the specified index of the input SIMD float vector (fvec) and returns it as a float
/// \param val Input value of type fvec to be converted
@@ -269,24 +289,6 @@ 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/kf/tools/CMakeLists.txt b/algo/kf/tools/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3a072d4a7390a40f9c69e21ecf66a206c2e931bb
--- /dev/null
+++ b/algo/kf/tools/CMakeLists.txt
@@ -0,0 +1,61 @@
+set(INCLUDE_DIRECTORIES
+ ${CMAKE_CURRENT_SOURCE_DIR}
+)
+
+set(SRCS
+ ${CMAKE_CURRENT_SOURCE_DIR}/KfMaterialMapCreator.cxx
+)
+
+SET_SOURCE_FILES_PROPERTIES(${SRCS} PROPERTIES COMPILE_FLAGS "-O3")
+
+
+If(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
+ ADD_DEFINITIONS(-Wall -Wextra -Wsign-promo -Wctor-dtor-privacy -Wreorder -Wno-deprecated -Wno-parentheses) # -Weffc++ -Wnon-virtual-dtor -Woverloaded-virtual -Wold-style-cast : wait for other parts of cbmroot\root.
+Else()
+ ADD_DEFINITIONS(-Wall -Wextra -Wsign-promo -Wno-pmf-conversions -Wctor-dtor-privacy -Wreorder -Wno-deprecated -Wstrict-null-sentinel -Wno-non-template-friend -Wno-parentheses -Wmissing-field-initializers) # -Weffc++ -Wnon-virtual-dtor -Woverloaded-virtual -Wold-style-cast : wait for other parts of cbmroot\root.
+EndIf()
+
+add_library(KfTools SHARED ${SRCS})
+
+#list(APPEND HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/utils/KfSimd.h)
+
+target_include_directories(KfTools
+ PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
+)
+
+#target_compile_definitions(KfTools PUBLIC NO_ROOT)
+
+target_link_libraries(KfTools
+ PUBLIC FairRoot::Base
+ FairLogger::FairLogger
+ ROOT::Core
+ ROOT::Gpad
+ ROOT::Graf
+ ROOT::Physics
+ fmt::fmt
+ PRIVATE KfCore
+ FairRoot::GeoBase
+ FairRoot::ParBase
+ ROOT::EG
+ ROOT::Geom
+ ROOT::Graf3d
+ ROOT::MathCore
+ ROOT::Matrix
+ ROOT::Minuit
+ ROOT::RIO
+ )
+
+install(TARGETS KfTools 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(
+ FILES
+ KfMaterialMapCreator.h
+ DESTINATION
+ include/
+)
diff --git a/algo/kf/tools/KfMaterialMapCreator.cxx b/algo/kf/tools/KfMaterialMapCreator.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b138863da49ce115fc504c84a900a4d9ec9fc171
--- /dev/null
+++ b/algo/kf/tools/KfMaterialMapCreator.cxx
@@ -0,0 +1,407 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file KfMaterialMapCreator.cxx
+/// \brief Utility to generate material budget map from the TGeoNavigator representation of the Setup (implementation)
+/// \author Sergey Gorbunov <se.gorbunov@gsi.de>
+/// \author Sergei Zharko <s.zharko@gsi.de>
+/// \date 29.08.2024
+
+#include "KfMaterialMapCreator.h"
+
+#include "KfMaterialMap.h"
+#include "TGeoManager.h"
+#include "TGeoMedium.h"
+#include "TGeoNavigator.h"
+#include "TGeoNode.h"
+#include "TGeoPhysicalNode.h"
+#include "TGeoVoxelFinder.h"
+
+#include <Logger.h>
+
+#include <sstream>
+#include <thread>
+
+using cbm::algo::kf::MaterialMap;
+using kf::tools::MaterialMapCreator;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+MaterialMapCreator::MaterialMapCreator(int verbose) : fVerbose(verbose)
+{
+ // constructor
+
+ // save the current number of threads to restore it when the work is done
+ fNthreadsOld = gGeoManager->GetMaxThreads();
+
+ // set defaut n threads to 1
+
+ fNthreads = 1;
+
+ std::stringstream msg;
+ msg << "\n------------------------------------------------------------------------------------";
+ msg << "\n Material budget map creator";
+ msg << '\n';
+ msg << "\n !! To run it with the full speed, set: \"export OMP_NUM_THREADS=<n CPUs>\" before running !!";
+ msg << '\n';
+ // if possible, read the allowed n threads from the environment variable
+ {
+ const char* env = std::getenv("OMP_NUM_THREADS");
+ if (env) {
+ fNthreads = TString(env).Atoi();
+ msg << "\n found environment variable OMP_NUM_THREADS = \"" << env << "\", read as integer: " << fNthreads;
+ }
+ }
+
+ // ensure that at least one thread is set
+
+ if (fNthreads < 1) {
+ fNthreads = 1;
+ }
+
+ msg << "\n Maps will be created using " << fNthreads << " CPU threads";
+ msg << '\n';
+ msg << "\n It might crash because of a ROOT bug. In this case do one of the following: ";
+ msg << '\n';
+ msg << "\n - fix illegal overlaps in the geometry via gGeoManager->CheckOverlaps()";
+ msg << "\n - run with CbmL1::SetSafeMaterialInitialization() option. ";
+ msg << "\n It fixes the crash, but slows down the initialization significantly.";
+ msg << "\n - manually disable voxel finders for problematic volumes in CaToolsMaterialMapCreator.cxx";
+ msg << '\n';
+ msg << "\n------------------------------------------------------------------------------------";
+ LOG_IF(info, fVerbose > 0) << msg.str();
+
+ InitThreads();
+
+ //
+ // It looks like TGeoVoxelFinder contains a bug and therefore sometimes produces a crash.
+ //
+ // Currently, the crash only happens on some of TGeoVolumeAssembly volumes
+ // and only when an (inproper?) alignment matrices are applied.
+ //
+ // Here we try to catch this case and disable those finders.
+ //
+ // Disabling all the voxel finders leads to a very long initialization, try to avoid it!
+ //
+ // TODO: make a bug report to the ROOT team; remove the following code once the bug is fixed.
+ //
+
+ if (fDoSafeInitialization) { // disable all voxel finders
+ TObjArray* volumes = gGeoManager->GetListOfVolumes();
+ if (volumes) {
+ for (int iv = 0; iv < volumes->GetEntriesFast(); iv++) {
+ TGeoVolume* vol = dynamic_cast<TGeoVolume*>(volumes->At(iv));
+ if (!vol) {
+ continue;
+ }
+ TGeoVoxelFinder* finder = vol->GetVoxels();
+ if (finder) {
+ finder->SetInvalid();
+ }
+ }
+ }
+ }
+ else { // disable some voxel finders in some cases
+
+ // check if any alignment was applied
+
+ bool isAlignmentApplied = false;
+ {
+ TObjArray* nodes = gGeoManager->GetListOfPhysicalNodes();
+ if (nodes) {
+ for (int in = 0; in < nodes->GetEntriesFast(); in++) {
+ TGeoPhysicalNode* node = dynamic_cast<TGeoPhysicalNode*>(nodes->At(in));
+ if (!node) continue;
+ if (node->IsAligned()) {
+ isAlignmentApplied = true;
+ break;
+ }
+ }
+ }
+ }
+
+ // disable potentially problematic voxel finders
+
+ if (isAlignmentApplied) {
+ TObjArray* volumes = gGeoManager->GetListOfVolumes();
+ if (volumes) {
+ for (int iv = 0; iv < volumes->GetEntriesFast(); iv++) {
+ TGeoVolumeAssembly* vol = dynamic_cast<TGeoVolumeAssembly*>(volumes->At(iv));
+ if (!vol) {
+ continue;
+ }
+ TGeoVoxelFinder* finder = vol->GetVoxels();
+ if (finder) {
+ finder->SetInvalid();
+ }
+ }
+ }
+ }
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+MaterialMapCreator::~MaterialMapCreator()
+{
+ // destructor
+
+ CleanUpThreads();
+
+ // delete the navigators, created in this class
+
+ for (auto i = fNavigators.begin(); i != fNavigators.end(); i++) {
+ gGeoManager->RemoveNavigator(*i);
+ }
+
+ // once TGeoManager is swithched in multithreaded mode, there is no way to make it non-mltithreaded again
+ // therefore we should set SetMaxThreads( >=0 )
+
+ if (fNthreadsOld > 0) {
+ gGeoManager->SetMaxThreads(fNthreadsOld);
+ }
+ else {
+ gGeoManager->SetMaxThreads(1);
+ }
+
+ // ensure that the default navigator exists
+
+ GetCurrentNavigator(-1);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+TGeoNavigator* MaterialMapCreator::GetCurrentNavigator(int iThread)
+{
+ // Get navigator for current thread, create it if it doesnt exist.
+ // Produce an error when anything goes wrong
+
+ TGeoNavigator* navi = gGeoManager->GetCurrentNavigator();
+ if (!navi) {
+
+ navi = gGeoManager->AddNavigator();
+
+ if (iThread >= 0) {
+ fNavigators.push_back(navi);
+ }
+
+ if (navi != gGeoManager->GetCurrentNavigator()) {
+ LOG(fatal) << "ca::tools::MaterialMapCreator: unexpected behavior of TGeoManager::AddNavigator() !!";
+ }
+ }
+
+ if (!navi) {
+ LOG(fatal) << "ca::tools::MaterialMapCreator: can not find / create TGeoNavigator for thread " << iThread;
+ }
+
+ return navi;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void MaterialMapCreator::CleanUpThreads()
+{
+ // clean up thread data in TGeoManager
+
+ gGeoManager->ClearThreadsMap();
+
+ // create a default navigator for multithreaded mode
+ // otherwise gGeoManager->GetCurrentNavigator() returns nullptr that might produces crashes in other code
+
+ GetCurrentNavigator(-1);
+}
+
+void MaterialMapCreator::InitThreads()
+{
+ // (re)initialise the number of threads in TGeoManager
+
+ // reserve enough threads in TGeoManager
+ if (fNthreads > gGeoManager->GetMaxThreads()) {
+ gGeoManager->SetMaxThreads(fNthreads);
+ // in case the number of threads is truncated by TGeoManager (must not happen)
+ fNthreads = gGeoManager->GetMaxThreads();
+ if (fNthreads < 1) {
+ fNthreads = 1;
+ }
+ }
+ CleanUpThreads();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+MaterialMap MaterialMapCreator::GenerateMaterialMap(double zRef, double zMin, double zMax, double xyMax, int nBinsDim)
+{
+ if (fDoRadialProjection) {
+ if (fTargetZ >= zRef - 0.1) {
+ LOG(warn)
+ << "kf::tools::MaterialMapCreator: the material reference z-coordinate must be at least 0.1 cm downstream "
+ "the target. Shifting the reference z-coordinate to the lower limit: target Z = "
+ << fTargetZ << ", material reference z = " << zRef;
+ }
+ zRef = std::max(zRef, fTargetZ + 1.);
+ zMin = std::max(fTargetZ, zMin);
+ }
+
+ if (!(zMin <= zRef && zRef <= zMax)) {
+ LOG(fatal) << "kf::tools::MaterialMapCreator: the material parameters are inconsistent. It must be: zMin " << zMin
+ << " <= zRef " << zRef << " <= zMax " << zMax;
+ }
+
+ MaterialMap matBudget(nBinsDim, xyMax, zRef, zMin, zMax);
+ double binSizeX = 2. * xyMax / nBinsDim;
+ double binSizeY = 2. * xyMax / nBinsDim;
+
+ // call it every time. for the case the number of threads was reset in some other code
+ InitThreads();
+
+ long nThreadCrosses[fNthreads];
+ long nThreadRays[fNthreads];
+ long nThreadRaysExpected = nBinsDim * nBinsDim * fNraysBinPerDim * fNraysBinPerDim / fNthreads;
+ for (int i = 0; i < fNthreads; i++) {
+ nThreadCrosses[i] = 0;
+ nThreadRays[i] = 0;
+ }
+
+ LOG_IF(info, fVerbose > 1) << "kf::tools::MaterialMapCreator: Generate material map using " << fNthreads
+ << " threads..";
+
+ auto naviThread = [&](int iThread) {
+ // create a navigator for this thread
+ TGeoNavigator* navi = GetCurrentNavigator(iThread);
+
+ for (int iBinX = iThread; iBinX < nBinsDim; iBinX += fNthreads) {
+ for (int iBinY = 0; iBinY < nBinsDim; ++iBinY) {
+ double radThick = 0.;
+ int nRays = 0;
+ double d = 1. / (fNraysBinPerDim);
+ for (int iStepX = 0; iStepX < fNraysBinPerDim; iStepX++) {
+ for (int iStepY = 0; iStepY < fNraysBinPerDim; iStepY++) {
+
+ // ray position at zRef
+ double rayX = -xyMax + (iBinX + d * (iStepX + 0.5)) * binSizeX;
+ double rayY = -xyMax + (iBinY + d * (iStepY + 0.5)) * binSizeY;
+
+ // ray slopes
+ double tx = 0.;
+ double ty = 0.;
+ double t = 1.;
+ if (fDoRadialProjection) {
+ tx = rayX / (zRef - fTargetZ);
+ ty = rayY / (zRef - fTargetZ);
+ t = sqrt(1. + tx * tx + ty * ty);
+ }
+
+ // ray position at zMin
+ double z = zMin;
+ double x = rayX + tx * (z - zRef);
+ double y = rayY + ty * (z - zRef);
+ LOG_IF(info, fVerbose > 2) << "ray at " << x << " " << y << " " << z << " zMin " << zMin << " zMax "
+ << zMax;
+
+ TGeoNode* node = navi->InitTrack(x, y, z, tx / t, ty / t, 1. / t);
+
+ bool doContinue = 1;
+ for (int iStep = 0; doContinue; iStep++) {
+
+ if (!node) {
+ // may happen when tracing outside of the CBM volume -> produce a warning
+ LOG(warning) << "kf::tools::MaterialMapCreator: TGeoNavigator can not find the geo node";
+ break;
+ }
+
+ TGeoMedium* medium = node->GetMedium();
+ if (!medium) {
+ LOG(fatal) << "kf::tools::MaterialMapCreator: TGeoNavigator can not find the geo medium";
+ }
+
+ TGeoMaterial* material = medium->GetMaterial();
+ if (!material) {
+ LOG(fatal) << "kf::tools::MaterialMapCreator: TGeoNavigator can not find the geo material";
+ }
+
+ double radLen = material->GetRadLen();
+
+ if (radLen < kMinRadLength) { // 0.5612 is rad. length of Lead at normal density
+ LOG(fatal) << "kf::tools::MaterialMapCreator: failed assertion! "
+ << " TGeoNavigator returns a suspicious material with an unrealistic "
+ << "radiation length of " << radLen << " cm. "
+ << " The allowed minimum is 0.56 cm (Lead has 0.5612 cm). Check your material! "
+ "Modify this assertion if necessary."
+ << " TGeoNode \"" << node->GetName() << "\", TGeoMedium \"" << medium->GetName()
+ << "\", TGeoMaterial \"" << material->GetName() << "\"";
+ }
+
+ // move navi to the next material border
+ node = navi->FindNextBoundaryAndStep(); //5., kTRUE);
+ nThreadCrosses[iThread]++;
+
+ LOG_IF(info, fVerbose > 2) << " RL " << radLen << " next pt " << navi->GetCurrentPoint()[0] << " "
+ << navi->GetCurrentPoint()[1] << " " << navi->GetCurrentPoint()[2];
+
+ double zNew = navi->GetCurrentPoint()[2];
+ if (zNew >= zMax) {
+ zNew = zMax;
+ doContinue = 0;
+ }
+
+ //if (zNew < z) {
+ //LOG(info) << " z curr " << z << " z new " << zNew << " dz " << zNew - z ; }
+ double dz = zNew - z;
+ if (dz < 0.) {
+ if (iStep == 0 && dz > -1.e-8) {
+ // when the ray starts exactly at the volume border, first dz might become negative,
+ // probably due to some roundoff errors. So we let it be a bit negative for the first intersection.
+ dz = 0.;
+ }
+ else {
+ LOG(fatal)
+ << "kf::tools::MaterialMapCreator: TGeoNavigator propagates the ray upstream. Something is wrong."
+ << " z old " << z << " z new " << zNew << ", dz " << dz;
+ }
+ }
+ radThick += dz / radLen;
+ z = zNew;
+ }
+
+ nRays++;
+ nThreadRays[iThread]++;
+ if (fVerbose > 2 && nThreadRays[iThread] % 1000000 == 0) {
+ LOG(info) << "kf::tools::MaterialMapCreator: report from thread " << iThread << ": material map is "
+ << 100. * nThreadRays[iThread] / nThreadRaysExpected << " \% done";
+ }
+ }
+ }
+ radThick /= nRays;
+ LOG_IF(info, fVerbose > 2) << " radThick " << radThick;
+ //doPrint = (radThick > 0.01);
+ matBudget.SetRadThickBin(iBinX, iBinY, radThick);
+ } // iBinX
+ } // iBinY
+ };
+
+ std::vector<std::thread> threads(fNthreads);
+
+ // run n threads
+ for (int i = 0; i < fNthreads; i++) {
+ threads[i] = std::thread(naviThread, i);
+ }
+
+ // wait for the threads to finish
+ for (auto& th : threads) {
+ th.join();
+ }
+
+ CleanUpThreads();
+ long nRays = 0;
+ long nCrosses = 0;
+ for (int i = 0; i < fNthreads; i++) {
+ nRays += nThreadRays[i];
+ nCrosses += nThreadCrosses[i];
+ }
+
+ LOG_IF(info, fVerbose > 0) << "kf::tools::MaterialMapCreator: shooted " << nRays << " rays. Each ray crossed "
+ << nCrosses * 1. / nRays << " material boundaries in average";
+
+ return matBudget;
+}
diff --git a/algo/kf/tools/KfMaterialMapCreator.h b/algo/kf/tools/KfMaterialMapCreator.h
new file mode 100644
index 0000000000000000000000000000000000000000..4814228e60328e7c59d11b975496450740378ac0
--- /dev/null
+++ b/algo/kf/tools/KfMaterialMapCreator.h
@@ -0,0 +1,99 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file KfMaterialMapCreator.h
+/// \brief Utility to generate material budget map from the TGeoNavigator representation of the Setup (implementation)
+/// \author Sergey Gorbunov <se.gorbunov@gsi.de>
+/// \author Sergei Zharko <s.zharko@gsi.de>
+/// \date 29.08.2024
+
+#pragma once
+
+#include "Rtypes.h"
+#include "TObject.h"
+
+class TGeoNavigator;
+namespace cbm::algo::kf
+{
+ class MaterialMap;
+}
+
+namespace kf::tools
+{
+ /// \class MaterialMapCreator
+ /// \brief An utility class to create a material budget map from the TGeo
+ class MaterialMapCreator {
+ public:
+ /// \brief Constructor from parameters
+ /// \param verbose Verbosity level (0 - only warning/error output)
+ MaterialMapCreator(int verbose = 0);
+
+ /// \brief Destructor
+ ~MaterialMapCreator();
+
+ /// \brief Project rays radially from the targetZ througth the XY-bins at a reference z.
+ /// \param targetZ z-coordinate of the target [cm]
+ void SetDoRadialProjection(double targetZ)
+ {
+ fDoRadialProjection = true;
+ fTargetZ = targetZ;
+ }
+
+ /// \brief Project rays horisontally along the Z axis (default)
+ void SetDoHorisontalProjection() { fDoRadialProjection = false; }
+
+ /// \brief Shoots nRaysDim * nRaysDim rays for each bin in the map
+ void SetNraysPerDim(int nRaysDim) { fNraysBinPerDim = (nRaysDim > 0) ? nRaysDim : 1; }
+
+
+ /// \brief Generates a material budget map
+ /// \param zRef Reference z-coordinate of the material layer [cm]
+ /// \param zMin z-coordinate of the material layer lower boundary [cm]
+ /// \param zMax z-coordinate of the material layer upper boundary [cm]
+ /// \param xyMax Transverse size of the material layer [cm]
+ /// \param nBinsDim Number of bins in the x(y) axis
+ ///
+ /// generate a material map, collecting the material between [zMin, zMax]
+ /// with radial rays from (0,0,targetZ) through the XY-bins at z == zRef.
+ ///
+ /// It creates a map with [nBinsDim x nBinsDim] bins, of a size of [+-xyMax, +-xyMax]
+ /// shooting fNraysBinPerDim x fNraysBinPerDim through each bin
+ ///
+ /// The calculated radiation thickness is a projection of the rad.thickness along the ray onto the Z axis.
+ /// RadThick = sum of (dZ / radLength) over ray trajectory pieces
+ ///
+ /// When doRadialProjection==false the rays are shoot horizontally along the Z axis
+ ///
+ cbm::algo::kf::MaterialMap GenerateMaterialMap(double zRef, double zMin, double zMax, double xyMax, int nBinsDim);
+
+ /// \brief Enables safe mode of the material initialization
+ void SetSafeMaterialInitialization(bool val = true) { fDoSafeInitialization = val; }
+
+ private:
+ /// \brief Initializes the necessary amount of threads in TGeoManager
+ void InitThreads();
+
+ /// \brief Cleans up the TGeoManager: threadIds, create a default navigator
+ void CleanUpThreads();
+
+ /// \brief Gets navigator for current thread, creates it if it does not exist
+ /// \param iThread Index of the thread
+ /// \throw std::runtime_error If the TGeoNavigator is not found
+ TGeoNavigator* GetCurrentNavigator(int iThread);
+
+ private:
+ static constexpr double kMinRadLength = 0.3; ///< Minimal radiational length allowed [cm]
+
+ std::vector<TGeoNavigator*> fNavigators{}; ///< list of created navigators
+ double fTargetZ{0.}; ///< z of the target for the radial projection
+ int fNthreadsOld{0}; ///< number of threads in TGeoManager before the helper was created
+ int fNthreads{0}; ///< number of threads
+ int fNraysBinPerDim{3}; ///< shoot fNraysBinPerDim * fNraysBinPerDim rays in each map bin
+ int fVerbose{0}; ///< verbosity level
+ bool fDoRadialProjection{false}; ///< if project rays horizontally along the Z axis (special mode)
+ bool fDoSafeInitialization{false}; ///< performs slow but safe initialization
+ ///< to get around the crashes in TGeoVoxelFinder
+ };
+
+} // namespace kf::tools
diff --git a/algo/kf/tools/README.md b/algo/kf/tools/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..0d02858a4918eb687e43d43c9b24bfd00531bf0c
--- /dev/null
+++ b/algo/kf/tools/README.md
@@ -0,0 +1,10 @@
+# KfTools Library
+
+## Introduction
+
+The KfTools contains additional utilities for the KF-framework, which cannot be included into
+the cbm::algo namespace, but in the same time are experiment-independent.
+
+## Library contents
+
+### **Material map generator**
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index f2cdc2a69df550c0461db31a41a437ced696ca0a..9b02fa975f02e4d404c1577c5c1e46539e8e5ecc 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -168,6 +168,7 @@ set(PRIVATE_DEPENDENCIES
CbmRecoSts
CbmSimSteer
CbmRecoBase
+ KfCbm
KFParticle
external::yaml-cpp
FairRoot::GeoBase
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 203dd07412599892a7f7a4afd1f02c360013a029..91a5616ade754a5a71d974bc880970e053b37bd3 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -37,6 +37,7 @@
#include "CaToolsDebugger.h"
#include "CaToolsField.h"
#include "CbmEvent.h"
+#include "CbmKfTrackingSetupInitializer.h"
#include "CbmMCDataObject.h"
#include "CbmStsFindTracks.h"
#include "CbmStsHit.h"
@@ -72,6 +73,7 @@ using cbm::algo::ca::Parameters;
using cbm::ca::MCModule;
using cbm::ca::TimeSliceReader;
using cbm::ca::tools::MaterialHelper;
+using cbm::kf::TrackingSetupInitializer;
ClassImp(CbmL1);
@@ -239,6 +241,7 @@ try {
CheckDetectorPresence();
+
// *****************************
// ** **
// ** GEOMETRY INITIALIZATION **
@@ -547,6 +550,81 @@ try {
}
fpAlgo->Init(fTrackingMode);
+ if constexpr (0) {
+ // **** FEATURE: KF-SETUP INITIALIZATION *****
+ // Creating a setup
+ auto pSetupInitializer = std::make_unique<TrackingSetupInitializer>();
+ pSetupInitializer->Use(fUseMVD, fUseSTS, fUseMUCH, fUseTRD, fUseTOF);
+ pSetupInitializer->InitProperties();
+ auto trackerSetup = pSetupInitializer->MakeSetup<double>(/*provide original field*/ true);
+ LOG(info) << "!!!!\n!!!!\n!!!!\n!!!! KF SETUP:\n" << trackerSetup.ToString(0) << "\n!!!!\n!!!!\n!!!!";
+
+ // Storing setup
+ TrackingSetupInitializer::Store(trackerSetup, "./trackerSetup.geo.kf.bin");
+
+ // Reading setup (now in fvec):
+ auto trackerSetupVec = TrackingSetupInitializer::Load<fvec>("./trackerSetup.geo.kf.bin");
+ LOG(info) << "!!!!\n!!!!\n!!!!\n!!!! KF SETUP (from file):\n"
+ << trackerSetupVec.ToString(1) << "\n!!!!\n!!!!\n!!!!";
+
+
+ {
+ // TEST: magnetic field
+ const auto& fld = trackerSetupVec.GetField<cbm::algo::kf::EFieldMode::Intrpl>();
+ 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);
+ auto caFldValue = fpParameters->GetStation(n.iSt).fieldSlice.GetFieldValue(n.x, n.y);
+ LOG(info) << "A: " << kfFldValue.ToString() << " --- " << caFldValue.ToString();
+ }
+ {
+ Node n{.iSt = 2, .x = -2.1, .y = 2.2};
+ auto kfFldValue = fld.GetFieldValue(n.iSt, n.x, n.y);
+ auto caFldValue = fpParameters->GetStation(n.iSt).fieldSlice.GetFieldValue(n.x, n.y);
+ LOG(info) << "B: " << kfFldValue.ToString() << " --- " << caFldValue.ToString();
+ }
+ {
+ Node n{.iSt = 5, .x = -2.1, .y = 2.2};
+ auto kfFldValue = fld.GetFieldValue(n.iSt, n.x, n.y);
+ auto caFldValue = fpParameters->GetStation(n.iSt).fieldSlice.GetFieldValue(n.x, n.y);
+ LOG(info) << "B: " << kfFldValue.ToString() << " --- " << caFldValue.ToString();
+ }
+
+ // TEST: material maps
+ {
+ Node n{.iSt = 1, .x = 2.3, .y = 1.2};
+ LOG(info) << trackerSetupVec.GetMaterial(n.iSt).GetThickness(n.x, n.y) << " -- "
+ << fpParameters->GetMaterialThickness(n.iSt, n.x, n.y);
+ }
+ {
+ Node n{.iSt = 1, .x = -2.3, .y = -1.2};
+ LOG(info) << trackerSetupVec.GetMaterial(n.iSt).GetThickness(n.x, n.y) << " -- "
+ << fpParameters->GetMaterialThickness(n.iSt, n.x, n.y);
+ }
+ {
+ Node n{.iSt = 1, .x = 2.5, .y = 0.2};
+ LOG(info) << trackerSetupVec.GetMaterial(n.iSt).GetThickness(n.x, n.y) << " -- "
+ << fpParameters->GetMaterialThickness(n.iSt, n.x, n.y);
+ }
+ {
+ Node n{.iSt = 4, .x = 5.3, .y = -2.2};
+ LOG(info) << trackerSetupVec.GetMaterial(n.iSt).GetThickness(n.x, n.y) << " -- "
+ << fpParameters->GetMaterialThickness(n.iSt, n.x, n.y);
+ }
+ {
+ Node n{.iSt = 4, .x = 0, .y = -1.2};
+ LOG(info) << trackerSetupVec.GetMaterial(n.iSt).GetThickness(n.x, n.y) << " -- "
+ << fpParameters->GetMaterialThickness(n.iSt, n.x, n.y);
+ }
+ }
+ }
+
+
// Initialize time-slice reader
fpTSReader = std::make_unique<TimeSliceReader>();
fpTSReader->SetDetector(ca::EDetectorID::kMvd, fUseMVD);
diff --git a/reco/kf/CMakeLists.txt b/reco/kf/CMakeLists.txt
index 9a5e8b40448fcddabfdd3c1ebce03f8acd373d32..ff911e51c2e38dad008f201720ebde1e36e9ff53 100644
--- a/reco/kf/CMakeLists.txt
+++ b/reco/kf/CMakeLists.txt
@@ -16,17 +16,17 @@ Else()
ADD_DEFINITIONS(-Wall -Wextra -Wsign-promo -Wno-pmf-conversions -Wctor-dtor-privacy -Wreorder -Wno-deprecated -Wstrict-null-sentinel -Wno-non-template-friend -Wno-parentheses -Wmissing-field-initializers) # -Weffc++ -Wnon-virtual-dtor -Woverloaded-virtual -Wold-style-cast : wait for other parts of cbmroot\root.
EndIf()
-add_library(CbmKf SHARED ${SRCS})
+add_library(KfCbm SHARED ${SRCS})
#list(APPEND HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/utils/KfSimd.h)
-target_include_directories(CbmKf
+target_include_directories(KfCbm
PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
)
#target_compile_definitions(KfCore PUBLIC NO_ROOT)
-target_link_libraries(CbmKf
+target_link_libraries(KfCbm
PUBLIC CbmMuchBase
CbmMvdBase
CbmStsBase
@@ -59,7 +59,7 @@ target_link_libraries(CbmKf
ROOT::RIO
)
-install(TARGETS CbmKf DESTINATION lib)
+install(TARGETS KfCbm 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")
@@ -69,6 +69,7 @@ install(TARGETS CbmKf DESTINATION lib)
install(
FILES
+ CbmKfOriginalField.h
CbmKfTarget.h
CbmKfTrackingSetupInitializer.h
DESTINATION
diff --git a/reco/kf/CbmKfOriginalField.h b/reco/kf/CbmKfOriginalField.h
new file mode 100644
index 0000000000000000000000000000000000000000..f87a18c1cf551cc476362285f581c57ee2ad4ff0
--- /dev/null
+++ b/reco/kf/CbmKfOriginalField.h
@@ -0,0 +1,55 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov[committer] */
+
+/// \file CbmKfOriginalField.h
+/// \brief Thread-safe representation of the magnetic field (header)
+/// \since 17.10.2023
+/// \author S.Gorbunov
+
+#pragma once // include this header only once per compilation unit
+
+#include "FairField.h"
+#include "FairRunAna.h"
+
+#include <mutex>
+#include <tuple>
+
+namespace cbm::kf
+{
+ /// \class OriginalField
+ /// \brief Thread-safe representation of the magnetic field in CBM
+ class OriginalField {
+ public:
+ /// \brief Returns magnetic field flux density in a spatial point
+ /// \param x x-coordinate of the point [cm]
+ /// \param y y-coordinate of the point [cm]
+ /// \param z z-coordinate of the point [cm]
+ /// \return tuple(Bx, By, Bz) of the magnetic field flux density components [kG]
+ std::tuple<double, double, double> operator()(double x, double y, double z) const
+ {
+ assert(FairRunAna::Instance());
+ assert(FairRunAna::Instance()->GetField());
+ double pos[3] = {x, y, z};
+ double B[3] = {0., 0., 0.};
+ // protect the field access
+ // TODO: make CbmField thread-safe
+ static std::mutex mymutex;
+ mymutex.lock();
+ FairRunAna::Instance()->GetField()->GetFieldValue(pos, B);
+ mymutex.unlock();
+ return std::tuple(B[0], B[1], B[2]);
+ }
+ };
+
+ /// \class ZeroField
+ class ZeroField {
+ public:
+ /// \brief Returns magnetic field flux density in a spatial point
+ /// \param x x-coordinate of the point [cm]
+ /// \param y y-coordinate of the point [cm]
+ /// \param z z-coordinate of the point [cm]
+ /// \return tuple(Bx, By, Bz) of the magnetic field flux density components [kG]
+ std::tuple<double, double, double> operator()(double, double, double) const { return std::tuple(0., 0., 0.); }
+ };
+} // namespace cbm::kf
diff --git a/reco/kf/CbmKfTarget.cxx b/reco/kf/CbmKfTarget.cxx
index a572edec2a9abc90a5807ed8c76c1b6a2e7c44b1..431e91d46aad433d316a754b234580e097f53f39 100644
--- a/reco/kf/CbmKfTarget.cxx
+++ b/reco/kf/CbmKfTarget.cxx
@@ -38,36 +38,39 @@ Target* Target::Instance()
// ---------------------------------------------------------------------------------------------------------------------
//
-void Target::Init()
+void Target::FindTargetNode(TString& targetPath, TGeoNode*& targetNode)
{
- if (!gGeoManager) {
- throw std::logic_error("cbm::kf::Target: the TGeoManager instance is not initialized on this step. Please be "
- "ensured to call the Target::Init() after the TGeoManager initialization.");
+ if (!targetNode) { // init at the top of the tree
+ targetNode = gGeoManager->GetTopNode();
+ targetPath = "/" + TString(targetNode->GetName());
}
- std::function<void(TString&, TGeoNode*)> FindTargetNode = [&](TString& targetPath, TGeoNode* targetNode) -> void {
- if (!targetNode) { // init at the top of the tree
- targetNode = gGeoManager->GetTopNode();
- targetPath = "/" + TString(targetNode->GetName());
- }
+ if (TString(targetNode->GetName()).Contains("target")) {
+ return;
+ }
- if (TString(targetNode->GetName()).Contains("target")) {
+ for (Int_t iNode = 0; iNode < targetNode->GetNdaughters(); ++iNode) {
+ TGeoNode* newNode = targetNode->GetDaughter(iNode);
+ TString newPath = targetPath + "/" + newNode->GetName();
+ FindTargetNode(newPath, newNode);
+ if (newNode) {
+ targetPath = newPath;
+ targetNode = newNode;
return;
}
+ }
+ targetPath = "";
+ targetNode = nullptr;
+}
- for (Int_t iNode = 0; iNode < targetNode->GetNdaughters(); ++iNode) {
- TGeoNode* newNode = targetNode->GetDaughter(iNode);
- TString newPath = targetPath + "/" + newNode->GetName();
- FindTargetNode(newPath, newNode);
- if (newNode) {
- targetPath = newPath;
- targetNode = newNode;
- return;
- }
- }
- targetPath = "";
- targetNode = nullptr;
- };
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void Target::Init()
+{
+ if (!gGeoManager) {
+ throw std::logic_error("cbm::kf::Target: the TGeoManager instance is not initialized on this step. Please be "
+ "ensured to call the Target::Init() after the TGeoManager initialization.");
+ }
TString targetPath;
TGeoNode* pTargetNode{nullptr};
diff --git a/reco/kf/CbmKfTarget.h b/reco/kf/CbmKfTarget.h
index 2440ad513e264b7b2bbd25075174e40661e65ad8..f9106f4c2c9af3651539ebe2ce57925f21827359 100644
--- a/reco/kf/CbmKfTarget.h
+++ b/reco/kf/CbmKfTarget.h
@@ -15,6 +15,9 @@
#include <mutex>
+class TString;
+class TGeoNode;
+
namespace cbm::kf
{
/// \class Target
@@ -55,6 +58,9 @@ namespace cbm::kf
/// \brief Target initializer
void Init();
+ /// \brief Finds a target
+ static void FindTargetNode(TString& targetPath, TGeoNode*& targetNode);
+
double fX{cbm::algo::kf::defs::Undef<double>}; ///< reference x-coordinate of the target position [cm]
double fY{cbm::algo::kf::defs::Undef<double>}; ///< reference y-coordinate of the target position [cm]
double fZ{cbm::algo::kf::defs::Undef<double>}; ///< reference z-coordinate of the target position [cm]
diff --git a/reco/kf/CbmKfTrackingSetupInitializer.cxx b/reco/kf/CbmKfTrackingSetupInitializer.cxx
index 4a418301bca9f8fe4d37442764546ae914f4a346..9650ac15884b69c2976eda0f5f5553d516acaa75 100644
--- a/reco/kf/CbmKfTrackingSetupInitializer.cxx
+++ b/reco/kf/CbmKfTrackingSetupInitializer.cxx
@@ -9,6 +9,7 @@
#include "CbmKfTrackingSetupInitializer.h"
+#include "CbmKfOriginalField.h"
#include "CbmKfTarget.h"
#include "CbmMuchTrackingInterface.h"
#include "CbmMvdTrackingInterface.h"
@@ -20,6 +21,7 @@
#include "FairRunAna.h"
#include "KfDefs.h"
#include "KfMaterialMapCreator.h"
+#include "Logger.h"
#include <functional>
@@ -41,12 +43,11 @@ try {
// Magnetic field initialization
if (FairRunAna::Instance()->GetField()) {
- this->SetFieldFunction(
- [](const double(&r)[3], double* B) { FairRunAna::Instance()->GetField()->GetFieldValue(r, B); },
- EFieldType::Normal);
+ this->SetFieldFunction(OriginalField(), EFieldType::Normal);
}
else {
- this->SetFieldFunction(cbm::algo::kf::defs::ZeroFieldFn, EFieldType::Null);
+ //this->SetFieldFunction(cbm::algo::kf::defs::ZeroFieldFn, EFieldType::Null);
+ this->SetFieldFunction(ZeroField(), EFieldType::Null);
}
// Tracking station property initialization
@@ -98,13 +99,18 @@ try {
double xyMax{1.3 * pTarget->GetRmax()};
int nBins{std::min(static_cast<int>(std::ceil(2. * xyMax / kMatCreatorPitch)), kMatCreatorMaxNbins)};
assert(nBins > 0);
- double zFirst{pTarget->GetZ() + pTarget->GetDz()}; // z-coordinate of the target material lower limit
+ double zFirst{pTarget->GetZ() - pTarget->GetDz()}; // z-coordinate of the target material lower limit
LOG(info) << "- collecting target material (z = " << zFirst << " - " << zLast << ", size = " << xyMax << ")\n";
this->SetTargetProperties(pTarget->GetX(), pTarget->GetY(), pTarget->GetZ(), kTargFieldInitStep,
materialCreator.GenerateMaterialMap(pTarget->GetZ(), zFirst, zLast, xyMax, nBins));
}
// Material and field layers initialization
+
+ for (const auto& layer : geoStations) {
+ LOG(info) << " >>>>> LAYER >>>>> " << layer.zRef;
+ }
+
for (auto layerIt = geoStations.begin(); layerIt != geoStations.end(); ++layerIt) {
double z1 = layerIt->zMax;
double z2 = z1;
@@ -121,7 +127,55 @@ try {
this->AddFieldSliceReference(layerIt->xMax, layerIt->yMax, layerIt->zRef);
zLast = zNew;
}
+
+ if constexpr (0) {
+ // Estimate acceptance
+ std::vector<cbm::algo::kf::MaterialMap> mapsTest;
+ mapsTest.reserve(geoStations.size());
+ double slope{0.};
+ zLast =
+ pTarget->GetZ() + pTarget->GetDz() + kTargMaterialOffset; // z-coordinate of the target material upper limit
+ for (auto layerIt = geoStations.begin(); layerIt != geoStations.end(); ++layerIt) {
+ double xyMax = 1.3 * std::max(std::fabs(layerIt->xMax), std::fabs(layerIt->yMax));
+ slope = std::max(slope, xyMax / (layerIt->zRef - pTarget->GetZ()));
+ }
+ for (auto layerIt = geoStations.begin(); layerIt != geoStations.end(); ++layerIt) {
+ double z1 = layerIt->zMax;
+ double z2 = z1;
+ if (std::next(layerIt) != geoStations.end()) {
+ z2 = std::next(layerIt)->zMin;
+ }
+ double zNew = 0.5 * (z1 + z2);
+ double xyMax = slope * (layerIt->zRef - pTarget->GetZ());
+ //int nBins = std::min(static_cast<int>(std::ceil(2. * xyMax / kMatCreatorPitch)), kMatCreatorMaxNbins);
+ int nBins = 100;
+ assert(nBins > 0);
+
+ LOG(info) << "* collecting material between z = " << zLast << " and " << zNew << " [cm]";
+ mapsTest.push_back(materialCreator.GenerateMaterialMap(layerIt->zRef, zLast, zNew, xyMax, nBins));
+ zLast = zNew;
+ }
+ for (const auto& map : mapsTest) {
+ LOG(info) << "TEST: " << map.ToString(0, 2);
+ }
+ {
+ // Addition test
+ int iInactive{2};
+ auto& M1{mapsTest[iInactive]};
+ auto& M2{mapsTest[iInactive + 1]};
+ auto M12gen{
+ materialCreator.GenerateMaterialMap(M2.GetZref(), M1.GetZmin(), M2.GetZmax(), M2.GetXYmax(), M2.GetNbins())};
+ auto M12add{M2};
+ M12add.Add(M1);
+ LOG(info) << "MATERIAL MAP ADDITION TEST:";
+ LOG(info) << "M1=\n" << M1.ToString(0, 2);
+ LOG(info) << "M2=\n" << M2.ToString(0, 2);
+ LOG(info) << "M12gen=\n" << M2.ToString(0, 2);
+ LOG(info) << "M12add=\n" << M2.ToString(0, 2);
+ }
+ }
}
+ LOG(info) << "kf::TrackingSetupInitializer: Tracking setup was initialized successfully";
return true;
}
catch (const std::exception& err) {
diff --git a/reco/kf/CbmKfTrackingSetupInitializer.h b/reco/kf/CbmKfTrackingSetupInitializer.h
index 4148c1ff636c7410c302217dbf94e6a21fb17185..50e643bb1589b47647162b80d4acb0443edfb221 100644
--- a/reco/kf/CbmKfTrackingSetupInitializer.h
+++ b/reco/kf/CbmKfTrackingSetupInitializer.h
@@ -9,7 +9,7 @@
#pragma once
-#include "KfSetup.h"
+#include "KfSetupInitializer.h"
#include <tuple>
@@ -37,12 +37,12 @@ namespace cbm::kf
private:
// Material map creator properties (TODO: Provide setters, if needed)
- static constexpr double kMatCreatorPitch{0.1}; ///< Material budget map minimal bin size [cm]
- static constexpr int kMatCreatorMaxNbins{100}; ///< Max number of bins in the material budget map in x(y) axis
- static constexpr int kMatCreatorNrays{3}; ///< Number of rays per dimension for the material budget
- static constexpr bool kMatCreatorSafeMode{true}; ///< Safe mode of the material map creation
- static constexpr double kTargFieldInitStep{2.5}; ///< Step between nodes in the target field initialization [cm]
- static constexpr double kTargMaterialOffset{0.2}; ///< Offset between target upper limit and its material zMax [cm]
+ static constexpr double kMatCreatorPitch{0.1}; ///< Material budget map minimal bin size [cm]
+ static constexpr int kMatCreatorMaxNbins{100}; ///< Max number of bins in the material budget map in x(y) axis
+ static constexpr int kMatCreatorNrays{3}; ///< Number of rays per dimension for the material budget
+ static constexpr bool kMatCreatorSafeMode{true}; ///< Safe mode of the material map creation
+ static constexpr double kTargFieldInitStep{2.5}; ///< Step between nodes in the target field initialization [cm]
+ static constexpr double kTargMaterialOffset{1}; ///< Offset between target upper limit and its material zMax [cm]
bool fbUseMvd{false}; ///< Are MVD stations included in the tracking setup
bool fbUseSts{false}; ///< Are STS stations included in the tracking setup