diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 2d5393ec9cef995525add61f2bd49b7207d9c7e4..00f4b16b20d1807c72a6779121f2d1ccd4b23d0e 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -41,12 +41,15 @@
#include "TClonesArray.h"
#include "TH1.h"
+#include <boost/functional/hash.hpp>
+
#include <algorithm>
#include <fstream>
#include <iostream>
#include <map>
#include <set>
#include <string_view>
+#include <unordered_map>
#include <utility>
#include "L1Algo/L1Algo.h"
@@ -67,6 +70,7 @@ class CbmMCDataObject;
class CbmEvent;
class TProfile2D;
+/// TODO: SZh 21.09.2022: Replace instances of this class with L1Hit
class CbmL1HitStore {
public:
int ExtIndex;
@@ -81,6 +85,36 @@ public:
int Det;
};
+/// Internal structure to handle link keys
+struct CbmL1LinkKey {
+ /// Constructor from links
+ CbmL1LinkKey(int32_t index, int32_t entry, int32_t file) : fIndex(index), fEntry(entry), fFile(file) {}
+ bool operator==(const CbmL1LinkKey& other) const
+ {
+ return fFile == other.fFile && fEntry == other.fEntry && fIndex == other.fIndex;
+ }
+
+ int32_t fIndex = -1; ///< index of point/track, saved to link
+ int32_t fEntry = -1; ///< index of link entry
+ int32_t fFile = -1; ///< index of link file
+};
+
+/// Hash for CbmL1LinkKey
+namespace std
+{
+ template<>
+ struct hash<CbmL1LinkKey> {
+ std::size_t operator()(const CbmL1LinkKey& key) const
+ {
+ std::size_t res = 0;
+ boost::hash_combine(res, key.fFile);
+ boost::hash_combine(res, key.fEntry);
+ boost::hash_combine(res, key.fIndex);
+ return res;
+ }
+ };
+} // namespace std
+
/// Enumeration for the detector subsystems used in L1 tracking
/// It is important for the subsystems to be specified in the actual order. The order is used
@@ -109,7 +143,6 @@ public:
// **********************
using DFSET = std::set<std::pair<int, int>>; // Why std::set<std::pair<>> instead of std::map<>?
- using DFEI2I = std::map<Double_t, Int_t>;
// **************************
// ** Friends classes list **
@@ -336,7 +369,7 @@ private:
// static bool compareZ(const int &a, const int &b );
// bool compareZ(const int &a, const int &b );
- /// Fills the fvMCTracks vector and the dFEI2vMCTracks set
+ /// Fills the fvMCTracks vector and the fmMCTracksMap
void Fill_vMCTracks();
/*
@@ -408,8 +441,6 @@ private:
static std::istream& eatwhite(std::istream& is); // skip spaces
static void writedir2current(TObject* obj); // help procedure
- inline Double_t dFEI(int file, int event, int idx) { return (1000 * idx) + file + (0.0001 * event); }
-
private:
std::string fInputDataFilename = ""; ///< File name to read/write input hits
@@ -569,8 +600,9 @@ public:
// L1Vector<int> vHitMCRef1;
// CbmMatch HitMatch;
private:
- DFEI2I dFEI2vMCPoints = {};
- DFEI2I dFEI2vMCTracks = {};
+ std::unordered_map<CbmL1LinkKey, int> fmMCPointsLinksMap = {}; /// Internal MC point index vs. link
+ std::unordered_map<CbmL1LinkKey, int> fmMCTracksLinksMap = {}; /// Internal MC track index vs. link
+
// *****************************
// ** Tracking performance QA **
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 8e5311961f2c202af4589d51960e6a33cc09d71c..a015b73cf0b565bc5ab3fcbb122d91ae3ea73254 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -188,8 +188,9 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
fvRecoTracks.clear(); /* <CbmL1Track> */ // TODO: Move from this function to more suitable place (S.Zharko)
fvHitPointIndexes.clear(); /* <int>: indexes of MC-points in fvMCPoints (by index of fpAlgo->vHits) */
fvHitStore.clear(); /* <CbmL1HitStore> */
- dFEI2vMCPoints.clear(); /* dFEI vs MC-point index: dFEI = index * 10000 + fileID + eventID * 0.0001 */
- dFEI2vMCTracks.clear(); /* dFEI vs MC-track index: dFEI = index * 10000 + fileID + eventID * 0.0001 */
+
+ fmMCPointsLinksMap.clear();
+ fmMCTracksLinksMap.clear();
if (fVerbose >= 10) cout << "ReadEvent: clear is done." << endl;
@@ -241,10 +242,10 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
*/
if (fPerformance) {
- // Fill fvMCTracks and dFEI2vMCTracks
+ // Fill fvMCTracks and fmMCTracksLinksMap
Fill_vMCTracks();
- // Fill fvMCPoints and dFEI2vMCPoints
+ // Fill fvMCPoints and fmMCPointsLinksMap
fvMCPoints.clear();
fvMCPoints.reserve(5 * fvMCTracks.size() * fNStations);
fvMCPointIndexesTs.clear();
@@ -276,12 +277,11 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
}
assert(MC.iStation >= 0);
if (fabs(maxDeviation) < fabs(bestDist)) { maxDeviation = bestDist; }
- Double_t dtrck = dFEI(iFile, iEvent, MC.ID);
- DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
- assert(trk_it != dFEI2vMCTracks.end());
- MC.ID = trk_it->second;
+ auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(MC.ID, iEvent, iFile));
+ assert(itTrack != fmMCTracksLinksMap.cend());
+ MC.ID = itTrack->second;
fvMCTracks[MC.ID].Points.push_back_no_warning(fvMCPoints.size());
- dFEI2vMCPoints.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMC), fvMCPoints.size()));
+ fmMCPointsLinksMap[CbmL1LinkKey(iMC, iEvent, iFile)] = fvMCPoints.size();
fvMCPoints.push_back(MC);
fvMCPointIndexesTs.push_back(0);
nMvdPoints++;
@@ -315,12 +315,11 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
}
assert(MC.iStation >= 0);
if (fabs(maxDeviation) < fabs(bestDist)) { maxDeviation = bestDist; }
- Double_t dtrck = dFEI(iFile, iEvent, MC.ID);
- DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
- assert(trk_it != dFEI2vMCTracks.end());
- MC.ID = trk_it->second;
+ auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(MC.ID, iEvent, iFile));
+ assert(itTrack != fmMCTracksLinksMap.cend());
+ MC.ID = itTrack->second;
fvMCTracks[MC.ID].Points.push_back_no_warning(fvMCPoints.size());
- dFEI2vMCPoints.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMC + nMvdPoints), fvMCPoints.size()));
+ fmMCPointsLinksMap[CbmL1LinkKey(iMC + nMvdPoints, iEvent, iFile)] = fvMCPoints.size();
fvMCPoints.push_back(MC);
fvMCPointIndexesTs.push_back(0);
nStsPoints++;
@@ -345,13 +344,11 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
if (MC.z > sta[iStActive].z[0] - 2.5) { MC.iStation = iStActive; }
}
assert(MC.iStation >= 0);
- Double_t dtrck = dFEI(iFile, iEvent, MC.ID);
- DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
- assert(trk_it != dFEI2vMCTracks.end());
- MC.ID = trk_it->second;
+ auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(MC.ID, iEvent, iFile));
+ assert(itTrack != fmMCTracksLinksMap.cend());
+ MC.ID = itTrack->second;
fvMCTracks[MC.ID].Points.push_back_no_warning(fvMCPoints.size());
- dFEI2vMCPoints.insert(
- DFEI2I::value_type(dFEI(iFile, iEvent, iMC + nMvdPoints + nStsPoints), fvMCPoints.size()));
+ fmMCPointsLinksMap[CbmL1LinkKey(iMC + nMvdPoints + nStsPoints, iEvent, iFile)] = fvMCPoints.size();
fvMCPoints.push_back(MC);
fvMCPointIndexesTs.push_back(0);
nMuchPoints++;
@@ -374,13 +371,12 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
}
if (MC.iStation < 0) continue;
assert(MC.iStation >= 0);
- Double_t dtrck = dFEI(iFile, iEvent, MC.ID);
- DFEI2I::iterator trk_it = dFEI2vMCTracks.find(dtrck);
- assert(trk_it != dFEI2vMCTracks.end());
- MC.ID = trk_it->second;
+ auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(MC.ID, iEvent, iFile));
+ assert(itTrack != fmMCTracksLinksMap.cend());
+ MC.ID = itTrack->second;
fvMCTracks[MC.ID].Points.push_back_no_warning(fvMCPoints.size());
- dFEI2vMCPoints.insert(
- DFEI2I::value_type(dFEI(iFile, iEvent, iMC + nMvdPoints + nStsPoints + nMuchPoints), fvMCPoints.size()));
+ fmMCPointsLinksMap[CbmL1LinkKey(iMC + nMvdPoints + nStsPoints + nMuchPoints, iEvent, iFile)] =
+ fvMCPoints.size();
fvMCPoints.push_back(MC);
fvMCPointIndexesTs.push_back(0);
nTrdPoints++;
@@ -424,10 +420,9 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
if (isTofPointMatched[iMC] == 0) continue;
CbmL1MCPoint MC;
if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 4)) {
- Double_t dtrck = dFEI(iFile, iEvent, MC.ID);
- auto trk_it = dFEI2vMCTracks.find(dtrck);
- if (trk_it == dFEI2vMCTracks.end()) continue;
- Int_t iTrack = trk_it->second;
+ auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(MC.ID, iEvent, iFile));
+ if (itTrack == fmMCTracksLinksMap.cend()) { continue; }
+ int iTrack = itTrack->second;
MC.iStation = -1;
const L1Station* sta = fpAlgo->GetParameters()->GetStations().begin();
@@ -469,13 +464,9 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
MC.ID = iTrack;
- dFEI2vMCPoints.insert(DFEI2I::value_type(
- dFEI(iFile, iEvent,
- fTofPointToTrack[iTofSta][iTrack] + nMvdPoints + nStsPoints + nMuchPoints + nTrdPoints),
- fvMCPoints.size()));
-
+ int iMC = fTofPointToTrack[iTofSta][iTrack] + nMvdPoints + nStsPoints + nMuchPoints + nTrdPoints;
+ fmMCPointsLinksMap[CbmL1LinkKey(iMC, iEvent, iFile)] = fvMCPoints.size();
fvMCPoints.push_back(MC);
-
fvMCPointIndexesTs.push_back(0);
nTofPoints++;
}
@@ -488,11 +479,11 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
sort(fvMCTracks[iTr].Points.begin(), fvMCTracks[iTr].Points.end(), compareZ);
if (fvMCTracks[iTr].mother_ID >= 0) {
- Double_t dtrck = dFEI(fvMCTracks[iTr].iFile, fvMCTracks[iTr].iEvent, fvMCTracks[iTr].mother_ID);
- DFEI2I::iterator map_it = dFEI2vMCTracks.find(dtrck);
- if (map_it == dFEI2vMCTracks.end()) fvMCTracks[iTr].mother_ID = -2;
- else
- fvMCTracks[iTr].mother_ID = map_it->second;
+ auto iFile = fvMCTracks[iTr].iFile;
+ auto iEvent = fvMCTracks[iTr].iEvent;
+ auto iMother = fvMCTracks[iTr].mother_ID;
+ auto itTrack = fmMCTracksLinksMap.find(CbmL1LinkKey(iMother, iEvent, iFile));
+ fvMCTracks[iTr].mother_ID = ((itTrack != fmMCTracksLinksMap.cend()) ? itTrack->second : -2);
}
} //iTr
if (fVerbose >= 10) cout << "Points in fvMCTracks are sorted." << endl;
@@ -772,10 +763,9 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
Int_t iFile = matchHitMatch->GetLink(0).GetFile();
Int_t iEvent = matchHitMatch->GetLink(0).GetEntry();
- Double_t dtrck = dFEI(iFile, iEvent, iIndex);
- DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck);
- if (trk_it == dFEI2vMCPoints.end()) continue;
- th.iMC = trk_it->second;
+ auto itPoint = fmMCPointsLinksMap.find(CbmL1LinkKey(iIndex, iEvent, iFile));
+ if (itPoint == fmMCPointsLinksMap.cend()) continue;
+ th.iMC = itPoint->second;
if ((1 == fMuchUseMcHit) && (th.iMC > -1))
th.SetHitFromPoint(fvMCPoints[th.iMC], fpAlgo->GetParameters()->GetStation(th.iStation));
}
@@ -1006,11 +996,9 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
//CbmTofPoint* pt = L1_DYNAMIC_CAST<CbmTofPoint*>(fTofPoints->Get(iFile, iEvent, iMC));
// pt->GetTrackID();
- Double_t dtrck = dFEI(iFile, iEvent, iIndex);
- DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck);
- if (trk_it == dFEI2vMCPoints.end()) continue;
-
- th.iMC = trk_it->second;
+ auto itPoint = fmMCPointsLinksMap.find(CbmL1LinkKey(iIndex, iEvent, iFile));
+ if (itPoint == fmMCPointsLinksMap.cend()) continue;
+ th.iMC = itPoint->second;
if ((1 == fTofUseMcHit) && (th.iMC > -1))
th.SetHitFromPoint(fvMCPoints[th.iMC], fpAlgo->GetParameters()->GetStation(th.iStation));
}
@@ -1196,8 +1184,7 @@ void CbmL1::Fill_vMCTracks()
T.isSignal = T.IsPrimary() && (T.z > header->GetZ() + 1.e-10);
fvMCTracks.push_back(T);
- // Double_t dtrck =dFEI(iFile,iEvent,iMCTrack);
- dFEI2vMCTracks.insert(DFEI2I::value_type(dFEI(iFile, iEvent, iMCTrack), fvMCTracks.size() - 1));
+ fmMCTracksLinksMap[CbmL1LinkKey(iMCTrack, iEvent, iFile)] = fvMCTracks.size() - 1;
} //iTracks
} //Links
@@ -1378,8 +1365,6 @@ void CbmL1::HitMatch()
CbmL1Hit& hit = fvExternalHits[iH];
if ((hit.Det == 1) && (2 != fStsUseMcHit)) {
- L1_SHOW(iH);
-
CbmStsHit* sh = L1_DYNAMIC_CAST<CbmStsHit*>(fpStsHits->At(fvExternalHits[iH].extIndex));
int iP = -1;
@@ -1416,15 +1401,12 @@ void CbmL1::HitMatch()
iEvent = fvFileEvent.begin()->second;
}
- Double_t dpnt = dFEI(iFile, iEvent, nMvdPoints + iIndex);
- DFEI2I::iterator pnt_it = dFEI2vMCPoints.find(dpnt);
-
-
- assert(pnt_it != dFEI2vMCPoints.end());
+ auto itPoint = fmMCPointsLinksMap.find(CbmL1LinkKey(iIndex + nMvdPoints, iEvent, iFile));
+ assert(itPoint != fmMCPointsLinksMap.cend());
if (link.GetWeight() > bestWeight) {
bestWeight = link.GetWeight();
- iP = pnt_it->second;
+ iP = itPoint->second;
}
}
} //mach cluster
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index df30f45a903d2364a3158b5945fa216ac6a52b4e..872c344ca481c0dbfdb408b8b35bb8efb0c2d5b6 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -216,8 +216,12 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
fld1.Set(b10, fld1Sta0.z, b11, fld1Sta1.z, b12, fld1Sta2.z);
T.chi2 = 0.;
- T.NDF = 2.;
+ T.NDF = 2.; /// Iterations -> Number of parameters - number of measurements,
if ((isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter)) T.NDF = fvec(0.);
+ // TODO: iteration parameter: "Starting NDF of track parameters"
+ // NDF = number of track parameters (6: x, y, tx, ty, qp, time) - number of measured parameters (3: x, y, time) on station or (2: x, y) on target
+ // Alternative: Iteration can find tracks starting from target or from station: -> use a FLAG
+
T.tx = tx;
T.ty = ty;
T.t = time;
@@ -1648,6 +1652,13 @@ inline void L1Algo::TripletsStaPort( /// creates triplets:
void L1Algo::CATrackFinder()
{
+#ifdef TRACKS_FROM_TRIPLETS
+ // TODO investigate kAllPrimJumpIter & kAllSecJumpIter
+ fNFindIterations = TRACKS_FROM_TRIPLETS_ITERATION + 1;
+#else
+ fNFindIterations = fParameters.GetNcaIterations();
+#endif
+
#ifdef _OPENMP
omp_set_num_threads(fNThreads);
@@ -1847,12 +1858,6 @@ void L1Algo::CATrackFinder()
// ---- Loop over Track Finder iterations ----------------------------------------------------------------//
-#ifdef TRACKS_FROM_TRIPLETS
- // TODO investigate kAllPrimJumpIter & kAllSecJumpIter
- fNFindIterations = TRACKS_FROM_TRIPLETS_ITERATION + 1;
-#else
- fNFindIterations = fParameters.GetNcaIterations();
-#endif
L1ASSERT(0, fNFindIterations == (int) fParameters.GetCAIterations().size());
isec = 0; // TODO: temporary! (S.Zharko)
@@ -2075,6 +2080,7 @@ void L1Algo::CATrackFinder()
);
if ((isec == kFastPrimJumpIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecJumpIter) || (fMissingHits)) {
+ // All iterations are "jump"!
Tindex nG_2;
hitsmG_2.clear();
i1G_2.clear();
@@ -2141,6 +2147,8 @@ void L1Algo::CATrackFinder()
if (isec == TRACKS_FROM_TRIPLETS_ITERATION) min_level = 0;
#endif
+ // TODO: Just remove it
+ // // min_level: lower then this triplets would never start
// int min_level = 1; // min level for start triplet. So min track length = min_level+3.
// if (isec == kAllPrimJumpIter) min_level = 1;
// if ( (isec == kAllSecIter) || (isec == kAllSecJumpIter) ) min_level = 2; // only the long low momentum tracks
diff --git a/reco/L1/L1Algo/L1Utils.h b/reco/L1/L1Algo/L1Utils.h
index 305d5025806f6d4d0e223940cd604d7bb4d41429..1fb299955b0a714fcc025b65e6c2d2ba9c1745c1 100644
--- a/reco/L1/L1Algo/L1Utils.h
+++ b/reco/L1/L1Algo/L1Utils.h
@@ -23,8 +23,10 @@
#include <cmath>
+#include "CaSimd.h"
#include "L1Def.h"
+
/// Class contains some utility functions for L1Algo
namespace L1Utils
{
@@ -42,10 +44,10 @@ namespace L1Utils
|| fabs(lhs - rhs) < std::numeric_limits<T>::min();
}
- [[gnu::always_inline]] inline void CheckSimdVectorEquality(fvec v, const char* name)
+ [[gnu::always_inline]] inline void CheckSimdVectorEquality(cbm::algo::ca::fvec v, const char* name)
{
bool ok = true;
- for (size_t i = 1; i < fvec::size(); i++) {
+ for (size_t i = 1; i < cbm::algo::ca::fvec::size(); i++) {
ok = ok && (v[i] == v[0]);
}
if (!ok) {