KfCore: Introduction of material budget map addition

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

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
 {

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;

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>

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