diff --git a/algo/ca/core/CMakeLists.txt b/algo/ca/core/CMakeLists.txt index 755fcaa22227efee64b9d565c9d460bde39574b5..92b404c128f5e14cafd3c4b3617c4fce6d27d176 100644 --- a/algo/ca/core/CMakeLists.txt +++ b/algo/ca/core/CMakeLists.txt @@ -10,6 +10,7 @@ set(SRCS ${CMAKE_CURRENT_SOURCE_DIR}/data/CaInputData.cxx ${CMAKE_CURRENT_SOURCE_DIR}/data/CaTrack.cxx ${CMAKE_CURRENT_SOURCE_DIR}/data/CaTrackParam.cxx + ${CMAKE_CURRENT_SOURCE_DIR}/data/CaGrid.cxx ${CMAKE_CURRENT_SOURCE_DIR}/pars/CaConfigReader.cxx ${CMAKE_CURRENT_SOURCE_DIR}/pars/CaField.cxx ${CMAKE_CURRENT_SOURCE_DIR}/pars/CaInitManager.cxx @@ -54,6 +55,9 @@ install( data/CaTrackParam.h data/CaTrack.h pars/CaConfigReader.h + data/CaGridEntry.h + data/CaGrid.h + data/CaGridArea.h pars/CaConstants.h pars/CaField.h pars/CaInitManager.h diff --git a/algo/ca/core/data/CaGrid.cxx b/algo/ca/core/data/CaGrid.cxx new file mode 100644 index 0000000000000000000000000000000000000000..4f06c5ee6d20c12fc4bd5a6337466dc27575b125 --- /dev/null +++ b/algo/ca/core/data/CaGrid.cxx @@ -0,0 +1,121 @@ +/* Copyright (C) 2017-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Maksym Zyzak [committer], Valentina Akishina */ + +/// \file L1Grid.cxx +/// \brief Class for storing 2d objects in a grid + + +#include "CaGrid.h" + +#include <algorithm> + +#include <string.h> + +#include "CaHit.h" + + +using namespace cbm::algo::ca; +using namespace cbm::algo; + + +void Grid::BuildBins(fscal xMin, fscal xMax, fscal yMin, fscal yMax, fscal binWidthX, fscal binWidthY) +{ + fMinX = std::min(xMin, xMax); + fMinY = std::min(yMin, yMax); + + xMax = std::max(xMin, xMax); + yMax = std::max(yMin, yMax); + + fBinWidthX = binWidthX; + fBinWidthY = binWidthY; + + // some sanity checks + if (fBinWidthX < 0.001) { fBinWidthX = 0.001; } + if (fBinWidthY < 0.001) { fBinWidthY = 0.001; } + + fBinWidthXinv = 1. / fBinWidthX; + fBinWidthYinv = 1. / fBinWidthY; + + fNx = static_cast<int>(std::ceil((xMax - fMinX) / fBinWidthX)); + fNy = static_cast<int>(std::ceil((yMax - fMinY) / fBinWidthY)); + + // some sanity checks + if (fNx < 1) fNx = 1; + if (fNy < 1) fNy = 1; + + fN = fNx * fNy; + + fEntries.clear(); + fFirstBinEntryIndex.reset(fN + 1, 0); + fNofBinEntries.reset(fN + 1, 0); +} + + +void Grid::StoreHits(const Vector<ca::Hit>& hits, ca::HitIndex_t hitStartIndex, ca::HitIndex_t nHits, + const Vector<unsigned char>& hitKeyFlags) +{ + fFirstBinEntryIndex.reset(fN + 1, 0); + fNofBinEntries.reset(fN + 1, 0); + + int nEntries = 0; + for (ca::HitIndex_t ih = 0; ih < nHits; ih++) { + const ca::Hit& hit = hits[hitStartIndex + ih]; + if (!(hitKeyFlags[hit.f] || hitKeyFlags[hit.b])) { + fNofBinEntries[GetBin(hit.x, hit.y)]++; + nEntries++; + } + } + + fEntries.reset(nEntries); + + for (int bin = 0; bin < fN; bin++) { + fFirstBinEntryIndex[bin + 1] = fFirstBinEntryIndex[bin] + fNofBinEntries[bin]; + fNofBinEntries[bin] = 0; + } + fNofBinEntries[fN] = 0; + + fMaxRangeX = 0.; + fMaxRangeY = 0.; + fMaxRangeT = 0.; + + for (ca::HitIndex_t ih = 0; ih < nHits; ih++) { + const ca::Hit& hit = hits[hitStartIndex + ih]; + if (!(hitKeyFlags[hit.f] || hitKeyFlags[hit.b])) { + int bin = GetBin(hit.x, hit.y); + fEntries[fFirstBinEntryIndex[bin] + fNofBinEntries[bin]].Set(hit, hitStartIndex + ih); + fNofBinEntries[bin]++; + fMaxRangeX = std::max(fMaxRangeX, hit.rangeX); + fMaxRangeY = std::max(fMaxRangeY, hit.rangeY); + fMaxRangeT = std::max(fMaxRangeT, hit.rangeT); + } + } +} + +void Grid::RemoveUsedHits(const Vector<ca::Hit>& hits, const Vector<unsigned char>& hitKeyFlags) +{ + int nEntries = 0; + fMaxRangeX = 0.; + fMaxRangeY = 0.; + fMaxRangeT = 0.; + + for (int bin = 0; bin < fN; bin++) { + ca::HitIndex_t firstEntryOld = fFirstBinEntryIndex[bin]; + fFirstBinEntryIndex[bin] = nEntries; + fNofBinEntries[bin] = 0; + for (ca::HitIndex_t i = firstEntryOld; i < fFirstBinEntryIndex[bin + 1]; i++) { + const ca::Hit& hit = hits[fEntries[i].GetObjectId()]; + if (!(hitKeyFlags[hit.f] || hitKeyFlags[hit.b])) { + fEntries[nEntries] = fEntries[i]; + nEntries++; + fNofBinEntries[bin]++; + fMaxRangeX = std::max(fMaxRangeX, hit.rangeX); + fMaxRangeY = std::max(fMaxRangeY, hit.rangeY); + fMaxRangeT = std::max(fMaxRangeT, hit.rangeT); + } + } + } + fFirstBinEntryIndex[fN] = nEntries; + fNofBinEntries[fN] = 0; + fEntries.reduce(nEntries); +} diff --git a/algo/ca/core/data/CaGrid.h b/algo/ca/core/data/CaGrid.h new file mode 100644 index 0000000000000000000000000000000000000000..01f193c1da9a5e667f1996fe6def8f31816b16a3 --- /dev/null +++ b/algo/ca/core/data/CaGrid.h @@ -0,0 +1,196 @@ +/* Copyright (C) 2010-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Igor Kulakov [committer], Sergey Gorbunov */ + +/// \file CaGrid.h +/// \brief A class to store hit information in a backet-sorted way on 2D grid +/// +/// This code is based on the code of the ALICE HLT Project +/// + +#ifndef CaGrid_h +#define CaGrid_h 1 + +#include "CaGridEntry.h" +#include "CaHit.h" +#include "CaSimd.h" +#include "CaVector.h" + +namespace cbm::algo::ca +{ + + namespace + { + using namespace cbm::algo; + } + + /// \class Grid + /// \brief Class for storing 2d objects in a grid + + /// It creates 2-dimensional backet-sorted grid of pointers to 2-dimensional objects. + /// + /// The class provides an access to the objects in selected 2D bin without touching the rest of the data. + /// To loop over the objects in arbitrary XY-area one can use the ca::GridArea class. + /// + /// The class is used by CaTracker to speed-up the hit search operations + /// The grid axis are named X and Y + /// + class Grid { + public: + /// \brief Default constructor + Grid() = default; + + /// \brief Destructor + ~Grid() = default; + + /// \brief Build the grid + /// \param xMin - minimal X value + /// \param xMax - maximal X value + /// \param yMin - minimal Y value + /// \param yMax - maximal Y value + /// \param binWidthX - bin width in X + /// \param binWidthY - bin width in Y + void BuildBins(fscal xMin, fscal xMax, fscal yMin, fscal yMax, fscal binWidthX, fscal binWidthY); + + /// \brief Store objects in the grid + // TODO: write the method with a template input + /// void StoreData(...); + + /// \brief Store ca::Hits in the grid + /// \param hits - vector of hits to store + /// \param hitStartIndex - index of the first hit to store + /// \param nHits - number of hits to store starting from the hitStartIndex + /// \param hitKeyFlags - vector of flags to recognise used hits and skip them + void StoreHits(const ca::Vector<ca::Hit>& hits, ca::HitIndex_t hitStartIndex, ca::HitIndex_t nHits, + const ca::Vector<unsigned char>& hitKeyFlags); + + /// \brief Remove grid entries that correspond to the used hits + void RemoveUsedHits(const ca::Vector<ca::Hit>& hits, const ca::Vector<unsigned char>& hitKeyFlags); + + /// \brief Get bin index for (X,Y) point with boundary check + /// \param X - point x coordinate + /// \param Y - point y coordinate + /// \return bin index in the range [0, fN-1] + int GetBin(fscal X, fscal Y) const; + + /// \brief Get bin X index with boundary check + /// \param X - x coordinate + /// \return binX index + int GetBinX(fscal X) const; + + /// \brief Get bin Y index with boundary check + /// \param Y - y coordinate + /// \return binY index + int GetBinY(fscal Y) const; + + /// \brief Get bin bounds along X + /// \param iBin - bin index + /// \return pair of (Xmin, Xmax) bounds + std::tuple<fscal, fscal> GetBinBoundsX(int iBin) const; + + /// \brief Get bin bounds along Y + /// \param iBin - bin index + /// \return pair of (Ymin, Ymax) bounds + std::tuple<fscal, fscal> GetBinBoundsY(int iBin) const; + + /// Get number of bins + int GetNofBins() const { return fN; } + + /// Get number of bins along X + int GetNofBinsX() const { return fNx; } + + /// Get number of bins along Y + int GetNofBinsY() const { return fNy; } + + /// Get index of the first bin entry in fHitsInBin array + ca::HitIndex_t GetFirstBinEntryIndex(int bin) const { return fFirstBinEntryIndex[(bin < fN) ? bin : fN]; } + + /// Get number of entries in the bin + //ca::HitIndex_t GetNofBinEntries(int bin) const { return fNofBinEntries[bin < fN ? bin : fN]; } + + /// Get entries + const ca::Vector<ca::GridEntry>& GetEntries() const { return fEntries; } + + /// Get minimal X value + fscal GetMinX() const { return fMinX; } + + /// Get minimal Y value + fscal GetMinY() const { return fMinY; } + + /// Get bin width in X + fscal GetBinWidthX() const { return fBinWidthX; } + + /// Get bin width in Y + fscal GetBinWidthY() const { return fBinWidthY; } + + /// Get maximal entry range in X + fscal GetMaxRangeX() const { return fMaxRangeX; } + + /// Get maximal entry range in Y + fscal GetMaxRangeY() const { return fMaxRangeY; } + + /// Get maximal entry range in T + fscal GetMaxRangeT() const { return fMaxRangeT; } + + private: + /// --- Data members --- + + 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 + + fscal fMaxRangeX {0.}; ///< maximal entry range in X + fscal fMaxRangeY {0.}; ///< maximal entry range in Y + fscal fMaxRangeT {0.}; ///< maximal entry range in T + + ca::Vector<ca::HitIndex_t> fFirstBinEntryIndex {"Grid::fFirstBinEntryIndex", 1, + 0}; ///< index of the first entry in the bin + + ca::Vector<ca::HitIndex_t> fNofBinEntries {"Grid::fNofBinEntries", 1, 0}; ///< number of hits in the bin + + ca::Vector<ca::GridEntry> fEntries { + "Ca::Grid::fEntries"}; ///< grid entries with references to the hit index in fWindowHits + }; + + /// --- Inline methods --- + + inline int Grid::GetBin(fscal X, fscal Y) const + { + // + return GetBinY(Y) * fNx + GetBinX(X); + } + + inline int Grid::GetBinX(fscal X) const + { + int binX = static_cast<int>((X - fMinX) * fBinWidthXinv); + return std::max(0, std::min(fNx - 1, binX)); + } + + inline int Grid::GetBinY(fscal Y) const + { + int binY = static_cast<int>((Y - fMinY) * fBinWidthYinv); + return std::max(0, std::min(fNy - 1, binY)); + } + + inline std::tuple<fscal, fscal> Grid::GetBinBoundsX(int iBin) const + { + fscal Xmin = fMinX + (iBin % fNx) * fBinWidthX; + return std::make_tuple(Xmin, Xmin + fBinWidthX); + } + + inline std::tuple<fscal, fscal> Grid::GetBinBoundsY(int iBin) const + { + fscal Ymin = fMinY + (iBin / fNx) * fBinWidthY; + return std::make_tuple(Ymin, Ymin + fBinWidthY); + } + +} // namespace cbm::algo::ca + +#endif diff --git a/algo/ca/core/data/CaGridArea.h b/algo/ca/core/data/CaGridArea.h new file mode 100644 index 0000000000000000000000000000000000000000..39632d9e4cb2980c16966802d41694c48d202293 --- /dev/null +++ b/algo/ca/core/data/CaGridArea.h @@ -0,0 +1,124 @@ +/* Copyright (C) 2012-2020 Frankfurt Institute for Advanced Studies, Goethe-Universität Frankfurt, Frankfurt + SPDX-License-Identifier: GPL-3.0-only + Authors: Maksym Zyzak, Igor Kulakov [committer] */ + +/// \file CaGrid.h + +#ifndef CaGridArea_H +#define CaGridArea_H 1 + +#include "CaGrid.h" +#include "CaHit.h" +#include "CaSimd.h" + +namespace cbm::algo::ca +{ + + namespace + { + using namespace cbm::algo; // to keep 'ca::' prefices in the code + } + + /// \brief Class for accessing objects in the 2D area that are stored in ca::Grid + /// + class GridArea { + public: + /// \brief Constructor + /// \param grid - the grid to work with + /// \param x - X coordinate of the center of the area + /// \param y - Y coordinate of the center of the area + /// \param dx - half-width of the area in X + /// \param dy - half-width of the area in Y + GridArea(const ca::Grid& grid, fscal x, fscal y, fscal dx, fscal dy); + + + /// \brief look up the next grid entry in the requested area + /// \return ind - the entry index in the grid.GetEntries() vector + /// \return true if the entry is found, false if there are no more entries in the area + bool GetNextGridEntry(ca::HitIndex_t& ind); + + /// \brief look up the next object id in the requested area + /// \return objectId - the id of the object + /// \return true if the object is found, false if there are no more pbjects in the area + bool GetNextObjectId(ca::HitIndex_t& objectId); + + /// \brief debug mode: loop over the entire GetEntries() vector ignoring the search area + void DoLoopOverEntireGrid(); + + private: + const ca::Grid& fGrid; + + int fAreaLastBinY {0}; // last Y bin of the area + int fAreaNbinsX {0}; // number of area bins in X + int fAreaFirstBin {0}; // first bin of the area (left-down corner of the area) + int fAreaCurrentBinY {0}; // current Y bin (incremented while iterating) + ca::HitIndex_t fCurentEntry {0}; // index of the current entry in fGrid.GetEntries() + ca::HitIndex_t fEntriesXend {0}; // end of the hit indices in current y-row + int fGridNbinsX {0}; // Number of grid bins in X (copy of fGrid::GetNofBinsX()) + }; + + inline GridArea::GridArea(const ca::Grid& grid, fscal x, fscal y, fscal dx, fscal dy) + : fGrid(grid) + , fGridNbinsX(fGrid.GetNofBinsX()) + { + int binXmin = fGrid.GetBinX(x - dx); + int binXmax = fGrid.GetBinX(x + dx); + int binYmin = fGrid.GetBinY(y - dy); + int binYmax = fGrid.GetBinY(y + dy); + + fAreaLastBinY = binYmax; + fAreaNbinsX = (binXmax - binXmin + 1); // bin index span in x direction + + fAreaFirstBin = (binYmin * fGridNbinsX + binXmin); + + fAreaCurrentBinY = binYmin; + fCurentEntry = fGrid.GetFirstBinEntryIndex(fAreaFirstBin); + fEntriesXend = fGrid.GetFirstBinEntryIndex(fAreaFirstBin + fAreaNbinsX); + } + + inline bool GridArea::GetNextGridEntry(ca::HitIndex_t& ind) + { + bool xIndexOutOfRange = (fCurentEntry >= fEntriesXend); // current entry is not in the area + + // jump to the next y row if fCurentEntry is outside of the X area + // TODO use ISUNLIKELY() macro here + while (xIndexOutOfRange) { + if (fAreaCurrentBinY >= fAreaLastBinY) { return false; } + fAreaCurrentBinY++; // get new y-line + fAreaFirstBin += fGridNbinsX; // move the left-down corner of the area to the next y-line + fCurentEntry = fGrid.GetFirstBinEntryIndex(fAreaFirstBin); // get first hit in cell, if y-line is new + fEntriesXend = fGrid.GetFirstBinEntryIndex(fAreaFirstBin + fAreaNbinsX); + xIndexOutOfRange = (fCurentEntry >= fEntriesXend); + } + + // TODO:: include L1_ASSERT and uncomment + //L1_ASSERT(fCurentEntry < fGrid.FirstHitInBin(fGrid.N()) || xIndexOutOfRange, + // fCurentEntry << " < " << fGrid.FirstHitInBin(fGrid.N()) << " || " << xIndexOutOfRange); + + ind = fCurentEntry; // return value + fCurentEntry++; // go to next + return true; + } + + inline bool GridArea::GetNextObjectId(ca::HitIndex_t& objectId) + { + ca::HitIndex_t entryIndex = 0; + + bool ok = GetNextGridEntry(entryIndex); + if (ok) { objectId = fGrid.GetEntries()[entryIndex].GetObjectId(); } + return ok; + } + + inline void GridArea::DoLoopOverEntireGrid() + { + fCurentEntry = 0; + fEntriesXend = fGrid.GetEntries().size(); + fAreaLastBinY = 0; + fAreaNbinsX = 0; + fAreaFirstBin = 0; + fAreaCurrentBinY = 0; + } + + +} // namespace cbm::algo::ca +#endif diff --git a/algo/ca/core/data/CaGridEntry.h b/algo/ca/core/data/CaGridEntry.h index 9a9b50d3bbe8fbc0b263c5a436a47d194d009cfb..5fa81820fbc8ee0d65ce08cd9d0ab079ba06149b 100644 --- a/algo/ca/core/data/CaGridEntry.h +++ b/algo/ca/core/data/CaGridEntry.h @@ -22,17 +22,17 @@ namespace cbm::algo::ca void Set(const ca::Hit& hit, ca::HitIndex_t id) { - fWindowHitId = id; - x = hit.x; - y = hit.y; - z = hit.z; - t = hit.t; - rangeX = hit.rangeX; - rangeY = hit.rangeY; - rangeT = hit.rangeT; + fObjectId = id; + x = hit.x; + y = hit.y; + z = hit.z; + t = hit.t; + rangeX = hit.rangeX; + rangeY = hit.rangeY; + rangeT = hit.rangeT; } - ca::HitIndex_t WindowHitId() const { return fWindowHitId; } + ca::HitIndex_t GetObjectId() const { return fObjectId; } fscal Z() const { return z; } fscal T() const { return t; } @@ -45,7 +45,7 @@ namespace cbm::algo::ca private: constexpr static fscal kUndef = 0.; // TODO: test with constants::Undef<fscal> - ca::HitIndex_t fWindowHitId {0}; ///< hit id in L1Algo::fWindowHits array + ca::HitIndex_t fObjectId {0}; ///< hit id in L1Algo::fWindowHits array fscal x {kUndef}; // x coordinate of the hit fscal y {kUndef}; // y coordinate of the hit diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt index 75c05e0b1edfcb40158df62cda0083770c2c1b7f..75be45045929bc256d9dd11aae0fbfa224d098e3 100644 --- a/reco/L1/CMakeLists.txt +++ b/reco/L1/CMakeLists.txt @@ -36,7 +36,6 @@ set(SRCS L1Algo/L1CaTrackFinderSlice.cxx L1Algo/L1BranchExtender.cxx L1Algo/L1TrackFitter.cxx - L1Algo/L1Grid.cxx CbmL1Util.cxx CbmL1Performance.cxx CbmL1ReadEvent.cxx diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx index 890470d89ea361003d1c76fe0f4a94b22c4b587b..82c52410de8aa337af0eaea0839af489bf9b7fb8 100644 --- a/reco/L1/L1Algo/L1Algo.cxx +++ b/reco/L1/L1Algo/L1Algo.cxx @@ -8,7 +8,6 @@ #include "CaGridEntry.h" #include "CaToolsDebugger.h" -#include "L1Grid.h" L1Algo::L1Algo() { diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h index ea83358281f6c02d817c515be6a3141194628f7d..535dcf0a29a2e33f1ec162beafc20bdbc702a5da 100644 --- a/reco/L1/L1Algo/L1Algo.h +++ b/reco/L1/L1Algo/L1Algo.h @@ -30,6 +30,7 @@ class L1AlgoDraw; #include "CaConstants.h" #include "CaField.h" +#include "CaGrid.h" #include "CaGridEntry.h" #include "CaHit.h" #include "CaInputData.h" @@ -41,7 +42,6 @@ class L1AlgoDraw; #include "L1Branch.h" #include "L1CloneMerger.h" #include "L1Fit.h" -#include "L1Grid.h" #include "L1Triplet.h" #include "L1Utils.h" // ? DEPRECATED ? @@ -288,11 +288,7 @@ private: public: Vector<L1HitTimeInfo> fHitTimeInfo; - L1Grid vGrid[constants::size::MaxNstations]; ///< - - fscal fMaxRangeX[constants::size::MaxNstations]; - fscal fMaxRangeY[constants::size::MaxNstations]; - fscal fMaxRangeT[constants::size::MaxNstations]; + ca::Grid vGrid[constants::size::MaxNstations]; ///< /// ----- Data used during track finding ----- /// diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/reco/L1/L1Algo/L1BranchExtender.cxx index 7f516a8190067cecf8eceadf0d31f7ed33dd3c66..6f524e3b53393ebb44e7cf282b6a6924cb212523 100644 --- a/reco/L1/L1Algo/L1BranchExtender.cxx +++ b/reco/L1/L1Algo/L1BranchExtender.cxx @@ -4,13 +4,12 @@ #include <iostream> -#include "CaGridEntry.h" +#include "CaGridArea.h" #include "CaTrack.h" #include "CaTrackParam.h" #include "CaVector.h" #include "L1Algo.h" #include "L1Branch.h" -#include "L1HitArea.h" // using namespace std; using cbm::algo::ca::Vector; // TMP!! @@ -200,24 +199,18 @@ void L1Algo::FindMoreHits(L1Branch& t, TrackParamV& Tout, const bool upstream, c TrackParamV& tr = fit.Tr(); - L1HitArea area(vGrid[ista], tr.X()[0], tr.Y()[0], - (sqrt(fPickGather * tr.C00()) + fMaxRangeX[ista] + fMaxDZ * abs(tr.Tx()))[0], - (sqrt(fPickGather * tr.C11()) + fMaxRangeY[ista] + fMaxDZ * abs(tr.Ty()))[0]); + ca::GridArea area(vGrid[ista], tr.X()[0], tr.Y()[0], + (sqrt(fPickGather * tr.C00()) + vGrid[ista].GetMaxRangeX() + fMaxDZ * abs(tr.Tx()))[0], + (sqrt(fPickGather * tr.C11()) + vGrid[ista].GetMaxRangeY() + fMaxDZ * abs(tr.Ty()))[0]); - for (ca::HitIndex_t ih = -1; true;) { // loop over the hits in the area + if (fParameters.DevIsIgnoreHitSearchAreas()) { area.DoLoopOverEntireGrid(); } - if (fParameters.DevIsIgnoreHitSearchAreas()) { - ih++; - if ((ca::HitIndex_t) ih >= vGrid[ista].Entries().size()) { break; } - } - else { - if (!area.GetNext(ih)) { break; } - } + ca::HitIndex_t ih = 0; - const ca::GridEntry& gridEntry = vGrid[ista].Entries()[ih]; + while (area.GetNextObjectId(ih)) { // loop over the hits in the area - if (fIsWindowHitSuppressed[gridEntry.WindowHitId()]) { continue; } - const ca::Hit& hit = fWindowHits[gridEntry.WindowHitId()]; + if (fIsWindowHitSuppressed[ih]) { continue; } + const ca::Hit& hit = fWindowHits[ih]; if (sta.timeInfo && tr.NdfTime()[0] > -2.) { fscal dt = hit.t - tr.Time()[0]; @@ -245,11 +238,11 @@ void L1Algo::FindMoreHits(L1Branch& t, TrackParamV& Tout, const bool upstream, c fscal d2 = d_x * d_x + d_y * d_y; if (d2 > r2_best) continue; - fscal dxm_est = sqrt(pickGather2 * C00)[0] + fMaxRangeX[ista]; + fscal dxm_est = sqrt(pickGather2 * C00)[0] + vGrid[ista].GetMaxRangeX(); if (fabs(d_x) > dxm_est) continue; r2_best = d2; - iHit_best = gridEntry.WindowHitId(); + iHit_best = ih; } if (iHit_best < 0) break; diff --git a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx index 4c9f6060b540ce383bc251d3705569dfccf49b8c..a184e58eaca57c2ba81805c6f5f70b04791443d5 100644 --- a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx +++ b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx @@ -22,6 +22,7 @@ #include "CbmL1.h" #include "CbmL1MCTrack.h" +#include "CaGrid.h" #include "CaGridEntry.h" #include "CaTrack.h" #include "CaTrackParam.h" @@ -29,8 +30,6 @@ #include "L1Assert.h" #include "L1Branch.h" #include "L1Fit.h" -#include "L1Grid.h" -#include "L1HitArea.h" #include "L1TripletConstructor.h" #ifdef DRAW #include "utils/L1AlgoDraw.h" @@ -226,9 +225,6 @@ void L1Algo::CaTrackFinderSlice() for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) { - fMaxRangeX[iS] = 0.; - fMaxRangeY[iS] = 0.; - fMaxRangeT[iS] = 0.; bool timeUninitialised = 1; fscal lasttime = 0; fscal starttime = 0; @@ -242,10 +238,6 @@ void L1Algo::CaTrackFinderSlice() for (ca::HitIndex_t ih = 0; ih < fSliceHitIds[iS].size(); ++ih) { const ca::Hit& h = fInputData.GetHit(fSliceHitIds[iS][ih]); - if (fMaxRangeX[iS] < h.rangeX) { fMaxRangeX[iS] = h.rangeX; } - if (fMaxRangeY[iS] < h.rangeY) { fMaxRangeY[iS] = h.rangeY; } - if (fMaxRangeT[iS] < h.rangeT) { fMaxRangeT[iS] = h.rangeT; } - if (h.x < gridMinX) { gridMinX = h.x; } if (h.x > gridMaxX) { gridMaxX = h.x; } if (h.y < gridMinY) { gridMinY = h.y; } @@ -363,10 +355,10 @@ void L1Algo::CaTrackFinderSlice() for (int istal = fParameters.GetNstationsActive() - 2; istal >= fFirstCAstation; istal--) { // start with downstream chambers - Tindex nGridEntriesL = vGrid[istal].Entries().size(); + Tindex nGridEntriesL = vGrid[istal].GetEntries().size(); for (Tindex iel = 0; iel < nGridEntriesL; ++iel) { - ca::HitIndex_t ihitl = vGrid[istal].Entries()[iel].WindowHitId(); + ca::HitIndex_t ihitl = vGrid[istal].GetEntries()[iel].GetObjectId(); constructor1.CreateTripletsForHit(istal, istal + 1, istal + 2, ihitl); if (fpCurrentIteration->GetJumpedFlag()) { diff --git a/reco/L1/L1Algo/L1Grid.cxx b/reco/L1/L1Algo/L1Grid.cxx deleted file mode 100644 index a31c1f8a03046f358918c045dd29ff0c9a488c98..0000000000000000000000000000000000000000 --- a/reco/L1/L1Algo/L1Grid.cxx +++ /dev/null @@ -1,119 +0,0 @@ -/* Copyright (C) 2017-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Maksym Zyzak [committer], Valentina Akishina */ - -#include "L1Grid.h" - -#include <algorithm> - -#include <string.h> - -#include "CaHit.h" - -using namespace std; // !! REMOVE - -namespace -{ - using cbm::algo::ca::Vector; // TMP!! - using namespace cbm::algo::ca; //TODO: remove - using namespace cbm::algo; //TODO: remove -} // namespace - - -/// Copy to memory block [@dest, @dest+@num] num number of times the value of i of type @T with size @typesize. -/// uses binary expansion of copied volume for speed up -//TODO: move it to ca::Vector -template<typename T> - -inline void memset(T* dest, T i, size_t num) -{ - 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; - memcpy(dest + lastBin + 1, dest, l * tsize); - lastBin += l; - } -} - - -void L1Grid::BuildBins(fscal xMin, fscal xMax, fscal yMin, fscal yMax, fscal sx, fscal sy) -{ - fMinX = std::min(xMin, xMax); - fMinY = std::min(yMin, yMax); - - xMax = std::max(xMin, xMax); - yMax = std::max(yMin, yMax); - - fBinWidthX = sx; - fBinWidthY = sy; - - fBinWidthXinv = 1. / fBinWidthX; - fBinWidthYinv = 1. / fBinWidthY; - - 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; - - fN = fNx * fNy; - - fEntries.clear(); - fBinFirstEntryIndex.reset(fN + 2, 0); - fNbinHits.reset(fN + 2, 0); -} - - -void L1Grid::StoreHits(const Vector<ca::Hit>& hits, ca::HitIndex_t hitStartIndex, ca::HitIndex_t nHits, - const Vector<unsigned char>& hitKeyFlags) -{ - fBinFirstEntryIndex.reset(fN + 2, 0); - fNbinHits.reset(fN + 2, 0); - - int nEntries = 0; - for (ca::HitIndex_t ih = 0; ih < nHits; ih++) { - const ca::Hit& hit = hits[hitStartIndex + ih]; - if (!(hitKeyFlags[hit.f] || hitKeyFlags[hit.b])) { - fNbinHits[GetBinBounded(hit.x, hit.y)]++; - nEntries++; - } - } - - fEntries.reset(nEntries); - - for (int bin = 0; bin < fN + 1; bin++) { - fBinFirstEntryIndex[bin + 1] = fBinFirstEntryIndex[bin] + fNbinHits[bin]; - fNbinHits[bin] = 0; - } - - for (ca::HitIndex_t ih = 0; ih < nHits; ih++) { - const ca::Hit& hit = hits[hitStartIndex + ih]; - if (!(hitKeyFlags[hit.f] || hitKeyFlags[hit.b])) { - int bin = GetBinBounded(hit.x, hit.y); - fEntries[fBinFirstEntryIndex[bin] + fNbinHits[bin]].Set(hit, hitStartIndex + ih); - fNbinHits[bin]++; - } - } -} - -void L1Grid::RemoveUsedHits(const Vector<ca::Hit>& hits, const Vector<unsigned char>& hitKeyFlags) -{ - int nEntries = 0; - for (int bin = 0; bin < fN; bin++) { - ca::HitIndex_t firstEntry = fBinFirstEntryIndex[bin]; - fBinFirstEntryIndex[bin] = nEntries; - fNbinHits[bin] = 0; - for (ca::HitIndex_t i = firstEntry; i < fBinFirstEntryIndex[bin + 1]; i++) { - const ca::Hit& hit = hits[fEntries[i].WindowHitId()]; - if (!(hitKeyFlags[hit.f] || hitKeyFlags[hit.b])) { - fEntries[nEntries] = fEntries[i]; - nEntries++; - fNbinHits[bin]++; - } - } - } - fBinFirstEntryIndex[fN] = nEntries; - fEntries.reduce(nEntries); -} diff --git a/reco/L1/L1Algo/L1Grid.h b/reco/L1/L1Algo/L1Grid.h deleted file mode 100644 index 3aab7ad51bea3050880ede31edf522c616cc3f1d..0000000000000000000000000000000000000000 --- a/reco/L1/L1Algo/L1Grid.h +++ /dev/null @@ -1,153 +0,0 @@ -/* Copyright (C) 2010 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt - SPDX-License-Identifier: GPL-3.0-only - Authors: Igor Kulakov [committer] */ - -/// \file CaGrid.h -/// \brief A class to store hit information in a backet-sorted way on 2D grid -/// -/// This code is based on the code of the ALICE HLT Project -/// - -#ifndef L1Grid_h -#define L1Grid_h - -#include "CaGridEntry.h" -#include "CaHit.h" -#include "CaSimd.h" -#include "CaVector.h" - -namespace -{ - using cbm::algo::ca::Vector; // TMP!! - using namespace cbm::algo::ca; //TODO: remove - using namespace cbm::algo; //TODO: remove -} // 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 BuildBins(fscal xMin, fscal xMax, fscal yMin, fscal yMax, fscal sx, fscal sy); - - void StoreHits(const Vector<ca::Hit>& hits, ca::HitIndex_t hitStartIndex, ca::HitIndex_t nHits, - const Vector<unsigned char>& hitKeyFlags); - - void RemoveUsedHits(const Vector<ca::Hit>& hits, const Vector<unsigned char>& hitKeyFlags); - - /// get bin index - int GetBin(fscal X, fscal Y) const; - - /// 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; - - /// get bin bounds - void GetBinBounds(int iBin, fscal& Xmin, fscal& Xmax, fscal& Ymin, fscal& Ymax) const; - - /// get number of bins - int N() const { return fN; } - - /// get number of bins in X - int Nx() const { return fNx; } - - /// get number of bins in Y - int Ny() const { return fNy; } - - /// get index of the first bin entry in fHitsInBin array - ca::HitIndex_t FirstHitInBin(int i) const - { - if (i < (fN + 1)) return fBinFirstEntryIndex[i]; - else - return fBinFirstEntryIndex[fN + 1]; - } - - /// get number of hits in the bin - ca::HitIndex_t NhitsInBin(int i) const - { - if (i < (fN + 1)) return fNbinHits[i]; - else - return 0; - } - - /// get entries - Vector<ca::GridEntry>& Entries() { return fEntries; } - -private: - 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"}; ///< index of the first entry in the bin - Vector<ca::HitIndex_t> fNbinHits {"L1Grid::fNbinHits"}; ///< number of hits in the bin - Vector<ca::GridEntry> fEntries { - "Ca::Grid::fEntries"}; ///< grid entries with references to the hit index in fWindowHits -}; - - -inline int L1Grid::GetBin(fscal X, fscal Y) const -{ - /// 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 < fNx); - assert(yBin < fNy); - return yBin * fNx + xBin; -} - - -inline void L1Grid::GetBinBounds(int iBin, fscal& Xmin, fscal& Xmax, fscal& Ymin, fscal& Ymax) const -{ - int yBin = iBin / fNx; - int xBin = iBin % fNx; - Xmin = fMinX + xBin * fBinWidthX; - Ymin = fMinY + yBin * fBinWidthY; - Xmax = Xmin + fBinWidthX; - Ymax = Ymin + fBinWidthY; -} - - -inline int L1Grid::GetBinBounded(fscal X, fscal Y) const -{ - /// get bin index - int xBin, yBin; - GetBinBounded(X, Y, xBin, yBin); - return yBin * fNx + xBin; -} - - -inline void L1Grid::GetBinBounded(fscal X, fscal Y, int& bX, int& bY) const -{ - 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)); -} - - -#endif diff --git a/reco/L1/L1Algo/L1HitArea.h b/reco/L1/L1Algo/L1HitArea.h deleted file mode 100644 index 206fe60ae1536fe3a156933120cf2e452892d910..0000000000000000000000000000000000000000 --- a/reco/L1/L1Algo/L1HitArea.h +++ /dev/null @@ -1,104 +0,0 @@ -/* Copyright (C) 2012-2020 Frankfurt Institute for Advanced Studies, Goethe-Universität Frankfurt, Frankfurt - SPDX-License-Identifier: GPL-3.0-only - Authors: Maksym Zyzak, Igor Kulakov [committer] */ - -#ifndef L1HitArea_H -#define L1HitArea_H - -#include "CaHit.h" -#include "L1Def.h" -#include "L1Grid.h" - -class L1Row; -class L1SliceData; - -using namespace cbm::algo::ca; //TODO: remove -using namespace cbm::algo; //TODO: remove - -class L1HitArea { -public: - L1HitArea(const L1Grid& grid, float x, float y, float dx, float dy); - - /** - * look up the next hit in the requested area. - * Sets h to the coordinates and returns the index for the hit data - */ - bool GetNext(ca::HitIndex_t& i); - -protected: - const L1Grid& fGrid; - - 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 -}; - -inline L1HitArea::L1HitArea(const L1Grid& grid, float x, float y, float dx, float dy) - : fGrid(grid) - , fBYmax(0) - , fBDX(0) - , fIndXmin(0) - , fIy(0) - , fHitXlst(0) - , fIh(0) - , fNx(fGrid.Nx()) -{ - const float minX = x - dx; - const float maxX = x + dx; - const float minY = y - dy; - const float maxY = y + dy; - int bXmin, bYmin, bXmax; // boundary bin indexes - fGrid.GetBinBounded(minX, minY, bXmin, bYmin); - fGrid.GetBinBounded(maxX, maxY, bXmax, fBYmax); - - fBDX = (bXmax - bXmin + 1); // bin index span in x direction - - fIndXmin = (bYmin * fNx + bXmin); // same as grid.GetBin(fminX, fMinY), i.e. the smallest bin index of interest - // fIndXmin + fBDX - 1 then is the largest bin index of interest with the same Y - - fGrid.GetBinBounds(fIndXmin, x, dx, y, dy); - - fIy = bYmin; - - fIh = fGrid.FirstHitInBin(fIndXmin); - fHitXlst = fGrid.FirstHitInBin(fIndXmin + fBDX); -} - -inline bool L1HitArea::GetNext(ca::HitIndex_t& i) -{ - bool xIndexOutOfRange = fIh >= fHitXlst; // current x is not in the area - bool nextYIndexOutOfRange = (fIy >= fBYmax); // there isn't any new y-line - - if (xIndexOutOfRange && nextYIndexOutOfRange) { // all iterators are over the end - return false; - } - - // at least one entry in the vector has (fIh >= fHitXlst && fIy < fBYmax) - bool needNextY = xIndexOutOfRange && !nextYIndexOutOfRange; - - // skip as long as fIh is outside of the interesting bin x-index - while (ISLIKELY(needNextY)) { //ISLIKELY to speed the programm and optimise the use of registers - fIy++; // get new y-line - // get next hit - fIndXmin += fNx; - fIh = fGrid.FirstHitInBin(fIndXmin); // get first hit in cell, if y-line is new - fHitXlst = fGrid.FirstHitInBin(fIndXmin + fBDX); - - xIndexOutOfRange = fIh >= fHitXlst; - nextYIndexOutOfRange = (fIy >= fBYmax); - needNextY = xIndexOutOfRange && !nextYIndexOutOfRange; - } - - L1_ASSERT(fIh < fGrid.FirstHitInBin(fGrid.N()) || xIndexOutOfRange, - fIh << " < " << fGrid.FirstHitInBin(fGrid.N()) << " || " << xIndexOutOfRange); - i = fIh; // return - - fIh++; // go to next - return !xIndexOutOfRange; -} - -#endif diff --git a/reco/L1/L1Algo/L1TripletConstructor.cxx b/reco/L1/L1Algo/L1TripletConstructor.cxx index ee418f1679f2c3a03458042aa341fef12d1e37a0..c44e486717e80c83f679ce1b13bd539fe2c210c9 100644 --- a/reco/L1/L1Algo/L1TripletConstructor.cxx +++ b/reco/L1/L1Algo/L1TripletConstructor.cxx @@ -10,14 +10,12 @@ #include <algorithm> #include <iostream> -#include "CaGridEntry.h" +#include "CaGridArea.h" #include "CaToolsDebugger.h" #include "CaTrackParam.h" #include "L1Algo.h" #include "L1Assert.h" #include "L1Fit.h" -#include "L1Grid.h" -#include "L1HitArea.h" using cbm::ca::tools::Debugger; @@ -606,11 +604,7 @@ void L1TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, co collectedHits.clear(); collectedHits.reserve(maxNhits); - if (fAlgo->vGrid[iSta].Entries().size() == 0) { return; } - - const ca::Station& sta = fAlgo->fParameters.GetStation(iSta); - const ca::GridEntry* entries = &(fAlgo->vGrid[iSta].Entries()[0]); - const ca::HitIndex_t nHits = fAlgo->vGrid[iSta].Entries().size(); + const ca::Station& sta = fAlgo->fParameters.GetStation(iSta); fFit.SetTrack(Tr); TrackParamV& T = fFit.Tr(); @@ -622,24 +616,24 @@ void L1TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, co const fscal timeError2 = T.C55()[0]; const fscal time = T.Time()[0]; - L1HitArea area(fAlgo->vGrid[iSta], T.X()[0], T.Y()[0], - (sqrt(Pick_m22 * T.C00()) + fAlgo->fMaxRangeX[iSta] + fAlgo->fMaxDZ * abs(T.Tx()))[0], - (sqrt(Pick_m22 * T.C11()) + fAlgo->fMaxRangeY[iSta] + fAlgo->fMaxDZ * abs(T.Ty()))[0]); + const auto& grid = fAlgo->vGrid[iSta]; + ca::GridArea area(grid, T.X()[0], T.Y()[0], + (sqrt(Pick_m22 * T.C00()) + grid.GetMaxRangeX() + fAlgo->fMaxDZ * abs(T.Tx()))[0], + (sqrt(Pick_m22 * T.C11()) + grid.GetMaxRangeY() + fAlgo->fMaxDZ * abs(T.Ty()))[0]); - for (ca::HitIndex_t ih = 0; (ih < nHits) && ((int) collectedHits.size() < maxNhits); - ih++) { // loop over all station hits + if (fAlgo->fParameters.DevIsIgnoreHitSearchAreas()) { area.DoLoopOverEntireGrid(); } - if (!fAlgo->fParameters.DevIsIgnoreHitSearchAreas()) { // only loop over the hits in the area - if (!area.GetNext(ih)) { break; } - } + ca::HitIndex_t ih = 0; + + while (area.GetNextObjectId(ih) && ((int) collectedHits.size() < maxNhits)) { // loop over station hits - const ca::GridEntry& entry = entries[ih]; - if (fAlgo->fIsWindowHitSuppressed[entry.WindowHitId()]) continue; - const ca::Hit& hit = fAlgo->fWindowHits[entry.WindowHitId()]; + if (fAlgo->fIsWindowHitSuppressed[ih]) continue; + + const ca::Hit& hit = fAlgo->fWindowHits[ih]; if (iMC >= 0) { // match via MC - if (iMC != fAlgo->GetMcTrackIdForWindowHit(fAlgo->vGrid[iSta].Entries()[ih].WindowHitId())) { continue; } + if (iMC != fAlgo->GetMcTrackIdForWindowHit(ih)) { continue; } } // check y-boundaries @@ -687,7 +681,7 @@ void L1TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, co if (chi2x[0] + chi2u[0] > chi2Cut) continue; } - collectedHits.push_back(entry.WindowHitId()); + collectedHits.push_back(ih); } // loop over the hits in the area -} \ No newline at end of file +}