From 3cf3f5240031e533268f4bc6c54891cefdc4cfb7 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Wed, 17 Jul 2024 10:52:55 +0200
Subject: [PATCH] KfCore: Introduction of material budget map addition
---
algo/ca/core/pars/CaConstants.h | 28 -----------
algo/kf/core/pars/KfMaterialMap.cxx | 75 ++++++++++++++++++++++-------
algo/kf/core/pars/KfMaterialMap.h | 17 ++++---
algo/kf/core/pars/KfSetup.h | 27 +++++++++++
algo/kf/core/utils/KfSimdVc.h | 14 +++---
5 files changed, 102 insertions(+), 59 deletions(-)
diff --git a/algo/ca/core/pars/CaConstants.h b/algo/ca/core/pars/CaConstants.h
index c858914bc1..dba34ef9e1 100644
--- a/algo/ca/core/pars/CaConstants.h
+++ b/algo/ca/core/pars/CaConstants.h
@@ -226,32 +226,4 @@ namespace cbm::algo::ca::constants
constexpr char GYbr[] = "\e[1;7;37m"; ///< bold-reverse grey
constexpr char WTbr[] = "\e[1;7;38m"; ///< bold-reverse white
} // namespace clrs
-
- // clang-format: off
- constexpr char Logo[] = "\n\n\n\n"
- " [][][] [][][] [][][][][] [][] [] "
- " \n"
- " [] [] [] [] [] [] [] [] [] "
- " \n"
- " [] [] [] [] [] "
- " \n"
- " [] [] [] [] [][][][] [][][] [][][] [] [][] [][] "
- "[][][][] [][][][]\n"
- " [] [] [] [] [] [] [] [] [] [] [] [] [] [] "
- " [] [] [] [] \n"
- " [] [] [] [] [] [][][] [] [] [] [] "
- " [] [] [] [] \n"
- " [] [] [] [] [] [] [] [] [] [] [] [] [] "
- " [] [] [] [] \n"
- " [][][] [][] [][] [][] [][] [][][][] [][][] [][] [][] [][] "
- "[][] [][] [][] [] \n"
- " "
- " [] \n"
- " By CBM Collaboration "
- " [] [] \n"
- " "
- " [][][]\n\n\n";
- // clang-format: on
-
-
} // namespace cbm::algo::ca::constants
diff --git a/algo/kf/core/pars/KfMaterialMap.cxx b/algo/kf/core/pars/KfMaterialMap.cxx
index 809e22753e..3490bd3d38 100644
--- a/algo/kf/core/pars/KfMaterialMap.cxx
+++ b/algo/kf/core/pars/KfMaterialMap.cxx
@@ -3,20 +3,23 @@
Authors: Sergey Gorbunov, Sergei Zharko [committer] */
#include "KfMaterialMap.h"
+#include "KfSimd.h"
#include "AlgoFairloggerCompat.h"
#include <iomanip>
#include <sstream>
#include <vector>
+#include <cmath>
using cbm::algo::kf::MaterialMap;
+using cbm::algo::kf::fvec;
-//------------------------------------------------------------------------------------------------------------------------------------
+// ---------------------------------------------------------------------------------------------------------------------
//
MaterialMap::MaterialMap(MaterialMap&& other) noexcept { this->Swap(other); }
-//------------------------------------------------------------------------------------------------------------------------------------
+// ---------------------------------------------------------------------------------------------------------------------
//
MaterialMap& MaterialMap::operator=(MaterialMap&& other) noexcept
{
@@ -27,7 +30,43 @@ MaterialMap& MaterialMap::operator=(MaterialMap&& other) noexcept
return *this;
}
-//------------------------------------------------------------------------------------------------------------------------------------
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void MaterialMap::Add(const MaterialMap& other, float zTarg)
+{
+ // The function allows to add a material layer either from the left or from the right to the station
+ // NOTE: A symmetry of x-y is assumed
+ constexpr int nRays = 3; // Number of rays in a bin in a dimension
+ const bool bRadialRays = !std::isnan(zTarg); // Are rays radial (true, coming from target) or parallel (false)?
+ const auto scaleFactor = bRadialRays ? ((other.fZref - zTarg) / (this->fZref - zTarg)) : 1.;
+ const auto binSize = 2. * scaleFactor * this->fXYmax / this->fNbins; // Size of each bin dimension [cm]
+ const auto stepSize = binSize / static_cast<float>(nRays); // Step between two neighboring rays [cm]
+
+ // The coordinates of the first ray intersection with the other material layer [cm]
+ float xBinOther = -this->fXYmax * scaleFactor + stepSize * 0.5;
+ float yBinOther = xBinOther;
+
+ // Loop over bins of the active (this)
+ for (int iBinY = 0; iBinY < this->fNbins; ++iBinY) {
+ for (int iBinX = 0; iBinX < this->fNbins; ++iBinX) {
+ // Collect material using ray shooting
+ float avgThickness = 0; // Collected average thickness
+ for (int iRayY = 0; iRayY < nRays; ++iRayY) {
+ for (int iRayX = 0; iRayX < nRays; ++iRayX) {
+ avgThickness += other.GetRadThickScal(xBinOther + iBinX * stepSize, yBinOther + iBinY * stepSize);
+ }
+ }
+ this->fTable[iBinX + this->fNbins * iBinY] += avgThickness / (nRays * nRays);
+ xBinOther += binSize;
+ }
+ yBinOther += binSize;
+ }
+ this->fZmin = std::min(this->fZmin, other.fZmin);
+ this->fZmax = std::max(this->fZmax, other.fZmax);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
int MaterialMap::GetBin(float x, float y) const
{
int i = static_cast<int>((x + fXYmax) * fFactor);
@@ -38,7 +77,7 @@ int MaterialMap::GetBin(float x, float y) const
return i + j * fNbins;
}
-//------------------------------------------------------------------------------------------------------------------------------------
+// ---------------------------------------------------------------------------------------------------------------------
//
float MaterialMap::GetRadThickScal(float x, float y) const
{
@@ -51,17 +90,17 @@ float MaterialMap::GetRadThickScal(float x, float y) const
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;
-//}
-
-//------------------------------------------------------------------------------------------------------------------------------------
+// ---------------------------------------------------------------------------------------------------------------------
+//
+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
{
@@ -85,7 +124,7 @@ void MaterialMap::CheckConsistency() const
}
}
-//------------------------------------------------------------------------------------------------------------------------------------
+// ---------------------------------------------------------------------------------------------------------------------
//
void MaterialMap::Initialize(int nBins, float xyMax, float zRef, float zMin, float zMax)
{
@@ -118,7 +157,7 @@ void MaterialMap::Initialize(int nBins, float xyMax, float zRef, float zMin, flo
fTable.resize(fNbins * fNbins);
}
-//------------------------------------------------------------------------------------------------------------------------------------
+// ---------------------------------------------------------------------------------------------------------------------
//
void MaterialMap::Swap(MaterialMap& other) noexcept
{
@@ -131,7 +170,7 @@ void MaterialMap::Swap(MaterialMap& other) noexcept
std::swap(fTable, other.fTable); // Probably can cause segmentation violation (did not understand)
}
-//------------------------------------------------------------------------------------------------------------------------------------
+// ---------------------------------------------------------------------------------------------------------------------
//
std::string MaterialMap::ToString() const
{
diff --git a/algo/kf/core/pars/KfMaterialMap.h b/algo/kf/core/pars/KfMaterialMap.h
index b97fde1f90..e8ea76f508 100644
--- a/algo/kf/core/pars/KfMaterialMap.h
+++ b/algo/kf/core/pars/KfMaterialMap.h
@@ -41,6 +41,11 @@ namespace cbm::algo::kf
/// \brief Destructor
~MaterialMap() noexcept = default;
+ /// \brief Adds material layer
+ /// \param other Other material layer
+ /// \param zTarg z-coordinate of the target
+ void Add(const MaterialMap& other, float zTarg = defs::Undef<float>);
+
/// \brief Gets number of bins (rows or columns) of the material table
int GetNbins() const { return fNbins; }
@@ -75,14 +80,14 @@ namespace cbm::algo::kf
/// \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;
+ 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);
- //}
+ bool IsNaN() const
+ {
+ return defs::IsUndefined(fNbins) || defs::IsUndefined(fXYmax) || defs::IsUndefined(fFactor)
+ || defs::IsUndefined(fZref) || defs::IsUndefined(fZmin) || defs::IsUndefined(fZmax);
+ }
/// \brief Verifies class invariant consistency
void CheckConsistency() const;
diff --git a/algo/kf/core/pars/KfSetup.h b/algo/kf/core/pars/KfSetup.h
index dcd211771e..1a5e57f598 100644
--- a/algo/kf/core/pars/KfSetup.h
+++ b/algo/kf/core/pars/KfSetup.h
@@ -22,6 +22,33 @@ namespace cbm::algo::kf
{
using MaterialContainer_t = Vector<MaterialMap>;
+ /// \class CollisionPoint
+ /// \brief KF-framework representation of the nominal collision point region
+ /// \tparam DataT Underlying data type
+ template<typename DataT>
+ class alignas(defs::MemAlign) CollisionPoint {
+ /// \brief Default constructor
+ CollisionPoint() = default;
+
+ /// \brief Copy constructor
+ CollisionPoint(const CollisionPoint&) = default;
+
+ /// \brief Move constructor
+ CollisionPoint(CollisionPoint&&) = default;
+
+ /// \brief Destructor
+ ~CollisionPoint() = default;
+
+ /// \brief Copy assignment operator
+ CollisionPoint& operator=(const CollisionPoint&) = default;
+
+ /// \brief Move assignment operator
+ CollisionPoint& operator=(CollisionPoint&&) = default;
+
+
+ };
+
+ /// \class Setup
/// \brief KF-framework representation of the detector setup
/// \tparam DataT Underlying data type
template<typename DataT>
diff --git a/algo/kf/core/utils/KfSimdVc.h b/algo/kf/core/utils/KfSimdVc.h
index 6162119184..20394e3fa4 100644
--- a/algo/kf/core/utils/KfSimdVc.h
+++ b/algo/kf/core/utils/KfSimdVc.h
@@ -12,13 +12,13 @@
namespace cbm::algo::kf
{
- typedef Vc::float_v fvec;
- //typedef Vc::double_v fvec;
- //typedef Vc::Vector<float, Vc::VectorAbi::Scalar> fvec;
- //typedef Vc::SimdArray<float, 4> fvec;
-
- typedef fvec::EntryType fscal;
- typedef fvec::MaskType fmask;
+ using fvec = Vc::float_v;
+ //using fvec = Vc::double_v;
+ //using fvec = Vc::Vector<float, Vc::VectorAbi::Scalar>;
+ //using fvec = Vc::SimdArray<float, 4>;
+
+ using fscal = fvec::EntryType;
+ using fmask = fvec::MaskType;
#define _fvecalignment __attribute__((aligned(Vc::VectorAlignment)))
} // namespace cbm::algo::kf
--
GitLab