Ca: clean up the grid class

algo/ca/core/pars/CaInitManager.h

@@ -295,6 +295,7 @@ namespace cbm::algo::ca
     FieldFunction_t fFieldFunction {[](const double (&)[3], double (&)[3]) {}};
     // NOTE: Stations of the detectors which are not assigned as active, are not included in the tracking!
 
+    // TODO: remove
     int fCAIterationsNumberCrosscheck {-1};  ///< Number of iterations to be passed (must be used for cross-checks)
 
     Parameters fParameters {};  ///< CA parameters object

reco/L1/L1Algo/L1Grid.cxx

@@ -30,9 +30,9 @@ template<typename T>
 
 inline void memset(T* dest, T i, size_t num)
 {
-  const size_t tsize   = sizeof(T);
-  unsigned int lastBin = 0;
-  dest[0]              = i;
+  const size_t tsize = sizeof(T);
+  size_t lastBin     = 0;
+  dest[0]            = i;
   while (lastBin + 1 < num) {
     unsigned int l = lastBin + 1;
     l              = lastBin + l + 1 > num ? num - lastBin - 1 : l;
@@ -41,82 +41,36 @@ inline void memset(T* dest, T i, size_t num)
   }
 }
 
-void L1Grid::UpdateIterGrid(unsigned int Nelements, ca::Hit* hits, Vector<ca::HitIndex_t>& indicesBuf,
-                            ca::HitIndex_t* indices, Vector<ca::Hit>& hitsBuf, Vector<ca::GridEntry>& pointsBuf,
-                            ca::GridEntry* points, int& NHitsOnStation, const Vector<unsigned char>& vSFlag)
-{
-  //L1_SHOW(vSFlag.size());
-  fFirstHitInBin.reset(fN + 2, 0);
-
-  for (ca::HitIndex_t x = 0; x < Nelements; x++) {
-
-    const ca::Hit& hit = hits[x];
-
-    if (!(vSFlag[hit.f] || vSFlag[hit.b])) {
-      const ca::HitIndex_t& bin = GetBinBounded(hit.x, hit.y);
-      fHitsInBin[x]             = fFirstHitInBin[bin + 1];
-      fFirstHitInBin[bin + 1]++;
-    }
-  }
-
-  int kk = 2;
-  /* Up-Sweep Phase */
-  for (unsigned int k = 1; k < fN + 2; k = kk) {
-    kk = k << 1;
-    for (unsigned int i = kk - 1; i < fN + 2; i += kk) {
-      fFirstHitInBin[i] = fFirstHitInBin[i - k] + fFirstHitInBin[i];
-    }
-  }
-
-  /* Down-Sweep Phase */
-  for (int k = kk >> 1; k > 1; k = kk) {
-    kk = k >> 1;
-    for (unsigned int i = k - 1; i < fN + 2 - kk; i += k) {
-      fFirstHitInBin[i + kk] = fFirstHitInBin[i] + fFirstHitInBin[i + kk];
-    }
-  }
-
-  for (ca::HitIndex_t x = 0; x < Nelements; x++) {
-
-    const ca::Hit& hit = hits[x];
-    if (!(vSFlag[hit.f] || vSFlag[hit.b])) {
-      const ca::HitIndex_t& bin           = GetBinBounded(hit.x, hit.y);
-      const ca::HitIndex_t& index1        = fHitsInBin[x] + fFirstHitInBin[bin];
-      hitsBuf[index1 + NHitsOnStation]    = hits[x];
-      indicesBuf[index1 + NHitsOnStation] = indices[x];
-      pointsBuf[index1 + NHitsOnStation]  = points[x];
-    }
-  }
-
-  NHitsOnStation += fFirstHitInBin[fN + 1];
-}
-
 
 void L1Grid::AllocateMemory()
 {
   int binsGrid = 600000;
 
-  fFirstHitInBin.reset(binsGrid, 0);
+  fBinFirstEntryIndex.reset(binsGrid, 0);
   fHitsInBin.reset(binsGrid, 0);
 }
 
 
