diff --git a/algo/ca/core/pars/CaInitManager.h b/algo/ca/core/pars/CaInitManager.h
index 826f0a353220b8a49ad3203ace27426778c4ce28..1b836d09fa0e48cf7edf8e6e9e8259cc64ab4f57 100644
--- a/algo/ca/core/pars/CaInitManager.h
+++ b/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
diff --git a/reco/L1/L1Algo/L1Grid.cxx b/reco/L1/L1Algo/L1Grid.cxx
index 85ff8dbe4b8489e26765ff5f0350ed7b7fd4aaf3..a07ed395ef40055d1d87741b9591b9a7f1077c5d 100644
--- a/reco/L1/L1Algo/L1Grid.cxx
+++ b/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];
+}
diff --git a/reco/L1/L1Algo/L1Grid.h b/reco/L1/L1Algo/L1Grid.h
index de031048555d388c4f742aff68a61317d797a962..f117a93b6c3a326711f6efc16ec374573d11cd9f 100644
--- a/reco/L1/L1Algo/L1Grid.h
+++ b/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));
}
diff --git a/reco/L1/L1Algo/L1HitArea.h b/reco/L1/L1Algo/L1HitArea.h
index bddd6868128c509e9055b4ca7c97e1fd70c0077c..206fe60ae1536fe3a156933120cf2e452892d910 100644
--- a/reco/L1/L1Algo/L1HitArea.h
+++ b/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);