diff --git a/algo/ca/core/data/CaGridEntry.h b/algo/ca/core/data/CaGridEntry.h
index cd44ddc65ad26ce0c59a9af2a8ecf1f190633d93..9a9b50d3bbe8fbc0b263c5a436a47d194d009cfb 100644
--- a/algo/ca/core/data/CaGridEntry.h
+++ b/algo/ca/core/data/CaGridEntry.h
@@ -18,59 +18,42 @@ namespace cbm::algo::ca
GridEntry() = default;
- GridEntry(const ca::Hit& hit) { Set(hit); };
+ GridEntry(const ca::Hit& hit, ca::HitIndex_t id) { Set(hit, id); };
- void Set(const ca::Hit& hit)
+ void Set(const ca::Hit& hit, ca::HitIndex_t id)
{
- x = hit.x;
- y = hit.y;
- z = hit.z;
- t = hit.t;
- dx2 = hit.dx2;
- dxy = hit.dxy;
- dy2 = hit.dy2;
- dt2 = hit.dt2;
- rangeX = hit.rangeX;
- rangeY = hit.rangeY;
- rangeT = hit.rangeT;
-
- fIsSuppressed = 0;
+ fWindowHitId = 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; }
+
fscal Z() const { return z; }
fscal T() const { return t; }
fscal X() const { return x; }
fscal Y() const { return y; }
- fscal dX2() const { return dx2; }
- fscal dXY() const { return dxy; }
- fscal dY2() const { return dy2; }
- fscal dT2() const { return dt2; }
fscal RangeX() const { return rangeX; }
fscal RangeY() const { return rangeY; }
fscal RangeT() const { return rangeT; }
- bool IsSuppressed() const { return fIsSuppressed; }
-
- void SetIsSuppresed(bool val) { fIsSuppressed = val; }
-
private:
constexpr static fscal kUndef = 0.; // TODO: test with constants::Undef<fscal>
+ ca::HitIndex_t fWindowHitId {0}; ///< hit id in L1Algo::fWindowHits array
+
fscal x {kUndef}; // x coordinate of the hit
fscal y {kUndef}; // y coordinate of the hit
fscal z {kUndef}; // z coordinate of the hit
fscal t {kUndef}; // t coordinate of the hit
- fscal dx2 {kUndef}; // x error^2 of the hit
- fscal dxy {kUndef}; // xy covariance of the hit
- fscal dy2 {kUndef}; // y error^2 of the hit
- fscal dt2 {kUndef}; // t error^2 of the hit
fscal rangeX {kUndef}; // x range of the hit
fscal rangeY {kUndef}; // y range of the hit
fscal rangeT {kUndef}; // t range of the hit
-
- /// fIsSuppressed flag is used to suppress duplicated hits at the module overlaps
- // TODO: collect those hits on the track instead of suppressing them
- bool fIsSuppressed {0}; // is the hit suppressed by another hit
};
} // namespace cbm::algo::ca
diff --git a/algo/ca/core/data/CaHit.h b/algo/ca/core/data/CaHit.h
index bd6a44cae062aef23d49657171001e3ba8c031d2..ff7fb114ca61db9526320285c745ffd9ec207888 100644
--- a/algo/ca/core/data/CaHit.h
+++ b/algo/ca/core/data/CaHit.h
@@ -15,20 +15,37 @@
namespace cbm::algo::ca
{
- using HitIndex_t = unsigned int; ///< Index of ca::Hit
- using StripIndex_t = unsigned int; ///< Index of the station strip
-
+ using HitIndex_t = unsigned int; ///< Index of ca::Hit
+ using HitKeyIndex_t = unsigned int; ///< Index of the hit key (e.g. front / back cluster id for STS)
/// \brief ca::Hit class describes a generic hit for the CA tracker
///
struct Hit {
friend class boost::serialization::access;
+ public:
+ Hit() = default;
+
+ fscal X() const { return x; }
+ fscal Y() const { return y; }
+ fscal Z() const { return z; }
+ fscal T() const { return t; }
+ fscal dX2() const { return dx2; }
+ fscal dY2() const { return dy2; }
+ fscal dXY() const { return dxy; }
+ fscal dT2() const { return dt2; }
+ fscal RangeX() const { return rangeX; }
+ fscal RangeY() const { return rangeY; }
+ fscal RangeT() const { return rangeT; }
+ HitIndex_t Id() const { return fId; }
+ int Ist() const { return iSt; }
+
+ public:
/// NOTE: For STS f and b correspond to the indexes of the front and back clusters of the hit in a dataset. For other
/// tracking detectors (MVD, MuCh, TRD, TOF) f == b and corresponds to the index of the hit. Indexes f and b
/// do not intersect between different detector stations.
- StripIndex_t f {0}; ///< front hit key index
- StripIndex_t b {0}; ///< back hit key index
+ HitKeyIndex_t f {0}; ///< front hit key index
+ HitKeyIndex_t b {0}; ///< back hit key index
fscal x {0.}; ///< measured X coordinate [cm]
fscal y {0.}; ///< measured Y coordinate [cm]
@@ -42,7 +59,7 @@ namespace cbm::algo::ca
fscal rangeY {0.}; ///< +/- range of uncertainty of Y coordinate [cm]
fscal rangeT {0.}; ///< +/- range of uncertainty of time [ns]
- HitIndex_t Id {0}; ///< id of the hit
+ HitIndex_t fId {0}; ///< id of the hit
int iSt {-1}; ///< index of station in the active stations array
private:
@@ -63,7 +80,7 @@ namespace cbm::algo::ca
ar& rangeX;
ar& rangeY;
ar& rangeT;
- ar& Id;
+ ar& fId;
ar& iSt;
}
};
diff --git a/reco/L1/CbmCaTimeSliceReader.cxx b/reco/L1/CbmCaTimeSliceReader.cxx
index cc76e90f96b05292d9cda63d01fee9b15be187dd..6d12d01060c46879dcdd6e16596663043169ad69 100644
--- a/reco/L1/CbmCaTimeSliceReader.cxx
+++ b/reco/L1/CbmCaTimeSliceReader.cxx
@@ -417,7 +417,7 @@ void TimeSliceReader::StoreHitRecord(const HitRecord& hitRecord)
aHit.rangeY = hitRecord.fRangeY;
aHit.rangeT = hitRecord.fRangeT;
- aHit.Id = static_cast<int>(fpIODataManager->GetNofHits());
+ aHit.fId = static_cast<int>(fpIODataManager->GetNofHits());
aHit.iSt = hitRecord.fStaId;
fpIODataManager->PushBackHit(aHit, hitRecord.fDataStream);
}
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 880fa8405ecb31f7344aa1b92b7659979f8e3ef6..62bb54989ae31113e105f5b9bdb5ebe6bd90e833 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -230,7 +230,8 @@ InitStatus CbmL1::Init()
fInitManager.SetMaxTripletPerDoublets(1000);
}
-
+ // uncomment for debug
+ //fInitManager.DevSetIgnoreHitSearchAreas(true);
//fInitManager.DevSetIsMatchDoubletsViaMc(true);
//fInitManager.DevSetIsMatchTripletsViaMc(true);
//fInitManager.DevSetIsExtendTracksViaMc(true);
@@ -801,7 +802,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
for (int i = 0; i < it->fNofHits; i++) {
int caHitId = fpAlgo->fRecoHits[trackFirstHit + i];
- int cbmHitID = fpAlgo->GetInputData().GetHit(caHitId).Id;
+ int cbmHitID = fpAlgo->GetInputData().GetHit(caHitId).Id();
t.Hits.push_back(cbmHitID);
}
fvRecoTracks.push_back(t);
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 82f19b74009ab33d1ff50c51c8c1d104c8c672c0..57d90763f40607fdb0590801902756622889d5c7 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -763,7 +763,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
h.rangeY = th.rangeY;
h.rangeT = th.rangeT;
- h.Id = iHit;
+ h.fId = iHit;
h.iSt = th.iStation;
// save hit
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 12d4b53fb315f3597b54f09ed3c930fd184aec81..890470d89ea361003d1c76fe0f4a94b22c4b587b 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -13,7 +13,6 @@
L1Algo::L1Algo()
{
for (unsigned int i = 0; i < constants::size::MaxNstations; i++) {
- vGrid[i].AllocateMemory();
fTriplets[i].SetName(std::stringstream() << "L1Algo::fTriplets[" << i << "]");
}
}
@@ -33,43 +32,6 @@ void L1Algo::ReceiveInputData(InputData&& inputData)
fInputData = std::move(inputData);
}
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void L1Algo::ResetSliceData()
-{
- int nHits = 0;
- for (int iS = 0; iS < fParameters.GetNstationsActive(); iS++) {
- nHits += fSliceHitIds[iS].size();
- }
-
- fGridHits.reset(nHits);
- fGridHitsBuf.reset(nHits);
- fGridHitIds.reset(nHits);
- fGridHitIdsBuf.reset(nHits);
- fGridPoints.reset(nHits);
- fGridPointsBuf.reset(nHits);
-
- fHitFirstTriplet.reset(nHits);
- fHitNtriplets.reset(nHits);
-
- // TODO: SG: reduce the array size to the number of keys in subslice
-
- fStripToTrack.reset(fInputData.GetNhitKeys(), -1);
-
- fSliceRecoTracks.clear();
- fSliceRecoTracks.reserve(2 * nHits / fParameters.GetNstationsActive());
-
- fSliceRecoHits.clear();
- fSliceRecoHits.reserve(2 * nHits);
-
- fTrackCandidates.clear();
- fTrackCandidates.reserve(nHits / 10);
- for (unsigned int iS = 0; iS < constants::size::MaxNstations; iS++) {
- int nHitsStation = fSliceHitIds[iS].size();
- fTriplets[iS].clear();
- fTriplets[iS].reserve(2 * nHitsStation);
- }
-}
// ---------------------------------------------------------------------------------------------------------------------
@@ -102,17 +64,17 @@ int L1Algo::GetMcTrackIdForCaHit(int iHit) const
return CbmL1::Instance()->GetMcPoints()[iMcPoint].ID;
}
-int L1Algo::GetMcTrackIdForGridHit(int iHit) const
+int L1Algo::GetMcTrackIdForWindowHit(int iHit) const
{
- int hitId = fGridHitIds[iHit];
+ int hitId = fWindowHits[iHit].Id();
int iMcPoint = CbmL1::Instance()->GetHitBestMcRefs()[hitId];
if (iMcPoint < 0) return -1;
return CbmL1::Instance()->GetMcPoints()[iMcPoint].ID;
}
-const CbmL1MCTrack* L1Algo::GetMcTrackForGridHit(int iHit) const
+const CbmL1MCTrack* L1Algo::GetMcTrackForWindowHit(int iHit) const
{
- int id = GetMcTrackIdForGridHit(iHit);
+ int id = GetMcTrackIdForWindowHit(iHit);
if (id < 0) return nullptr;
return &CbmL1::Instance()->GetMcTracks()[id];
}
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index d4dc73a350666eb435d27a1c697909e3dd86c5af..ea83358281f6c02d817c515be6a3141194628f7d 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -182,7 +182,7 @@ public:
/// ----- Subroutines used by L1Algo::CaTrackFinder() ------
- void ResetSliceData();
+ void ReadWindowData();
void CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best_L, fscal& best_chi2, const L1Triplet* curr_trip,
L1Branch& curr_tr, unsigned char& curr_L, fscal& curr_chi2, unsigned char min_best_l,
@@ -261,9 +261,9 @@ public:
int GetMcTrackIdForCaHit(int iHit) const;
/// Get mc track ID for a hit (debug tool)
- int GetMcTrackIdForGridHit(int iGridHit) const;
+ int GetMcTrackIdForWindowHit(int iGridHit) const;
- const CbmL1MCTrack* GetMcTrackForGridHit(int iHit) const;
+ const CbmL1MCTrack* GetMcTrackForWindowHit(int iHit) const;
private:
bool checkTripletMatch(const L1Triplet& l, const L1Triplet& r, fscal& dchi2) const;
@@ -288,18 +288,36 @@ private:
public:
Vector<L1HitTimeInfo> fHitTimeInfo;
- L1Grid vGrid[constants::size::MaxNstations]; ///<
+ L1Grid vGrid[constants::size::MaxNstations]; ///<
fscal fMaxRangeX[constants::size::MaxNstations];
fscal fMaxRangeY[constants::size::MaxNstations];
fscal fMaxRangeT[constants::size::MaxNstations];
+ /// ----- Data used during track finding -----
+ ///
+
+ ///\brief hits of the current time window
+ /// it is a portion of fInputData to process in the current time window
+ /// hit.Id is replaced by the hit index in fInputData
+ Vector<ca::Hit> fWindowHits {"L1Algo::fWindowHits"};
+
+ Vector<unsigned char> fIsWindowHitSuppressed {"L1Algo::fIsWindowHitSuppressed"};
+
+ ///\brief first index of the station hits in the time window
+ ca::HitIndex_t fStationHitsStartIndex[constants::size::MaxNstations + 1] {0};
+
+ ///\brief number of station hits in the time window
+ ca::HitIndex_t fStationNhits[constants::size::MaxNstations + 1] {0};
+
+ /// ----------
+
double fCaRecoTime {0.}; // time of the track finder + fitter
- Vector<Track> fRecoTracks {"L1Algo::fRecoTracks"}; ///< reconstructed tracks
+ Vector<Track> fRecoTracks {"L1Algo::fRecoTracks"}; ///< reconstructed tracks
Vector<ca::HitIndex_t> fRecoHits {"L1Algo::fRecoHits"}; ///< packed hits of reconstructed tracks
- Vector<Track> fSliceRecoTracks {"L1Algo::fSliceRecoTracks"}; ///< reconstructed tracks in sub-timeslice
+ Vector<Track> fSliceRecoTracks {"L1Algo::fSliceRecoTracks"}; ///< reconstructed tracks in sub-timeslice
Vector<ca::HitIndex_t> fSliceRecoHits {"L1Algo::fSliceRecoHits"}; ///< packed hits of reconstructed tracks
/// Created triplets vs station index
@@ -312,18 +330,6 @@ public:
///< indices of the sub-slice hits
Vector<ca::HitIndex_t> fSliceHitIds[constants::size::MaxNstations] {"L1Algo::fSliceHitIds"};
- Vector<ca::Hit> fGridHits {"L1Algo::fGridHits"}; ///< hits, ordered with respect to their grid bins
- Vector<ca::Hit> fGridHitsBuf {"L1Algo::fGridHitsBuf"}; ///< hits, ordered with respect to their grid bins
-
- Vector<ca::HitIndex_t> fGridHitIds {"L1Algo::fGridHitIds"}; ///< indices of grid hits: iGridHit -> iCaHit
- Vector<ca::HitIndex_t> fGridHitIdsBuf {"L1Algo::fGridHitIdsBuf"}; ///< buffer for a new fGridHitIds
-
- Vector<ca::GridEntry> fGridPoints {"L1Algo::fGridPoints"}; ///< grid points parallel to fGridHits
- Vector<ca::GridEntry> fGridPointsBuf {"L1Algo::fGridPointsBuf"};
-
- ca::HitIndex_t fGridHitStartIndex[constants::size::MaxNstations + 1] {0};
- ca::HitIndex_t fGridHitStopIndex[constants::size::MaxNstations + 1] {0};
-
Vector<int> fStripToTrack {"L1Algo::fStripToTrack"}; // strip to track pointers
TrackingMode fTrackingMode {kSts};
@@ -332,8 +338,8 @@ public:
fvec Err {0.f};
/// --- data used during finding iterations
- int isec {0}; // iteration TODO: to be dispatched (S.Zharko, 21.06.2022)
- const Iteration* fpCurrentIteration = nullptr; ///< pointer to the current CA track finder iteration
+ unsigned int fCurrentIterationIndex {0}; ///< index of the corrent iteration, needed for debug output
+ const Iteration* fpCurrentIteration = nullptr; ///< pointer to the current CA track finder iteration
Vector<int> fHitFirstTriplet {"L1Algo::fHitFirstTriplet"}; /// link hit -> first triplet { hit, *, *}
Vector<int> fHitNtriplets {"L1Algo::fHitNtriplets"}; /// link hit ->n triplets { hit, *, *}
@@ -381,7 +387,7 @@ private:
ca::FieldValue fTargB _fvecalignment {}; // field in the target point (modifiable, do not touch!!)
L1XYMeasurementInfo TargetXYInfo _fvecalignment {}; // target constraint [cm]
- int fGhostSuppression {1}; // NOTE: Should be equal to 0 in TRACKS_FROM_TRIPLETS mode!
+ int fGhostSuppression {1}; // NOTE: Should be equal to 0 in TRACKS_FROM_TRIPLETS mode!
/// ----- Debug features -----
#ifdef PULLS
diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/reco/L1/L1Algo/L1BranchExtender.cxx
index 0216c7f96fb25dfa67652880084231e9ce07b443..7f516a8190067cecf8eceadf0d31f7ed33dd3c66 100644
--- a/reco/L1/L1Algo/L1BranchExtender.cxx
+++ b/reco/L1/L1Algo/L1BranchExtender.cxx
@@ -208,18 +208,16 @@ void L1Algo::FindMoreHits(L1Branch& t, TrackParamV& Tout, const bool upstream, c
if (fParameters.DevIsIgnoreHitSearchAreas()) {
ih++;
- if ((ca::HitIndex_t) ih >= (fGridHitStopIndex[ista] - fGridHitStartIndex[ista])) { break; }
+ if ((ca::HitIndex_t) ih >= vGrid[ista].Entries().size()) { break; }
}
else {
if (!area.GetNext(ih)) { break; }
}
- ca::HitIndex_t globalInd = fGridHitStartIndex[ista] + ih;
+ const ca::GridEntry& gridEntry = vGrid[ista].Entries()[ih];
- const ca::GridEntry& hitPoint = fGridPoints[globalInd];
- if (hitPoint.IsSuppressed()) { continue; }
-
- const ca::Hit& hit = fGridHits[globalInd];
+ if (fIsWindowHitSuppressed[gridEntry.WindowHitId()]) { continue; }
+ const ca::Hit& hit = fWindowHits[gridEntry.WindowHitId()];
if (sta.timeInfo && tr.NdfTime()[0] > -2.) {
fscal dt = hit.t - tr.Time()[0];
@@ -251,14 +249,14 @@ void L1Algo::FindMoreHits(L1Branch& t, TrackParamV& Tout, const bool upstream, c
if (fabs(d_x) > dxm_est) continue;
r2_best = d2;
- iHit_best = globalInd;
+ iHit_best = gridEntry.WindowHitId();
}
if (iHit_best < 0) break;
- newHits.push_back(fGridHitIds[iHit_best]);
+ newHits.push_back(fWindowHits[iHit_best].Id());
- const ca::Hit& hit = fGridHits[iHit_best];
+ const ca::Hit& hit = fWindowHits[iHit_best];
fit.Extrapolate(hit.z, fld);
fit.FilterHit(sta, hit);
diff --git a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
index 0225b41d9d4a256e0d36d96a37bd785f7d39a9b3..4c9f6060b540ce383bc251d3705569dfccf49b8c 100644
--- a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
+++ b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
@@ -119,6 +119,52 @@ bool L1Algo::checkTripletMatch(const L1Triplet& l, const L1Triplet& r, fscal& dc
// * *
// **************************************************************************************************
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void L1Algo::ReadWindowData()
+{
+ int nHits = 0;
+ for (int iS = 0; iS < fParameters.GetNstationsActive(); iS++) {
+ fStationHitsStartIndex[iS] = nHits;
+ fStationNhits[iS] = fSliceHitIds[iS].size();
+ nHits += fStationNhits[iS];
+ }
+ fStationHitsStartIndex[fParameters.GetNstationsActive()] = nHits;
+
+ fWindowHits.reset(nHits);
+ fIsWindowHitSuppressed.reset(nHits, 0);
+
+ for (int iS = 0; iS < fParameters.GetNstationsActive(); iS++) {
+ for (ca::HitIndex_t ih = 0; ih < fSliceHitIds[iS].size(); ++ih) {
+ ca::Hit h = fInputData.GetHit(fSliceHitIds[iS][ih]);
+ h.fId = fSliceHitIds[iS][ih];
+
+ fWindowHits[fStationHitsStartIndex[iS] + ih] = h;
+ }
+ }
+
+ fHitFirstTriplet.reset(nHits, 0);
+ fHitNtriplets.reset(nHits, 0);
+
+ // TODO: SG: reduce the array size to the number of keys in subslice
+
+ fStripToTrack.reset(fInputData.GetNhitKeys(), -1);
+
+ fSliceRecoTracks.clear();
+ fSliceRecoTracks.reserve(2 * nHits / fParameters.GetNstationsActive());
+
+ fSliceRecoHits.clear();
+ fSliceRecoHits.reserve(2 * nHits);
+
+ fTrackCandidates.clear();
+ fTrackCandidates.reserve(nHits / 10);
+ for (unsigned int iS = 0; iS < constants::size::MaxNstations; iS++) {
+ int nHitsStation = fSliceHitIds[iS].size();
+ fTriplets[iS].clear();
+ fTriplets[iS].reserve(2 * nHitsStation);
+ }
+}
+
void L1Algo::CaTrackFinderSlice()
{
fNFindIterations = fParameters.GetNcaIterations();
@@ -152,12 +198,10 @@ void L1Algo::CaTrackFinderSlice()
#endif
-#if defined(XXX) || defined(COUNTERS)
+#if defined(COUNTERS)
static unsigned int stat_N = 0; // number of events
stat_N++;
-#endif
-#ifdef COUNTERS
static Tindex stat_nStartHits = 0;
static Tindex stat_nHits[fNFindIterations] = {0};
@@ -178,10 +222,8 @@ void L1Algo::CaTrackFinderSlice()
* CATrackFinder routine setting
***********************************/
- ResetSliceData();
-
+ ReadWindowData();
- ca::HitIndex_t nGridHitsFilled = 0;
for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
fMaxRangeX[iS] = 0.;
@@ -219,16 +261,6 @@ void L1Algo::CaTrackFinderSlice()
// TODO: changing the grid also changes the result. Investigate why it happens.
// TODO: find the optimal grid size
- /* SG: old code from Valentina
- int nSliceHits = fSliceHitIds.size();
- int nSliceHits = fSliceHitIdsStopIndex[iS] - fSliceHitIdsStartIndex[iS];
- 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 = fSliceHitIds[iS].size();
fscal sizeY = gridMaxY - gridMinY;
fscal sizeX = gridMaxX - gridMinX;
@@ -244,62 +276,41 @@ void L1Algo::CaTrackFinderSlice()
if (xStep < 0.01 * scale) xStep = 0.01 * scale;
vGrid[iS].BuildBins(gridMinX, gridMaxX, gridMinY, gridMaxY, xStep, yStep);
- vGrid[iS].StoreHits(*this, iS, nGridHitsFilled);
+ vGrid[iS].StoreHits(fWindowHits, fStationHitsStartIndex[iS], fStationNhits[iS], fvHitKeyFlags);
}
- fGridHits.reduce(nGridHitsFilled);
- fGridHitsBuf.reduce(nGridHitsFilled);
- fGridHitIds.reduce(nGridHitsFilled);
- fGridHitIdsBuf.reduce(nGridHitsFilled);
- fGridPoints.reduce(nGridHitsFilled);
- fGridPointsBuf.reduce(nGridHitsFilled);
+ // ---- Loop over Track Finder iterations ----------------------------------------------------------------//
- for (unsigned int ih = 0; ih < fGridHits.size(); ++ih) {
- fGridPoints[ih].Set(fGridHits[ih]);
- }
- /********************************/ /**
- * Loop over tracking iterations
- ***********************************/
+ L1ASSERT(0, fNFindIterations == (int) fParameters.GetCAIterations().size());
+ for (fCurrentIterationIndex = 0; fCurrentIterationIndex < fParameters.GetCAIterations().size();
+ ++fCurrentIterationIndex) {
- // ---- Loop over Track Finder iterations ----------------------------------------------------------------//
+ // current iteration of the CA tracker
+ const auto& caIteration = fParameters.GetCAIterations()[fCurrentIterationIndex];
+ fpCurrentIteration = &caIteration;
+ if (fCurrentIterationIndex > 0) {
+ for (int ista = 0; ista < fParameters.GetNstationsActive(); ++ista) {
+ vGrid[ista].RemoveUsedHits(fWindowHits, fvHitKeyFlags);
+ }
+ }
- L1ASSERT(0, fNFindIterations == (int) fParameters.GetCAIterations().size());
- isec = 0; // TODO: temporary! (S.Zharko)
-
- for (const auto& caIteration : fParameters.GetCAIterations()) // all finder
- {
- fpCurrentIteration = &caIteration; // select current iteration
+ fIsWindowHitSuppressed.reset(fWindowHits.size(), 0);
for (int j = 0; j < fParameters.GetNstationsActive(); j++) {
fTriplets[j].clear();
}
- /// isec - number of current iterations, fNFindIterations - number of all iterations
+
#ifdef COUNTERS
Tindex nSinglets = 0;
#endif
- if (isec != 0) {
+ fHitFirstTriplet.reset(fWindowHits.size(), 0);
+ fHitNtriplets.reset(fWindowHits.size(), 0);
- fGridHitIds.swap(fGridHitIdsBuf);
- fGridHits.swap(fGridHitsBuf);
- fGridPoints.swap(fGridPointsBuf);
- }
- // TODO: Replace NStations with fInitManager.GetNstationsGeom() (S.Zharko)
- for (int ist = 0; ist < fParameters.GetNstationsActive(); ++ist) {
- for (ca::HitIndex_t ih = fGridHitStartIndex[ist]; ih < fGridHitStopIndex[ist]; ++ih) {
- //SG!!
- fHitFirstTriplet[ih] = 0;
- fHitNtriplets[ih] = 0;
- }
- }
- /*
- fHitFirstTriplet.assign(fGridHitStopIndex[fParameters.GetNstationsActive() - 1],0);
- fHitNtriplets.assign(fGridHitStopIndex[fParameters.GetNstationsActive() - 1],0);
-*/
{
// #pragma omp task
{
@@ -341,14 +352,6 @@ void L1Algo::CaTrackFinderSlice()
/// If sort by y then it is max diff between same station's modules (~0.4cm)
fMaxDZ = caIteration.GetMaxDZ(); //0;
}
-
-#ifndef L1_NO_ASSERT
- for (int istal = fParameters.GetNstationsActive() - 1; istal >= 0; istal--) {
- L1_ASSERT(fGridHitStopIndex[istal] >= fGridHitStartIndex[istal],
- fGridHitStopIndex[istal] << " >= " << fGridHitStartIndex[istal]);
- L1_ASSERT(fGridHitStopIndex[istal] <= fGridHits.size(), fGridHitStopIndex[istal] << " <= " << fGridHits.size());
- }
-#endif // L1_NO_ASSERT
}
@@ -360,12 +363,15 @@ void L1Algo::CaTrackFinderSlice()
for (int istal = fParameters.GetNstationsActive() - 2; istal >= fFirstCAstation;
istal--) { // start with downstream chambers
- Tindex NHits_l = fGridHitStopIndex[istal] - fGridHitStartIndex[istal];
- for (Tindex ih = 0; ih < NHits_l; ++ih) {
- constructor1.CreateTripletsForHit(istal, istal + 1, istal + 2, ih);
+ Tindex nGridEntriesL = vGrid[istal].Entries().size();
+
+ for (Tindex iel = 0; iel < nGridEntriesL; ++iel) {
+ ca::HitIndex_t ihitl = vGrid[istal].Entries()[iel].WindowHitId();
+ constructor1.CreateTripletsForHit(istal, istal + 1, istal + 2, ihitl);
+
if (fpCurrentIteration->GetJumpedFlag()) {
- constructor2.CreateTripletsForHit(istal, istal + 1, istal + 3, ih);
- constructor3.CreateTripletsForHit(istal, istal + 2, istal + 3, ih);
+ constructor2.CreateTripletsForHit(istal, istal + 1, istal + 3, ihitl);
+ constructor3.CreateTripletsForHit(istal, istal + 2, istal + 3, ihitl);
}
auto oldSize = fTriplets[istal].size();
@@ -377,9 +383,8 @@ void L1Algo::CaTrackFinderSlice()
fTriplets[istal].insert(fTriplets[istal].end(), constructor3.GetTriplets().begin(),
constructor3.GetTriplets().end());
- const ca::HitIndex_t ihitl = ih + fGridHitStartIndex[istal];
- fHitFirstTriplet[ihitl] = PackTripletId(istal, oldSize);
- fHitNtriplets[ihitl] = fTriplets[istal].size() - oldSize;
+ fHitFirstTriplet[ihitl] = PackTripletId(istal, oldSize);
+ fHitNtriplets[ihitl] = fTriplets[istal].size() - oldSize;
}
} // istal
@@ -478,7 +483,8 @@ void L1Algo::CaTrackFinderSlice()
for (Tindex itrip = 0; itrip < (Tindex) fTriplets[istaF].size(); ++itrip) {
L1Triplet& first_trip = (fTriplets[istaF][itrip]);
- if (fvHitKeyFlags[fGridHits[first_trip.GetLHit()].f] || fvHitKeyFlags[fGridHits[first_trip.GetLHit()].b]) {
+ if (fvHitKeyFlags[fWindowHits[first_trip.GetLHit()].f]
+ || fvHitKeyFlags[fWindowHits[first_trip.GetLHit()].b]) {
continue;
}
// ghost suppression !!!
@@ -486,7 +492,7 @@ void L1Algo::CaTrackFinderSlice()
if (!fpCurrentIteration->GetTrackFromTripletsFlag()) { // ghost suppression !!!
int nHits = 3 + first_trip.GetLevel();
- if (fGridHits[first_trip.GetLHit()].iSt == 0) {
+ if (fWindowHits[first_trip.GetLHit()].iSt == 0) {
if (nHits < fpCurrentIteration->GetMinNhitsStation0()) { continue; }
}
else {
@@ -503,7 +509,7 @@ void L1Algo::CaTrackFinderSlice()
// curr_tr.Lengtha = 0;
curr_tr.ista = 0;
curr_tr.fHits.clear();
- curr_tr.fHits.push_back(fGridHitIds[first_trip.GetLHit()]);
+ curr_tr.fHits.push_back(fWindowHits[first_trip.GetLHit()].Id());
curr_tr.NHits = 1;
curr_L = 1;
@@ -675,45 +681,21 @@ void L1Algo::CaTrackFinderSlice()
// suppress strips of suppressed hits
- for (ca::HitIndex_t ip = 0; ip < fGridPoints.size(); ip++) {
- const ca::GridEntry& hp = fGridPoints[ip];
- if (hp.IsSuppressed()) {
- const ca::Hit& hit = fGridHits[ip];
+ for (unsigned int ih = 0; ih < fWindowHits.size(); ih++) {
+ if (fIsWindowHitSuppressed[ih]) {
+ const ca::Hit& hit = fWindowHits[ih];
fvHitKeyFlags[hit.f] = 1;
fvHitKeyFlags[hit.b] = 1;
}
}
- if (isec < (fNFindIterations - 1)) {
- int NHitsOnStation = 0;
-
- for (int ista = 0; ista < fParameters.GetNstationsActive(); ++ista) {
- int start = fGridHitStartIndex[ista];
- int Nelements = fGridHitStopIndex[ista] - start;
- ca::Hit* staHits = nullptr; // to avoid out-of-range error in ..[start]
- ca::HitIndex_t* staHitIndices = nullptr;
- ca::GridEntry* staHitPoints = nullptr;
- if (Nelements > 0) {
- staHits = &(fGridHits[start]);
- staHitIndices = &(fGridHitIds[start]);
- staHitPoints = &(fGridPoints[start]);
- }
-
- int NHitsOnStationTmp = NHitsOnStation;
-
- vGrid[ista].UpdateIterGrid(Nelements, staHits, fGridHitIdsBuf, staHitIndices, fGridHitsBuf, fGridPointsBuf,
- staHitPoints, NHitsOnStation, fvHitKeyFlags);
-
- fGridHitStartIndex[ista] = NHitsOnStationTmp;
- fGridHitStopIndex[ista] = NHitsOnStation;
- }
-
+ {
// #ifdef DRAW
// draw->ClearVeiw();
// // draw->DrawInfo();
// draw->DrawRestHits(fGridHitStartIndex, fGridHitStopIndex, RealIHit);
// draw->DrawRecoTracks();
- // draw->SaveCanvas("Reco_"+isec+"_");
+ // draw->SaveCanvas("Reco_"+iCaIteration+"_");
// draw->DrawAsk();
// #endif
@@ -721,19 +703,16 @@ void L1Algo::CaTrackFinderSlice()
// fL1Pulls->Build(1);
// #endif
#ifdef COUNTERS
- // stat_nHits[isec] += (fGridHits*)->Size();
-
- cout << "iter = " << isec << endl;
- cout << " NSinglets = " << stat_nSinglets[isec] / stat_N << endl;
- // cout << " NDoublets = " << stat_nDoublets[isec]/stat_N << endl;
- cout << " NTriplets = " << stat_nTriplets[isec] / stat_N << endl;
- cout << " NfGridHit = " << stat_nHits[isec] / stat_N << endl;
-
+ // stat_nHits[iCaIteration] += (fGridHits*)->Size();
+ cout << "iter = " << iCaIteration << endl;
+ cout << " NSinglets = " << stat_nSinglets[iCaIteration] / stat_N << endl;
+ // cout << " NDoublets = " << stat_nDoublets[iCaIteration]/stat_N << endl;
+ cout << " NTriplets = " << stat_nTriplets[iCaIteration] / stat_N << endl;
+ cout << " NfGridHit = " << stat_nHits[iCaIteration] / stat_N << endl;
#endif // COUNTERS
}
- ++isec; // TODO: temporary, to be removed! ()
- } // for (int isec
- // ---- Loop over Track Finder iterations: END -----------------------------------------------------------//
+
+ } // ---- Loop over Track Finder iterations: END -----------------------------------------------------------//
// Fit tracks
this->L1KFTrackFitter();
@@ -799,11 +778,11 @@ void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best_L, fsc
//L1_SHOW(fvHitKeyFlags.size());
//L1_SHOW(fGridHits[ihitm].f);
//L1_SHOW(fGridHits[ihitm].b);
- if (!(fvHitKeyFlags[fGridHits[ihitm].f] || fvHitKeyFlags[fGridHits[ihitm].b])) {
+ if (!(fvHitKeyFlags[fWindowHits[ihitm].f] || fvHitKeyFlags[fWindowHits[ihitm].b])) {
// curr_tr.Hits.push_back(fGridHitIds[ihitm]);
- curr_tr.fHits.push_back(fGridHitIds[ihitm]);
+ curr_tr.fHits.push_back(fWindowHits[ihitm].Id());
curr_tr.NHits++;
@@ -815,10 +794,10 @@ void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best_L, fsc
//L1_SHOW(fvHitKeyFlags.size());
//L1_SHOW(fGridHits[ihitr].f);
//L1_SHOW(fGridHits[ihitr].b);
- if (!(fvHitKeyFlags[fGridHits[ihitr].f] || fvHitKeyFlags[fGridHits[ihitr].b])) {
+ if (!(fvHitKeyFlags[fWindowHits[ihitr].f] || fvHitKeyFlags[fWindowHits[ihitr].b])) {
//curr_tr.Hits.push_back(fGridHitIds[ihitr]);
- curr_tr.fHits.push_back(fGridHitIds[ihitr]);
+ curr_tr.fHits.push_back(fWindowHits[ihitr].Id());
curr_tr.NHits++;
@@ -871,8 +850,8 @@ void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best_L, fsc
fscal dchi2 = 0.;
if (!checkTripletMatch(*curr_trip, new_trip, dchi2)) continue;
- if (fvHitKeyFlags[fGridHits[new_trip.GetLHit()].f]
- || fvHitKeyFlags[fGridHits[new_trip.GetLHit()].b]) { //hits are used
+ if (fvHitKeyFlags[fWindowHits[new_trip.GetLHit()].f]
+ || fvHitKeyFlags[fWindowHits[new_trip.GetLHit()].b]) { //hits are used
// no used hits allowed -> compare and store track
if ((curr_L > best_L) || ((curr_L == best_L) && (curr_chi2 < best_chi2))) {
best_tr = curr_tr;
@@ -888,12 +867,12 @@ void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best_L, fsc
if (0) { //SGtrd2d debug!!
- int mc01 = GetMcTrackIdForGridHit(curr_trip->GetLHit());
- int mc02 = GetMcTrackIdForGridHit(curr_trip->GetMHit());
- int mc03 = GetMcTrackIdForGridHit(curr_trip->GetRHit());
- int mc11 = GetMcTrackIdForGridHit(new_trip.GetLHit());
- int mc12 = GetMcTrackIdForGridHit(new_trip.GetMHit());
- int mc13 = GetMcTrackIdForGridHit(new_trip.GetRHit());
+ int mc01 = GetMcTrackIdForWindowHit(curr_trip->GetLHit());
+ int mc02 = GetMcTrackIdForWindowHit(curr_trip->GetMHit());
+ int mc03 = GetMcTrackIdForWindowHit(curr_trip->GetRHit());
+ int mc11 = GetMcTrackIdForWindowHit(new_trip.GetLHit());
+ int mc12 = GetMcTrackIdForWindowHit(new_trip.GetMHit());
+ int mc13 = GetMcTrackIdForWindowHit(new_trip.GetRHit());
if ((mc01 == mc02) && (mc02 == mc03)) {
cout << " sta " << ista << " mc0 " << mc01 << " " << mc02 << " " << mc03 << " mc1 " << mc11 << " " << mc12
<< " " << mc13 << " chi2 " << curr_chi2 / (2 * (curr_L + 2) - 4) << " new "
@@ -919,7 +898,7 @@ void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best_L, fsc
// add new hit
new_tr[ista] = curr_tr;
- new_tr[ista].fHits.push_back(fGridHitIds[new_trip.GetLHit()]);
+ new_tr[ista].fHits.push_back(fWindowHits[new_trip.GetLHit()].Id());
new_tr[ista].NHits++;
const int new_ista = ista + new_trip.GetMSta() - new_trip.GetLSta();
diff --git a/reco/L1/L1Algo/L1Grid.cxx b/reco/L1/L1Algo/L1Grid.cxx
index a07ed395ef40055d1d87741b9591b9a7f1077c5d..a31c1f8a03046f358918c045dd29ff0c9a488c98 100644
--- a/reco/L1/L1Algo/L1Grid.cxx
+++ b/reco/L1/L1Algo/L1Grid.cxx
@@ -5,14 +5,10 @@
#include "L1Grid.h"
#include <algorithm>
-#include <cstdio>
-#include <assert.h>
#include <string.h>
#include "CaHit.h"
-#include "L1Algo.h"
-#include "L1Def.h"
using namespace std; // !! REMOVE
@@ -26,6 +22,7 @@ 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)
@@ -42,15 +39,6 @@ inline void memset(T* dest, T i, size_t num)
}
-void L1Grid::AllocateMemory()
-{
- int binsGrid = 600000;
-
- fBinFirstEntryIndex.reset(binsGrid, 0);
- fHitsInBin.reset(binsGrid, 0);
-}
-
-
void L1Grid::BuildBins(fscal xMin, fscal xMax, fscal yMin, fscal yMax, fscal sx, fscal sy)
{
fMinX = std::min(xMin, xMax);
@@ -71,103 +59,61 @@ void L1Grid::BuildBins(fscal xMin, fscal xMax, fscal yMin, fscal yMax, fscal sx,
if (fNy < 1) fNy = 1;
fN = fNx * fNy;
-}
-
-
-void L1Grid::StoreHits(L1Algo& algo, int iS, ca::HitIndex_t& nGridHitsFilled)
-{
- ca::HitIndex_t nHits = algo.fSliceHitIds[iS].size();
-
- algo.fGridHitStartIndex[iS] = nGridHitsFilled;
+ fEntries.clear();
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] = 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 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);
-
- const ca::HitIndex_t& index1 = fBinFirstEntryIndex[bin] + fHitsInBin[ih];
- algo.fGridHits[nGridHitsFilled + index1] = hit;
- algo.fGridHitIds[nGridHitsFilled + index1] = caHitId;
- }
-
- nGridHitsFilled += nHits;
- algo.fGridHitStopIndex[iS] = nGridHitsFilled;
+ fNbinHits.reset(fN + 2, 0);
}
-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)
+void L1Grid::StoreHits(const Vector<ca::Hit>& hits, ca::HitIndex_t hitStartIndex, ca::HitIndex_t nHits,
+ const Vector<unsigned char>& hitKeyFlags)
{
- //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];
+ fNbinHits.reset(fN + 2, 0);
- 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 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++;
}
}
- 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];
- }
+ fEntries.reset(nEntries);
+
+ for (int bin = 0; bin < fN + 1; bin++) {
+ fBinFirstEntryIndex[bin + 1] = fBinFirstEntryIndex[bin] + fNbinHits[bin];
+ fNbinHits[bin] = 0;
}
- /* 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 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]++;
}
}
+}
- 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];
+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]++;
+ }
}
}
-
- NHitsOnStation += fBinFirstEntryIndex[fN + 1];
+ fBinFirstEntryIndex[fN] = nEntries;
+ fEntries.reduce(nEntries);
}
diff --git a/reco/L1/L1Algo/L1Grid.h b/reco/L1/L1Algo/L1Grid.h
index f117a93b6c3a326711f6efc16ec374573d11cd9f..3aab7ad51bea3050880ede31edf522c616cc3f1d 100644
--- a/reco/L1/L1Algo/L1Grid.h
+++ b/reco/L1/L1Algo/L1Grid.h
@@ -16,8 +16,6 @@
#include "CaSimd.h"
#include "CaVector.h"
-class L1Algo;
-
namespace
{
using cbm::algo::ca::Vector; // TMP!!
@@ -45,15 +43,12 @@ public:
/// 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 StoreHits(const Vector<ca::Hit>& hits, ca::HitIndex_t hitStartIndex, ca::HitIndex_t nHits,
+ const Vector<unsigned char>& hitKeyFlags);
- 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 RemoveUsedHits(const Vector<ca::Hit>& hits, const Vector<unsigned char>& hitKeyFlags);
/// get bin index
int GetBin(fscal X, fscal Y) const;
@@ -84,6 +79,17 @@ public:
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
@@ -95,8 +101,10 @@ private:
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<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
};
diff --git a/reco/L1/L1Algo/L1TripletConstructor.cxx b/reco/L1/L1Algo/L1TripletConstructor.cxx
index 8c4ee86c59fc88a022f5d69b2e76f0539d846173..ee418f1679f2c3a03458042aa341fef12d1e37a0 100644
--- a/reco/L1/L1Algo/L1TripletConstructor.cxx
+++ b/reco/L1/L1Algo/L1TripletConstructor.cxx
@@ -37,9 +37,6 @@ void L1TripletConstructor::InitStations(int istal, int istam, int istar)
fStaM = &fAlgo->fParameters.GetStation(fIstaM);
fStaR = &fAlgo->fParameters.GetStation(fIstaR);
- fHitsL = &(fAlgo->fGridPoints[0]) + fAlgo->fGridHitStartIndex[fIstaL];
- fHitsM = &(fAlgo->fGridPoints[0]) + fAlgo->fGridHitStartIndex[fIstaM];
- fHitsR = &(fAlgo->fGridPoints[0]) + fAlgo->fGridHitStartIndex[fIstaR];
{ // two stations for approximating the field between the target and the left hit
int sta1 = fIstaL;
@@ -79,16 +76,14 @@ void L1TripletConstructor::InitStations(int istal, int istam, int istar)
}
-void L1TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, ca::HitIndex_t ihl)
+void L1TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, ca::HitIndex_t iHitL)
{
-
InitStations(istal, istam, istar);
if (!fIsInitialized) return;
- fIhitL = ihl;
-
- ca::GridEntry& hitl = fHitsL[fIhitL];
+ fIhitL = iHitL;
+ ca::Hit& hitl = fAlgo->fWindowHits[fIhitL];
// fit a straight line through the target and the left hit.
{
@@ -182,8 +177,7 @@ void L1TripletConstructor::FindDoublets()
int iMC = -1;
if (fAlgo->fParameters.DevIsMatchDoubletsViaMc()) {
- int indL = fAlgo->fGridHitStartIndex[fIstaL] + fIhitL;
- iMC = fAlgo->GetMcTrackIdForGridHit(indL);
+ iMC = fAlgo->GetMcTrackIdForWindowHit(fIhitL);
if (iMC < 0) { return; }
}
@@ -206,10 +200,10 @@ void L1TripletConstructor::FitDoublets()
for (unsigned int i2 = 0; i2 < hitsMtmp.size(); i2++) {
- const ca::HitIndex_t imh = hitsMtmp[i2];
- const ca::GridEntry& hitm = fHitsM[imh];
+ const ca::HitIndex_t indM = hitsMtmp[i2];
+ const ca::Hit& hitm = fAlgo->fWindowHits[indM];
- if (hitm.IsSuppressed()) continue;
+ if (fAlgo->fIsWindowHitSuppressed[indM]) continue;
TrackParamV& T2 = fFit.Tr();
T2 = fTrL;
@@ -247,37 +241,35 @@ void L1TripletConstructor::FitDoublets()
fscal ty = T2.Ty()[0];
fscal tt = T2.Vi()[0] * sqrt(1. + tx * tx + ty * ty); // dt/dl * dl/dz
- for (unsigned int j2 = i2 + 1; j2 < hitsMtmp.size(); j2++) {
+ for (unsigned int iClone = i2 + 1; iClone < hitsMtmp.size(); iClone++) {
- ca::GridEntry& hitm1 = fHitsM[hitsMtmp[j2]];
+ int indClone = hitsMtmp[iClone];
+ const ca::Hit& hitClone = fAlgo->fWindowHits[indClone];
- fscal dz = hitm1.Z() - T2.Z()[0];
+ fscal dz = hitClone.Z() - T2.Z()[0];
if ((staM().timeInfo) && (T2.NdfTime()[0] > -2)) {
- fscal dt = T2.Time()[0] + tt * dz - hitm.T();
- if (fabs(dt) > sqrt(30. * T2.C55()[0]) + hitm1.RangeT()) { continue; }
+ fscal dt = T2.Time()[0] + tt * dz - hitClone.T();
+ if (fabs(dt) > sqrt(30. * T2.C55()[0]) + hitClone.RangeT()) { continue; }
}
- fscal dx = T2.GetX()[0] + tx * dz - hitm1.X();
- if (fabs(dx) > 1.3 * (3.5 * sqrt(T2.C00()[0]) + hitm1.RangeX())) { continue; }
+ fscal dx = T2.GetX()[0] + tx * dz - hitClone.X();
+ if (fabs(dx) > 1.3 * (3.5 * sqrt(T2.C00()[0]) + hitClone.RangeX())) { continue; }
- fscal dy = T2.Y()[0] + ty * dz - hitm1.Y();
- if (fabs(dy) > 1.6 * (3.5 * sqrt(T2.C11()[0]) + hitm1.RangeY())) { continue; }
+ fscal dy = T2.Y()[0] + ty * dz - hitClone.Y();
+ if (fabs(dy) > 1.6 * (3.5 * sqrt(T2.C11()[0]) + hitClone.RangeY())) { continue; }
if (fAlgo->fParameters.DevIsSuppressOverlapHitsViaMc()) {
- int indL = fAlgo->fGridHitStartIndex[fIstaL] + fIhitL;
- int indM = fAlgo->fGridHitStartIndex[fIstaM] + imh;
- int indM1 = fAlgo->fGridHitStartIndex[fIstaM] + hitsMtmp[j2];
- int iMC = fAlgo->GetMcTrackIdForGridHit(indL);
- if ((iMC != fAlgo->GetMcTrackIdForGridHit(indM)) || (iMC != fAlgo->GetMcTrackIdForGridHit(indM1))) {
+ int iMC = fAlgo->GetMcTrackIdForWindowHit(fIhitL);
+ if ((iMC != fAlgo->GetMcTrackIdForWindowHit(indM)) || (iMC != fAlgo->GetMcTrackIdForWindowHit(indClone))) {
continue;
}
}
- hitm1.SetIsSuppresed(1);
+ fAlgo->fIsWindowHitSuppressed[indClone] = 1;
}
}
- fHitsM_2.push_back(imh);
+ fHitsM_2.push_back(indM);
fTracks_2.push_back(T2);
} // i2
}
@@ -330,10 +322,9 @@ void L1TripletConstructor::FindRightHit()
int iMC = -1;
if (fAlgo->fParameters.DevIsMatchTripletsViaMc()) {
- int indL = fAlgo->fGridHitStartIndex[fIstaL] + fIhitL;
- int indM = fAlgo->fGridHitStartIndex[fIstaM] + fHitsM_2[i2];
- iMC = fAlgo->GetMcTrackIdForGridHit(indL);
- if (iMC < 0 || iMC != fAlgo->GetMcTrackIdForGridHit(indM)) { continue; }
+ int indM = fHitsM_2[i2];
+ iMC = fAlgo->GetMcTrackIdForWindowHit(fIhitL);
+ if (iMC < 0 || iMC != fAlgo->GetMcTrackIdForWindowHit(indM)) { continue; }
}
Vector<ca::HitIndex_t> collectedHits;
@@ -347,12 +338,9 @@ void L1TripletConstructor::FindRightHit()
}
for (unsigned int ih = 0; ih < collectedHits.size(); ih++) {
- ca::HitIndex_t irh = collectedHits[ih];
- const ca::GridEntry& hitr = fHitsR[irh];
- if (hitr.IsSuppressed()) continue;
-
+ ca::HitIndex_t irh = collectedHits[ih];
+ if (fAlgo->fIsWindowHitSuppressed[irh]) continue;
// pack the triplet
-
fHitsM_3.push_back(fHitsM_2[i2]);
fHitsR_3.push_back(irh);
@@ -404,9 +392,10 @@ void L1TripletConstructor::FitTriplets()
for (int i3 = 0; i3 < n3; ++i3) {
// prepare data
- ca::HitIndex_t ihit[NHits] = {fAlgo->fGridHitIds[fIhitL + fAlgo->fGridHitStartIndex[ista[0]]],
- fAlgo->fGridHitIds[fHitsM_3[i3] + fAlgo->fGridHitStartIndex[ista[1]]],
- fAlgo->fGridHitIds[fHitsR_3[i3] + fAlgo->fGridHitStartIndex[ista[2]]]};
+ ca::HitIndex_t iwhit[NHits] = {fIhitL, fHitsM_3[i3], fHitsR_3[i3]};
+
+ ca::HitIndex_t ihit[NHits] = {fAlgo->fWindowHits[iwhit[0]].Id(), fAlgo->fWindowHits[iwhit[1]].Id(),
+ fAlgo->fWindowHits[iwhit[2]].Id()};
if (fAlgo->fParameters.DevIsMatchTripletsViaMc()) {
int mc1 = fAlgo->GetMcTrackIdForCaHit(ihit[0]);
@@ -547,8 +536,8 @@ void L1TripletConstructor::FitTriplets()
if ((mc1 >= 0) && (mc1 == mc2) && (mc1 == mc3)) {
const CbmL1MCTrack& mctr = CbmL1::Instance()->GetMcTracks()[mc1];
cbm::ca::tools::Debugger::Instance().FillNtuple(
- "triplets", mctr.iEvent, fAlgo->isec, ih0, h0.x, h0.y, h0.z, ih1, h1.x, h1.y, h1.z, ih2, h2.x, h2.y, h2.z,
- mc1, fIstaL, mctr.p, mctr.x, mctr.y, mctr.z, (float) T.GetChiSq()[0], (float) T.Ndf()[0],
+ "triplets", mctr.iEvent, fAlgo->fCurrentIterationIndex, ih0, h0.x, h0.y, h0.z, ih1, h1.x, h1.y, h1.z, ih2,
+ h2.x, h2.y, h2.z, mc1, fIstaL, mctr.p, mctr.x, mctr.y, mctr.z, (float) T.GetChiSq()[0], (float) T.Ndf()[0],
(float) T.ChiSqTime()[0], (float) T.NdfTime()[0]);
}
}
@@ -575,13 +564,16 @@ void L1TripletConstructor::StoreTriplets()
fscal chi2 = T3.GetChiSq()[0]; // / T3.NDF[0];
- const ca::HitIndex_t ihitl = fIhitL + fAlgo->fGridHitStartIndex[fIstaL];
- const ca::HitIndex_t ihitm = fHitsM_3[i3] + fAlgo->fGridHitStartIndex[fIstaM];
- const ca::HitIndex_t ihitr = fHitsR_3[i3] + fAlgo->fGridHitStartIndex[fIstaR];
+ const ca::HitIndex_t ihitl = fIhitL;
+ const ca::HitIndex_t ihitm = fHitsM_3[i3];
+ const ca::HitIndex_t ihitr = fHitsR_3[i3];
- L1_ASSERT(ihitl < fAlgo->fGridHitStopIndex[fIstaL], ihitl << " < " << fAlgo->fGridHitStopIndex[fIstaL]);
- L1_ASSERT(ihitm < fAlgo->fGridHitStopIndex[fIstaM], ihitm << " < " << fAlgo->fGridHitStopIndex[fIstaM]);
- L1_ASSERT(ihitr < fAlgo->fGridHitStopIndex[fIstaR], ihitr << " < " << fAlgo->fGridHitStopIndex[fIstaR]);
+ L1_ASSERT(ihitl < fAlgo->fStationHitsStartIndex[fIstaL + 1],
+ ihitl << " < " << fAlgo->fStationHitsStartIndex[fIstaL + 1]);
+ L1_ASSERT(ihitm < fAlgo->fStationHitsStartIndex[fIstaM + 1],
+ ihitm << " < " << fAlgo->fStationHitsStartIndex[fIstaM + 1]);
+ L1_ASSERT(ihitr < fAlgo->fStationHitsStartIndex[fIstaR + 1],
+ ihitr << " < " << fAlgo->fStationHitsStartIndex[fIstaR + 1]);
if (!fAlgo->fpCurrentIteration->GetTrackFromTripletsFlag()) {
if (chi2 > fAlgo->fTripletFinalChi2Cut) { continue; }
@@ -614,9 +606,11 @@ void L1TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, co
collectedHits.clear();
collectedHits.reserve(maxNhits);
- const ca::Station& sta = fAlgo->fParameters.GetStation(iSta);
- const ca::GridEntry* hits = &(fAlgo->fGridPoints[0]) + fAlgo->fGridHitStartIndex[iSta];
- const ca::HitIndex_t nHits = fAlgo->fGridHitStopIndex[iSta] - fAlgo->fGridHitStartIndex[iSta];
+ 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();
fFit.SetTrack(Tr);
TrackParamV& T = fFit.Tr();
@@ -639,11 +633,13 @@ void L1TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, co
if (!area.GetNext(ih)) { break; }
}
- const ca::GridEntry& hit = hits[ih];
- if (hit.IsSuppressed()) continue;
+ const ca::GridEntry& entry = entries[ih];
+ if (fAlgo->fIsWindowHitSuppressed[entry.WindowHitId()]) continue;
+
+ const ca::Hit& hit = fAlgo->fWindowHits[entry.WindowHitId()];
if (iMC >= 0) { // match via MC
- if (iMC != fAlgo->GetMcTrackIdForGridHit(fAlgo->fGridHitStartIndex[iSta] + ih)) { continue; }
+ if (iMC != fAlgo->GetMcTrackIdForWindowHit(fAlgo->vGrid[iSta].Entries()[ih].WindowHitId())) { continue; }
}
// check y-boundaries
@@ -691,7 +687,7 @@ void L1TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, co
if (chi2x[0] + chi2u[0] > chi2Cut) continue;
}
- collectedHits.push_back(ih);
+ collectedHits.push_back(entry.WindowHitId());
} // loop over the hits in the area
-}
+}
\ No newline at end of file
diff --git a/reco/L1/L1Algo/L1TripletConstructor.h b/reco/L1/L1Algo/L1TripletConstructor.h
index 9d127457580a67cce6a5a740a498fbf559d03bae..3a3eba9fab6cad917a138637d66a2f1bb5bd9fce 100644
--- a/reco/L1/L1Algo/L1TripletConstructor.h
+++ b/reco/L1/L1Algo/L1TripletConstructor.h
@@ -94,14 +94,10 @@ private:
const ca::Station* fStaM {nullptr}; ///< mid station
const ca::Station* fStaR {nullptr}; ///< right station
- ca::GridEntry* fHitsL {nullptr}; ///< hits on the left station
- ca::GridEntry* fHitsM {nullptr}; ///< hits on the middle station
- ca::GridEntry* fHitsR {nullptr}; ///< hits on the right station
-
const ca::Station* fFld0Sta[2]; // two stations for approximating the field between the target and the left hit
const ca::Station* fFld1Sta[3]; // three stations for approximating the field between the left and the right hit
- ca::HitIndex_t fIhitL;
+ ca::HitIndex_t fIhitL; ///< index of the left hit in fAlgo->fWindowHits
TrackParamV fTrL;
ca::FieldRegion fFldL;
diff --git a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
index c381d686c9d39817dde380e01f39f9e582be540b..18865838859cba9643163e4b7b413ed0068fc857 100644
--- a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
+++ b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
@@ -77,14 +77,14 @@ void L1AlgoDraw::InitL1Draw(L1Algo* algo_)
algo = algo_;
vHits.clear();
- vHits.reserve(algo->fGridHits.size());
- for (unsigned int i = 0; i < algo->GetInputData().GetNhits(); i++) {
- vHits.push_back(algo->fGridHits[i]);
+ vHits.reserve(algo->fWindowHits.size());
+ for (unsigned int i = 0; i < algo->fWindowHits.size(); i++) {
+ vHits.push_back(algo->fWindowHits[i]);
}
NStations = algo->GetParameters()->GetNstationsActive();
for (int i = 0; i < NStations; i++) {
- HitsStartIndex[i] = algo->fGridHitStartIndex[i];
- HitsStopIndex[i] = algo->fGridHitStopIndex[i];
+ HitsStartIndex[i] = algo->fStationHitsStartIndex[i];
+ HitsStopIndex[i] = algo->fStationHitsStartIndex[i] + algo->fStationNhits[i];
vStations[i] = algo->GetParameters()->GetStation(i);
}
}