-void L1Grid::BuildBins(float xMin, float xMax, float yMin, float yMax, float sx, float sy)
+void L1Grid::BuildBins(fscal xMin, fscal xMax, fscal yMin, fscal yMax, fscal sx, fscal sy)
 {
-  fStepXInv = 1.f / sx;
-  fStepYInv = 1.f / sy;
+  fMinX = std::min(xMin, xMax);
+  fMinY = std::min(yMin, yMax);
 
-  fXMinOverStep = xMin * fStepXInv;
-  fYMinOverStep = yMin * fStepYInv;
+  xMax = std::max(xMin, xMax);
+  yMax = std::max(yMin, yMax);
 
-  fNx = (xMax * fStepXInv - fXMinOverStep + 1.f);
-  fNy = (yMax * fStepYInv - fYMinOverStep + 1.f);
+  fBinWidthX = sx;
+  fBinWidthY = sy;
 
-  // unsigned int Nold = fN;
+  fBinWidthXinv = 1. / fBinWidthX;
+  fBinWidthYinv = 1. / fBinWidthY;
 
-  fN = fNx * fNy;
+  fNx = static_cast<int>(std::ceil((xMax - fMinX) / fBinWidthX));
+  fNy = static_cast<int>(std::ceil((yMax - fMinY) / fBinWidthY));
+  if (fNx < 1) fNx = 1;
+  if (fNy < 1) fNy = 1;
 
-  fBinInGrid = fN + 1;
+  fN = fNx * fNy;
 }
 
 
@@ -126,30 +80,30 @@ void L1Grid::StoreHits(L1Algo& algo, int iS, ca::HitIndex_t& nGridHitsFilled)
 
   algo.fGridHitStartIndex[iS] = nGridHitsFilled;
 
-  fFirstHitInBin.reset(fN + 2, 0);
-
+  fBinFirstEntryIndex.reset(fN + 2, 0);
+  fHitsInBin.reset(nHits, 0);
   for (ca::HitIndex_t ih = 0; ih < nHits; ih++) {
     ca::HitIndex_t caHitId = algo.fSliceHitIds[iS][ih];
     const ca::Hit& hit     = algo.GetInputData().GetHit(caHitId);
     auto bin               = GetBinBounded(hit.x, hit.y);
-    fHitsInBin[ih]         = fFirstHitInBin[bin + 1];
-    fFirstHitInBin[bin + 1]++;
+    fHitsInBin[ih]         = fBinFirstEntryIndex[bin + 1];
+    fBinFirstEntryIndex[bin + 1]++;
   }
 
   int kk = 2;
