From a6af70f5a347cda9f51971bc04e2e230c77c1da7 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Wed, 18 Sep 2024 01:25:04 +0200
Subject: [PATCH] KF-setup initialization improvement
---
algo/kf/core/CMakeLists.txt | 6 +-
algo/kf/core/KfDefs.h | 4 +-
algo/kf/core/geo/KfField.h | 2 +-
algo/kf/core/geo/KfIMaterialMapFactory.h | 43 ++++
algo/kf/core/geo/KfModuleIndexMap.cxx | 119 +++++++++
algo/kf/core/geo/KfModuleIndexMap.h | 193 +++++++++++++++
algo/kf/core/geo/KfSetup.cxx | 2 +-
algo/kf/core/geo/KfSetup.h | 62 ++---
algo/kf/core/geo/KfSetupBuilder.cxx | 108 ++++++++
algo/kf/core/geo/KfSetupBuilder.h | 233 ++++++++++++++++++
algo/kf/core/geo/KfSetupInitializer.cxx | 63 -----
algo/kf/core/geo/KfSetupInitializer.h | 166 -------------
algo/kf/core/geo/KfTarget.cxx | 8 +-
algo/kf/core/geo/KfTarget.h | 29 ++-
reco/L1/CMakeLists.txt | 1 -
reco/L1/CbmKfTrackingSetupInitializer.cxx | 195 ---------------
reco/L1/CbmKfTrackingSetupInitializer.h | 53 ----
reco/L1/CbmL1.cxx | 11 +-
reco/L1/CbmL1DetectorID.h | 2 +-
reco/L1/CbmL1Util.cxx | 3 +-
reco/kfnew/CbmKfTrackingSetupInitializer.cxx | 172 +++----------
reco/kfnew/CbmKfTrackingSetupInitializer.h | 32 ++-
reco/kfnew/tools/CMakeLists.txt | 4 +-
...apCreator.cxx => KfMaterialMapFactory.cxx} | 51 ++--
...ialMapCreator.h => KfMaterialMapFactory.h} | 14 +-
25 files changed, 862 insertions(+), 714 deletions(-)
create mode 100644 algo/kf/core/geo/KfIMaterialMapFactory.h
create mode 100644 algo/kf/core/geo/KfModuleIndexMap.cxx
create mode 100644 algo/kf/core/geo/KfModuleIndexMap.h
create mode 100644 algo/kf/core/geo/KfSetupBuilder.cxx
create mode 100644 algo/kf/core/geo/KfSetupBuilder.h
delete mode 100644 algo/kf/core/geo/KfSetupInitializer.cxx
delete mode 100644 algo/kf/core/geo/KfSetupInitializer.h
delete mode 100644 reco/L1/CbmKfTrackingSetupInitializer.cxx
delete mode 100644 reco/L1/CbmKfTrackingSetupInitializer.h
rename reco/kfnew/tools/{KfMaterialMapCreator.cxx => KfMaterialMapFactory.cxx} (89%)
rename reco/kfnew/tools/{KfMaterialMapCreator.h => KfMaterialMapFactory.h} (91%)
diff --git a/algo/kf/core/CMakeLists.txt b/algo/kf/core/CMakeLists.txt
index 7ee3882d8f..a39e049f2f 100644
--- a/algo/kf/core/CMakeLists.txt
+++ b/algo/kf/core/CMakeLists.txt
@@ -21,7 +21,8 @@ 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}/geo/KfSetupBuilder.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/geo/KfModuleIndexMap.cxx
${CMAKE_CURRENT_SOURCE_DIR}/utils/KfUtils.cxx
)
@@ -110,7 +111,8 @@ install(
geo/KfMaterialMap.h
geo/KfModuleIndexMap.h
geo/KfSetup.h
- geo/KfSetupInitializer.h
+ geo/KfSetupBuilder.h
+ geo/KfIMaterialMapFactory.h
geo/KfTarget.h
geo/KfField.h
diff --git a/algo/kf/core/KfDefs.h b/algo/kf/core/KfDefs.h
index efd12d4c76..59268c388c 100644
--- a/algo/kf/core/KfDefs.h
+++ b/algo/kf/core/KfDefs.h
@@ -72,8 +72,8 @@ namespace cbm::algo::kf::defs
//constexpr int MaxNofFieldSlices = 64; ///< Max number of field slices
//constexpr int MaxNofMaterialLayers = 32; ///< Max number of the material layers
// NOTE: max uint8_t is assumed in order to satisfy fles::Subsystem
- constexpr int MaxNofDetSubsystems = 256; ///< Max number of detector types (STS, TRD, RICH,...)
- constexpr int MaxNofDetComponents = 128; ///< Max number of detector components (stations, layers, ...)
+ constexpr int MaxNofDetSubsystems = 256; ///< Max number of detector types (STS, TRD, RICH,...)
+ constexpr int MaxNofDetComponents = 128; ///< Max number of detector components (stations, layers, ...)
// ----- Control -----------------------------------------------------------------------------------------------------
constexpr int DebugLvl = 0; ///< Level of debug output
diff --git a/algo/kf/core/geo/KfField.h b/algo/kf/core/geo/KfField.h
index 680f2587c7..218c7ed1ea 100644
--- a/algo/kf/core/geo/KfField.h
+++ b/algo/kf/core/geo/KfField.h
@@ -456,7 +456,7 @@ namespace cbm::algo::kf
fld[iNode] = GlobalField::GetFieldValue<double>(fFieldFn, fTarget[0], fTarget[1], z);
z += fTargetStep;
}
- field.fPrimVertexField = FieldRegion<T>(fFieldType, fld[0], pos[0], fld[1], pos[1], fld[2], pos[2]);
+ field.fPrimVertexField = FieldRegion<T>(fld[0], pos[0], fld[1], pos[1], fld[2], pos[2]);
// TODO: Provide alternative method for Orig
}
}
diff --git a/algo/kf/core/geo/KfIMaterialMapFactory.h b/algo/kf/core/geo/KfIMaterialMapFactory.h
new file mode 100644
index 0000000000..403a4019aa
--- /dev/null
+++ b/algo/kf/core/geo/KfIMaterialMapFactory.h
@@ -0,0 +1,43 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file KfIMaterialMapFactory.h
+/// \brief Interface to the material map creator
+/// \since 28.08.2024
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+namespace cbm::algo::kf
+{
+ class MaterialMap;
+
+ /// \class IMaterialMapFactory
+ /// \brief Interface to the material map creator
+ class IMaterialMapFactory {
+ public:
+ /// \brief Destructor
+ virtual ~IMaterialMapFactory() = default;
+
+ /// \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
+ ///
+ virtual MaterialMap GenerateMaterialMap(double zRef, double zMin, double zMax, double xyMax, int nBinsDim) = 0;
+ };
+} // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfModuleIndexMap.cxx b/algo/kf/core/geo/KfModuleIndexMap.cxx
new file mode 100644
index 0000000000..bdb3a1d724
--- /dev/null
+++ b/algo/kf/core/geo/KfModuleIndexMap.cxx
@@ -0,0 +1,119 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file KfModuleIndexMap.cxx
+/// \brief A helper class to map external indices with the ones of KF-setup (implementation)
+/// \author Sergei Zharko <s.zharko@gsi.de>
+/// \since 16.09.2024
+
+#include "KfModuleIndexMap.h"
+
+namespace cbm::algo::kf
+{
+ // -------------------------------------------------------------------------------------------------------------------
+ //
+ void ModuleIndexMap::Reset()
+ {
+ fvLocToGlb.clear();
+ fvGlbToLoc.clear();
+ fvDetLocOffset.clear();
+ fvDetIntToExt.clear();
+ fvDetExtToInt.clear();
+ }
+
+ // -------------------------------------------------------------------------------------------------------------------
+ //
+ std::string ModuleIndexMap::ToString(int indentLevel) const
+ {
+ using fmt::format;
+ constexpr char IndentCh = '\t';
+ std::string indent(indentLevel, IndentCh);
+ std::stringstream msg;
+ msg << '\n' << indent << format("{:<50} ", "Detector internal to external index map: ");
+ for (auto id : fvDetIntToExt) {
+ msg << format("{:>8} ", id);
+ }
+ msg << '\n' << indent << format("{:<50} ", "Detector external to internal index map: ");
+ for (auto id : fvDetExtToInt) {
+ msg << format("{:>8} ", id);
+ }
+ msg << '\n' << indent << format("{:<50} ", "First global ID vs. internal detector: ");
+ for (auto id : fvDetLocOffset) {
+ msg << format("{:>8} ", id);
+ }
+ msg << '\n' << indent << format("{:<50} ", "Local ID to global ID: ");
+ for (auto id : fvLocToGlb) {
+ msg << format("{:>8} ", id);
+ }
+ msg << '\n' << indent << format("{:<50} ", "Geo ID to (detID:locID): ");
+ for (const auto& [iDetExt, iLoc] : fvGlbToLoc) {
+ msg << format("{:>3}:{:<3}", iDetExt, iLoc);
+ }
+
+ return msg.str();
+ }
+
+ // -------------------------------------------------------------------------------------------------------------------
+ //
+ ModuleIndexMap ModuleIndexMapFactory::MakeIndexMap() const
+ {
+ ModuleIndexMap indexMap;
+
+ auto& vDetIntToExt = indexMap.fvDetIntToExt;
+ auto& vDetExtToInt = indexMap.fvDetExtToInt;
+ auto& vDetLocOffset = indexMap.fvDetLocOffset;
+ auto& vLocToGlb = indexMap.fvLocToGlb;
+ auto& vGlbToLoc = indexMap.fvGlbToLoc;
+
+ // ------- Component layout
+ // Initialization of internal and external detector subsystem mapping
+ for (const auto& component : fvComponentLayers) {
+ auto iDetExt = component.fDetId;
+ if (iDetExt < 0) {
+ std::stringstream msg;
+ msg << "ModuleIndexMapFactory::MakeIndexMap: illegal detector subsystem enumeration: " << iDetExt
+ << " is out of range [0, " << defs::MaxNofDetSubsystems << "). Enties of the EDetID enumeration are "
+ << "required to be non-negative";
+ throw std::out_of_range(msg.str());
+ }
+ if (iDetExt >= static_cast<int>(vDetExtToInt.size())) {
+ vDetExtToInt.resize(iDetExt + 1, -1);
+ }
+
+ if (vDetExtToInt[iDetExt] == -1) { // Detector subsystem is not registered
+ // Assign an internal index to the detID:
+ vDetExtToInt[iDetExt] = static_cast<int>(vDetIntToExt.size());
+ vDetIntToExt.push_back(iDetExt);
+ }
+ }
+ { // Estimate offsets
+ int nDet = vDetIntToExt.size();
+ vDetLocOffset.resize(nDet + 1, 0);
+ for (const auto& component : fvComponentLayers) {
+ int pick = vDetExtToInt[component.fDetId] + 1;
+ vDetLocOffset[pick] = std::max(vDetLocOffset[pick], component.fLocId + 1);
+ }
+ for (int iDet{0}; iDet < nDet; ++iDet) {
+ vDetLocOffset[iDet + 1] += vDetLocOffset[iDet];
+ }
+ }
+
+ vLocToGlb.resize(vDetLocOffset.back(), -1);
+ vGlbToLoc.resize(fvComponentLayers.size());
+ for (auto itComp = fvComponentLayers.cbegin(); itComp != fvComponentLayers.cend(); ++itComp) {
+ const auto& component{*itComp};
+ int iDetInt = vDetExtToInt[component.fDetId];
+ int iLoc = component.fLocId;
+ int iGlb = std::distance(fvComponentLayers.begin(), itComp);
+ vLocToGlb[vDetLocOffset[iDetInt] + iLoc] = iGlb;
+ vGlbToLoc[iGlb] = std::pair(component.fDetId, iLoc);
+ }
+
+ return indexMap;
+ }
+
+ // -------------------------------------------------------------------------------------------------------------------
+ //
+ void ModuleIndexMapFactory::Reset() { fvComponentLayers.clear(); }
+} // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfModuleIndexMap.h b/algo/kf/core/geo/KfModuleIndexMap.h
new file mode 100644
index 0000000000..78c186de40
--- /dev/null
+++ b/algo/kf/core/geo/KfModuleIndexMap.h
@@ -0,0 +1,193 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file KfModuleIndexMap.h
+/// \brief A helper class to map external indices with the ones of KF-setup
+/// \author Sergei Zharko <s.zharko@gsi.de>
+/// \since 10.09.2024
+
+#pragma once
+
+#include "KfDefs.h"
+
+#include <boost/serialization/access.hpp>
+#include <boost/serialization/utility.hpp>
+
+#include <set>
+#include <sstream>
+#include <string>
+#include <vector>
+
+#include <fmt/format.h>
+
+namespace cbm::algo::kf
+{
+
+ /// \class ModuleIndexMap
+ /// \brief Maps local detector and station indices to the material maps and field slices
+ ///
+ /// The class maps external indices of components combined into subsets to a plain global index.
+ /// The indices both of the components inside a subset (iLoc) and the indices of the
+ /// subsets themselves (iDetExt) must be unique, but can be unsorted and can have holes in their
+ /// ranges. The global index (iGlob) is assigned automatically depending on the element ordering.
+ ///
+ /// The mapping is organized with the sequential contiguous containers (std::vector). The containers
+ /// used are:
+ /// (i) fvDetExtToInt [ iDetExt ] => iDetInt
+ /// size: max(iDetExt)
+ /// (ii) fvDetIntToExt [ iDetInt ] => iDetExt (int)
+ /// size: number of registered subsystems
+ /// (iii) fvDetLocOffset [ iDetInt ] => max number of elements for iDetInt - 1
+ /// size: max possible number of elements (including inactive/unused)
+ /// fvDetLocOffset[0] = 0, fvDetLocOffset.back() = n of max possible elements
+ /// (iv) fvGlbToLoc [ iGlob ] => pair(iDetInt, iLoc)
+ /// size: number of registered components
+ /// (v) fvLocToGlb [ fvDetLocOffset[iDetInt] + iLoc] => iGlob
+ /// size: max possible number of components (including inactive/unused)
+ ///
+ class alignas(VcMemAlign) ModuleIndexMap {
+ friend class ModuleIndexMapFactory;
+ friend class boost::serialization::access;
+
+ public:
+ /// \brief Gets total number of components
+ int GetNofLayers() const { return fvDetLocOffset.back(); }
+
+ /// \brief Gets number of layers for a particular detector subsystem
+ /// \tparam EDetID concrete index type of the det ID (can be either an enum, or an integral type)
+ /// \param detId Detector ID of the component
+ ///
+ /// \note The maximum possible number of components in the subset is returned. It accounts for
+ /// inactive/unused components.
+ template<class EDetID>
+ int GetNofLayers(EDetID detId) const
+ {
+ return fvDetLocOffset[fvDetExtToInt[static_cast<int>(detId)] + 1]
+ - fvDetLocOffset[fvDetExtToInt[static_cast<int>(detId)]];
+ }
+
+ /// \brief Converts internal layer index to pair (detID, locID)
+ /// \tparam EDetID concrete index type of the det ID (can be either an enum, or an integral type)
+ /// \param globId Internal layer index
+ /// \return pair(detector ID, local ID) of the layer
+ template<class EDetID>
+ std::pair<EDetID, int> GlobalToLocal(int globId) const;
+
+ /// \brief Converts external pair (detID, locID) to internal layer index
+ /// \tparam EDetID concrete index type of the det ID (can be either an enum, or an integral type)
+ /// \param locId Local ID of the layer
+ /// \param detId Detector ID of the component
+ /// \note if the (detId, locId) pair was not registered, returns -1
+ /// \return global index of the layer
+ template<class EDetID>
+ int LocalToGlobal(EDetID detId, int locId) const
+ {
+ return fvLocToGlb[fvDetLocOffset[fvDetExtToInt[static_cast<int>(detId)]] + locId];
+ }
+
+ /// \brief Disables a component
+ ///
+ ///
+ template<class EDetID>
+ void Disable(int locId, EDetID detId);
+
+ /// \brief Resets the instance
+ void Reset();
+
+ /// \brief String representation of the instance
+ std::string ToString(int indentLevel = 0) const;
+
+ private:
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int /*version*/)
+ {
+ ar& fvLocToGlb;
+ ar& fvGlbToLoc;
+ ar& fvDetLocOffset;
+ ar& fvDetIntToExt;
+ ar& fvDetExtToInt;
+ }
+
+ /// \brief Map of global component index to local (localId, intDetId)
+
+ std::vector<int> fvLocToGlb{};
+ std::vector<std::pair<int, int>> fvGlbToLoc{};
+ std::vector<int> fvDetLocOffset{}; ///< First index of component for det
+ std::vector<int> fvDetIntToExt{}; ///< Maps external detID to internal
+ std::vector<int> fvDetExtToInt{}; ///< Maps internal detID to external
+ };
+
+ /// \class ModuleIndexMapFactory
+ /// \brief Creates a valid module mapper
+ class ModuleIndexMapFactory {
+ /// \struct Component
+ /// \brief Structure to keep information on layers
+ struct Component {
+ Component(int detId, int locId, double z) : fZ(z), fLocId(locId), fDetId(detId) {}
+ double fZ; ///< Reference z-coordinate of the component
+ int fLocId; ///< Local index of component
+ int fDetId; ///< External index of detector subsystem
+ bool operator<(const Component& rhs) const { return fZ < rhs.fZ; }
+ };
+
+ public:
+ /// \brief Adds component info
+ /// \param detId Detector subsytem index (external)
+ /// \param locId Detector component local index (external)
+ /// \param z Reference z-component of the component [cm]
+ template<class EDetID>
+ void AddComponent(EDetID detId, int locId, double z);
+
+ /// \brief Creates a module index map
+ ModuleIndexMap MakeIndexMap() const;
+
+ /// \brief Resets the instance
+ void Reset();
+
+ private:
+ std::set<Component> fvComponentLayers;
+ };
+
+ // -------------------------------------------------------------------------------------------------------------------
+ //
+ template<class EDetID>
+ void ModuleIndexMap::Disable(int locIdDisable, EDetID detIdDisable)
+ {
+ // Check, if the detector id is there
+ auto iDetExtDsbl{static_cast<int>(detIdDisable)};
+ if (iDetExtDsbl >= fvDetExtToInt.size()) {
+ return; // Nothing to disable, detector is already not in the map
+ }
+
+ //auto iGeoDsbl{LocalToGeo(locIdDisable, detIdDisable)};
+ }
+
+ // -------------------------------------------------------------------------------------------------------------------
+ //
+ template<class EDetID>
+ inline std::pair<EDetID, int> ModuleIndexMap::GlobalToLocal(int globId) const
+ {
+ const auto& [iDetExt, iLoc] = fvGlbToLoc[globId];
+ return std::pair(static_cast<EDetID>(iDetExt), iLoc);
+ }
+
+ // -------------------------------------------------------------------------------------------------------------------
+ //
+ template<class EDetID>
+ void ModuleIndexMapFactory::AddComponent(EDetID detId, int locId, double z)
+ {
+ if (!fvComponentLayers.emplace(static_cast<int>(detId), locId, z).second) {
+ std::stringstream msg;
+ msg << "ModuleIndexMapFactory: attempt of adding a duplicating component with z = " << z
+ << ", detID = " << static_cast<int>(detId) << " and locId = " << locId
+ << ".\n The next components were already added:";
+ for (const auto& c : fvComponentLayers) {
+ msg << "\n\t- {z, detId, locId} = {" << c.fZ << ", " << c.fDetId << ", " << c.fLocId << "}";
+ }
+ throw std::logic_error(msg.str());
+ }
+ }
+
+
+} // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfSetup.cxx b/algo/kf/core/geo/KfSetup.cxx
index dd4fa54a17..3854fd48dc 100644
--- a/algo/kf/core/geo/KfSetup.cxx
+++ b/algo/kf/core/geo/KfSetup.cxx
@@ -25,8 +25,8 @@ std::string Setup<T>::ToString(int verbosity, int indentLevel) const
std::string indent(indentLevel, indentCh);
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 << "\nMODULE INDEXING SCHEME:\n" << fModuleIndexMap.ToString(indentLevel + 1);
msg << indent << "\nMATERIAL LAYERS:\n";
for (const auto& layer : fvMaterialLayers) {
msg << layer.ToString(indentLevel + 1, verbosity) << '\n';
diff --git a/algo/kf/core/geo/KfSetup.h b/algo/kf/core/geo/KfSetup.h
index 2f1a8a3d84..39ce0ca775 100644
--- a/algo/kf/core/geo/KfSetup.h
+++ b/algo/kf/core/geo/KfSetup.h
@@ -12,11 +12,11 @@
#include "KfDefs.h"
#include "KfField.h"
#include "KfMaterialMap.h"
+#include "KfModuleIndexMap.h"
#include "KfTarget.h"
#include "KfVector.h"
#include <boost/serialization/access.hpp>
-#include <boost/serialization/split_free.hpp>
#include <fstream>
#include <optional>
@@ -32,10 +32,10 @@ namespace cbm::algo::kf
template<typename T>
class alignas(VcMemAlign) Setup {
friend class boost::serialization::access; // Boost serializer methods
+ friend class SetupBuilder;
template<typename>
friend class Setup;
-
public:
/// \brief Constructor
/// \param fldMode Field mode
@@ -53,7 +53,8 @@ namespace cbm::algo::kf
/// \tparam I Underlying floating-point data type of the source
template<typename I>
Setup(const Setup<I>& other)
- : fvMaterialLayers(other.fvMaterialLayers)
+ : fModuleIndexMap(other.fModuleIndexMap)
+ , fvMaterialLayers(other.fvMaterialLayers)
, fField(other.fField)
, fTarget(other.fTarget)
{
@@ -68,10 +69,6 @@ namespace cbm::algo::kf
/// \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
/// \throw std::logic_error If the particular field member is undefined (nullopt)
const Field<T>& GetField() const { return fField; }
@@ -80,6 +77,19 @@ namespace cbm::algo::kf
/// \param iLayer Index of layer
const MaterialMap& GetMaterial(int iLayer) const { return fvMaterialLayers[iLayer]; }
+ /// \brief Gets material layer from external indices
+ /// \tparam EDetID Index of the detector subsystem
+ /// \param iDet DetectorID
+ /// \param iLoc Local index of the module
+ template<class EDetID>
+ const MaterialMap& GetMaterial(EDetID iDet, int iLoc) const
+ {
+ return fvMaterialLayers[fModuleIndexMap.LocalToGlobal(iDet, iLoc)];
+ }
+
+ /// \brief Gets module index map
+ const ModuleIndexMap& GetModuleIndexMap() const { return fModuleIndexMap; }
+
/// \brief Gets number of defined material layers
int GetNofMaterialLayers() const { return static_cast<int>(fvMaterialLayers.size()); }
@@ -91,40 +101,20 @@ namespace cbm::algo::kf
/// \param indentLevel Indent level of the string output
std::string ToString(int verbosity = 1, int indentLevel = 0) const;
- /// \brief Sets field
- void SetField(const Field<double>& field) { fField = Field<T>(field); }
-
- /// \brief Sets target (moves from the source)
- /// \param target Target instance
- void SetTarget(const Target<double>& target) { fTarget = target; }
-
private:
- /// \brief BOOST serialization load method
- template<class Archive>
- void load(Archive& ar, const unsigned int /*version*/)
- {
- ar >> fvMaterialLayers;
- ar >> fTarget;
- ar >> fField;
- }
-
- template<class Archive>
- void save(Archive& ar, const unsigned int /*version*/) const
- {
- ar << fvMaterialLayers;
- ar << fTarget;
- ar << fField;
- }
-
template<class Archive>
- void serialize(Archive& ar, const unsigned int version)
+ void serialize(Archive& ar, const unsigned int /*version*/)
{
- boost::serialization::split_member(ar, *this, version);
+ ar& fModuleIndexMap;
+ ar& fvMaterialLayers;
+ ar& fTarget;
+ ar& fField;
}
- std::vector<MaterialMap> fvMaterialLayers = {}; ///< Container of the material maps
- Field<T> fField; ///< Interpolated field (NOTE: maybe make optional)
- Target<T> fTarget; ///< Target layer
+ ModuleIndexMap fModuleIndexMap{}; ///< Index map (internal->external)
+ std::vector<MaterialMap> fvMaterialLayers{}; ///< Container of the material maps
+ Field<T> fField; ///< Interpolated field (NOTE: maybe make optional)
+ Target<T> fTarget; ///< Target layer
};
} // namespace cbm::algo::kf
diff --git a/algo/kf/core/geo/KfSetupBuilder.cxx b/algo/kf/core/geo/KfSetupBuilder.cxx
new file mode 100644
index 0000000000..e8cff76b4c
--- /dev/null
+++ b/algo/kf/core/geo/KfSetupBuilder.cxx
@@ -0,0 +1,108 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file KfSetupBuilder.h
+/// @brief A base KF-Setup initialization class (source)
+/// @since 28.03.2024
+/// @author Sergei Zharko <s.zharko@gsi.de>
+
+#include "KfSetupBuilder.h"
+
+#include <sstream>
+
+using cbm::algo::kf::SetupBuilder;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetupBuilder::Init()
+{
+ // Check setters
+ if (fGeoLayers.size() == 0) {
+ throw std::runtime_error("kf::SetupBuilder: no geometry layers initialized");
+ }
+ if (!fpMaterialMapFactory.get()) {
+ throw std::runtime_error("kf::SetupBuilder: no material map factory provided");
+ }
+
+ // Init target properties
+ {
+ // Material budget
+ double zMin{fTarget.GetZ() + kTargetCenterOffset};
+ double zMax{fTarget.GetZ() + fTarget.GetDz() * kTargetMaterialOffset};
+ double rMax{fTarget.GetR() * kTargetMaterialTransverseSizeMargin};
+ fTarget.SetMaterial(
+ fpMaterialMapFactory->GenerateMaterialMap(fTarget.GetZ(), zMin, zMax, rMax, kTargetMaterialMapNofBins));
+
+ // Field region
+ fFieldFactory.SetStep(kTargetFieldInitStep);
+ }
+
+ // Init geometry layers
+ {
+ // Estimate acceptance slope
+ double acceptanceSlope = -1.; // arbitrary negative value
+ for (const auto& layer : fGeoLayers) {
+ double tCurr = std::max(std::fabs(layer.fXmax), std::fabs(layer.fYmax)) / (layer.fZref - fTarget.GetZ());
+ assert(tCurr > 0);
+ acceptanceSlope = std::max(acceptanceSlope, tCurr);
+ }
+
+ fvMaterial.reserve(fGeoLayers.size());
+ double zLast{fTarget.GetZ() + fTarget.GetDz() * kTargetMaterialOffset};
+ for (auto layerIt = fGeoLayers.cbegin(); layerIt != fGeoLayers.cend(); ++layerIt) {
+ double z1 = layerIt->fZmax;
+ double z2 = z1;
+ if (std::next(layerIt) != fGeoLayers.cend()) {
+ z2 = std::next(layerIt)->fZmin;
+ }
+ double zNew{0.5 * (z1 + z2)};
+ double xyMax{acceptanceSlope * (layerIt->fZref - fTarget.GetZ())};
+ auto material{fpMaterialMapFactory->GenerateMaterialMap(layerIt->fZref, zLast, zNew, xyMax, fMatMapNofBins)};
+ fvMaterial.push_back(std::move(material));
+ // Note: square material maps to follow the material budget map regions
+ fFieldFactory.AddSliceReference(xyMax, xyMax, layerIt->fZref);
+ fModuleIndexFactory.AddComponent<int>(layerIt->fDetID, layerIt->fLocID, layerIt->fZref);
+ zLast = zNew;
+ }
+ }
+ LOG(info) << "Initialized";
+ fbPropertiesInitialized = true;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetupBuilder::Reset()
+{
+ fbPropertiesInitialized = false;
+ fFieldFactory.Reset();
+ fGeoLayers.clear();
+ fpMaterialMapFactory = nullptr;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetupBuilder::SetTargetProperty(double x, double y, double z, double dz, double r)
+{
+ fbPropertiesInitialized = false;
+ fTarget.SetX(x);
+ fTarget.SetY(y);
+ fTarget.SetZ(z);
+ fTarget.SetDz(dz);
+ fTarget.SetR(r);
+ fFieldFactory.SetTarget(x, y, z);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void SetupBuilder::Store(const Setup<double>& setup, const std::string& fileName)
+{
+ std::ofstream ofs(fileName, std::ios::binary);
+ if (!ofs) {
+ std::stringstream msg;
+ msg << "kf::SetupBuilder::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/KfSetupBuilder.h b/algo/kf/core/geo/KfSetupBuilder.h
new file mode 100644
index 0000000000..b67338ea5b
--- /dev/null
+++ b/algo/kf/core/geo/KfSetupBuilder.h
@@ -0,0 +1,233 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// @file KfSetupBuilder.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 "KfIMaterialMapFactory.h"
+#include "KfSetup.h"
+
+#include <boost/archive/binary_iarchive.hpp>
+#include <boost/archive/binary_oarchive.hpp>
+
+#include <fstream>
+#include <memory>
+#include <optional>
+#include <sstream>
+#include <string>
+#include <vector>
+
+namespace cbm::algo::kf
+{
+ /// \struct GeoLayer
+ /// \brief A helper structure to store geometrical information of the layers
+ template<class EDetID>
+ struct GeoLayer {
+ GeoLayer(EDetID detID, int iLoc, double zRef, double zMin, double zMax, double xMax, double yMax)
+ : fZref(zRef)
+ , fZmin(zMin)
+ , fZmax(zMax)
+ , fXmax(xMax)
+ , fYmax(yMax)
+ , fDetID(detID)
+ , fLocID(iLoc)
+ {
+ }
+
+ double fZref; ///< Reference z-coordinate [cm]
+ double fZmin; ///< Lower z-coordinate boundary [cm]
+ double fZmax; ///< Upper z-coordinate boundary [cm]
+ double fXmax; ///< Half size of the layer in the x-direction [cm]
+ double fYmax; ///< Half size of the layer in the y-direction [cm]
+ EDetID fDetID; ///< Index of detector subsystem
+ int fLocID; ///< Local index of the detector module
+ bool operator<(const GeoLayer& rhs) const { return fZref < rhs.fZref; }
+ };
+
+ /// \class SetupBuilder
+ /// \brief Creates a valid initialized Setup instance
+ class SetupBuilder {
+ public:
+ /// \brief Default constructor
+ SetupBuilder() = default;
+
+ /// \brief Copy constructor
+ SetupBuilder(const SetupBuilder&) = delete;
+
+ /// \brief Move constructor
+ SetupBuilder(SetupBuilder&&) = delete;
+
+ /// \brief Destructor
+ virtual ~SetupBuilder() = default;
+
+ /// \brief Copy assignment operator
+ SetupBuilder& operator=(const SetupBuilder&) = delete;
+
+ /// \brief Move assignment operator
+ SetupBuilder& operator=(SetupBuilder&&) = delete;
+
+ /// \brief Adds material layer
+ /// \tparam EDetID detector ID index type
+ /// \param geoLayer A geometrical layer
+ template<class EDetID>
+ void AddLayer(const GeoLayer<EDetID>& geoLayer);
+
+ /// \brief Initializes, validates and cashes the parameters
+ /// \throw std::runtime_error If pre-initialization was incomplete
+ void Init();
+
+ /// \brief Creates a setup instance
+ template<typename T>
+ Setup<T> MakeSetup() const;
+
+ /// \brief Resets the instance
+ void Reset();
+
+ /// \brief Sets magnetic field function
+ /// \param fieldMode Magnetic field mode:
+ /// \param fieldFn Magnetic field function
+ /// \param fieldType Magnetic field type
+ void SetFieldProperty(EFieldMode fieldMode, const FieldFn_t& fieldFn, EFieldType fieldType)
+ {
+ fbPropertiesInitialized = false;
+ fFieldFactory.SetFieldFunction(fieldFn, fieldType);
+ fFieldFactory.SetFieldMode(fieldMode);
+ }
+
+ /// \brief Sets material map creator
+ /// \param pMaterialFactory Pointer to the actual material map creator instance
+ void SetMaterialMapFactory(const std::shared_ptr<IMaterialMapFactory>& pMaterialFactory)
+ {
+ fbPropertiesInitialized = false;
+ fpMaterialMapFactory = pMaterialFactory;
+ }
+
+ /// \brief Sets target initialization properties
+ /// \param x Target x-coordinate [cm]
+ /// \param y Target y-coordinate [cm]
+ /// \param z Target z-coordinate [cm]
+ /// \param dz Target half-thickness [cm]
+ /// \param r Target transverse size (XYmax)
+ void SetTargetProperty(double x, double y, double z, double dz, double r);
+
+
+ // ********************
+ // ** Static methods **
+ // ********************
+
+ /// \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:
+ // TODO: Define target material more precisely
+ static constexpr double kTargetCenterOffset{0.05}; // Offset from target center [cm]
+ static constexpr double kTargetMaterialOffset{2.}; // Offset of the target material [in dz]
+ static constexpr double kTargetMaterialTransverseSizeMargin{1.3}; // Margin of target transverse size
+ static constexpr double kTargetFieldInitStep{2.5}; // Step between nodal points to determine field near target [cm]
+ static constexpr int kTargetMaterialMapNofBins{20};
+
+ std::set<GeoLayer<int>> fGeoLayers{}; ///< Set of geo layers
+ std::vector<MaterialMap> fvMaterial{}; ///< Material map container
+ std::shared_ptr<IMaterialMapFactory> fpMaterialMapFactory{nullptr}; ///< Material map creator
+ ModuleIndexMapFactory fModuleIndexFactory; ///< Module index factory
+ FieldFactory fFieldFactory; ///< Instance of field factory
+ Target<double> fTarget; ///< Target properties
+ int fMatMapNofBins{100}; ///< Number of bins in material maps
+ bool fbPropertiesInitialized{false}; ///< Flag: if all the necessary fields are initialized
+ };
+
+
+ // ********************************
+ // ** Template method definition **
+ // ********************************
+
+
+ // ---------------------------------------------------------------------------------------------------------------------
+ //
+ template<class EDetID>
+ void SetupBuilder::AddLayer(const GeoLayer<EDetID>& geoLayer)
+ {
+ fbPropertiesInitialized = false;
+ auto layer = fGeoLayers.emplace(static_cast<int>(geoLayer.fDetID), geoLayer.fLocID, geoLayer.fZref, geoLayer.fZmin,
+ geoLayer.fZmax, geoLayer.fXmax, geoLayer.fYmax);
+ if (!layer.second) {
+ std::stringstream msg;
+ msg << "SetupBuilder::AddLayer: attempt of adding a duplicating geometry layer: ";
+ msg << "DetID = " << static_cast<int>(geoLayer.fDetID) << ", localID = " << geoLayer.fLocID
+ << "fZref = " << geoLayer.fZref << ", fZmin = " << geoLayer.fZmin << ", fZmax = " << geoLayer.fZmax
+ << ", fXmax = " << geoLayer.fXmax << ", fYmax = " << geoLayer.fYmax
+ << ".\nThe next layers were already added:";
+ for (const auto& el : fGeoLayers) {
+ msg << "\n\t- DetID = " << static_cast<int>(el.fDetID) << ", localID = " << el.fLocID << "fZref = " << el.fZref
+ << ", fZmin = " << el.fZmin << ", fZmax = " << el.fZmax << ", fXmax = " << el.fXmax
+ << ", fYmax = " << el.fYmax;
+ }
+ throw std::logic_error(msg.str());
+ }
+ }
+
+ // ---------------------------------------------------------------------------------------------------------------------
+ //
+ template<typename T>
+ Setup<T> SetupBuilder::MakeSetup() const
+ {
+ Setup<T> setup(fFieldFactory.GetFieldMode());
+
+ if (!fbPropertiesInitialized) {
+ throw std::logic_error("kf::SetupBuilder: attempt to create an uninitialized setup");
+ }
+
+ setup.fTarget = this->fTarget;
+ for (const auto& material : this->fvMaterial) {
+ setup.fvMaterialLayers.push_back(material);
+ }
+ setup.fField = fFieldFactory.MakeField<T>();
+ setup.fModuleIndexMap = fModuleIndexFactory.MakeIndexMap();
+ return setup;
+ }
+
+ // -------------------------------------------------------------------------------------------------------------------
+ //
+ template<typename T>
+ Setup<T> SetupBuilder::Load(const std::string& fileName)
+ {
+ Setup<double> setup(EFieldMode::Intrpl);
+ std::ifstream ifs(fileName, std::ios::binary);
+ if (!ifs) {
+ std::stringstream msg;
+ msg << "kf::SetupBuilder::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::SetupBuilder::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
deleted file mode 100644
index 3148dc428d..0000000000
--- a/algo/kf/core/geo/KfSetupInitializer.cxx
+++ /dev/null
@@ -1,63 +0,0 @@
-/* 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
deleted file mode 100644
index 1f2f51e51c..0000000000
--- a/algo/kf/core/geo/KfSetupInitializer.h
+++ /dev/null
@@ -1,166 +0,0 @@
-/* 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
- template<typename T>
- Setup<T> MakeSetup() const;
-
- /// \brief Resets the instance
- void Reset();
-
- /// \brief Sets magnetic field function
- /// \param fieldMode Magnetic field mode:
- /// \param fieldFn Magnetic field function
- /// \param fieldType Magnetic field type
- void SetFieldProperty(EFieldMode fieldMode, const FieldFn_t& fieldFn, EFieldType fieldType)
- {
- fFieldFactory.SetFieldFunction(fieldFn, fieldType);
- fFieldFactory.SetFieldMode(fieldMode);
- }
-
- /// \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() const
- {
- Setup<T> setup(fFieldFactory.GetFieldMode());
- // Target initialization
- setup.SetTarget(fTarget);
-
- // Field initialization
- setup.SetField(fFieldFactory.MakeField<double>());
-
- // 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(EFieldMode::Intrpl);
- 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 6f4864d125..bc8b00a188 100644
--- a/algo/kf/core/geo/KfTarget.cxx
+++ b/algo/kf/core/geo/KfTarget.cxx
@@ -17,11 +17,12 @@ using cbm::algo::kf::Target;
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename T>
-Target<T>::Target(const MaterialMap& material, double x, double y, double z)
- : fMaterial(material)
- , fX(utils::simd::Cast<double, T>(x))
+Target<T>::Target(double x, double y, double z, double dz, double r)
+ : fX(utils::simd::Cast<double, T>(x))
, fY(utils::simd::Cast<double, T>(y))
, fZ(utils::simd::Cast<double, T>(z))
+ , fDz(utils::simd::Cast<double, T>(dz))
+ , fR(utils::simd::Cast<double, T>(r))
{
}
@@ -46,6 +47,7 @@ std::string Target<T>::ToString(int indentLevel, int verbose) const
std::string indent(indentLevel, IndentChar);
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 << "half-thickness: " << Cnv(fDz) << " [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 51ab42ed1c..e901c6c58f 100644
--- a/algo/kf/core/geo/KfTarget.h
+++ b/algo/kf/core/geo/KfTarget.h
@@ -31,11 +31,12 @@ namespace cbm::algo::kf
Target() = default;
/// \brief Constructor from parameters
- /// \param material Instance of the material budget map near the target region
/// \param x x-coordinate of the nominal target center [cm]
/// \param y y-coordinate of the nominal target center [cm]
/// \param z z-coordinate of the nominal target center [cm]
- Target(const MaterialMap& material, double x, double y, double z);
+ /// \param dz half-thickness of the target [cm]
+ /// \param r size of target
+ Target(double x, double y, double z, double dz, double r);
/// \brief Copy constructor
/// \tparam I Underlying floating-point type of the source
@@ -45,6 +46,8 @@ namespace cbm::algo::kf
, fX(utils::simd::Cast<I, T>(other.fX))
, fY(utils::simd::Cast<I, T>(other.fY))
, fZ(utils::simd::Cast<I, T>(other.fZ))
+ , fDz(utils::simd::Cast<I, T>(other.fDz))
+ , fR(utils::simd::Cast<I, T>(other.fR))
{
}
@@ -57,12 +60,18 @@ namespace cbm::algo::kf
/// \brief Gets x-coordinate of the nominal target center
const T& GetX() const { return fX; }
- /// \brief Gets x-coordinate of the nominal target center
+ /// \brief Gets y-coordinate of the nominal target center
const T& GetY() const { return fY; }
- /// \brief Gets x-coordinate of the nominal target center
+ /// \brief Gets z-coordinate of the nominal target center
const T& GetZ() const { return fZ; }
+ /// \brief Gets target half-thickness
+ const T& GetDz() const { return fDz; }
+
+ /// \brief Gets transverse size of target
+ const T& GetR() const { return fR; }
+
/// \brief Gets material map
const MaterialMap& GetMaterial() const { return fMaterial; }
@@ -82,6 +91,14 @@ namespace cbm::algo::kf
/// \param z x-coordinate [cm]
void SetZ(const T& z) { fZ = z; }
+ /// \brief Sets target half-thickness
+ /// \param dz half-thickness [cm]
+ void SetDz(const T& dz) { fDz = dz; }
+
+ /// \brief Sets target transverse size
+ /// \brief Transverse size of target [cm]
+ void SetR(const T& r) { fR = r; }
+
/// \brief String representation of the class
/// \param indentLevel Indent level of the string output
/// \param verbose Verbosity level
@@ -97,11 +114,15 @@ namespace cbm::algo::kf
ar& fX;
ar& fY;
ar& fZ;
+ ar& fDz;
+ ar& fR;
}
MaterialMap fMaterial{}; ///< Material map in the target region
T fX{defs::Undef<T>}; ///< x-coordinate of the nominal target center [cm]
T fY{defs::Undef<T>}; ///< y-coordinate of the nominal target center [cm]
T fZ{defs::Undef<T>}; ///< z-coordinate of the nominal target center [cm]
+ T fDz{defs::Undef<T>}; ///< Half-thickness of the target [cm]
+ T fR{defs::Undef<T>}; ///< Half-thickness of the target [cm]
};
} // namespace cbm::algo::kf
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 8201ea78ec..8d8b26c1ae 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -162,7 +162,6 @@ set(PRIVATE_DEPENDENCIES
CbmRecoSts # <==== !!!! was bringing dependency on the "online" version of the Algo library which brought the NO_ROOT flag !!!!!
CbmSimSteer
CbmRecoBase
- KfCbm
KFParticle
external::yaml-cpp
FairRoot::GeoBase
diff --git a/reco/L1/CbmKfTrackingSetupInitializer.cxx b/reco/L1/CbmKfTrackingSetupInitializer.cxx
deleted file mode 100644
index 9650ac1588..0000000000
--- a/reco/L1/CbmKfTrackingSetupInitializer.cxx
+++ /dev/null
@@ -1,195 +0,0 @@
-/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergei Zharko [committer] */
-
-/// \file CbmKfTrackingSetupInitializer.cxx
-/// \brief Tracking setup initializer in CBM (source)
-/// \since 28.08.2024
-/// \author Sergei Zharko <s.zharko@gsi.de>
-
-#include "CbmKfTrackingSetupInitializer.h"
-
-#include "CbmKfOriginalField.h"
-#include "CbmKfTarget.h"
-#include "CbmMuchTrackingInterface.h"
-#include "CbmMvdTrackingInterface.h"
-#include "CbmSetup.h"
-#include "CbmStsTrackingInterface.h"
-#include "CbmTofTrackingInterface.h"
-#include "CbmTrdTrackingInterface.h"
-#include "FairField.h"
-#include "FairRunAna.h"
-#include "KfDefs.h"
-#include "KfMaterialMapCreator.h"
-#include "Logger.h"
-
-#include <functional>
-
-using cbm::kf::TrackingSetupInitializer;
-using kf::tools::MaterialMapCreator;
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-bool TrackingSetupInitializer::InitProperties()
-try {
- using cbm::algo::kf::EFieldType;
-
- // Check, if a subsystem exists in the setup
- fbUseMvd &= CbmSetup::Instance()->IsActive(ECbmModuleId::kMvd);
- fbUseSts &= CbmSetup::Instance()->IsActive(ECbmModuleId::kSts);
- fbUseMuch &= CbmSetup::Instance()->IsActive(ECbmModuleId::kMuch);
- fbUseTrd &= CbmSetup::Instance()->IsActive(ECbmModuleId::kTrd);
- fbUseTof &= CbmSetup::Instance()->IsActive(ECbmModuleId::kTof);
-
- // Magnetic field initialization
- if (FairRunAna::Instance()->GetField()) {
- this->SetFieldFunction(OriginalField(), EFieldType::Normal);
- }
- else {
- //this->SetFieldFunction(cbm::algo::kf::defs::ZeroFieldFn, EFieldType::Null);
- this->SetFieldFunction(ZeroField(), EFieldType::Null);
- }
-
- // Tracking station property initialization
- {
- // Geometrical layer (contains information on the station sizes)
- //
- struct GeoLayer {
- double zRef; // Reference z-coordinate [cm]
- double zMin; // Lower z-coordinate boundary [cm]
- double zMax; // Upper z-coordinate boundary [cm]
- double xMax; // Half size of the layer in the x-direction [cm]
- double yMax; // Half size of the layer in the y-direction [cm]
- bool operator<(const GeoLayer& rhs) const { return zRef < rhs.zRef; }
- };
- std::set<GeoLayer> geoStations;
-
- auto CollectStations = [&](const auto* pIfs) -> void {
- if (!pIfs) {
- return;
- }
- for (int iSt = 0; iSt < pIfs->GetNtrackingStations(); ++iSt) {
- geoStations.emplace(GeoLayer{.zRef = pIfs->GetZref(iSt),
- .zMin = pIfs->GetZmin(iSt),
- .zMax = pIfs->GetZmax(iSt),
- .xMax = std::max(std::fabs(pIfs->GetXmin(iSt)), std::fabs(pIfs->GetXmax(iSt))),
- .yMax = std::max(std::fabs(pIfs->GetYmin(iSt)), std::fabs(pIfs->GetYmax(iSt)))});
- }
- };
- CollectStations(fbUseMvd ? CbmMvdTrackingInterface::Instance() : nullptr);
- CollectStations(fbUseSts ? CbmStsTrackingInterface::Instance() : nullptr);
- CollectStations(fbUseMuch ? CbmMuchTrackingInterface::Instance() : nullptr);
- CollectStations(fbUseTrd ? CbmTrdTrackingInterface::Instance() : nullptr);
- CollectStations(fbUseTof ? CbmTofTrackingInterface::Instance() : nullptr);
-
- // Retriev target properties
- const auto* pTarget = Target::Instance();
-
- // Init material map creator
- MaterialMapCreator materialCreator{};
- materialCreator.SetSafeMaterialInitialization(kMatCreatorSafeMode);
- materialCreator.SetDoRadialProjection(pTarget->GetZ());
- materialCreator.SetNraysPerDim(kMatCreatorNrays);
-
- // Target initialization
- double zLast =
- pTarget->GetZ() + pTarget->GetDz() + kTargMaterialOffset; // z-coordinate of the target material upper limit
- {
- // Material of the target
- 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
- 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;
- if (std::next(layerIt) != geoStations.end()) {
- z2 = std::next(layerIt)->zMin;
- }
- double zNew = 0.5 * (z1 + z2);
- double xyMax = 1.3 * std::max(std::fabs(layerIt->xMax), std::fabs(layerIt->yMax));
- int nBins = std::min(static_cast<int>(std::ceil(2. * xyMax / kMatCreatorPitch)), kMatCreatorMaxNbins);
- assert(nBins > 0);
-
- LOG(info) << "- collecting material between z = " << zLast << " and " << zNew << " [cm]";
- this->AddMaterial(materialCreator.GenerateMaterialMap(layerIt->zRef, zLast, zNew, xyMax, nBins));
- 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) {
- LOG(error) << "kf::TrackingSetupInitializer: Tracking setup was not initialized. Reason: " << err.what();
- return false;
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void TrackingSetupInitializer::Use(bool mvd, bool sts, bool much, bool trd, bool tof)
-{
- fbUseMvd = mvd;
- fbUseSts = sts;
- fbUseMuch = much;
- fbUseTrd = trd;
- fbUseTof = tof;
-}
diff --git a/reco/L1/CbmKfTrackingSetupInitializer.h b/reco/L1/CbmKfTrackingSetupInitializer.h
deleted file mode 100644
index 50e643bb15..0000000000
--- a/reco/L1/CbmKfTrackingSetupInitializer.h
+++ /dev/null
@@ -1,53 +0,0 @@
-/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergei Zharko [committer] */
-
-/// \file CbmKfTrackingSetupInitializer.h
-/// \brief Tracking setup initializer in CBM (source)
-/// \since 28.08.2024
-/// \author Sergei Zharko <s.zharko@gsi.de>
-
-#pragma once
-
-#include "KfSetupInitializer.h"
-
-#include <tuple>
-
-namespace cbm::kf
-{
- /// \class TrackingSetupInitializer
- /// \brief Encapsulation of the kf::Setup initialization routines for CBM
- class TrackingSetupInitializer : public cbm::algo::kf::SetupInitializer {
- public:
- /// \brief Initializes the instance
- /// \param mvd Is MVD used
- /// \param sts Is STS used
- /// \param much Is MuCh used
- /// \param trd Is TRD used
- /// \param tof Is TOF used
- bool InitProperties() override;
-
- /// \brief Enables/disables detector subsystems in the setup
- /// \param mvd Is MVD used
- /// \param sts Is STS used
- /// \param much Is MuCh used
- /// \param trd Is TRD used
- /// \param tof Is TOF used
- void Use(bool mvd, bool sts, bool much, bool trd, bool tof);
-
- 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{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
- bool fbUseMuch{false}; ///< Are MuCh stations included in the tracking setup
- bool fbUseTrd{false}; ///< Are TRD stations included in the tracking setup
- bool fbUseTof{false}; ///< Are TOF stations included in the tracking setup
- };
-} // namespace cbm::kf
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index ca54da6acc..55488de212 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -544,17 +544,18 @@ try {
// Creating a setup
auto pSetupInitializer = std::make_unique<TrackingSetupInitializer>();
pSetupInitializer->Use(fUseMVD, fUseSTS, fUseMUCH, fUseTRD, fUseTOF);
- pSetupInitializer->InitProperties();
+ pSetupInitializer->Init(cbm::algo::kf::EFieldMode::Intrpl);
auto trackerSetup = pSetupInitializer->MakeSetup<double>();
- LOG(info) << "!!!!\n!!!!\n!!!!\n!!!! KF SETUP:\n" << trackerSetup.ToString(0) << "\n!!!!\n!!!!\n!!!!";
+ LOG(info) << "!!!!\n!!!!\n!!!!\n!!!! KF SETUP:\n"
+ << trackerSetup.ToString(0) << "\nsize: " << sizeof(trackerSetup) << "\n!!!!\n!!!!\n!!!!";
// Storing setup
- TrackingSetupInitializer::Store(trackerSetup, "./trackerSetup.geo.kf.bin");
+ cbm::algo::kf::SetupBuilder::Store(trackerSetup, "./trackerSetup.geo.kf.bin");
// Reading setup (now in fvec):
- auto trackerSetupVec = TrackingSetupInitializer::Load<fvec>("./trackerSetup.geo.kf.bin");
+ auto trackerSetupVec = cbm::algo::kf::SetupBuilder::Load<fvec>("./trackerSetup.geo.kf.bin");
LOG(info) << "!!!!\n!!!!\n!!!!\n!!!! KF SETUP (from file):\n"
- << trackerSetupVec.ToString(1) << "\n!!!!\n!!!!\n!!!!";
+ << trackerSetupVec.ToString(1) << "\nsize:" << sizeof(trackerSetupVec) << "\n!!!!\n!!!!\n!!!!";
{
diff --git a/reco/L1/CbmL1DetectorID.h b/reco/L1/CbmL1DetectorID.h
index 5ecdda9522..72258838ac 100644
--- a/reco/L1/CbmL1DetectorID.h
+++ b/reco/L1/CbmL1DetectorID.h
@@ -10,8 +10,8 @@
#pragma once
//#include "CaDefs.h"
-#include "CbmEnumArray.h"
#include "CbmDefs.h"
+#include "CbmEnumArray.h"
#include <string>
diff --git a/reco/L1/CbmL1Util.cxx b/reco/L1/CbmL1Util.cxx
index 393b966555..745f6db483 100644
--- a/reco/L1/CbmL1Util.cxx
+++ b/reco/L1/CbmL1Util.cxx
@@ -3,8 +3,9 @@
Authors: Sergey Gorbunov [committer] */
#include "CbmL1Util.h"
-#include "FairTrackParam.h"
+
#include "CaDefs.h"
+#include "FairTrackParam.h"
namespace cbm::L1Util
{
diff --git a/reco/kfnew/CbmKfTrackingSetupInitializer.cxx b/reco/kfnew/CbmKfTrackingSetupInitializer.cxx
index 107431d132..bfb13b1c3f 100644
--- a/reco/kfnew/CbmKfTrackingSetupInitializer.cxx
+++ b/reco/kfnew/CbmKfTrackingSetupInitializer.cxx
@@ -20,22 +20,21 @@
#include "FairField.h"
#include "FairRunAna.h"
#include "KfDefs.h"
-#include "KfMaterialMapCreator.h"
+#include "KfMaterialMapFactory.h"
#include "Logger.h"
-#include "KfModuleIndexMap.h"
-
#include <functional>
using cbm::kf::TrackingSetupInitializer;
-using kf::tools::MaterialMapCreator;
+using kf::tools::MaterialMapFactory;
// ---------------------------------------------------------------------------------------------------------------------
//
-bool TrackingSetupInitializer::InitProperties()
+bool TrackingSetupInitializer::Init(cbm::algo::kf::EFieldMode fldMode)
try {
- using cbm::algo::kf::EFieldMode;
using cbm::algo::kf::EFieldType;
+ using cbm::algo::kf::GeoLayer;
+ fBuilder.Reset();
// Check, if a subsystem exists in the setup
fbUseMvd &= CbmSetup::Instance()->IsActive(ECbmModuleId::kMvd);
@@ -46,146 +45,51 @@ try {
// Magnetic field initialization
if (FairRunAna::Instance()->GetField()) {
- this->SetFieldProperty(EFieldMode::Intrpl, OriginalField(), EFieldType::Normal);
+ fBuilder.SetFieldProperty(fldMode, OriginalField(), EFieldType::Normal);
}
else {
- this->SetFieldProperty(EFieldMode::Intrpl, ZeroField(), EFieldType::Null);
+ fBuilder.SetFieldProperty(fldMode, ZeroField(), EFieldType::Null);
}
// Tracking station property initialization
- {
- // Geometrical layer (contains information on the station sizes)
- //
- struct GeoLayer {
- double zRef; // Reference z-coordinate [cm]
- double zMin; // Lower z-coordinate boundary [cm]
- double zMax; // Upper z-coordinate boundary [cm]
- double xMax; // Half size of the layer in the x-direction [cm]
- double yMax; // Half size of the layer in the y-direction [cm]
- ETrackingDetID detID; ///< Index of tracking detector
- int iStLoc; ///< Local station of the layer
-
- bool operator<(const GeoLayer& rhs) const { return zRef < rhs.zRef; }
- };
- std::set<GeoLayer> geoStations;
-
- auto CollectStations = [&](ETrackingDetID detID, const auto* pIfs) -> void {
- if (!pIfs) {
- return;
- }
- for (int iSt = 0; iSt < pIfs->GetNtrackingStations(); ++iSt) {
- geoStations.emplace(GeoLayer{.zRef = pIfs->GetZref(iSt),
- .zMin = pIfs->GetZmin(iSt),
- .zMax = pIfs->GetZmax(iSt),
- .xMax = std::max(std::fabs(pIfs->GetXmin(iSt)), std::fabs(pIfs->GetXmax(iSt))),
- .yMax = std::max(std::fabs(pIfs->GetYmin(iSt)), std::fabs(pIfs->GetYmax(iSt))),
- .detID = detID,
- .iStLoc = iSt});
- }
- };
- CollectStations(ETrackingDetID::Mvd, fbUseMvd ? CbmMvdTrackingInterface::Instance() : nullptr);
- CollectStations(ETrackingDetID::Sts, fbUseSts ? CbmStsTrackingInterface::Instance() : nullptr);
- CollectStations(ETrackingDetID::Much, fbUseMuch ? CbmMuchTrackingInterface::Instance() : nullptr);
- CollectStations(ETrackingDetID::Trd, fbUseTrd ? CbmTrdTrackingInterface::Instance() : nullptr);
- CollectStations(ETrackingDetID::Tof, fbUseTof ? CbmTofTrackingInterface::Instance() : nullptr);
-
- // Retriev target properties
- const auto* pTarget = Target::Instance();
-
- // Init material map creator
- MaterialMapCreator materialCreator{};
- materialCreator.SetSafeMaterialInitialization(kMatCreatorSafeMode);
- materialCreator.SetDoRadialProjection(pTarget->GetZ());
- materialCreator.SetNraysPerDim(kMatCreatorNrays);
-
- // Target initialization
- double zLast =
- pTarget->GetZ() + pTarget->GetDz() + kTargMaterialOffset; // z-coordinate of the target material upper limit
- {
- // Material of the target
- 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
- 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));
+ auto CollectStations = [&](const auto* pIfs, ETrackingDetID detID) -> void {
+ if (!pIfs) {
+ return;
}
-
-
-
- // Material and field layers initialization
-
- cbm::algo::kf::ModuleIndexMapInitializer<ETrackingDetID> indexMapInitializer;
- 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;
+ for (int iSt = 0; iSt < pIfs->GetNtrackingStations(); ++iSt) {
+ if (iSt == 1 && detID == ETrackingDetID::Mvd || iSt == 2 && detID == ETrackingDetID::Sts) {
+ continue; // !!!!!!!!!!!!! TEST
}
- double zNew = 0.5 * (z1 + z2);
- double xyMax = 1.3 * std::max(std::fabs(layerIt->xMax), std::fabs(layerIt->yMax));
- int nBins = std::min(static_cast<int>(std::ceil(2. * xyMax / kMatCreatorPitch)), kMatCreatorMaxNbins);
- assert(nBins > 0);
-
- LOG(info) << "- collecting material between z = " << zLast << " and " << zNew << " [cm]";
- this->AddMaterial(materialCreator.GenerateMaterialMap(layerIt->zRef, zLast, zNew, xyMax, nBins));
- this->AddFieldSliceReference(layerIt->xMax, layerIt->yMax, layerIt->zRef);
- indexMapInitializer.AddComponent(layerIt->detID, layerIt->iStLoc, layerIt->zRef);
- zLast = zNew;
+ fBuilder.AddLayer(
+ GeoLayer<ETrackingDetID>{.fDetID = detID, // ca::Tracking detector id scheme
+ .fLocID = iSt, // ca::Tracking station indexing
+ .fZref = pIfs->GetZref(iSt),
+ .fZmin = pIfs->GetZmin(iSt),
+ .fZmax = pIfs->GetZmax(iSt),
+ .fXmax = std::max(std::fabs(pIfs->GetXmin(iSt)), std::fabs(pIfs->GetXmax(iSt))),
+ .fYmax = std::max(std::fabs(pIfs->GetYmin(iSt)), std::fabs(pIfs->GetYmax(iSt)))});
}
+ };
+ CollectStations(fbUseMvd ? CbmMvdTrackingInterface::Instance() : nullptr, ETrackingDetID::Mvd);
+ CollectStations(fbUseSts ? CbmStsTrackingInterface::Instance() : nullptr, ETrackingDetID::Sts);
+ CollectStations(fbUseMuch ? CbmMuchTrackingInterface::Instance() : nullptr, ETrackingDetID::Much);
+ CollectStations(fbUseTrd ? CbmTrdTrackingInterface::Instance() : nullptr, ETrackingDetID::Trd);
+ CollectStations(fbUseTof ? CbmTofTrackingInterface::Instance() : nullptr, ETrackingDetID::Tof);
- auto indexMap = indexMapInitializer.MakeIndexMap();
- LOG(info) << "---> Index map: " << indexMap.ToString();
+ // Retriev target properties
+ const auto* pTarget = Target::Instance();
+ fBuilder.SetTargetProperty(pTarget->GetX(), pTarget->GetY(), pTarget->GetZ(), pTarget->GetDz(), pTarget->GetRmax());
+ // Init material map creator
+ auto pMaterialFactory = std::make_shared<MaterialMapFactory>();
+ pMaterialFactory->SetSafeMaterialInitialization(kMatCreatorSafeMode);
+ pMaterialFactory->SetDoRadialProjection(pTarget->GetZ());
+ pMaterialFactory->SetNraysPerDim(kMatCreatorNrays);
+ fBuilder.SetMaterialMapFactory(pMaterialFactory);
- 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);
- }
- }
- }
+ fBuilder.Init();
+
LOG(info) << "kf::TrackingSetupInitializer: Tracking setup was initialized successfully";
return true;
}
diff --git a/reco/kfnew/CbmKfTrackingSetupInitializer.h b/reco/kfnew/CbmKfTrackingSetupInitializer.h
index 42ff65f9f0..73fc04f2a3 100644
--- a/reco/kfnew/CbmKfTrackingSetupInitializer.h
+++ b/reco/kfnew/CbmKfTrackingSetupInitializer.h
@@ -9,7 +9,7 @@
#pragma once
-#include "KfSetupInitializer.h"
+#include "KfSetupBuilder.h"
#include <tuple>
@@ -21,26 +21,23 @@ namespace cbm::kf
// TODO: SZh 12.09.2024: For me it's unclear, if we should move this entire file to CbmL1, or move the tracking
// detector IDs declaration here. For now this enum should fully follow the definition
// of cbm::algo::ca::EDetectorID.
- enum class ETrackingDetID {
- Mvd = 0x10,
- Sts = 0x20,
+ enum class ETrackingDetID
+ {
+ Mvd = 0,
+ Sts,
Much,
- Trd = 0x40,
- Tof = 0x30,
+ Trd,
+ Tof,
END
};
/// \class TrackingSetupInitializer
/// \brief Encapsulation of the kf::Setup initialization routines for CBM
- class TrackingSetupInitializer : public cbm::algo::kf::SetupInitializer {
+ class TrackingSetupInitializer {
public:
/// \brief Initializes the instance
- /// \param mvd Is MVD used
- /// \param sts Is STS used
- /// \param much Is MuCh used
- /// \param trd Is TRD used
- /// \param tof Is TOF used
- bool InitProperties() override;
+ /// \param fldMode Field mode
+ bool Init(cbm::algo::kf::EFieldMode fldMode);
/// \brief Enables/disables detector subsystems in the setup
/// \param mvd Is MVD used
@@ -50,7 +47,16 @@ namespace cbm::kf
/// \param tof Is TOF used
void Use(bool mvd, bool sts, bool much, bool trd, bool tof);
+ /// \brief Makes setup object
+ template<typename T>
+ cbm::algo::kf::Setup<T> MakeSetup()
+ {
+ return fBuilder.MakeSetup<T>();
+ }
+
private:
+ cbm::algo::kf::SetupBuilder fBuilder{};
+
// 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
diff --git a/reco/kfnew/tools/CMakeLists.txt b/reco/kfnew/tools/CMakeLists.txt
index 3a072d4a73..2f4a9baf52 100644
--- a/reco/kfnew/tools/CMakeLists.txt
+++ b/reco/kfnew/tools/CMakeLists.txt
@@ -3,7 +3,7 @@ set(INCLUDE_DIRECTORIES
)
set(SRCS
- ${CMAKE_CURRENT_SOURCE_DIR}/KfMaterialMapCreator.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/KfMaterialMapFactory.cxx
)
SET_SOURCE_FILES_PROPERTIES(${SRCS} PROPERTIES COMPILE_FLAGS "-O3")
@@ -55,7 +55,7 @@ install(TARGETS KfTools DESTINATION lib)
install(
FILES
- KfMaterialMapCreator.h
+ KfMaterialMapFactory.h
DESTINATION
include/
)
diff --git a/reco/kfnew/tools/KfMaterialMapCreator.cxx b/reco/kfnew/tools/KfMaterialMapFactory.cxx
similarity index 89%
rename from reco/kfnew/tools/KfMaterialMapCreator.cxx
rename to reco/kfnew/tools/KfMaterialMapFactory.cxx
index b138863da4..49437007ee 100644
--- a/reco/kfnew/tools/KfMaterialMapCreator.cxx
+++ b/reco/kfnew/tools/KfMaterialMapFactory.cxx
@@ -2,13 +2,13 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Sergei Zharko [committer] */
-/// \file KfMaterialMapCreator.cxx
+/// \file KfMaterialMapFactory.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 "KfMaterialMapFactory.h"
#include "KfMaterialMap.h"
#include "TGeoManager.h"
@@ -24,11 +24,11 @@
#include <thread>
using cbm::algo::kf::MaterialMap;
-using kf::tools::MaterialMapCreator;
+using kf::tools::MaterialMapFactory;
// ---------------------------------------------------------------------------------------------------------------------
//
-MaterialMapCreator::MaterialMapCreator(int verbose) : fVerbose(verbose)
+MaterialMapFactory::MaterialMapFactory(int verbose) : fVerbose(verbose)
{
// constructor
@@ -67,7 +67,7 @@ MaterialMapCreator::MaterialMapCreator(int verbose) : fVerbose(verbose)
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 - manually disable voxel finders for problematic volumes in CaToolsMaterialMapFactory.cxx";
msg << '\n';
msg << "\n------------------------------------------------------------------------------------";
LOG_IF(info, fVerbose > 0) << msg.str();
@@ -143,7 +143,7 @@ MaterialMapCreator::MaterialMapCreator(int verbose) : fVerbose(verbose)
// ---------------------------------------------------------------------------------------------------------------------
//
-MaterialMapCreator::~MaterialMapCreator()
+MaterialMapFactory::~MaterialMapFactory()
{
// destructor
@@ -172,7 +172,7 @@ MaterialMapCreator::~MaterialMapCreator()
// ---------------------------------------------------------------------------------------------------------------------
//
-TGeoNavigator* MaterialMapCreator::GetCurrentNavigator(int iThread)
+TGeoNavigator* MaterialMapFactory::GetCurrentNavigator(int iThread)
{
// Get navigator for current thread, create it if it doesnt exist.
// Produce an error when anything goes wrong
@@ -187,12 +187,12 @@ TGeoNavigator* MaterialMapCreator::GetCurrentNavigator(int iThread)
}
if (navi != gGeoManager->GetCurrentNavigator()) {
- LOG(fatal) << "ca::tools::MaterialMapCreator: unexpected behavior of TGeoManager::AddNavigator() !!";
+ LOG(fatal) << "ca::tools::MaterialMapFactory: unexpected behavior of TGeoManager::AddNavigator() !!";
}
}
if (!navi) {
- LOG(fatal) << "ca::tools::MaterialMapCreator: can not find / create TGeoNavigator for thread " << iThread;
+ LOG(fatal) << "ca::tools::MaterialMapFactory: can not find / create TGeoNavigator for thread " << iThread;
}
return navi;
@@ -200,7 +200,7 @@ TGeoNavigator* MaterialMapCreator::GetCurrentNavigator(int iThread)
// ---------------------------------------------------------------------------------------------------------------------
//
-void MaterialMapCreator::CleanUpThreads()
+void MaterialMapFactory::CleanUpThreads()
{
// clean up thread data in TGeoManager
@@ -212,7 +212,7 @@ void MaterialMapCreator::CleanUpThreads()
GetCurrentNavigator(-1);
}
-void MaterialMapCreator::InitThreads()
+void MaterialMapFactory::InitThreads()
{
// (re)initialise the number of threads in TGeoManager
@@ -230,21 +230,22 @@ void MaterialMapCreator::InitThreads()
// ---------------------------------------------------------------------------------------------------------------------
//
-MaterialMap MaterialMapCreator::GenerateMaterialMap(double zRef, double zMin, double zMax, double xyMax, int nBinsDim)
+MaterialMap MaterialMapFactory::GenerateMaterialMap(double zRef, double zMin, double zMax, double xyMax, int nBinsDim)
{
if (fDoRadialProjection) {
- if (fTargetZ >= zRef - 0.1) {
+ if (fTargetZ + 0.05 >= zRef) {
LOG(warn)
- << "kf::tools::MaterialMapCreator: the material reference z-coordinate must be at least 0.1 cm downstream "
+ << "kf::tools::MaterialMapFactory: 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.);
+ zRef = std::max(zRef, fTargetZ + 0.05);
zMin = std::max(fTargetZ, zMin);
+ zMax = std::max(zMax, zRef + 0.05);
}
if (!(zMin <= zRef && zRef <= zMax)) {
- LOG(fatal) << "kf::tools::MaterialMapCreator: the material parameters are inconsistent. It must be: zMin " << zMin
+ LOG(fatal) << "kf::tools::MaterialMapFactory: the material parameters are inconsistent. It must be: zMin " << zMin
<< " <= zRef " << zRef << " <= zMax " << zMax;
}
@@ -263,8 +264,8 @@ MaterialMap MaterialMapCreator::GenerateMaterialMap(double zRef, double zMin, do
nThreadRays[i] = 0;
}
- LOG_IF(info, fVerbose > 1) << "kf::tools::MaterialMapCreator: Generate material map using " << fNthreads
- << " threads..";
+ LOG_IF(info, fVerbose > -1) << "kf::tools::MaterialMapFactory: Generate material map using " << fNthreads
+ << " threads..";
auto naviThread = [&](int iThread) {
// create a navigator for this thread
@@ -306,24 +307,24 @@ MaterialMap MaterialMapCreator::GenerateMaterialMap(double zRef, double zMin, do
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";
+ LOG(warning) << "kf::tools::MaterialMapFactory: 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";
+ LOG(fatal) << "kf::tools::MaterialMapFactory: 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";
+ LOG(fatal) << "kf::tools::MaterialMapFactory: 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! "
+ LOG(fatal) << "kf::tools::MaterialMapFactory: 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! "
@@ -356,7 +357,7 @@ MaterialMap MaterialMapCreator::GenerateMaterialMap(double zRef, double zMin, do
}
else {
LOG(fatal)
- << "kf::tools::MaterialMapCreator: TGeoNavigator propagates the ray upstream. Something is wrong."
+ << "kf::tools::MaterialMapFactory: TGeoNavigator propagates the ray upstream. Something is wrong."
<< " z old " << z << " z new " << zNew << ", dz " << dz;
}
}
@@ -367,7 +368,7 @@ MaterialMap MaterialMapCreator::GenerateMaterialMap(double zRef, double zMin, do
nRays++;
nThreadRays[iThread]++;
if (fVerbose > 2 && nThreadRays[iThread] % 1000000 == 0) {
- LOG(info) << "kf::tools::MaterialMapCreator: report from thread " << iThread << ": material map is "
+ LOG(info) << "kf::tools::MaterialMapFactory: report from thread " << iThread << ": material map is "
<< 100. * nThreadRays[iThread] / nThreadRaysExpected << " \% done";
}
}
@@ -400,7 +401,7 @@ MaterialMap MaterialMapCreator::GenerateMaterialMap(double zRef, double zMin, do
nCrosses += nThreadCrosses[i];
}
- LOG_IF(info, fVerbose > 0) << "kf::tools::MaterialMapCreator: shooted " << nRays << " rays. Each ray crossed "
+ LOG_IF(info, fVerbose > 0) << "kf::tools::MaterialMapFactory: shooted " << nRays << " rays. Each ray crossed "
<< nCrosses * 1. / nRays << " material boundaries in average";
return matBudget;
diff --git a/reco/kfnew/tools/KfMaterialMapCreator.h b/reco/kfnew/tools/KfMaterialMapFactory.h
similarity index 91%
rename from reco/kfnew/tools/KfMaterialMapCreator.h
rename to reco/kfnew/tools/KfMaterialMapFactory.h
index 4814228e60..60de48e78f 100644
--- a/reco/kfnew/tools/KfMaterialMapCreator.h
+++ b/reco/kfnew/tools/KfMaterialMapFactory.h
@@ -2,7 +2,7 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Sergei Zharko [committer] */
-/// \file KfMaterialMapCreator.h
+/// \file KfMaterialMapFactory.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>
@@ -10,6 +10,7 @@
#pragma once
+#include "KfIMaterialMapFactory.h"
#include "Rtypes.h"
#include "TObject.h"
@@ -21,16 +22,16 @@ namespace cbm::algo::kf
namespace kf::tools
{
- /// \class MaterialMapCreator
+ /// \class MaterialMapFactory
/// \brief An utility class to create a material budget map from the TGeo
- class MaterialMapCreator {
+ class MaterialMapFactory : public cbm::algo::kf::IMaterialMapFactory {
public:
/// \brief Constructor from parameters
/// \param verbose Verbosity level (0 - only warning/error output)
- MaterialMapCreator(int verbose = 0);
+ MaterialMapFactory(int verbose = 0);
/// \brief Destructor
- ~MaterialMapCreator();
+ ~MaterialMapFactory();
/// \brief Project rays radially from the targetZ througth the XY-bins at a reference z.
/// \param targetZ z-coordinate of the target [cm]
@@ -65,7 +66,8 @@ namespace kf::tools
///
/// 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);
+ cbm::algo::kf::MaterialMap GenerateMaterialMap(double zRef, double zMin, double zMax, double xyMax,
+ int nBinsDim) override;
/// \brief Enables safe mode of the material initialization
void SetSafeMaterialInitialization(bool val = true) { fDoSafeInitialization = val; }
--
GitLab