From 4dfdfe093bdde8c13109561320c8d8fea52f3c78 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Sat, 30 Mar 2024 23:58:45 +0100
Subject: [PATCH] kf: KfCore library
---
algo/CMakeLists.txt | 1 +
algo/ca/TrackingChain.cxx | 2 +-
algo/ca/core/CMakeLists.txt | 1 +
algo/ca/core/pars/CaConstants.h | 8 ++
algo/ca/core/pars/CaInitManager.cxx | 8 ++
algo/ca/core/pars/CaInitManager.h | 16 +++
algo/ca/core/tracking/CaFramework.cxx | 14 +++
algo/ca/core/tracking/CaFramework.h | 14 +++
algo/kf/CMakeLists.txt | 1 +
algo/kf/core/CMakeLists.txt | 59 +++++++++++
algo/kf/core/KfDefs.h | 77 ++++++++++++++
algo/kf/core/KfFramework.cxx | 48 +++++++++
algo/kf/core/KfFramework.h | 65 ++++++++++++
algo/kf/core/pars/KfMaterialMap.cxx | 143 ++++++++++++++++++++++++++
algo/kf/core/pars/KfMaterialMap.h | 137 ++++++++++++++++++++++++
algo/kf/core/pars/KfParameters.cxx | 59 +++++++++++
algo/kf/core/pars/KfParameters.h | 83 +++++++++++++++
algo/kf/core/pars/KfSetup.cxx | 61 +++++++++++
algo/kf/core/pars/KfSetup.h | 92 +++++++++++++++++
algo/kf/core/utils/KfVector.h | 95 +++++++++++++++++
reco/L1/CbmL1.cxx | 4 +-
21 files changed, 985 insertions(+), 3 deletions(-)
create mode 100644 algo/kf/CMakeLists.txt
create mode 100644 algo/kf/core/CMakeLists.txt
create mode 100644 algo/kf/core/KfDefs.h
create mode 100644 algo/kf/core/KfFramework.cxx
create mode 100644 algo/kf/core/KfFramework.h
create mode 100644 algo/kf/core/pars/KfMaterialMap.cxx
create mode 100644 algo/kf/core/pars/KfMaterialMap.h
create mode 100644 algo/kf/core/pars/KfParameters.cxx
create mode 100644 algo/kf/core/pars/KfParameters.h
create mode 100644 algo/kf/core/pars/KfSetup.cxx
create mode 100644 algo/kf/core/pars/KfSetup.h
create mode 100644 algo/kf/core/utils/KfVector.h
diff --git a/algo/CMakeLists.txt b/algo/CMakeLists.txt
index b04d23c36f..85b74016e4 100644
--- a/algo/CMakeLists.txt
+++ b/algo/CMakeLists.txt
@@ -53,6 +53,7 @@ endif()
add_subdirectory(log)
add_subdirectory(data)
+add_subdirectory(kf)
add_subdirectory(ca)
# exclude unittests from being build inside the container
diff --git a/algo/ca/TrackingChain.cxx b/algo/ca/TrackingChain.cxx
index 1599bf88fe..93347a0b08 100644
--- a/algo/ca/TrackingChain.cxx
+++ b/algo/ca/TrackingChain.cxx
@@ -63,8 +63,8 @@ void TrackingChain::Init()
// ------ Initialize CA framework
fCaMonitor.Reset();
fCaFramework.SetNofThreads(fConfig.fNofThreads);
- fCaFramework.Init(ca::Framework::TrackingMode::kMcbm);
fCaFramework.ReceiveParameters(std::move(parameters));
+ fCaFramework.Init(ca::Framework::TrackingMode::kMcbm);
// ------ Initialize QA modules
if (fInputQa.IsSenderDefined()) {
diff --git a/algo/ca/core/CMakeLists.txt b/algo/ca/core/CMakeLists.txt
index 90a19d9f52..b255a4bfdb 100644
--- a/algo/ca/core/CMakeLists.txt
+++ b/algo/ca/core/CMakeLists.txt
@@ -62,6 +62,7 @@ target_include_directories(CaCore
target_compile_definitions(CaCore PUBLIC NO_ROOT)
target_link_libraries(CaCore
+ KfCore
Vc::Vc
Boost::serialization
OnlineDataLog # needed for the logger
diff --git a/algo/ca/core/pars/CaConstants.h b/algo/ca/core/pars/CaConstants.h
index 4fd8d9bf84..c858914bc1 100644
--- a/algo/ca/core/pars/CaConstants.h
+++ b/algo/ca/core/pars/CaConstants.h
@@ -10,9 +10,17 @@
#pragma once // include this header only once per compilation unit
#include "CaSimd.h"
+#include "KfFramework.h"
#include <limits>
+namespace cbm::algo::ca
+{
+ using KfFramework_t = cbm::algo::kf::Framework<float>;
+ using KfParameters_t = cbm::algo::kf::Parameters<float>;
+ using KfSetup_t = cbm::algo::kf::Setup<float>;
+} // namespace cbm::algo::ca
+
/// Namespace contains compile-time constants definition for the CA tracking algorithm
///
namespace cbm::algo::ca::constants
diff --git a/algo/ca/core/pars/CaInitManager.cxx b/algo/ca/core/pars/CaInitManager.cxx
index b2318db200..3ab0d39044 100644
--- a/algo/ca/core/pars/CaInitManager.cxx
+++ b/algo/ca/core/pars/CaInitManager.cxx
@@ -504,6 +504,14 @@ void InitManager::SetTargetPosition(double x, double y, double z)
//
Parameters<fvec>&& InitManager::TakeParameters() { return std::move(fParameters); }
+// ----------------------------------------------------------------------------------------------------------------------
+//
+KfParameters_t&& InitManager::TakeKfParameters() { return std::move(fKfParameters); }
+
+// ----------------------------------------------------------------------------------------------------------------------
+//
+KfSetup_t&& InitManager::TakeKfSetup() { return std::move(fKfSetup); }
+
// ----------------------------------------------------------------------------------------------------------------------
//
void InitManager::WriteParametersObject(const std::string& fileName) const
diff --git a/algo/ca/core/pars/CaInitManager.h b/algo/ca/core/pars/CaInitManager.h
index 1667d07b61..7cb30eea68 100644
--- a/algo/ca/core/pars/CaInitManager.h
+++ b/algo/ca/core/pars/CaInitManager.h
@@ -17,6 +17,8 @@
#include "CaSimd.h"
#include "CaStationInitializer.h"
#include "CaUtils.h"
+#include "KfParameters.h"
+#include "KfSetup.h"
#include <array>
#include <bitset>
@@ -247,10 +249,20 @@ namespace cbm::algo::ca
fParameters.fMisalignmentT[static_cast<int>(detectorId)] = t;
}
+ /// \brief Sets default fitter mass
+ /// \param mass Particle mass [GeV/c2]
+ void SetDefaultMass(double mass) { fKfParameters.SetDefaultMass(mass); }
+
/// \brief Takes parameters object from the init-manager instance
/// \return A parameter object
Parameters<fvec>&& TakeParameters();
+ /// \brief Takes KF parameters object
+ KfParameters_t&& TakeKfParameters();
+
+ /// \brief Takes KF setup object
+ KfSetup_t&& TakeKfSetup();
+
/// \brief Writes parameters object from boost-serialized binary file
/// \param fileName Name of input file
void WriteParametersObject(const std::string& fileName) const;
@@ -310,6 +322,10 @@ namespace cbm::algo::ca
int fCAIterationsNumberCrosscheck{-1}; ///< Number of iterations to be passed (must be used for cross-checks)
Parameters<fvec> fParameters{}; ///< CA parameters object
+ KfParameters_t fKfParameters{}; ///< Kf-parameters object
+ KfSetup_t fKfSetup{}; ///< Kf setup object
+ // TODO: With a separate KF-framework instance we need to figure it out, how to store and read the corresponding
+ // parameters (essential for the online reconstruction!!!)
std::string fsConfigInputMain = ""; ///< name for the input configuration file
std::string fsConfigInputUser = ""; ///< name for the input configuration file
diff --git a/algo/ca/core/tracking/CaFramework.cxx b/algo/ca/core/tracking/CaFramework.cxx
index 82326d6755..43a7085dbe 100644
--- a/algo/ca/core/tracking/CaFramework.cxx
+++ b/algo/ca/core/tracking/CaFramework.cxx
@@ -22,6 +22,8 @@ Framework::Framework()
using cbm::algo::ca::ECounter; // monitor counter key type
using cbm::algo::ca::EDetectorID;
using cbm::algo::ca::ETimer; // monitor timer key type
+using cbm::algo::ca::KfFramework_t;
+using cbm::algo::ca::KfParameters_t;
//using cbm::ca::tools::Debugger;
// ---------------------------------------------------------------------------------------------------------------------
@@ -37,6 +39,7 @@ void Framework::Init(const TrackingMode mode)
fvTrackFinderWindow.emplace_back(*this, fvWData[iThread], fvMonitorDataThread[iThread]);
}
fTrackFinder.Init();
+ fKf.Init();
}
// ---------------------------------------------------------------------------------------------------------------------
@@ -76,6 +79,17 @@ void Framework::ReceiveParameters(Parameters<fvec>&& parameters)
ca::FieldRegion<fvec>::ForceUseOfOriginalField(fParameters.DevIsUseOfOriginalField());
}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+
+void Framework::ReceiveKfSettings(KfParameters_t&& pars, KfSetup_t&& setup)
+{
+ fKf.Pars() = std::move(pars);
+ fKf.Setup() = std::move(setup);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
int Framework::GetMcTrackIdForCaHit(int /*iHit*/) const
{
return -1;
diff --git a/algo/ca/core/tracking/CaFramework.h b/algo/ca/core/tracking/CaFramework.h
index 7239170b0a..409b165cb9 100644
--- a/algo/ca/core/tracking/CaFramework.h
+++ b/algo/ca/core/tracking/CaFramework.h
@@ -27,6 +27,7 @@
#include "CaTriplet.h"
#include "CaVector.h"
#include "CaWindowData.h"
+#include "KfFramework.h"
#include <array>
#include <iomanip>
@@ -140,6 +141,11 @@ namespace cbm::algo::ca
/// Receives tracking parameters
void ReceiveParameters(Parameters<fvec>&& parameters);
+ /// \brief Receives KF-framework parameters and setup
+ /// \param pars KF parameters object
+ /// \param setup KF setup object
+ void ReceiveKfSettings(KfParameters_t&& pars, KfSetup_t&& setup);
+
/// Gets pointer to input data object for external access
const InputData& GetInputData() const { return fInputData; }
@@ -207,7 +213,15 @@ namespace cbm::algo::ca
/// \brief Gets number of threads
int GetNofThreads() const { return fNofThreads; }
+ /// \brief Accesses KF framework instance
+ KfFramework_t& Kf() { return fKf; }
+
+ /// \brief Accesses KF framework instance (immutable)
+ const KfFramework_t& Kf() const { return fKf; }
+
private:
+ KfFramework_t fKf; ///< KF framework instance
+
int fNstationsBeforePipe{0}; ///< number of stations before pipe (MVD stations in CBM)
int fNfieldStations{0}; ///< number of stations in the field region
float fDefaultMass{constants::phys::MuonMass}; ///< mass of the propagated particle [GeV/c2]
diff --git a/algo/kf/CMakeLists.txt b/algo/kf/CMakeLists.txt
new file mode 100644
index 0000000000..ad6d4787c2
--- /dev/null
+++ b/algo/kf/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(core)
diff --git a/algo/kf/core/CMakeLists.txt b/algo/kf/core/CMakeLists.txt
new file mode 100644
index 0000000000..f89e9015c6
--- /dev/null
+++ b/algo/kf/core/CMakeLists.txt
@@ -0,0 +1,59 @@
+set(INCLUDE_DIRECTORIES
+ ${CMAKE_CURRENT_SOURCE_DIR}
+ ${CMAKE_CURRENT_SOURCE_DIR}/data
+ ${CMAKE_CURRENT_SOURCE_DIR}/pars
+ ${CMAKE_CURRENT_SOURCE_DIR}/utils
+)
+
+set(SRCS
+ ${CMAKE_CURRENT_SOURCE_DIR}/KfFramework.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/pars/KfMaterialMap.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/pars/KfParameters.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/pars/KfSetup.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(KfCore SHARED ${SRCS})
+
+target_include_directories(KfCore
+ PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/data
+ ${CMAKE_CURRENT_SOURCE_DIR}/pars
+ ${CMAKE_CURRENT_SOURCE_DIR}/utils
+ ${CMAKE_CURRENT_SOURCE_DIR}
+)
+
+target_compile_definitions(KfCore PUBLIC NO_ROOT)
+
+target_link_libraries(KfCore
+ Vc::Vc
+ Boost::serialization
+ OnlineDataLog # needed for the logger
+ external::fles_logging # needed for the logger
+ external::fles_ipc # needed for the logger
+ external::yaml-cpp
+ )
+
+install(TARGETS KfCore DESTINATION lib)
+install(DIRECTORY kf TYPE INCLUDE FILES_MATCHING PATTERN "*.h")
+install(DIRECTORY kf/utils TYPE INCLUDE FILES_MATCHING PATTERN "*.h")
+install(DIRECTORY kf/data TYPE INCLUDE FILES_MATCHING PATTERN "*.h")
+install(DIRECTORY kf/pars TYPE INCLUDE FILES_MATCHING PATTERN "*.h")
+
+install(
+ FILES
+ KfFramework.h
+ KfDefs.h
+ pars/KfMaterialMap.h
+ pars/KfParameters.h
+ pars/KfSetup.h
+ utils/KfVector.h
+ DESTINATION
+ include/
+)
diff --git a/algo/kf/core/KfDefs.h b/algo/kf/core/KfDefs.h
new file mode 100644
index 0000000000..3fafcc560e
--- /dev/null
+++ b/algo/kf/core/KfDefs.h
@@ -0,0 +1,77 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file KfDefs.h
+/// @brief Common constant definitions for the Kalman Filter library
+/// @since 28.03.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include <limits>
+
+namespace cbm::algo::kf::defs
+{
+ // ----- Control -----------------------------------------------------------------------------------------------------
+ constexpr int DebugLvl = 0; ///< Level of debug output
+ constexpr bool GetterCheck = true; ///< Bound check in getters
+
+
+ // ----- Math --------------------------------------------------------------------------------------------------------
+ template<class T = double>
+ constexpr T Pi = T(3.14159265358979323846L); ///< Value of PI, used in ROOT TMath
+
+
+ // ----- Physics -----------------------------------------------------------------------------------------------------
+ template<class T = double>
+ constexpr T MuonMass = T(0.105658375523L); ///< Muon mass [GeV/c2]
+
+ template<class T = double>
+ constexpr T PionMass = T(0.1395703918L); ///< Pion mass [GeV/c2]
+
+ template<class T = double>
+ constexpr T KaonMass = T(0.493677L); ///< Kaon mass [GeV/c2] (PDG 22.08.2023)
+
+ template<class T = double>
+ constexpr T ElectronMass = T(0.0005109989500015L); ///< Electron mass [GeV/c2]
+
+ template<class T = double>
+ constexpr T ProtonMass = T(0.93827208816L); ///< Proton mass [GeV/c2] (PDG 11.08.2022)
+
+ template<class T = double>
+ constexpr T SpeedOfLight = T(29.9792458L); ///< Speed of light [cm/ns]
+
+ template<class T = double>
+ constexpr T SpeedOfLightInv = T(1.L) / SpeedOfLight<T>; ///< Inverse speed of light [ns/cm]
+
+ template<class T = double>
+ constexpr T MinFeild = T(1.e-4L); ///< Minimal (negligible) magnetic field value [kG]
+
+
+ // ----- Misc --------------------------------------------------------------------------------------------------------
+ constexpr int MemAlign = 16; ///< Default alignment of memory (bytes)
+
+
+ // ----- Undefined values --------------------------------------------------------------------------------------------
+ /// \brief Undefined values
+ template<typename T1, typename T2 = T1>
+ constexpr T2 Undef;
+
+ template<>
+ inline constexpr int Undef<int> = std::numeric_limits<int>::min();
+
+ template<>
+ inline constexpr unsigned Undef<unsigned> = std::numeric_limits<unsigned>::max();
+
+ template<>
+ inline constexpr float Undef<float> = std::numeric_limits<float>::signaling_NaN();
+
+ template<>
+ inline constexpr double Undef<double> = std::numeric_limits<double>::signaling_NaN();
+
+ // TODO: provide SIMD
+ //template<>
+ //inline constexpr fscal Undef<fvec> = std::numeric_limits<fscal>::signaling_NaN();
+
+} // namespace cbm::algo::kf::defs
diff --git a/algo/kf/core/KfFramework.cxx b/algo/kf/core/KfFramework.cxx
new file mode 100644
index 0000000000..a11fe973de
--- /dev/null
+++ b/algo/kf/core/KfFramework.cxx
@@ -0,0 +1,48 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file KfFramework.cxx
+/// @brief The Kalman-filter framework main class (header)
+/// @since 28.03.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#include "KfFramework.h"
+
+using cbm::algo::kf::Framework;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename DataT>
+bool Framework<DataT>::Init()
+try {
+ std::stringstream errMsg;
+ try {
+ fPars.Init();
+ }
+ catch (const std::runtime_error& err) {
+ errMsg << "\tParameters initialization errors: " << err.what() << '\n';
+ }
+
+ try {
+ fSetup.Init();
+ }
+ catch (const std::runtime_error& err) {
+ errMsg << "\tSetup initialization errors: " << err.what() << '\n';
+ }
+
+ if (!errMsg.str().empty()) {
+ throw std::runtime_error(errMsg.str());
+ }
+
+ L_(info) << '\n' << fPars.ToString(3) << '\n' << fSetup.ToString(3);
+ return true;
+}
+catch (const std::exception& err) {
+ L_(error) << "kf::Framework: initialization failed. Reason: " << err.what();
+ return false;
+}
+
+template class cbm::algo::kf::Framework<float>;
+template class cbm::algo::kf::Framework<double>;
+// template class cbm::algo::kf::Framework<fvec>;
diff --git a/algo/kf/core/KfFramework.h b/algo/kf/core/KfFramework.h
new file mode 100644
index 0000000000..902126523d
--- /dev/null
+++ b/algo/kf/core/KfFramework.h
@@ -0,0 +1,65 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file KfFramework.h
+/// @brief The Kalman-filter framework main class (header)
+/// @since 28.03.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include "KfDefs.h"
+#include "KfParameters.h"
+#include "KfSetup.h"
+
+namespace cbm::algo::kf
+{
+ /// \class Framework
+ /// \brief Main class of the KfCore library
+ /// \tparam DataT Underlying data-type
+ template<typename DataT>
+ class Framework {
+ public:
+ /// \brief Default constructor
+ Framework() = default;
+
+ /// \brief Copy constructor
+ Framework(const Framework&) = delete;
+
+ /// \brief Move constructor
+ Framework(Framework&&) = delete;
+
+ /// \brief Destructor
+ ~Framework() = default;
+
+ /// \brief Copy assignment operator
+ Framework& operator=(const Framework&) = delete;
+
+ /// \brief Move assignment operator
+ Framework& operator=(Framework&&) = delete;
+
+ /// \brief Initialization method
+ bool Init();
+
+ /// \brief Initialization satus
+ bool IsInitialized() const { return fbInitialized; }
+
+ /// \brief Parameters access
+ Parameters<DataT>& Pars() { return fPars; }
+
+ /// \brief Parameters access
+ const kf::Parameters<DataT>& Pars() const { return fPars; }
+
+ /// \brief Setup access (mutable)
+ kf::Setup<DataT>& Setup() { return fSetup; }
+
+ /// \brief Setup access
+ const kf::Setup<DataT>& Setup() const { return fSetup; }
+
+ private:
+ kf::Parameters<DataT> fPars; ///< KF parameters
+ kf::Setup<DataT> fSetup; ///< KF setup
+ bool fbInitialized = false; ///< Initialization status
+ };
+} // namespace cbm::algo::kf
diff --git a/algo/kf/core/pars/KfMaterialMap.cxx b/algo/kf/core/pars/KfMaterialMap.cxx
new file mode 100644
index 0000000000..809e22753e
--- /dev/null
+++ b/algo/kf/core/pars/KfMaterialMap.cxx
@@ -0,0 +1,143 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+#include "KfMaterialMap.h"
+
+#include "AlgoFairloggerCompat.h"
+
+#include <iomanip>
+#include <sstream>
+#include <vector>
+
+using cbm::algo::kf::MaterialMap;
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+MaterialMap::MaterialMap(MaterialMap&& other) noexcept { this->Swap(other); }
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+MaterialMap& MaterialMap::operator=(MaterialMap&& other) noexcept
+{
+ if (this != &other) {
+ MaterialMap tmp(std::move(other));
+ this->Swap(tmp);
+ }
+ return *this;
+}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+int MaterialMap::GetBin(float x, float y) const
+{
+ 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;
+}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+float MaterialMap::GetRadThickScal(float x, float y) const
+{
+ x = (x < fXYmax && x >= -fXYmax) ? x : 0;
+ y = (y < fXYmax && y >= -fXYmax) ? y : 0;
+ int i = static_cast<int>((x + fXYmax) * fFactor);
+ int j = static_cast<int>((y + fXYmax) * fFactor);
+ i = (i < fNbins && i >= 0) ? i : fNbins / 2;
+ j = (j < fNbins && j >= 0) ? j : fNbins / 2;
+ return fTable[i + j * fNbins];
+}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+// TODO: provide SIMD
+//fvec MaterialMap::GetRadThickVec(fvec x, fvec y) const
+//{
+// fvec r;
+// for (size_t i = 0; i < fvec::size(); i++)
+// r[i] = GetRadThickScal(x[i], y[i]);
+// return r;
+//}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+void MaterialMap::CheckConsistency() const
+{
+ /* (i) Check, if the object was initialized */
+ // TODO: provide undef checking utility
+ //if (IsNaN()) {
+ // throw std::logic_error("MaterialMap: class was not initialized");
+ //}
+
+ /* (ii) Check if the thickness values correct (non-negative) */
+ bool isThicknessOk = true;
+ for (int i = 0; i < int(fTable.size()); ++i) {
+ if (fTable[i] < 0) {
+ isThicknessOk = false;
+ LOG(error) << "MaterialMap: found illegal negative station thickness value " << fTable[i]
+ << " at iBinX = " << (i % fNbins) << ", iBin = " << (i / fNbins);
+ }
+ }
+ if (!isThicknessOk) {
+ throw std::logic_error("MaterialMap: incorrect station thickness values found in the thickness map");
+ }
+}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+void MaterialMap::Initialize(int nBins, float xyMax, float zRef, float zMin, float zMax)
+{
+ fNbins = nBins;
+ fXYmax = xyMax;
+ fZref = zRef;
+ fZmin = zMin;
+ fZmax = zMax;
+
+ if (fNbins < 1) {
+ std::stringstream aStream;
+ aStream << "MaterialMap: object cannot be initialized with non-positive nBins = " << fNbins;
+ throw std::logic_error(aStream.str());
+ }
+
+ if (fXYmax < 0.) {
+ std::stringstream aStream;
+ aStream << "MaterialMap: object cannot be initialized with non-positive XYmax = " << fXYmax << " [cm]";
+ throw std::logic_error(aStream.str());
+ }
+
+ if (!((fZmin <= fZref) && (fZref <= zMax))) {
+ std::stringstream aStream;
+ aStream << "MaterialMap: object cannot be initialized with inconsistent Z: min " << fZmin << " ref " << fZref
+ << " max " << fZmax << " [cm]";
+ throw std::logic_error(aStream.str());
+ }
+
+ fFactor = 0.5 * fNbins / fXYmax;
+ fTable.resize(fNbins * fNbins);
+}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+void MaterialMap::Swap(MaterialMap& other) noexcept
+{
+ std::swap(fNbins, other.fNbins);
+ std::swap(fXYmax, other.fXYmax);
+ std::swap(fFactor, other.fFactor);
+ std::swap(fZref, other.fZref);
+ std::swap(fZmin, other.fZmin);
+ std::swap(fZmax, other.fZmax);
+ std::swap(fTable, other.fTable); // Probably can cause segmentation violation (did not understand)
+}
+
+//------------------------------------------------------------------------------------------------------------------------------------
+//
+std::string MaterialMap::ToString() const
+{
+ std::stringstream msg;
+ 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;
+ return msg.str();
+}
diff --git a/algo/kf/core/pars/KfMaterialMap.h b/algo/kf/core/pars/KfMaterialMap.h
new file mode 100644
index 0000000000..b97fde1f90
--- /dev/null
+++ b/algo/kf/core/pars/KfMaterialMap.h
@@ -0,0 +1,137 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Igor Kulakov, Sergey Gorbunov, Andrey Lebedev, Sergei Zharko [committer] */
+
+#pragma once // include this header only once per compilation unit
+
+// NOTE: No dependency from CaCore is allowed
+//#include "CaSimd.h"
+//#include "CaUtils.h"
+#include "AlgoFairloggerCompat.h"
+#include "KfDefs.h"
+
+#include <boost/serialization/vector.hpp>
+
+#include <iomanip>
+#include <sstream>
+#include <string>
+#include <vector>
+
+namespace cbm::algo::kf
+{
+ /// \class MaterialMap
+ /// \brief A map of station thickness in units of radiation length (X0) to the specific point in XY plane
+ class MaterialMap {
+ public:
+ /// \brief Default constructor
+ MaterialMap() = default;
+
+ /// \brief Copy constructor
+ MaterialMap(const MaterialMap& other) = default;
+
+ /// \brief Copy assignment operator
+ MaterialMap& operator=(const MaterialMap& other) = default;
+
+ /// \brief Move constructor
+ MaterialMap(MaterialMap&& other) noexcept;
+
+ /// \brief Move assignment operator
+ MaterialMap& operator=(MaterialMap&& other) noexcept;
+
+ /// \brief Destructor
+ ~MaterialMap() noexcept = default;
+
+ /// \brief Gets number of bins (rows or columns) of the material table
+ int GetNbins() const { return fNbins; }
+
+ /// \brief Gets radius in cm of the material table
+ float GetXYmax() const { return fXYmax; }
+
+ /// \brief Gets reference Z of the material in cm
+ float GetZref() const { return fZref; }
+
+ /// \brief Gets minimal Z of the collected material in cm
+ float GetZmin() const { return fZmin; }
+
+ /// \brief Gets maximal Z of the collected material in 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
+ /// \param x X coordinate of the point [cm]
+ /// \param y Y coordinate of the point [cm]
+ float GetRadThickScal(float x, float y) const;
+
+ /// \brief Gets material thickness in units of X0 in (x,y) point of the station
+ /// cbm::algo::ca::fvec type can be float, that is why "Vec" and "Scal" specifications
+ /// \param x X coordinate of the point [cm] (SIMDized vector)
+ /// \param y Y coordinate of the point [cm] (SIMDized veotor)
+ // TODO: provide SIMD
+ //fvec GetRadThickVec(fvec x, fvec y) const;
+
+ /// \brief Checks, if the fields are NaN
+ //bool IsNaN() const
+ //{
+ // return utils::IsUndefined(fNbins) || utils::IsUndefined(fXYmax) || utils::IsUndefined(fFactor)
+ // || utils::IsUndefined(fZref) || utils::IsUndefined(fZmin) || utils::IsUndefined(fZmax);
+ //}
+
+ /// \brief Verifies class invariant consistency
+ void CheckConsistency() const;
+
+ /// \brief Sets value of material thickness in units of X0 for a given cell of the material table
+ /// \param iBinX Index of table column
+ /// \param iBinY Index of table row
+ /// \param thickness Thickness of the material in units of X0
+ /// \note Indices of rows and columns in the table runs from 0 to nBins-1 inclusively, where nBins is the number
+ /// both of rows and columns. One should be careful while reading and storing the table from ROOT-file,
+ /// because iBinX = 0 and iBinY = 0 in the TH1::SetBinContent method of usually defines the underflow bin.
+ void SetRadThickBin(int iBinX, int iBinY, float thickness) { fTable[iBinX + fNbins * iBinY] = thickness; }
+
+ /// \brief Sets properties of the material table -- number of rows or columnts and the size of station in XY plane
+ /// \param nBins Number of rows or columns
+ /// \param xyMax Size of station in x and y dimensions [cm]
+ void Initialize(int nBins, float xyMax, float zRef, float zMin, float zMax);
+
+ /// \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;
+
+ private:
+ 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
+ float fZref = defs::Undef<float>; ///< Reference Z of the collected material [cm]
+ float fZmin = defs::Undef<float>; ///< Minimal Z of the collected material [cm]
+ float fZmax = defs::Undef<float>; ///< Minimal Z of the collected material [cm]
+ std::vector<float> fTable{}; ///< Material budget table
+
+ /// \brief Serialization function
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int)
+ {
+ ar& fNbins;
+ ar& fXYmax;
+ ar& fFactor;
+ ar& fZref;
+ ar& fZmin;
+ ar& fZmax;
+ ar& fTable;
+ }
+ };
+
+} // namespace cbm::algo::kf
diff --git a/algo/kf/core/pars/KfParameters.cxx b/algo/kf/core/pars/KfParameters.cxx
new file mode 100644
index 0000000000..999aec0b9a
--- /dev/null
+++ b/algo/kf/core/pars/KfParameters.cxx
@@ -0,0 +1,59 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file KfParameters.cxx
+/// @brief Parameters representation for the Kalman-filter framework (source)
+/// @since 28.03.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#include "KfParameters.h"
+
+#include <sstream>
+
+using cbm::algo::kf::Parameters;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename DataT>
+void Parameters<DataT>::CheckConsistency() const
+{
+ std::stringstream errMsg;
+
+ // checks go here... the errMsg is to be filled on each checking stage
+
+ if (!errMsg.str().empty()) {
+ throw std::runtime_error(errMsg.str());
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename DataT>
+void Parameters<DataT>::Init()
+{
+ fbInitialized = true;
+ this->CheckConsistency();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename DataT>
+std::string Parameters<DataT>::ToString(int verbosity, int indentLevel) const
+{
+ std::stringstream msg;
+
+ if (verbosity > 0) {
+ constexpr char indentCh = '\t';
+ std::string indent(indentLevel, indentCh);
+ msg << indent << " ----- KF Parameters list ----- \n";
+ msg << indent << "RUNTIME CONSTANTS:\n";
+ msg << indent << indentCh << "Default particle mass: " << fDefaultMass << " [GeV/c2] \n";
+ }
+
+ return msg.str();
+}
+
+template class cbm::algo::kf::Parameters<float>;
+template class cbm::algo::kf::Parameters<double>;
+// template class cbm::algo::kf::Parameters<fvec>;
diff --git a/algo/kf/core/pars/KfParameters.h b/algo/kf/core/pars/KfParameters.h
new file mode 100644
index 0000000000..922f4454a4
--- /dev/null
+++ b/algo/kf/core/pars/KfParameters.h
@@ -0,0 +1,83 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file KfParameters.cxx
+/// @brief Parameter set for the Kalman-filter framework (header)
+/// @since 28.03.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include "KfDefs.h"
+
+#include <boost/serialization/access.hpp>
+
+#include <string>
+
+namespace cbm::algo::kf
+{
+ /// \class Parameters
+ /// \brief Parameter set of the Kalman-filter framework
+ /// \tparam DataT Underlying data type
+ template<typename DataT>
+ class Parameters {
+ public:
+ /// \brief Default constructor
+ Parameters() = default;
+
+ /// \brief Copy constructor
+ Parameters(const Parameters&) = default;
+
+ /// \brief Move constructor
+ Parameters(Parameters&&) = default;
+
+ /// \brief Destructor
+ ~Parameters() = default;
+
+ /// \brief Copy assignment operator
+ Parameters& operator=(const Parameters&) = default;
+
+ /// \brief Move assignment operator
+ Parameters& operator=(Parameters&&) = default;
+
+ /// \brief Checks parameters instance
+ /// \note Throws std::runtime_error, if the class inconsistent
+ void CheckConsistency() const;
+
+ /// \brief Gets default mass
+ const DataT GetDefaultMass() const { return fDefaultMass; }
+
+ /// \brief Initialization method
+ void Init();
+
+ /// \brief Returns initialization status
+ bool IsInitialized() const { return fbInitialized; }
+
+ /// \brief Resets the contents of the parameters
+ void Reset() { *this = std::move(Parameters<DataT>{}); }
+
+ /// \brief Sets default mass
+ void SetDefaultMass(DataT mass) { fDefaultMass = mass; }
+
+ /// \brief String representation of the 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;
+
+ private:
+ /// \brief Serialization method
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int /*version*/)
+ {
+ ar& fDefaultMass;
+ ar& fbInitialized;
+ }
+
+ DataT fDefaultMass = defs::PionMass<DataT>; ///< Default particle mass [GeV/c2]
+
+ bool fbInitialized = false; ///< Is instance initialized
+ };
+
+} // namespace cbm::algo::kf
diff --git a/algo/kf/core/pars/KfSetup.cxx b/algo/kf/core/pars/KfSetup.cxx
new file mode 100644
index 0000000000..bce84eece3
--- /dev/null
+++ b/algo/kf/core/pars/KfSetup.cxx
@@ -0,0 +1,61 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file KfSetup.cxx
+/// @brief Setup representation for the Kalman-filter framework (source)
+/// @since 28.03.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#include "KfSetup.h"
+
+#include <sstream>
+
+using cbm::algo::kf::Setup;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename DataT>
+void Setup<DataT>::CheckConsistency() const
+{
+ std::stringstream errMsg;
+
+ // checks go here... the errMsg is to be filled on each checking stage
+
+ if (!errMsg.str().empty()) {
+ throw std::runtime_error(errMsg.str());
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename DataT>
+void Setup<DataT>::Init()
+{
+ fbInitialized = true;
+ CheckConsistency();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+template<typename DataT>
+std::string Setup<DataT>::ToString(int verbosity, int indentLevel) const
+{
+ std::stringstream msg;
+
+ if (verbosity > 0) {
+ constexpr char indentCh = '\t';
+ std::string indent(indentLevel, indentCh);
+ msg << indent << " ----- KF setup -----\n";
+ msg << indent << "MATERIAL LAYERS:\n";
+ for (const auto& layer : fvMatLayers) {
+ msg << indent << indentCh << layer.ToString() << '\n';
+ }
+ }
+
+ return msg.str();
+}
+
+template class cbm::algo::kf::Setup<float>;
+template class cbm::algo::kf::Setup<double>;
+//template class cbm::algo::kf::Setup<fvec>;
diff --git a/algo/kf/core/pars/KfSetup.h b/algo/kf/core/pars/KfSetup.h
new file mode 100644
index 0000000000..dcd211771e
--- /dev/null
+++ b/algo/kf/core/pars/KfSetup.h
@@ -0,0 +1,92 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file KfSetup.h
+/// @brief Setup representation for the Kalman-filter framework (header)
+/// @since 28.03.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once // include this header only once per compilation unit
+
+#include "KfDefs.h"
+#include "KfMaterialMap.h"
+#include "KfVector.h"
+
+#include <boost/serialization/access.hpp>
+
+#include <string>
+#include <vector>
+
+namespace cbm::algo::kf
+{
+ using MaterialContainer_t = Vector<MaterialMap>;
+
+ /// \brief KF-framework representation of the detector setup
+ /// \tparam DataT Underlying data type
+ template<typename DataT>
+ class alignas(defs::MemAlign) Setup {
+ public:
+ /// \brief Default constructor
+ Setup() = default;
+
+ /// \brief Destructor
+ ~Setup() = default;
+
+ /// \brief Copy constructor
+ Setup(const Setup&) = default;
+
+ /// \brief Move constructor
+ Setup(Setup&&) = default;
+
+ /// \brief Copy assignment operator
+ Setup& operator=(const Setup&) = default;
+
+ /// \brief Move assignment operator
+ Setup& operator=(Setup&&) = default;
+
+ /// \brief Adds material layer
+ void AddMaterialLayer(const MaterialMap& layer) { fvMatLayers.push_back(layer); }
+
+ /// \brief Adds material layer, moving it from the source
+ void AddMaterialLayer(MaterialMap&& layer) { fvMatLayers.push_back(std::move(layer)); }
+
+ /// \brief Check method
+ /// \note Throws std::runtime_error, if the class inconsistent
+ void CheckConsistency() const;
+
+ /// \brief Gets material layer
+ [[gnu::always_inline]] const MaterialMap& GetMaterialLayer(int iLayer) const
+ {
+ return fvMatLayers.at<defs::GetterCheck>(iLayer);
+ }
+
+ /// \brief Initializes setup
+ void Init();
+
+ /// \brief Returns initialization status
+ bool IsInitialized() const { return fbInitialized; }
+
+ /// \brief Resets initialization
+ void Reset() { *this = std::move(Setup{}); }
+
+ /// \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;
+
+ private:
+ /// \brief Serialization method
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int /*version*/)
+ {
+ ar& fvMatLayers;
+ ar& fbInitialized;
+ }
+
+ MaterialContainer_t fvMatLayers = {}; ///< Material layers map
+
+ bool fbInitialized = false; ///< Is instance initialized
+ };
+} // namespace cbm::algo::kf
diff --git a/algo/kf/core/utils/KfVector.h b/algo/kf/core/utils/KfVector.h
new file mode 100644
index 0000000000..9074b71a78
--- /dev/null
+++ b/algo/kf/core/utils/KfVector.h
@@ -0,0 +1,95 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file KfVector.h
+/// \brief std::vector with an additional utility set
+/// \since 30.03.2024
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include "AlgoFairloggerCompat.h"
+
+#include <boost/serialization/access.hpp>
+#include <boost/serialization/base_object.hpp>
+#include <boost/serialization/vector.hpp>
+//#include <boost/stacktrace.hpp>
+
+#include <exception>
+#include <sstream>
+#include <vector>
+
+namespace cbm::algo::kf
+{
+ /// \class Vector
+ /// \brief A std::vector with additional debugging utility set
+ template<class T, class Allocator = std::allocator<T>>
+ class Vector : public std::vector<T, Allocator> {
+ public:
+ using BaseVector_t = std::vector<T, Allocator>;
+
+ /// \brief Generic constructor from parameters
+ template<typename... Args>
+ Vector(Args... values) : BaseVector_t(values...)
+ {
+ }
+
+ /// \brief Copy constructor
+ Vector(const Vector& other) : BaseVector_t(other) {}
+
+ /// \brief Move constructor
+ Vector(Vector&& other) : BaseVector_t(std::move(other)) {}
+
+ /// \brief Copy assignment operator
+ //Vector& operator=(const Vector& other) { return BaseVector_t::operator=(other); }
+ Vector& operator=(const Vector& other) = default;
+
+ /// \brief Move assignment operator
+ //Vector& operator=(Vector&& other) { return BaseVector_t::operator=(std::move(other)); }
+ Vector& operator=(Vector&& other) = default;
+
+ /// \brief Access operator overload
+ template<bool CheckBoundaries = false>
+ T& at(std::size_t pos)
+ {
+ if constexpr (CheckBoundaries) {
+ if (pos >= this->size()) {
+ std::stringstream msg;
+ msg << "kf::Vector::operator[]: accessing element with index out of boundaries (pos = " << pos
+ << ", size = " << this->size() << "). "
+ << "Stack trace:\n";
+ //<< boost::stacktrace::stacktrace();
+ throw std::runtime_error(msg.str());
+ }
+ }
+ return BaseVector_t::operator[](pos);
+ }
+
+ /// \brief Constant access operator overload
+ template<bool CheckBoundaries = false>
+ const T& at(std::size_t pos) const
+ {
+ if constexpr (CheckBoundaries) {
+ if (pos >= this->size()) {
+ std::stringstream msg;
+ msg << "kf::Vector::operator[]: accessing element with index out of boundaries (pos = " << pos
+ << ", size = " << this->size() << "). "
+ << "Stack trace:\n";
+ //<< boost::stacktrace::stacktrace();
+ throw std::runtime_error(msg.str());
+ }
+ }
+ return BaseVector_t::operator[](pos);
+ }
+
+ private:
+ /// \brief Serialization method
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int /*version*/)
+ {
+ ar& boost::serialization::base_object<BaseVector_t>(*this);
+ }
+ };
+} // namespace cbm::algo::kf
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 0e38eeb956..493e306ea0 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -526,8 +526,6 @@ try {
// Re-organize the the relation between CbmL1 and ca::Framework. I believe, we don't need a global pointer to the ca::Framework
// instance.
fpAlgo = &gAlgo;
- fpAlgo->Init(fTrackingMode);
-
//
// ** Send formed parameters object to ca::Framework instance **
@@ -536,7 +534,9 @@ try {
auto parameters = fInitManager.TakeParameters();
fpParameters = std::make_shared<ca::Parameters<ca::fvec>>(parameters);
fpAlgo->ReceiveParameters(std::move(parameters));
+ fpAlgo->ReceiveKfSettings(fInitManager.TakeKfParameters(), fInitManager.TakeKfSetup());
}
+ fpAlgo->Init(fTrackingMode);
// Initialize time-slice reader
fpTSReader = std::make_unique<TimeSliceReader>();
--
GitLab