-  /* Up-Sweep Phase */
-  for (unsigned int k = 1; k < fN + 2; k = kk) {
+  // Up-Sweep Phase
+  for (int k = 1; k < fN + 2; k = kk) {
     kk = k << 1;
-    for (unsigned int i = kk - 1; i < fN + 2; i += kk) {
-      fFirstHitInBin[i] = fFirstHitInBin[i - k] + fFirstHitInBin[i];
+    for (int i = kk - 1; i < fN + 2; i += kk) {
+      fBinFirstEntryIndex[i] = fBinFirstEntryIndex[i - k] + fBinFirstEntryIndex[i];
     }
   }
 
-  /* Down-Sweep Phase */
-  for (unsigned int k = kk >> 1; k > 1; k = kk) {
+  // Down-Sweep Phase
+  for (int k = kk >> 1; k > 1; k = kk) {
     kk = k >> 1;
-    for (unsigned int i = k - 1; i < fN + 2 - kk; i += k) {
-      fFirstHitInBin[i + kk] = fFirstHitInBin[i] + fFirstHitInBin[i + kk];
+    for (int i = k - 1; i < fN + 2 - kk; i += k) {
+      fBinFirstEntryIndex[i + kk] = fBinFirstEntryIndex[i] + fBinFirstEntryIndex[i + kk];
     }
   }
 
@@ -158,7 +112,7 @@ void L1Grid::StoreHits(L1Algo& algo, int iS, ca::HitIndex_t& nGridHitsFilled)
     const ca::Hit& hit     = algo.GetInputData().GetHit(caHitId);
     auto bin               = GetBinBounded(hit.x, hit.y);
 
-    const ca::HitIndex_t& index1               = fFirstHitInBin[bin] + fHitsInBin[ih];
+    const ca::HitIndex_t& index1               = fBinFirstEntryIndex[bin] + fHitsInBin[ih];
     algo.fGridHits[nGridHitsFilled + index1]   = hit;
     algo.fGridHitIds[nGridHitsFilled + index1] = caHitId;
   }
@@ -166,3 +120,54 @@ void L1Grid::StoreHits(L1Algo& algo, int iS, ca::HitIndex_t& nGridHitsFilled)
   nGridHitsFilled += nHits;
   algo.fGridHitStopIndex[iS] = nGridHitsFilled;
 }
+
+
+void L1Grid::UpdateIterGrid(ca::HitIndex_t Nelements, ca::Hit* hits, Vector<ca::HitIndex_t>& indicesBuf,
+                            ca::HitIndex_t* indices, Vector<ca::Hit>& hitsBuf, Vector<ca::GridEntry>& pointsBuf,
+                            ca::GridEntry* points, int& NHitsOnStation, const Vector<unsigned char>& vSFlag)
+{
+  //L1_SHOW(vSFlag.size());
+  fBinFirstEntryIndex.reset(fN + 2, 0);
+  fHitsInBin.reset(Nelements, 0);
+  for (ca::HitIndex_t x = 0; x < Nelements; x++) {
+
+    const ca::Hit& hit = hits[x];
+
+    if (!(vSFlag[hit.f] || vSFlag[hit.b])) {
+      const ca::HitIndex_t& bin = GetBinBounded(hit.x, hit.y);
+      fHitsInBin[x]             = fBinFirstEntryIndex[bin + 1];
+      fBinFirstEntryIndex[bin + 1]++;
+    }
+  }
+
+  int kk = 2;
+  /* Up-Sweep Phase */
+  for (int k = 1; k < fN + 2; k = kk) {
+    kk = k << 1;
+    for (int i = kk - 1; i < fN + 2; i += kk) {
+      fBinFirstEntryIndex[i] = fBinFirstEntryIndex[i - k] + fBinFirstEntryIndex[i];
+    }
+  }
+
+  /* Down-Sweep Phase */
+  for (int k = kk >> 1; k > 1; k = kk) {
+    kk = k >> 1;
+    for (int i = k - 1; i < fN + 2 - kk; i += k) {
+      fBinFirstEntryIndex[i + kk] = fBinFirstEntryIndex[i] + fBinFirstEntryIndex[i + kk];
+    }
+  }
+
+  for (ca::HitIndex_t x = 0; x < Nelements; x++) {
+
+    const ca::Hit& hit = hits[x];
+    if (!(vSFlag[hit.f] || vSFlag[hit.b])) {
+      const ca::HitIndex_t& bin           = GetBinBounded(hit.x, hit.y);
+      const ca::HitIndex_t& index1        = fHitsInBin[x] + fBinFirstEntryIndex[bin];
+      hitsBuf[index1 + NHitsOnStation]    = hits[x];
+      indicesBuf[index1 + NHitsOnStation] = indices[x];
+      pointsBuf[index1 + NHitsOnStation]  = points[x];
+    }
+  }
+
+  NHitsOnStation += fBinFirstEntryIndex[fN + 1];
+}

reco/L1/L1Algo/L1Grid.h

@@ -13,6 +13,7 @@
 
 #include "CaGridEntry.h"
 #include "CaHit.h"
+#include "CaSimd.h"
 #include "CaVector.h"
 
 class L1Algo;
@@ -28,112 +29,116 @@ namespace
 /// \class L1Grid
 ///
 /// \brief 2-dimensional grid of pointers.
+///
 /// pointers to (x,y)-like objects are assigned to the corresponding grid bin
 /// used by L1Tracker to speed-up the hit operations
 /// grid axis are named X,Y
 ///
 class L1Grid {
 public:
+  /// Constructor
   L1Grid() = default;
 
+  /// Constructor
   L1Grid(const L1Grid& grid) : fN(grid.N()), fNx(grid.Nx()), fNy(grid.Ny()) {}
 
+  /// Destructor
   ~L1Grid() = default;
 
+  void AllocateMemory();
+
+  void BuildBins(fscal xMin, fscal xMax, fscal yMin, fscal yMax, fscal sx, fscal sy);
+
   void StoreHits(L1Algo& Algo, int iS, ca::HitIndex_t& nGridHitsFilled);
 
-  void BuildBins(float xMin, float xMax, float yMin, float yMax, float sx, float sy);
+  void UpdateIterGrid(ca::HitIndex_t Nelements, ca::Hit* hits, Vector<ca::HitIndex_t>& indicesBuf,
+                      ca::HitIndex_t* indices, Vector<ca::Hit>& hitsBuf, Vector<ca::GridEntry>& entriesBuf,
+                      ca::GridEntry* entries, int& NHitsOnStation, const Vector<unsigned char>& vSFlag);
 
-  void AllocateMemory();
-  void Create(float xMin, float xMax, float yMin, float yMax, float sx, float sy);
+  /// get bin index
+  int GetBin(fscal X, fscal Y) const;
 
-  void Fill(const ca::GridEntry* entries, ca::HitIndex_t n);  // call after sort
-  void FillPar(const ca::GridEntry* entries, ca::HitIndex_t n);
+  /// get bin index with boundary check
+  int GetBinBounded(fscal X, fscal Y) const;
 
+  /// get bin X,Y indices with boundary check
+  void GetBinBounded(fscal X, fscal Y, int& bX, int& bY) const;
 
-  int GetBin(float X, float Y) const;
+  /// get bin bounds
+  void GetBinBounds(int iBin, fscal& Xmin, fscal& Xmax, fscal& Ymin, fscal& Ymax) const;
 
-  unsigned int GetBinBounded(const float& X, const float& Y) const;
-  void GetBinBounded(const float& X, const float& Y, unsigned short& bX, unsigned short& bY) const;
+  /// get number of bins
+  int N() const { return fN; }
 
-  void GetBinBounds(unsigned int iBin, float& Xmin, float& Xmax, float& Ymin, float& Ymax) const;
+  /// get number of bins in X
+  int Nx() const { return fNx; }
 
-  unsigned int N() const { return fN; }
-  unsigned short Nx() const { return fNx; }
-  unsigned short Ny() const { return fNy; }
+  /// get number of bins in Y
+  int Ny() const { return fNy; }
 
-  ca::HitIndex_t FirstHitInBin(unsigned int i) const
+  /// get index of the first bin entry in fHitsInBin array
+  ca::HitIndex_t FirstHitInBin(int i) const
   {
-    if (i < (fN + 1)) return fFirstHitInBin[i];
+    if (i < (fN + 1)) return fBinFirstEntryIndex[i];
     else
-      return fFirstHitInBin[fN + 1];
+      return fBinFirstEntryIndex[fN + 1];
   }
 
-
-  void UpdateIterGrid(unsigned int Nelements, ca::Hit* hits, Vector<ca::HitIndex_t>& indicesBuf,
-                      ca::HitIndex_t* indices, Vector<ca::Hit>& hitsBuf, Vector<ca::GridEntry>& entriesBuf,
-                      ca::GridEntry* entries, int& NHitsOnStation, const Vector<unsigned char>& vSFlag);
-
-
 private:
-  unsigned int fN {0};        //* total N bins
-  unsigned short fNx {0};     //* N bins in X
-  unsigned short fNy {0};     //* N bins in Y
-  float fXMinOverStep {0.f};  //* minimal X value * fStepXInv
-  float fYMinOverStep {0.f};  //* minimal Y value * fStepYInv
-  float fStepXInv {0.f};      //* inverse bin size in X
-  float fStepYInv {0.f};      //* inverse bin size in Y
-  int fBinInGrid {0};
-
-  Vector<ca::HitIndex_t> fFirstHitInBin {"L1Grid::fFirstHitInBin"};
+  int fN {0};                ///< total N bins
+  int fNx {0};               ///< N bins in X
+  int fNy {0};               ///< N bins in Y
+  fscal fMinX {0.};          ///< minimal X value
+  fscal fMinY {0.};          ///< minimal Y value
+  fscal fBinWidthX {0.};     ///< bin width in X
+  fscal fBinWidthY {0.};     ///< bin width in Y
+  fscal fBinWidthXinv {0.};  ///< inverse bin width in X
+  fscal fBinWidthYinv {0.};  ///< inverse bin width in Y
+
+  Vector<ca::HitIndex_t> fBinFirstEntryIndex {"L1Grid::fBinFirstEntryIndex"};
   Vector<ca::HitIndex_t> fHitsInBin {"L1Grid::fHitsInBin"};
-
-  // vector <omp_lock_t> lock;
 };
 
 
-inline int L1Grid::GetBin(float X, float Y) const
+inline int L1Grid::GetBin(fscal X, fscal Y) const
 {
-  //* get the bin pointer
-  const int& xBin = static_cast<int>(X * fStepXInv - fXMinOverStep);
-  const int& yBin = static_cast<int>(Y * fStepYInv - fYMinOverStep);
+  /// get bin index
+  int xBin = static_cast<int>((X - fMinX) * fBinWidthXinv);
+  int yBin = static_cast<int>((Y - fMinY) * fBinWidthYinv);
   assert(xBin >= 0);
   assert(yBin >= 0);
-  assert(xBin < static_cast<int>(fNx));
-  assert(yBin < static_cast<int>(fNy));
-  const int& bin = yBin * fNx + xBin;
-  return bin;
+  assert(xBin < fNx);
+  assert(yBin < fNy);
+  return yBin * fNx + xBin;
 }
 
 
-inline void L1Grid::GetBinBounds(unsigned int iBin, float& Xmin, float& Xmax, float& Ymin, float& Ymax) const
+inline void L1Grid::GetBinBounds(int iBin, fscal& Xmin, fscal& Xmax, fscal& Ymin, fscal& Ymax) const
 {
   int yBin = iBin / fNx;
   int xBin = iBin % fNx;
-  Xmin     = (fXMinOverStep + xBin) / fStepXInv;
-  Ymin     = (fYMinOverStep + yBin) / fStepYInv;
-  Xmax     = Xmin + 1. / fStepXInv;
-  Ymax     = Ymin + 1. / fStepYInv;
+  Xmin     = fMinX + xBin * fBinWidthX;
+  Ymin     = fMinY + yBin * fBinWidthY;
+  Xmax     = Xmin + fBinWidthX;
+  Ymax     = Ymin + fBinWidthY;
 }
 
 
-inline unsigned int L1Grid::GetBinBounded(const float& X, const float& Y) const
+inline int L1Grid::GetBinBounded(fscal X, fscal Y) const
 {
-  //* get the bin pointer at
-
-  unsigned short xBin, yBin;
+  /// get bin index
+  int xBin, yBin;
   GetBinBounded(X, Y, xBin, yBin);
-  return (unsigned int) yBin * (unsigned int) fNx + (unsigned int) xBin;
+  return yBin * fNx + xBin;
 }
 
 
-inline void L1Grid::GetBinBounded(const float& X, const float& Y, unsigned short& bX, unsigned short& bY) const
+inline void L1Grid::GetBinBounded(fscal X, fscal Y, int& bX, int& bY) const
 {
-  const short& xBin = (X * fStepXInv - fXMinOverStep);
-  const short& yBin = (Y * fStepYInv - fYMinOverStep);
-
-  bX = std::max(short(0), std::min(short(fNx - 1), xBin));
-  bY = std::max(short(0), std::min(short(fNy - 1), yBin));
+  int xBin = (X - fMinX) * fBinWidthXinv;
+  int yBin = (Y - fMinY) * fBinWidthYinv;
+  bX       = std::max(0, std::min(fNx - 1, xBin));
+  bY       = std::max(0, std::min(fNy - 1, yBin));
 }
 
 

reco/L1/L1Algo/L1HitArea.h

@@ -28,10 +28,10 @@ public:
 protected:
   const L1Grid& fGrid;
 
-  unsigned short fBYmax;  // maximal Y bin index
-  unsigned short fBDX;    // X distance of bin indexes
-  unsigned int fIndXmin;  // minimum index for
-  unsigned short fIy;     // current Y bin index (incremented while iterating)
+  int fBYmax;               // maximal Y bin index
+  int fBDX;                 // X distance of bin indexes
+  int fIndXmin;             // minimum index for
+  int fIy;                  // current Y bin index (incremented while iterating)
   ca::HitIndex_t fHitXlst;  // last possible hit index in current y-line
   ca::HitIndex_t fIh;       // hit index iterating inside the bins
   int fNx;                // Number of bins in X direction
@@ -51,7 +51,7 @@ inline L1HitArea::L1HitArea(const L1Grid& grid, float x, float y, float dx, floa
   const float maxX = x + dx;
   const float minY = y - dy;
   const float maxY = y + dy;
-  unsigned short bXmin, bYmin, bXmax;  // boundary bin indexes
+  int bXmin, bYmin, bXmax;  // boundary bin indexes
   fGrid.GetBinBounded(minX, minY, bXmin, bYmin);
   fGrid.GetBinBounded(maxX, maxY, bXmax, fBYmax);