From 62e507d8aec2268e0c5b0b2ae3004bd9055fceb1 Mon Sep 17 00:00:00 2001 From: "se.gorbunov" <se.gorbunov@gsi.de> Date: Wed, 1 Mar 2023 02:38:41 +0000 Subject: [PATCH] L1: cleanup grid classes --- reco/L1/CbmL1.cxx | 2 + reco/L1/CbmL1Performance.cxx | 5 + reco/L1/L1Algo/L1Algo.cxx | 10 +- reco/L1/L1Algo/L1Algo.h | 2 +- reco/L1/L1Algo/L1CaTrackFinderSlice.cxx | 54 ++++--- reco/L1/L1Algo/L1Grid.cxx | 30 ++-- reco/L1/L1Algo/L1Grid.h | 160 +++++++++---------- reco/L1/L1Algo/L1HitArea.h | 198 ++++++++++++------------ 8 files changed, 241 insertions(+), 220 deletions(-) diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx index 0f956a441c..0a1dec0c40 100644 --- a/reco/L1/CbmL1.cxx +++ b/reco/L1/CbmL1.cxx @@ -899,6 +899,7 @@ void CbmL1::Reconstruct(CbmEvent* event) #endif // ----- Read data from branches and send data from IODataManager to L1Algo ---------------------------------------- + ReadEvent(event); @@ -910,6 +911,7 @@ void CbmL1::Reconstruct(CbmEvent* event) } } + if ((fPerformance) && (fSTAPDataMode < 2)) { InputPerformance(); } // FieldApproxCheck(); diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx index 9feb4cb4eb..547921ac25 100644 --- a/reco/L1/CbmL1Performance.cxx +++ b/reco/L1/CbmL1Performance.cxx @@ -1240,6 +1240,7 @@ void CbmL1::TrackFitPerformance() gDirectory = currentDir; } // if first call + for (vector<CbmL1Track>::iterator it = fvRecoTracks.begin(); it != fvRecoTracks.end(); ++it) { if (it->IsGhost()) continue; @@ -1482,6 +1483,7 @@ void CbmL1::TrackFitPerformance() // last TOF point do { + break; // only produce these plots in debug mode // only produce these plots in debug mode if (!fpTofPoints) break; const CbmL1MCTrack& mcTrack = *(it->GetMCTracks()[0]); @@ -1745,6 +1747,9 @@ void CbmL1::FieldIntegralCheck() void CbmL1::InputPerformance() { + // The input performance is currently evaluated by other CBM tasks + // TODO: obsolete code, remove it + // static TH1I *nStripFHits, *nStripBHits, *nStripFMC, *nStripBMC; static TH1F *resXsts, *resYsts, *resTsts, *resXmvd, *resYmvd /*, *pullZ*/; diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx index 20e47bdfbf..4980c5ef9a 100644 --- a/reco/L1/L1Algo/L1Algo.cxx +++ b/reco/L1/L1Algo/L1Algo.cxx @@ -169,13 +169,13 @@ void L1Algo::ReceiveParameters(L1Parameters&& parameters) } /// TODO: Move to L1Hit -void L1Algo::GetHitCoor(const L1Hit& _h, fscal& _x, fscal& _y, char iS) +std::pair<fscal, fscal> L1Algo::GetHitCoorOnGrid(const L1Hit& h, char iS) { const L1Station& sta = fParameters.GetStation(int(iS)); - std::tie(_x, _y) = sta.ConvUVtoXY<fscal>(_h.u, _h.v); - float dz = _h.z - fParameters.GetTargetPositionZ()[0]; - _x /= dz; - _y /= dz; + fscal x, y; + std::tie(x, y) = sta.ConvUVtoXY<fscal>(h.u, h.v); + float dz = h.z - fParameters.GetTargetPositionZ()[0]; + return {x / dz, y / dz}; } /// TODO: Move to L1Hit diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h index 71f0568377..4e8f2dcfbb 100644 --- a/reco/L1/L1Algo/L1Algo.h +++ b/reco/L1/L1Algo/L1Algo.h @@ -174,8 +174,8 @@ public: void GetHitCoor(const L1Hit& _h, fscal& _x, fscal& _y, fscal& _z, const L1Station& sta); + std::pair<fscal, fscal> GetHitCoorOnGrid(const L1Hit& h, char iS); - void GetHitCoor(const L1Hit& _h, fscal& _x, fscal& _y, char iS); L1HitPoint CreateHitPoint(const L1Hit& hit); // full the hit point by hit information. diff --git a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx index 1bbd894c68..6281c7f60e 100644 --- a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx +++ b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx @@ -1561,27 +1561,10 @@ void L1Algo::CaTrackFinderSlice() ResetSliceData(); - int nSliceHits = fSliceHitIds.size(); - #ifdef XXX c_timerG.Start(); #endif // XXX - fscal yStep = 1.5 / sqrt(nSliceHits); // empirics. 0.01*sqrt(2374) ~= 0.5 - - - // fscal yStep = 0.5 / sqrt(nSliceHits); // empirics. 0.01*sqrt(2374) ~= 0.5 - if (yStep > 0.3) yStep = 0.3; - fscal xStep = yStep * 3; - // fscal xStep = yStep * 3; - - // yStep = 0.0078; - // const fscal hitDensity = sqrt( nSliceHits ); - - // fscal yStep = 0.7*4/hitDensity; // empirics. 0.01*sqrt(2374) ~= 0.5 - // if (yStep > 0.3) - // yStep = 1.25; - // xStep = 2.05; L1HitIndex_t nGridHitsFilled = 0; @@ -1594,8 +1577,15 @@ void L1Algo::CaTrackFinderSlice() fscal lasttime = 0; fscal starttime = 0; + + fscal gridMinX = -0.1; + fscal gridMaxX = 0.1; + fscal gridMinY = -0.1; + fscal gridMaxY = 0.1; + for (L1HitIndex_t ih = fSliceHitIdsStartIndex[iS]; ih < fSliceHitIdsStopIndex[iS]; ++ih) { - const L1Hit& h = fInputData.GetHit(fSliceHitIds[ih]); + const L1Hit& h = fInputData.GetHit(fSliceHitIds[ih]); + auto [dxx, dxy, dyy] = st.FormXYCovarianceMatrix(h.du2, h.dv2); fscal dx = sqrt(dxx); fscal dy = sqrt(dyy); @@ -1605,14 +1595,40 @@ void L1Algo::CaTrackFinderSlice() if (fMaxDy[iS] < dy) { fMaxDy[iS] = dy; } if (fMaxDt[iS] < dt) { fMaxDt[iS] = dt; } + auto [x, y] = GetHitCoorOnGrid(h, iS); + if (x < gridMinX) { gridMinX = x; } + if (x > gridMaxX) { gridMaxX = x; } + if (y < gridMinY) { gridMinY = y; } + if (y > gridMaxY) { gridMaxY = y; } + const fscal time = h.t; assert(std::isfinite(time)); if (timeUninitialised || lasttime < time) { lasttime = time; } if (timeUninitialised || starttime > time) { starttime = time; } timeUninitialised = false; } - + /* + int nSliceHits = fSliceHitIds.size(); + fscal yStep = 1.5 / sqrt(nSliceHits + 1); // empirics. 0.01*sqrt(2374) ~= 0.5 + if (yStep > 0.3) yStep = 0.3; + if (yStep < 0.01) yStep = 0.01; + fscal xStep = yStep * 3; vGridTime[iS].BuildBins(-1, 1, -0.6, 0.6, starttime, lasttime, xStep, yStep, (lasttime - starttime + 1)); + */ + + int nSliceHits = fSliceHitIdsStopIndex[iS] - fSliceHitIdsStartIndex[iS]; + fscal sizeY = gridMaxY - gridMinY; + fscal sizeX = gridMaxX - gridMinX; + int nBins2D = 1 + nSliceHits; + fscal yStep = sizeY / sqrt(nBins2D); + fscal xStep = sizeX / sqrt(nBins2D); + if (yStep > 0.3) yStep = 0.3; + if (yStep < 0.01) yStep = 0.01; + if (xStep > 0.3) xStep = 0.3; + if (xStep < 0.01) xStep = 0.01; + + vGridTime[iS].BuildBins(gridMinX, gridMaxX, gridMinY, gridMaxY, starttime, lasttime, xStep, yStep, + (lasttime - starttime + 1)); vGridTime[iS].StoreHits(*this, iS, nGridHitsFilled); } diff --git a/reco/L1/L1Algo/L1Grid.cxx b/reco/L1/L1Algo/L1Grid.cxx index 5b80e0e806..835ed39f97 100644 --- a/reco/L1/L1Algo/L1Grid.cxx +++ b/reco/L1/L1Algo/L1Grid.cxx @@ -56,8 +56,7 @@ void L1Grid::UpdateIterGrid(unsigned int Nelements, L1Hit* hits, L1Vector<L1HitI const L1Hit& hit = hits[x]; if (!(vSFlag[hit.f] || vSFlag[hit.b])) { - Algo.GetHitCoor(hit, xs, ys, iS); - + std::tie(xs, ys) = Algo.GetHitCoorOnGrid(hit, iS); const L1HitIndex_t& bin = GetBinBounded(xs, ys, hit.t); fHitsInBin[x] = fFirstHitInBin[bin + 1]; @@ -98,7 +97,7 @@ void L1Grid::UpdateIterGrid(unsigned int Nelements, L1Hit* hits, L1Vector<L1HitI const L1Hit& hit = hits[x]; if (!(vSFlag[hit.f] || vSFlag[hit.b])) { - Algo.GetHitCoor(hit, xs, ys, iS); + std::tie(xs, ys) = Algo.GetHitCoorOnGrid(hit, iS); const L1HitIndex_t& bin = GetBinBounded(xs, ys, hit.t); @@ -143,24 +142,24 @@ void L1Grid::AllocateMemory(int NThreads) } -void L1Grid::BuildBins(float yMin, float yMax, float zMin, float zMax, float tMin, float tMax, float sy, float sz, +void L1Grid::BuildBins(float xMin, float xMax, float yMin, float yMax, float tMin, float tMax, float sx, float sy, float st) { + fStepXInv = 1.f / sx; fStepYInv = 1.f / sy; - fStepZInv = 1.f / sz; fStepTInv = 1.f / st; + fXMinOverStep = xMin * fStepXInv; fYMinOverStep = yMin * fStepYInv; - fZMinOverStep = zMin * fStepZInv; fTMinOverStep = tMin * fStepTInv; + fNx = (xMax * fStepXInv - fXMinOverStep + 1.f); fNy = (yMax * fStepYInv - fYMinOverStep + 1.f); - fNz = (zMax * fStepZInv - fZMinOverStep + 1.f); fNt = (tMax * fStepTInv - fTMinOverStep + 1.f); // unsigned int Nold = fN; - fN = fNy * fNz * fNt; + fN = fNx * fNy * fNt; fBinInGrid = (((fN) / fNThreads) + 1); } @@ -174,18 +173,16 @@ void L1Grid::StoreHits(L1Algo& algo, int iS, L1HitIndex_t& nGridHitsFilled) algo.fGridHitStartIndex[iS] = nGridHitsFilled; - fscal x = 0; - fscal y = 0; - fFirstHitInBin.reset(fN + 2, 0); #ifdef _OPENMP -#pragma omp parallel for firstprivate(x, y) +#pragma omp parallel for #endif for (L1HitIndex_t ih = 0; ih < nHits; ih++) { L1HitIndex_t caHitId = algo.fSliceHitIds[firstHitIndex + ih]; const L1Hit& h = algo.GetInputData().GetHit(caHitId); - algo.GetHitCoor(h, x, y, iS); + fscal x, y; + std::tie(x, y) = algo.GetHitCoorOnGrid(h, iS); auto bin = GetBinBounded(x, y, h.t); fHitsInBin[ih] = fFirstHitInBin[bin + 1]; #ifdef _OPENMP @@ -219,12 +216,13 @@ void L1Grid::StoreHits(L1Algo& algo, int iS, L1HitIndex_t& nGridHitsFilled) } -#pragma omp parallel for firstprivate(x, y) +#pragma omp parallel for for (L1HitIndex_t ih = 0; ih < nHits; ih++) { L1HitIndex_t caHitId = algo.fSliceHitIds[firstHitIndex + ih]; const L1Hit& h = algo.GetInputData().GetHit(caHitId); - algo.GetHitCoor(h, x, y, iS); - auto bin = GetBinBounded(x, y, h.t); + fscal x, y; + std::tie(x, y) = algo.GetHitCoorOnGrid(h, iS); + auto bin = GetBinBounded(x, y, h.t); { const L1HitIndex_t& index1 = fFirstHitInBin[bin] + fHitsInBin[ih]; algo.fGridHits[nGridHitsFilled + index1] = h; diff --git a/reco/L1/L1Algo/L1Grid.h b/reco/L1/L1Algo/L1Grid.h index 58fa660eb4..f6a8653fa0 100644 --- a/reco/L1/L1Algo/L1Grid.h +++ b/reco/L1/L1Algo/L1Grid.h @@ -31,53 +31,53 @@ class L1Algo; * @class L1Grid * * 2-dimensional grid of pointers. - * pointers to (y,z)-like objects are assigned to the corresponding grid bin + * 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 Z,Y to be similar to TPC row coordinates. + * grid axis are named Y,x to be similar to TPC row coordinates. */ class L1Grid { public: L1Grid() = default; - L1Grid(const L1Grid& grid) : fN(grid.N()), fNy(grid.Ny()), fNz(grid.Nz()), fNt(grid.Nt()) {} + L1Grid(const L1Grid& grid) : fN(grid.N()), fNx(grid.Nx()), fNy(grid.Ny()), fNt(grid.Nt()) {} ~L1Grid() = default; void StoreHits(L1Algo& Algo, int iS, L1HitIndex_t& nGridHitsFilled); - void CreatePar0(float yMin, float yMax, float zMin, float zMax, float sy, float sz); - void BuildBins(float yMin, float yMax, float zMin, float zMax, float tMin, float tMax, float sy, float sz, float st); + void CreatePar0(float xMin, float xMax, float yMin, float yMax, float sx, float sy); + void BuildBins(float xMin, float xMax, float yMin, float yMax, float tMin, float tMax, float sx, float sy, float st); void Initial1(int NThreads); void AllocateMemory(int NThreads); - void Create(float yMin, float yMax, float zMin, float zMax, float sy, float sz); + void Create(float xMin, float xMax, float yMin, float yMax, float sx, float sy); void Fill(const L1HitPoint* points, L1HitIndex_t n); // call after sort void FillPar(const L1HitPoint* points, L1HitIndex_t n); - int GetBin(float Y, float Z) const; + int GetBin(float X, float Y) const; - unsigned int GetBinBounded(const float& Y, const float& Z) const; - void GetBinBounded(const float& Y, const float& Z, unsigned short& bY, unsigned short& bZ) 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; - int GetBin(float Y, float Z, float T) const; + int GetBin(float X, float Y, float T) const; - // static unsigned short GetBinBounded( const L1Grid *array, const unsigned short &indexes, const float &Y, const float &Z ); - // static void GetBinBounded( const L1Grid *array, const unsigned short &indexes, const float &Y, const float &Z, unsigned short *bY, unsigned short *bZ ); - // static unsigned short Ny( const L1Grid *array, const unsigned short &indexes ) { return unsigned short( array, &L1Grid::fNy, indexes ); } + // static unsigned short GetBinBounded( const L1Grid *array, const unsigned short &indexes, const float &X, const float &Y ); + // static void GetBinBounded( const L1Grid *array, const unsigned short &indexes, const float &X, const float &Y, unsigned short *bX, unsigned short *bY ); + // static unsigned short Nx( const L1Grid *array, const unsigned short &indexes ) { return unsigned short( array, &L1Grid::fNx, indexes ); } - void GetBinBounds(unsigned int iBin, float& Ymin, float& Ymax, float& Zmin, float& Zmax) const; - unsigned int GetBinBounded(const float& Y, const float& Z, const float& T) const; - void GetBinBounded(const float& Y, const float& Z, const float& T, unsigned short& bY, unsigned short& bZ, + void GetBinBounds(unsigned int iBin, float& Xmin, float& Xmax, float& Ymin, float& Ymax) const; + unsigned int GetBinBounded(const float& X, const float& Y, const float& T) const; + void GetBinBounded(const float& X, const float& Y, const float& T, unsigned short& bX, unsigned short& bY, unsigned short& bT) const; - void GetBinBounds(unsigned int iBin, float& Ymin, float& Ymax, float& Zmin, float& Zmax, float& Tmin, + void GetBinBounds(unsigned int iBin, float& Xmin, float& Xmax, float& Ymin, float& Ymax, float& Tmin, float& Tmax) const; unsigned int N() const { return fN; } + unsigned short Nx() const { return fNx; } unsigned short Ny() const { return fNy; } - unsigned short Nz() const { return fNz; } unsigned short Nt() const { return fNt; } L1HitIndex_t FirstHitInBin(unsigned int i) const @@ -103,15 +103,15 @@ public: private: unsigned int fN {0}; //* total N bins + unsigned short fNx {0}; //* N bins in X unsigned short fNy {0}; //* N bins in Y - unsigned short fNz {0}; //* N bins in Z - unsigned short fNt {0}; //* N bins in Z + unsigned short fNt {0}; //* N bins in T + float fXMinOverStep {0.f}; //* minimal X value * fStepXInv float fYMinOverStep {0.f}; //* minimal Y value * fStepYInv - float fZMinOverStep {0.f}; //* minimal Z value * fStepZInv - float fTMinOverStep {0.f}; //* minimal Z value * fStepZInv + float fTMinOverStep {0.f}; //* minimal T value * fStepTInv + float fStepXInv {0.f}; //* inverse bin size in X float fStepYInv {0.f}; //* inverse bin size in Y - float fStepZInv {0.f}; //* inverse bin size in Z - float fStepTInv {0.f}; //* inverse bin size in Z + float fStepTInv {0.f}; //* inverse bin size in T int fBinInGrid {0}; unsigned short fNThreads {0}; @@ -121,136 +121,136 @@ private: // vector <omp_lock_t> lock; }; -// inline unsigned short L1Grid::GetBinBounded( const L1Grid *array, const unsigned short &indexes, const float &Y, const float &Z ) +// inline unsigned short L1Grid::GetBinBounded( const L1Grid *array, const unsigned short &indexes, const float &X, const float &Y ) // { -// const float fZMinOverStep( array, &L1Grid::fZMinOverStep, indexes ); -// const float fStepZInv( array, &L1Grid::fStepZInv, indexes ); -// const unsigned short fNz( array, &L1Grid::fNz, indexes ); -// short_v zBin = ( Z * fStepZInv - fZMinOverStep ).staticCast<short_v>(); -// unsigned short zBin2 = CAMath::Max( short_v( Vc::Zero ), CAMath::Min( short_v( fNz - 1 ), zBin ) ).staticCast<unsigned short>(); - // const float fYMinOverStep( array, &L1Grid::fYMinOverStep, indexes ); // const float fStepYInv( array, &L1Grid::fStepYInv, indexes ); // const unsigned short fNy( array, &L1Grid::fNy, indexes ); // short_v yBin = ( Y * fStepYInv - fYMinOverStep ).staticCast<short_v>(); // unsigned short yBin2 = CAMath::Max( short_v( Vc::Zero ), CAMath::Min( short_v( fNy - 1 ), yBin ) ).staticCast<unsigned short>(); -// return zBin2 * fNy + yBin2; + +// const float fXMinOverStep( array, &L1Grid::fXMinOverStep, indexes ); +// const float fStepXInv( array, &L1Grid::fStepXInv, indexes ); +// const unsigned short fNx( array, &L1Grid::fNx, indexes ); +// short_v xBin = ( X * fStepXInv - fXMinOverStep ).staticCast<short_v>(); +// unsigned short xBin2 = CAMath::Max( short_v( Vc::Zero ), CAMath::Min( short_v( fNx - 1 ), xBin ) ).staticCast<unsigned short>(); +// return yBin2 * fNx + xBin2; // } -// inline void L1Grid::GetBinBounded( const L1Grid *array, const unsigned short &indexes, const float &Y, const float &Z, unsigned short *bY, unsigned short *bZ ) +// inline void L1Grid::GetBinBounded( const L1Grid *array, const unsigned short &indexes, const float &X, const float &Y, unsigned short *bX, unsigned short *bY ) // { +// const float fXMinOverStep( array, &L1Grid::fXMinOverStep, indexes ); +// const float fStepXInv( array, &L1Grid::fStepXInv, indexes ); +// const unsigned short fNx( array, &L1Grid::fNx, indexes ); +// const short fNx2 = fNx.staticCast<short_v>(); +// const short &xBin = static_cast<short>( X * fStepXInv - fXMinOverStep ); +// *bX = static_cast<unsigned short>(max( 0, min( fNx2 - 1, xBin ) )); + // const float fYMinOverStep( array, &L1Grid::fYMinOverStep, indexes ); // const float fStepYInv( array, &L1Grid::fStepYInv, indexes ); // const unsigned short fNy( array, &L1Grid::fNy, indexes ); -// const short fNy2 = fNy.staticCast<short_v>(); -// const short &yBin = static_cast<short>( Y * fStepYInv - fYMinOverStep ); -// *bY = static_cast<unsigned short>(max( 0, min( fNy2 - 1, yBin ) )); - -// const float fZMinOverStep( array, &L1Grid::fZMinOverStep, indexes ); -// const float fStepZInv( array, &L1Grid::fStepZInv, indexes ); -// const unsigned short fNz( array, &L1Grid::fNz, indexes ); -// const short_v fNz2 = fNz.staticCast<short_v>(); -// const short_v &zBin = ( Z * fStepZInv - fZMinOverStep ).staticCast<short_v>(); -// *bZ = CAMath::Max( short_v( Vc::Zero ), CAMath::Min( fNz2 - 1, zBin ) ).staticCast<unsigned short>(); +// const short_v fNy2 = fNy.staticCast<short_v>(); +// const short_v &yBin = ( Y * fStepYInv - fYMinOverStep ).staticCast<short_v>(); +// *bY = CAMath::Max( short_v( Vc::Zero ), CAMath::Min( fNy2 - 1, yBin ) ).staticCast<unsigned short>(); // } -inline int L1Grid::GetBin(float Y, float Z) const +inline int L1Grid::GetBin(float X, float Y) const { //* get the bin pointer + const int& xBin = static_cast<int>(X * fStepXInv - fXMinOverStep); const int& yBin = static_cast<int>(Y * fStepYInv - fYMinOverStep); - const int& zBin = static_cast<int>(Z * fStepZInv - fZMinOverStep); + assert(xBin >= 0); assert(yBin >= 0); - assert(zBin >= 0); + assert(xBin < static_cast<int>(fNx)); assert(yBin < static_cast<int>(fNy)); - assert(zBin < static_cast<int>(fNz)); - const int& bin = zBin * fNy + yBin; + const int& bin = yBin * fNx + xBin; return bin; } -inline int L1Grid::GetBin(float Y, float Z, float T) const +inline int L1Grid::GetBin(float X, float Y, float T) const { //* get the bin pointer + const int& xBin = static_cast<int>(X * fStepXInv - fXMinOverStep); const int& yBin = static_cast<int>(Y * fStepYInv - fYMinOverStep); - const int& zBin = static_cast<int>(Z * fStepZInv - fZMinOverStep); const int& tBin = static_cast<int>(T * fStepTInv - fTMinOverStep); + assert(xBin >= 0); assert(yBin >= 0); - assert(zBin >= 0); assert(tBin >= 0); + assert(xBin < static_cast<int>(fNx)); assert(yBin < static_cast<int>(fNy)); - assert(zBin < static_cast<int>(fNz)); assert(tBin < static_cast<int>(fNt)); - const int& bin = zBin * fNy + yBin + tBin * (fNy * fNz); + const int& bin = yBin * fNx + xBin + tBin * (fNx * fNy); return bin; } -inline void L1Grid::GetBinBounds(unsigned int iBin, float& Ymin, float& Ymax, float& Zmin, float& Zmax) const +inline void L1Grid::GetBinBounds(unsigned int iBin, float& Xmin, float& Xmax, float& Ymin, float& Ymax) const { - int zBin = iBin / fNy; - int yBin = iBin % fNy; + int yBin = iBin / fNx; + int xBin = iBin % fNx; + Xmin = (fXMinOverStep + xBin) / fStepXInv; Ymin = (fYMinOverStep + yBin) / fStepYInv; - Zmin = (fZMinOverStep + zBin) / fStepZInv; + Xmax = Xmin + 1. / fStepXInv; Ymax = Ymin + 1. / fStepYInv; - Zmax = Zmin + 1. / fStepZInv; } -inline void L1Grid::GetBinBounds(unsigned int iBin, float& Ymin, float& Ymax, float& Zmin, float& Zmax, float& Tmin, +inline void L1Grid::GetBinBounds(unsigned int iBin, float& Xmin, float& Xmax, float& Ymin, float& Ymax, float& Tmin, float& Tmax) const { - int zBin = (iBin % (fNy * fNz)) / fNy; - int yBin = (iBin % (fNy * fNz)) % fNy; - int tBin = (iBin / (fNy * fNz)); + int yBin = (iBin % (fNx * fNy)) / fNx; + int xBin = (iBin % (fNx * fNy)) % fNx; + int tBin = (iBin / (fNx * fNy)); + Xmin = (fXMinOverStep + xBin) / fStepXInv; Ymin = (fYMinOverStep + yBin) / fStepYInv; - Zmin = (fZMinOverStep + zBin) / fStepZInv; Tmin = (fTMinOverStep + tBin) / fStepTInv; + Xmax = Xmin + 1. / fStepXInv; Ymax = Ymin + 1. / fStepYInv; - Zmax = Zmin + 1. / fStepZInv; Tmax = Tmin + 1. / fStepTInv; } -inline unsigned int L1Grid::GetBinBounded(const float& Y, const float& Z) const +inline unsigned int L1Grid::GetBinBounded(const float& X, const float& Y) const { //* get the bin pointer at - unsigned short yBin, zBin; - GetBinBounded(Y, Z, yBin, zBin); - return (unsigned int) zBin * (unsigned int) fNy + (unsigned int) yBin; + unsigned short xBin, yBin; + GetBinBounded(X, Y, xBin, yBin); + return (unsigned int) yBin * (unsigned int) fNx + (unsigned int) xBin; } -inline unsigned int L1Grid::GetBinBounded(const float& Y, const float& Z, const float& T) const +inline unsigned int L1Grid::GetBinBounded(const float& X, const float& Y, const float& T) const { //* get the bin pointer at - unsigned short yBin, zBin, tBin; - GetBinBounded(Y, Z, T, yBin, zBin, tBin); - return (unsigned int) zBin * (unsigned int) fNy + (unsigned int) yBin - + tBin * ((unsigned int) fNy * (unsigned int) fNz); + unsigned short xBin, yBin, tBin; + GetBinBounded(X, Y, T, xBin, yBin, tBin); + return (unsigned int) yBin * (unsigned int) fNx + (unsigned int) xBin + + tBin * ((unsigned int) fNx * (unsigned int) fNy); } -inline void L1Grid::GetBinBounded(const float& Y, const float& Z, unsigned short& bY, unsigned short& bZ) const +inline void L1Grid::GetBinBounded(const float& X, const float& Y, unsigned short& bX, unsigned short& bY) const { + const short& xBin = (X * fStepXInv - fXMinOverStep); const short& yBin = (Y * fStepYInv - fYMinOverStep); - const short& zBin = (Z * fStepZInv - fZMinOverStep); + bX = std::max(short(0), std::min(short(fNx - 1), xBin)); bY = std::max(short(0), std::min(short(fNy - 1), yBin)); - bZ = std::max(short(0), std::min(short(fNz - 1), zBin)); } -inline void L1Grid::GetBinBounded(const float& Y, const float& Z, const float& T, unsigned short& bY, - unsigned short& bZ, unsigned short& bT) const +inline void L1Grid::GetBinBounded(const float& X, const float& Y, const float& T, unsigned short& bX, + unsigned short& bY, unsigned short& bT) const { + const short& xBin = (X * fStepXInv - fXMinOverStep); const short& yBin = (Y * fStepYInv - fYMinOverStep); - const short& zBin = (Z * fStepZInv - fZMinOverStep); const short& tBin = (T * fStepTInv - fTMinOverStep); // cout<<fStepTInv<<" fStepTInv "<<fTMinOverStep<<" fTMinOverStep "<<T<<" T "<<endl; + bX = std::max(short(0), std::min(short(fNx - 1), xBin)); bY = std::max(short(0), std::min(short(fNy - 1), yBin)); - bZ = std::max(short(0), std::min(short(fNz - 1), zBin)); bT = std::max(short(0), std::min(short(fNt - 1), tBin)); // cout<<(fNt - 1)<<" (fNt - 1) "<<tBin<<" tBin "<<bT<<" bT "<<endl; diff --git a/reco/L1/L1Algo/L1HitArea.h b/reco/L1/L1Algo/L1HitArea.h index 0fcd83d699..64c40229c0 100644 --- a/reco/L1/L1Algo/L1HitArea.h +++ b/reco/L1/L1Algo/L1HitArea.h @@ -13,7 +13,7 @@ class L1SliceData; class L1HitArea { public: - L1HitArea(const L1Grid& grid, float y, float z, float dy, float dz); + L1HitArea(const L1Grid& grid, float x, float y, float dx, float dy); /** * look up the next hit in the requested area. @@ -24,82 +24,82 @@ public: protected: const L1Grid& fGrid; - unsigned short fBZmax; // maximal Z bin index - unsigned short fBDY; // Y distance of bin indexes - unsigned int fIndYmin; // minimum index for - unsigned short fIz; // current Z bin index (incremented while iterating) - L1HitIndex_t fHitYlst; // last possible hit index in current z-line + 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) + L1HitIndex_t fHitXlst; // last possible hit index in current y-line L1HitIndex_t fIh; // hit index iterating inside the bins - int fNy; // Number of bins in Y direction + int fNx; // Number of bins in X direction }; -inline L1HitArea::L1HitArea(const L1Grid& grid, float y, float z, float dy, float dz) +inline L1HitArea::L1HitArea(const L1Grid& grid, float x, float y, float dx, float dy) : fGrid(grid) - , fBZmax(0) - , fBDY(0) - , fIndYmin(0) - , fIz(0) - , fHitYlst(0) + , fBYmax(0) + , fBDX(0) + , fIndXmin(0) + , fIy(0) + , fHitXlst(0) , fIh(0) - , fNy(fGrid.Ny()) + , fNx(fGrid.Nx()) { + const float minX = x - dx; + const float maxX = x + dx; const float minY = y - dy; const float maxY = y + dy; - const float minZ = z - dz; - const float maxZ = z + dz; - unsigned short bYmin, bZmin, bYmax; // boundary bin indexes - fGrid.GetBinBounded(minY, minZ, bYmin, bZmin); - fGrid.GetBinBounded(maxY, maxZ, bYmax, fBZmax); + unsigned short bXmin, bYmin, bXmax; // boundary bin indexes + fGrid.GetBinBounded(minX, minY, bXmin, bYmin); + fGrid.GetBinBounded(maxX, maxY, bXmax, fBYmax); - fBDY = (bYmax - bYmin + 1); // bin index span in y direction + fBDX = (bXmax - bXmin + 1); // bin index span in x direction - fIndYmin = (bZmin * fNy + bYmin); // same as grid.GetBin(fMinY, fMinZ), i.e. the smallest bin index of interest - // fIndYmin + fBDY - 1 then is the largest bin index of interest with the same Z + 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(fIndYmin, y, dy, z, dz); + fGrid.GetBinBounds(fIndXmin, x, dx, y, dy); - fIz = bZmin; + fIy = bYmin; - fIh = fGrid.FirstHitInBin(fIndYmin); - fHitYlst = fGrid.FirstHitInBin(fIndYmin + fBDY); + fIh = fGrid.FirstHitInBin(fIndXmin); + fHitXlst = fGrid.FirstHitInBin(fIndXmin + fBDX); } inline bool L1HitArea::GetNext(L1HitIndex_t& i) { - bool yIndexOutOfRange = fIh >= fHitYlst; // current y is not in the area - bool nextZIndexOutOfRange = (fIz >= fBZmax); // there isn't any new z-line + bool xIndexOutOfRange = fIh >= fHitXlst; // current x is not in the area + bool nextYIndexOutOfRange = (fIy >= fBYmax); // there isn't any new y-line - if (yIndexOutOfRange && nextZIndexOutOfRange) { // all iterators are over the end + if (xIndexOutOfRange && nextYIndexOutOfRange) { // all iterators are over the end return false; } - // at least one entry in the vector has (fIh >= fHitYlst && fIz < fBZmax) - bool needNextZ = yIndexOutOfRange && !nextZIndexOutOfRange; + // 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 y-index - while (ISLIKELY(needNextZ)) { //ISLIKELY to speed the programm and optimise the use of registers - fIz++; // get new z-line + // 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 - fIndYmin += fNy; - fIh = fGrid.FirstHitInBin(fIndYmin); // get first hit in cell, if z-line is new - fHitYlst = fGrid.FirstHitInBin(fIndYmin + fBDY); + fIndXmin += fNx; + fIh = fGrid.FirstHitInBin(fIndXmin); // get first hit in cell, if y-line is new + fHitXlst = fGrid.FirstHitInBin(fIndXmin + fBDX); - yIndexOutOfRange = fIh >= fHitYlst; - nextZIndexOutOfRange = (fIz >= fBZmax); - needNextZ = yIndexOutOfRange && !nextZIndexOutOfRange; + xIndexOutOfRange = fIh >= fHitXlst; + nextYIndexOutOfRange = (fIy >= fBYmax); + needNextY = xIndexOutOfRange && !nextYIndexOutOfRange; } - L1_ASSERT(fIh < fGrid.FirstHitInBin(fGrid.N()) || yIndexOutOfRange, - fIh << " < " << fGrid.FirstHitInBin(fGrid.N()) << " || " << yIndexOutOfRange); + L1_ASSERT(fIh < fGrid.FirstHitInBin(fGrid.N()) || xIndexOutOfRange, + fIh << " < " << fGrid.FirstHitInBin(fGrid.N()) << " || " << xIndexOutOfRange); i = fIh; // return fIh++; // go to next - return !yIndexOutOfRange; + return !xIndexOutOfRange; } class L1HitAreaTime { public: - L1HitAreaTime(const L1Grid& grid, float y, float z, float dy, float dz, float t, float dt); + L1HitAreaTime(const L1Grid& grid, float x, float y, float dx, float dy, float t, float dt); /** * look up the next hit in the requested area. @@ -110,122 +110,122 @@ public: protected: const L1Grid& fGrid; - unsigned short fBZmax; // maximal Z bin index - unsigned short fBDY; // Y distance of bin indexes - unsigned int fIndYmin; // minimum index for - unsigned short fIz; // current Z bin index (incremented while iterating) - L1HitIndex_t fHitYlst; // last possible hit index in current z-line + 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) + L1HitIndex_t fHitXlst; // last possible hit index in current y-line L1HitIndex_t fIh; // hit index iterating inside the bins - int fNy; // Number of bins in Y direction - int fNz; + int fNx; // Number of bins in X direction + int fNy; - unsigned short fBTmax; // maximal Z bin index - unsigned short fIt; // current Z bin index (incremented while iterating) + unsigned short fBTmax; // maximal Y bin index + unsigned short fIt; // current Y bin index (incremented while iterating) - unsigned short fBZmin; unsigned short fBYmin; + unsigned short fBXmin; - unsigned int fIndZmin; // minimum index for + unsigned int fIndYmin; // minimum index for }; -inline L1HitAreaTime::L1HitAreaTime(const L1Grid& grid, float y, float z, float dy, float dz, float t, float dt) +inline L1HitAreaTime::L1HitAreaTime(const L1Grid& grid, float x, float y, float dx, float dy, float t, float dt) : fGrid(grid) - , fBZmax(0) - , fBDY(0) - , fIndYmin(0) - , fIz(0) - , fHitYlst(0) + , fBYmax(0) + , fBDX(0) + , fIndXmin(0) + , fIy(0) + , fHitXlst(0) , fIh(0) + , fNx(fGrid.Nx()) , fNy(fGrid.Ny()) - , fNz(fGrid.Nz()) , fBTmax(0) , fIt(0) - , fBZmin(0) , fBYmin(0) - , fIndZmin(0) + , fBXmin(0) + , fIndYmin(0) { + const float minX = x - dx; + const float maxX = x + dx; const float minY = y - dy; const float maxY = y + dy; - const float minZ = z - dz; - const float maxZ = z + dz; const float minT = t - dt; const float maxT = t + dt; - unsigned short bYmin, bZmin, bYmax, bTmin; // boundary bin indexes - fGrid.GetBinBounded(minY, minZ, minT, bYmin, bZmin, bTmin); - fGrid.GetBinBounded(maxY, maxZ, maxT, bYmax, fBZmax, fBTmax); - fBZmin = bZmin; + unsigned short bXmin, bYmin, bXmax, bTmin; // boundary bin indexes + fGrid.GetBinBounded(minX, minY, minT, bXmin, bYmin, bTmin); + fGrid.GetBinBounded(maxX, maxY, maxT, bXmax, fBYmax, fBTmax); fBYmin = bYmin; + fBXmin = bXmin; - fBDY = (bYmax - bYmin + 1); // bin index span in y direction + fBDX = (bXmax - bXmin + 1); // bin index span in x direction - fIndYmin = (bTmin * fNy * fNz + bZmin * fNy + bYmin); + fIndXmin = (bTmin * fNx * fNy + bYmin * fNx + bXmin); - //fGrid.GetBinBounds(fIndYmin, y, dy, z, dz, t, dt); + //fGrid.GetBinBounds(fIndXmin, x, dx, y, dy, t, dt); - fIz = bZmin; + fIy = bYmin; fIt = bTmin; - fIh = fGrid.FirstHitInBin(fIndYmin); + fIh = fGrid.FirstHitInBin(fIndXmin); - fHitYlst = fGrid.FirstHitInBin(fIndYmin + fBDY); + fHitXlst = fGrid.FirstHitInBin(fIndXmin + fBDX); } inline bool L1HitAreaTime::GetNext(L1HitIndex_t& i) { - bool yIndexOutOfRange = fIh >= fHitYlst; // current y is not in the area - bool nextZIndexOutOfRange = (fIz >= fBZmax); // there isn't any new z-line - bool nextTIndexOutOfRange = (fIt >= fBTmax); // there isn't any new z-line + bool xIndexOutOfRange = fIh >= fHitXlst; // current x is not in the area + bool nextYIndexOutOfRange = (fIy >= fBYmax); // there isn't any new y-line + bool nextTIndexOutOfRange = (fIt >= fBTmax); // there isn't any new y-line - if (yIndexOutOfRange && nextZIndexOutOfRange && nextTIndexOutOfRange) { // all iterators are over the end + if (xIndexOutOfRange && nextYIndexOutOfRange && nextTIndexOutOfRange) { // all iterators are over the end return false; } - // at least one entry in the vector has (fIh >= fHitYlst && fIz < fBZmax) - bool needNextZ = yIndexOutOfRange && !nextZIndexOutOfRange; - bool needNextT = yIndexOutOfRange && nextZIndexOutOfRange && !nextTIndexOutOfRange; + // at least one entry in the vector has (fIh >= fHitXlst && fIy < fBYmax) + bool needNextY = xIndexOutOfRange && !nextYIndexOutOfRange; + bool needNextT = xIndexOutOfRange && nextYIndexOutOfRange && !nextTIndexOutOfRange; - while (ISLIKELY((needNextZ) || (needNextT))) { //ISLIKELY to speed the programm and optimise the use of registers + while (ISLIKELY((needNextY) || (needNextT))) { //ISLIKELY to speed the programm and optimise the use of registers if (needNextT) { fIt++; - fIz = fBZmin; + fIy = fBYmin; - fIndYmin = (fIt * fNy * fNz + fBZmin * fNy + fBYmin); - fIh = fGrid.FirstHitInBin(fIndYmin); // get first hit in cell, if z-line is new - fHitYlst = fGrid.FirstHitInBin(fIndYmin + fBDY); + fIndXmin = (fIt * fNx * fNy + fBYmin * fNx + fBXmin); + fIh = fGrid.FirstHitInBin(fIndXmin); // get first hit in cell, if y-line is new + fHitXlst = fGrid.FirstHitInBin(fIndXmin + fBDX); } else { - fIz++; // get new z-line + fIy++; // get new y-line // get next hit - fIndYmin += fNy; - fIh = fGrid.FirstHitInBin(fIndYmin); // get first hit in cell, if z-line is new + fIndXmin += fNx; + fIh = fGrid.FirstHitInBin(fIndXmin); // get first hit in cell, if y-line is new - fHitYlst = fGrid.FirstHitInBin(fIndYmin + fBDY); + fHitXlst = fGrid.FirstHitInBin(fIndXmin + fBDX); } - yIndexOutOfRange = fIh >= fHitYlst; - nextZIndexOutOfRange = (fIz >= fBZmax); - needNextZ = yIndexOutOfRange && !nextZIndexOutOfRange; + xIndexOutOfRange = fIh >= fHitXlst; + nextYIndexOutOfRange = (fIy >= fBYmax); + needNextY = xIndexOutOfRange && !nextYIndexOutOfRange; nextTIndexOutOfRange = (fIt >= fBTmax); - needNextT = yIndexOutOfRange && nextZIndexOutOfRange && !nextTIndexOutOfRange; + needNextT = xIndexOutOfRange && nextYIndexOutOfRange && !nextTIndexOutOfRange; } - L1_ASSERT(fIh < fGrid.FirstHitInBin(fGrid.N() + 1) || yIndexOutOfRange, - fIh << " < " << fGrid.FirstHitInBin(fGrid.N() + 1) << " || " << yIndexOutOfRange); + L1_ASSERT(fIh < fGrid.FirstHitInBin(fGrid.N() + 1) || xIndexOutOfRange, + fIh << " < " << fGrid.FirstHitInBin(fGrid.N() + 1) << " || " << xIndexOutOfRange); i = fIh; // return fIh++; // go to next - return !yIndexOutOfRange; + return !xIndexOutOfRange; } #endif -- GitLab