diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 2e64b389e81482795a98e2d2745850466cf23bea..4d73329cc7853b18eb3eeb25ce4ad6c781684774 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -359,6 +359,14 @@ InitStatus CbmL1::Init()
cout << endl << endl;
}
+ if (fVerbose > 1) {
+#ifdef _OPENMP
+ LOG(info) << "L1 is running with OpenMP" << endl;
+#else
+ LOG(info) << "L1 is running without OpenMP" << endl;
+#endif
+ }
+
FairRootManager* fManger = FairRootManager::Instance();
FairRunAna* Run = FairRunAna::Instance();
@@ -1299,7 +1307,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
if (!fTimesliceMode) {
- fData->vSFlag.clear();
+ fData->fStripFlag.clear();
newTS = 0;
TsStart = 0;
@@ -1334,18 +1342,18 @@ void CbmL1::Reconstruct(CbmEvent* event)
#ifdef USE_EVENT_NUMBER
h.n = mcp.event;
#endif
- const int ista = (*algo->vSFlag)[h.f] / 4;
+ const int ista = (*algo->fStripFlag)[h.f] / 4;
const L1Station& sta = algo->vStations[ista];
if (std::find(strips.begin(), strips.end(), h.f) != strips.end()) { // separate strips
- const_cast<vector<unsigned char>*>(algo->vSFlag)->push_back((*algo->vSFlag)[h.f]);
+ (*algo->fStripFlag).push_back((*algo->fStripFlag)[h.f]);
h.f = algo->NStsStrips;
algo->NStsStrips++;
}
strips.push_back(h.f);
if (std::find(strips.begin(), strips.end(), h.b) != strips.end()) {
- const_cast<vector<unsigned char>*>(algo->vSFlag)->push_back((*algo->vSFlag)[h.b]);
+ (*algo->fStripFlag).push_back((*algo->fStripFlag)[h.b]);
h.b = algo->NStsStrips;
algo->NStsStrips++;
}
@@ -1468,8 +1476,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
float TsStart_new = TsStart + TsLength - TsOverlap;
- for (vector<L1Track>::iterator it = algo->vTracks.begin(); it != (algo->vTracks.begin() + algo->NTracksIsecAll);
- it++) {
+ for (L1Vector<L1Track>::iterator it = algo->fTracks.begin(); it != algo->fTracks.end(); it++) {
CbmL1Track t;
for (int i = 0; i < 7; i++)
@@ -1494,13 +1501,13 @@ void CbmL1::Reconstruct(CbmEvent* event)
for (int i = 0; i < it->NHits; i++) {
int start_hit1 = start_hit;
- if (algo->vRecoHits[start_hit1] > vStsHits.size() - 1) start_hit++;
+ if (algo->fRecoHits[start_hit1] > vStsHits.size() - 1) start_hit++;
else if (fTimesliceMode)
- t.StsHits.push_back(((*algo->vStsHits)[algo->vRecoHits[start_hit]]).ID);
+ t.StsHits.push_back(((*algo->vStsHits)[algo->fRecoHits[start_hit]]).ID);
else
- t.StsHits.push_back(algo->vRecoHits[start_hit]);
+ t.StsHits.push_back(algo->fRecoHits[start_hit]);
- StsHitsLocal.push_back(algo->vRecoHits[start_hit++]);
+ StsHitsLocal.push_back(algo->fRecoHits[start_hit++]);
}
t.mass = algo->fDefaultMass; // pion mass
@@ -1542,8 +1549,8 @@ void CbmL1::Reconstruct(CbmEvent* event)
/// set strips as unused
for (int i = 0; i < StsHitsLocal.size(); i++) {
- algo->SetFUnUsed(const_cast<unsigned char&>((*algo->vSFlag)[vStsHits[StsHitsLocal[i]].f]));
- algo->SetFUnUsed(const_cast<unsigned char&>((*algo->vSFlag)[vStsHits[StsHitsLocal[i]].b]));
+ algo->SetFUnUsed(const_cast<unsigned char&>((*algo->fStripFlag)[vStsHits[StsHitsLocal[i]].f]));
+ algo->SetFUnUsed(const_cast<unsigned char&>((*algo->fStripFlag)[vStsHits[StsHitsLocal[i]].b]));
}
}
vRTracksCur.push_back(t);
@@ -1702,8 +1709,8 @@ void CbmL1::writedir2current(TObject* obj)
void CbmL1::IdealTrackFinder()
{
- algo->vTracks.clear();
- algo->vRecoHits.clear();
+ algo->fTracks.clear();
+ algo->fRecoHits.clear();
for (vector<CbmL1MCTrack>::iterator i = vMCTracks.begin(); i != vMCTracks.end(); ++i) {
CbmL1MCTrack& MC = *i;
@@ -1732,7 +1739,7 @@ void CbmL1::IdealTrackFinder()
const int hitI = hitIndices[iH];
if (hitI < 0) continue;
- // algo->vRecoHits.push_back(hitI);
+ // algo->fRecoHits.push_back(hitI);
algoTr.NHits++;
}
@@ -1742,7 +1749,7 @@ void CbmL1::IdealTrackFinder()
for (int iH = 0; iH < algo->NStations; iH++) {
const int hitI = hitIndices[iH];
if (hitI < 0) continue;
- algo->vRecoHits.push_back(hitI);
+ algo->fRecoHits.push_back(hitI);
}
@@ -1754,10 +1761,9 @@ void CbmL1::IdealTrackFinder()
algoTr.TFirst[4] = MC.q / MC.p;
algoTr.TFirst[5] = MC.z;
- algo->vTracks.push_back(algoTr);
+ algo->fTracks.push_back(algoTr);
}
- algo->NTracksIsecAll = algo->vTracks.size();
-}; // void CbmL1::IdealTrackFinder()
+} // void CbmL1::IdealTrackFinder()
/// ----- STandAlone Package service-functions -----------------------------
@@ -1812,20 +1818,20 @@ void CbmL1::WriteSTAPAlgoData() // must be called after ReadEvent
cout << "vStsZPos[" << n << "]"
<< " have been written." << endl;
}
- // write vSFlag
- n = (*algo->vSFlag).size();
+ // write fStripFlag
+ n = (*algo->fStripFlag).size();
fadata << n << endl;
unsigned char element;
for (int i = 0; i < n; i++) {
- element = (*algo->vSFlag)[i];
+ element = (*algo->fStripFlag)[i];
fadata << static_cast<int>(element) << endl;
};
if (fVerbose >= 4) {
- cout << "vSFlag[" << n << "]"
+ cout << "fStripFlag[" << n << "]"
<< " have been written." << endl;
}
if (fVerbose >= 4) {
- cout << "vSFlagB[" << n << "]"
+ cout << "fStripFlagB[" << n << "]"
<< " have been written." << endl;
}
// write vStsHits
@@ -2062,7 +2068,7 @@ void CbmL1::ReadSTAPAlgoData()
if (algo->vStsHits) const_cast<std::vector<L1StsHit>*>(algo->vStsHits)->clear();
algo->NStsStrips = 0;
if (algo->vStsZPos) const_cast<std::vector<float>*>(algo->vStsZPos)->clear();
- if (algo->vSFlag) const_cast<vector<unsigned char>*>(algo->vSFlag)->clear();
+ if (algo->fStripFlag) algo->fStripFlag->clear();
// check correct position in file
char s[] = "Event: ";
@@ -2091,15 +2097,15 @@ void CbmL1::ReadSTAPAlgoData()
cout << "vStsZPos[" << n << "]"
<< " have been read." << endl;
}
- // read algo->vSFlag
+ // read algo->fStripFlag
fadata >> n;
for (int i = 0; i < n; i++) {
int element;
fadata >> element;
- const_cast<vector<unsigned char>*>(algo->vSFlag)->push_back(static_cast<unsigned char>(element));
+ algo->fStripFlag->push_back(static_cast<unsigned char>(element));
}
if (fVerbose >= 4) {
- cout << "vSFlag[" << n << "]"
+ cout << "fStripFlag[" << n << "]"
<< " have been read." << endl;
}
// read algo->vStsHits
diff --git a/reco/L1/CbmL1Def.h b/reco/L1/CbmL1Def.h
index fc483056f9cc12d3314832f37a3ca4943612ddc8..603f1e04f356ebff1fb8c8d53d7d496f5b9aeec0 100644
--- a/reco/L1/CbmL1Def.h
+++ b/reco/L1/CbmL1Def.h
@@ -3,7 +3,7 @@
Authors: Maksym Zyzak, Igor Kulakov [committer], Sergey Gorbunov */
#ifndef CbmL1Def_h
-#define CbmL1Def_h 1
+#define CbmL1Def_h
// #define FAST_CODE // FAST_CODE = more unsafe
@@ -67,77 +67,4 @@ using namespace std;
typedef int index_type;
-
-template<typename T>
-class L1Vector : public std::vector<T> {
-public:
- L1Vector(const char* name = "no name") : std::vector<T>(), fSize(0), fName(name) {};
- L1Vector(const char* name, const unsigned int n) : std::vector<T>(n), fSize(0), fName(name) {};
- L1Vector(const char* name, const unsigned int n, const unsigned int value)
- : std::vector<T>(n, value)
- , fSize(0)
- , fName(name) {};
-
-
- unsigned int Size() const { return fSize; } // Size() return number
- void Reset() { fSize = 0; }
-
- void Resize(const unsigned int n)
- {
- if (n > std::vector<T>::size()) {
-#ifdef _OPENMP
-#pragma omp critical
-#endif
- std::vector<T>::resize(n);
- }
-
- fSize = n;
- }
-
- void Store(const T& element)
- {
- if (fSize >= std::vector<T>::size()) {
-#ifdef _OPENMP
-#pragma omp critical
-#endif
- std::vector<T>::push_back(element);
- }
- else
- std::vector<T>::at(fSize) = element;
-
- fSize++;
- }
-
-
- T& operator[](const size_t index)
- {
-
- //assert(index <= fSize); // allow auto-resize by 1 element only
-
- if (index >= std::vector<T>::size()) {
-#ifdef _OPENMP
-#pragma omp critical
-#endif
- std::vector<T>::resize((index + 1) * 2);
- std::cout << "Warning: L1Vector " << fName << " autoresize to " << (index + 1) * 2 << std::endl;
- }
- if (index >= fSize) fSize = index + 1;
-
- return std::vector<T>::operator[](index);
- }
-
- const T& operator[](const size_t index) const { return std::vector<T>::operator[](index); }
-
- std::string getName()
- {
- std::string s = " L1Vector<";
- s += fName + "> ";
- return s;
- }
-
-private:
- unsigned int fSize;
- std::string fName;
-};
-
#endif // CbmL1Def_h
diff --git a/reco/L1/CbmL1MCTrack.cxx b/reco/L1/CbmL1MCTrack.cxx
index d9fcaad28103b7da56533d4c6d2130445a4a64e9..620dc497b11f3a4321fc44b691282d8dcc2e46f0 100644
--- a/reco/L1/CbmL1MCTrack.cxx
+++ b/reco/L1/CbmL1MCTrack.cxx
@@ -148,7 +148,7 @@ void CbmL1MCTrack::CalculateHitCont()
for (int ih = 0; ih < nhits; ih++) {
int jh = StsHits[ih];
const L1StsHit& h = (*algo->vStsHits)[jh];
- int ista = (*algo->vSFlag)[h.f] / 4;
+ int ista = (*algo->fStripFlag)[h.f] / 4;
if (ista - istaold == 1) ncont++;
else if (ista - istaold > 1) {
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 6c3bbe4c690fc9f8b6d4ffcbcffc34bac9822268..ffda23f659b58becf28d6a4d04dc0b5adbdf0656 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -910,12 +910,14 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
// save strips in L1Algo
fData_->NStsStrips = NStrips;
- fData_->vSFlag.resize(NStrips, 0);
+ fData_->fStripFlag.clear();
+ fData_->fStripFlag.reserve(NStrips);
+ fData_->fStripFlag.assign(NStrips, 0);
for (int ih = 0; ih < nHits; ih++) {
TmpHit& th = tmpHits[ih];
char flag = th.iStation * 4;
- fData_->vSFlag[th.iStripF] = flag;
- fData_->vSFlag[th.iStripB] = flag;
+ fData_->fStripFlag[th.iStripF] = flag;
+ fData_->fStripFlag[th.iStripB] = flag;
} // ih
if (fVerbose >= 10) { cout << "ReadEvent: strips are read." << endl; }
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index ed2dd7c5764bff2beaac305d3a444acf8233712b..7b9644eb1033d24fe4ab33c7a3795f2f3ea32006 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -7,6 +7,179 @@
#include "L1Grid.h"
#include "L1HitPoint.h"
+L1Algo::L1Algo(int nThreads, int ExpectedHits)
+ : fDupletPortionSize("L1Algo::fDupletPortionSize")
+ , fMergerTrackFirstStation("L1Algo::fMergerTrackFirstStation")
+ , fMergerTrackLastStation("L1Algo::fMergerTrackLastStation")
+ , fMergerTrackFirstHit("L1Algo::fMergerTrackFirstHit")
+ , fMergerTrackLastHit("L1Algo::fMergerTrackLastHit")
+ , fMergerTrackNeighbour("L1Algo::fMergerTrackNeighbour")
+ , fMergerTrackChi2("L1Algo::fMergerTrackChi2")
+ , fMergerTrackIsStored("L1Algo::fMergerTrackIsStored")
+ , fMergerTrackIsDownstreamNeighbour("L1Algo::fMergerTrackIsDownstreamNeighbour")
+ , fMergerTracksNew("L1Algo::fMergerTracksNew")
+ , fMergerRecoHitsNew("L1Algo::fMergerRecoHitsNew")
+ , NStations(0)
+ , // number of all detector stations
+ NMvdStations(0)
+ , // number of mvd stations
+ NStsStations(0)
+ , NFStations(0)
+ , fRadThick()
+ , NStsStrips(0) // strips positions created from hits
+ , vStsZPos(0)
+ , // all possible z-positions of hits
+ vStsHits(0) // hits as a combination of front-, backstrips and z-position
+ , fStripFlag(0) // information of hits station & using hits in tracks(),
+ , CATime(0)
+ , // time of trackfinding
+ fTracks("L1Algo::fTracks")
+ , // reconstructed tracks
+ fRecoHits("L1Algo::fRecoHits")
+ , // packed hits of reconstructed tracks
+ StsHitsStartIndex(nullptr)
+ , StsHitsStopIndex(nullptr)
+ , NHitsIsecAll(0)
+ , vStsDontUsedHits_A(ExpectedHits)
+ , vStsDontUsedHits_B(ExpectedHits)
+ , vStsDontUsedHits_Buf(ExpectedHits)
+ , vStsDontUsedHitsxy_A(ExpectedHits)
+ , vStsDontUsedHitsxy_buf(ExpectedHits)
+ , vStsDontUsedHitsxy_B(ExpectedHits)
+ , RealIHit_v(ExpectedHits)
+ , RealIHit_v_buf(ExpectedHits)
+ , RealIHit_v_buf2(ExpectedHits)
+ ,
+
+#ifdef _OPENMP
+ fHitToBestTrackF("L1Algo::fHitToBestTrackF")
+ , fHitToBestTrackB("L1Algo::fHitToBestTrackB")
+ ,
+#endif
+ fStripToTrack("L1Algo::fStripToTrack")
+ , fStripToTrackB("L1Algo::fStripToTrackB")
+ ,
+ //sh (),
+ fNThreads(nThreads)
+ , fUseHitErrors(0)
+ , fmCBMmode(0)
+ , fGlobal(0)
+ , isec(0)
+ , vStsHitsUnused()
+ , RealIHitP()
+ , RealIHitPBuf()
+ , vStsHitPointsUnused()
+ , RealIHit(nullptr)
+ , FIRSTCASTATION()
+ , threadNumberToCpuMap()
+ , TRACK_CHI2_CUT(10.)
+ , TRIPLET_CHI2_CUT(5.)
+ , DOUBLET_CHI2_CUT(5.)
+ , TIME_CUT1(0.)
+ , TIME_CUT2(0.)
+ , MaxDZ(0.)
+ ,
+#ifdef DRAW
+ draw(0)
+ ,
+#endif
+ Pick_gather(0)
+ , PickNeighbour(0)
+ , // (PickNeighbour < dp/dp_error) => triplets are neighbours
+ MaxInvMom(0)
+ , // max considered q/p for tracks
+ MaxSlope(0)
+ , targX(0)
+ , targY(0)
+ , targZ(0)
+ , // target coor
+ targB()
+ , // field in the target point
+ TargetXYInfo()
+ , // target constraint [cm]
+ vtxFieldRegion()
+ , // really doesn't used
+ vtxFieldValue()
+ , // field at teh vertex position.
+ //vTripletsP(), // container for triplets got in finding
+ fTrackingLevel(0)
+ , fGhostSuppression(0)
+ , // really doesn't used
+ fMomentumCutOff(0) // really doesn't used
+{
+
+ fDupletPortionSize.reserve(100000);
+ fTracks.reserve(40000);
+ fRecoHits.reserve(400000);
+
+ fStripToTrack.reserve(ExpectedHits * 4);
+ fStripToTrackB.reserve(ExpectedHits * 4);
+
+ for (int i = 0; i < fNThreads; i++) {
+
+ fTracks_local[i].SetName(std::stringstream() << "L1Algo::fTracks_local[" << i << "]");
+ fTracks_local[i].clear();
+ fTracks_local[i].reserve(100000);
+
+ fRecoHits_local[i].SetName(std::stringstream() << "L1Algo::fRecoHits_local[" << i << "]");
+ fRecoHits_local[i].clear();
+ fRecoHits_local[i].reserve(400000);
+
+ TripForHit[0].resize(ExpectedHits);
+ TripForHit[1].resize(ExpectedHits);
+
+ fTrackCandidates[i].SetName(std::stringstream() << "L1Algo::fTrackCandidates[" << i << "]");
+ fTrackCandidates[i].clear();
+ fTrackCandidates[i].reserve(10000);
+
+ fT_3[i].reserve(MaxPortionTriplets / fvecLen);
+ fhitsl_3[i].reserve(MaxPortionTriplets);
+ fhitsm_3[i].reserve(MaxPortionTriplets);
+ fhitsr_3[i].reserve(MaxPortionTriplets);
+ fu_front3[i].reserve(MaxPortionTriplets / fvecLen);
+ fu_back3[i].reserve(MaxPortionTriplets / fvecLen);
+ fz_pos3[i].reserve(MaxPortionTriplets / fvecLen);
+ fTimeR[i].reserve(MaxPortionTriplets / fvecLen);
+ fTimeER[i].reserve(MaxPortionTriplets / fvecLen);
+ dx[i].reserve(MaxPortionTriplets / fvecLen);
+ dy[i].reserve(MaxPortionTriplets / fvecLen);
+ du[i].reserve(MaxPortionTriplets / fvecLen);
+ dv[i].reserve(MaxPortionTriplets / fvecLen);
+
+ for (int j = 0; j < MaxNStations; j++) {
+ fTriplets[j][i].SetName(std::stringstream() << "L1Algo::fTriplets[" << i << "][" << j << "]");
+ fTriplets[j][i].reserve(ExpectedHits);
+ fTriplets[j][i].clear();
+ }
+ }
+
+ for (int i = 0; i < MaxNStations; i++) {
+ vGridTime[i].AllocateMemory(fNThreads);
+ }
+
+#ifdef _OPENMP
+ fHitToBestTrackF.reserve(ExpectedHits * 2);
+ fHitToBestTrackB.reserve(ExpectedHits * 2);
+#endif
+
+ NHitsIsecAll = ExpectedHits;
+ const int kExpectedTracks = ExpectedHits / 8;
+
+ fMergerTrackFirstStation.reserve(kExpectedTracks);
+ fMergerTrackLastStation.reserve(kExpectedTracks);
+ fMergerTrackFirstHit.reserve(kExpectedTracks);
+ fMergerTrackLastHit.reserve(kExpectedTracks);
+ fMergerTrackNeighbour.reserve(kExpectedTracks);
+ fMergerTrackChi2.reserve(kExpectedTracks);
+ fMergerTrackIsStored.reserve(kExpectedTracks);
+ fMergerTrackIsDownstreamNeighbour.reserve(kExpectedTracks);
+ fMergerTracksNew.reserve(kExpectedTracks);
+ fMergerRecoHitsNew.reserve(ExpectedHits);
+
+ // IsNext.resize(kExpectedTracks);
+}
+
+
void L1Algo::Init(const vector<fscal>& geo, const bool UseHitErrors, const bool mCBMmode)
{
@@ -172,15 +345,15 @@ void L1Algo::Init(const vector<fscal>& geo, const bool UseHitErrors, const bool
}
-void L1Algo::SetData(const vector<L1StsHit>& StsHits_, int nStsStrips_, const vector<fscal>& StsZPos_,
- const vector<unsigned char>& SFlag_, const THitI* StsHitsStartIndex_,
- const THitI* StsHitsStopIndex_, const int NhitsGlobal)
+void L1Algo::SetData(vector<L1StsHit>& StsHits_, int nStsStrips_, const vector<fscal>& StsZPos_,
+ L1Vector<unsigned char>& SFlag_, const THitI* StsHitsStartIndex_, const THitI* StsHitsStopIndex_,
+ const int NhitsGlobal)
{
vStsHits = &StsHits_;
NStsStrips = nStsStrips_;
vStsZPos = &StsZPos_;
- vSFlag = &SFlag_;
+ fStripFlag = &SFlag_;
StsHitsStartIndex = StsHitsStartIndex_;
StsHitsStopIndex = StsHitsStopIndex_;
@@ -200,23 +373,34 @@ void L1Algo::SetData(const vector<L1StsHit>& StsHits_, int nStsStrips_, const ve
RealIHit_v_buf2.resize(nHits);
#ifdef _OPENMP
- hitToBestTrackF.resize(NhitsGlobal * 2);
- hitToBestTrackB.resize(NhitsGlobal * 2);
+ fHitToBestTrackF.resize(NhitsGlobal * 2);
+ fHitToBestTrackB.resize(NhitsGlobal * 2);
+ for (unsigned int j = 0; j < fHitToBestTrackB.size(); j++) {
+ omp_init_lock(&fHitToBestTrackB[j]);
+ omp_init_lock(&fHitToBestTrackF[j]);
+ }
#endif
- vStripToTrack.resize(NStsStrips);
- vStripToTrackB.resize(NStsStrips);
+
+ fStripToTrack.reserve(NStsStrips);
+ fStripToTrackB.reserve(NStsStrips);
TripForHit[0].resize(nHits);
TripForHit[1].resize(nHits);
NHitsIsecAll = nHits;
- n_g1.resize(2 * nHits);
+ fDupletPortionSize.clear();
+ fDupletPortionSize.reserve(2 * nHits);
for (int i = 0; i < fNThreads; i++) {
- vTracks_local[i].resize(nHits / 10);
- vRecoHits_local[i].resize(nHits);
- CandidatesTrack[i].resize(nHits / 10);
- for (int j = 0; j < MaxNStations; j++)
- TripletsLocal1[j][i].resize(2 * nHits);
+ fTracks_local[i].clear();
+ fTracks_local[i].reserve(nHits / 10);
+ fRecoHits_local[i].clear();
+ fRecoHits_local[i].reserve(nHits);
+ fTrackCandidates[i].clear();
+ fTrackCandidates[i].reserve(nHits / 10);
+ for (int j = 0; j < MaxNStations; j++) {
+ fTriplets[j][i].clear();
+ fTriplets[j][i].reserve(2 * nHits);
+ }
}
/*
@@ -224,15 +408,15 @@ void L1Algo::SetData(const vector<L1StsHit>& StsHits_, int nStsStrips_, const ve
vStsStrips.resize(StsStrips_.size());
vStsStripsB.resize(StsStripsB_.size());
vStsZPos.resize(StsZPos_.size());
- vSFlag.resize(SFlag_.size());
- vSFlagB.resize(SFlagB_.size());
+ fStripFlag.resize(SFlag_.size());
+ fStripFlagB.resize(SFlagB_.size());
for(Tindex i=0; i< static_cast<Tindex>(StsHits_.size()); ++i ) vStsHits[i] = StsHits_[i];
for(Tindex i=0; i< static_cast<Tindex>(StsStrips_.size()); ++i ) vStsStrips[i] = StsStrips_[i];
for(Tindex i=0; i< static_cast<Tindex>(StsStripsB_.size()); ++i ) vStsStripsB[i] = StsStripsB_[i];
for(Tindex i=0; i< static_cast<Tindex>(StsZPos_.size()); ++i ) vStsZPos[i] = StsZPos_[i];
- for(Tindex i=0; i< static_cast<Tindex>(SFlag_.size()); ++i ) vSFlag[i] = SFlag_[i];
- for(Tindex i=0; i< static_cast<Tindex>(SFlagB_.size()); ++i ) vSFlagB[i] = SFlagB_[i];
+ for(Tindex i=0; i< static_cast<Tindex>(SFlag_.size()); ++i ) fStripFlag[i] = SFlag_[i];
+ for(Tindex i=0; i< static_cast<Tindex>(SFlagB_.size()); ++i ) fStripFlagB[i] = SFlagB_[i];
for(Tindex i=0; i<MaxNStations+1; ++i) StsHitsStartIndex[i] = StsHitsStartIndex_[i];
for(Tindex i=0; i<MaxNStations+1; ++i) StsHitsStopIndex[i] = StsHitsStopIndex_[i];*/
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index 135ef6a092db005683b4faa31b900ba7ad192a18..ef0cd8d234c8c773daac15e345a9a9e0c7c6865c 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -3,7 +3,7 @@
Authors: Maksym Zyzak, Valentina Akishina, Igor Kulakov */
#ifndef L1Algo_h
-#define L1Algo_h 1
+#define L1Algo_h
// #define TBB // TODO: Doesn't work now. Renew
@@ -54,6 +54,7 @@ class L1AlgoDraw;
#include "L1TrackPar.h"
#include "L1TrackParFit.h"
#include "L1Triplet.h"
+#include "L1Vector.h"
#ifdef _OPENMP
#include "omp.h"
@@ -80,170 +81,7 @@ typedef int Tindex;
class L1Algo {
public:
// L1Algo(int nThreads=7):
- L1Algo(int nThreads = 1, int TypicalSize = 200000)
- : n_g1("L1Algo::n_g1")
- , FirstHit("L1Algo::FirstHit")
- , LastHit("L1Algo::LastHit")
- , FirstHitIndex("L1Algo::FirstHitIndex")
- , LastHitIndex("L1Algo::LastHitIndex")
- , Neighbour("L1Algo::Neighbour")
- , TrackChi2("L1Algo::TrackChi2")
- , vRecoHitsNew("L1Algo::vRecoHitsNew")
- , vTracksNew("L1Algo::vTracksNew")
- , NStations(0)
- , // number of all detector stations
- NMvdStations(0)
- , // number of mvd stations
- NStsStations(0)
- , NFStations(0)
- , fRadThick()
- , NStsStrips(0) // strips positions created from hits
- , vStsZPos(0)
- , // all possible z-positions of hits
- vStsHits(0) // hits as a combination of front-, backstrips and z-position
- , vSFlag(0) // information of hits station & using hits in tracks(),
- , CATime(0)
- , // time of trackfinding
- vTracks("L1Algo::vTracks", 40000)
- , // reconstructed tracks
- vRecoHits("L1Algo::vRecoHits", 400000)
- , // packed hits of reconstructed tracks
- StsHitsStartIndex(nullptr)
- , StsHitsStopIndex(nullptr)
- , NHitsIsecAll(0)
- , NTracksIsecAll(0)
- , vStsDontUsedHits_A(TypicalSize)
- , vStsDontUsedHits_B(TypicalSize)
- , vStsDontUsedHits_Buf(TypicalSize)
- , vStsDontUsedHitsxy_A(TypicalSize)
- , vStsDontUsedHitsxy_buf(TypicalSize)
- , vStsDontUsedHitsxy_B(TypicalSize)
- , RealIHit_v(TypicalSize)
- , RealIHit_v_buf(TypicalSize)
- , RealIHit_v_buf2(TypicalSize)
- ,
-
-#ifdef _OPENMP
- hitToBestTrackF("L1Algo::hitToBestTrackF", TypicalSize * 2)
- , hitToBestTrackB("L1Algo::hitToBestTrackB", TypicalSize * 2)
- ,
-#endif
- vStripToTrack("L1Algo::vStripToTrack", TypicalSize * 4)
- , vStripToTrackB("L1Algo::vStripToTrackB", TypicalSize * 4)
- ,
- //sh (),
- fNThreads(nThreads)
- , fUseHitErrors(0)
- , fmCBMmode(0)
- , fGlobal(0)
- , isec(0)
- , vStsHitsUnused()
- , RealIHitP()
- , RealIHitPBuf()
- , vStsHitPointsUnused()
- , RealIHit(nullptr)
- , FIRSTCASTATION()
- , threadNumberToCpuMap()
- , TRACK_CHI2_CUT(10.)
- , TRIPLET_CHI2_CUT(5.)
- , DOUBLET_CHI2_CUT(5.)
- , TIME_CUT1(0.)
- , TIME_CUT2(0.)
- , MaxDZ(0.)
- ,
-#ifdef DRAW
- draw(0)
- ,
-#endif
- Pick_gather(0)
- , PickNeighbour(0)
- , // (PickNeighbour < dp/dp_error) => triplets are neighbours
- MaxInvMom(0)
- , // max considered q/p for tracks
- MaxSlope(0)
- , targX(0)
- , targY(0)
- , targZ(0)
- , // target coor
- targB()
- , // field in the target point
- TargetXYInfo()
- , // target constraint [cm]
- vtxFieldRegion()
- , // really doesn't used
- vtxFieldValue()
- , // field at teh vertex position.
- //vTripletsP(), // container for triplets got in finding
- fTrackingLevel(0)
- , fGhostSuppression(0)
- , // really doesn't used
- fMomentumCutOff(0) // really doesn't used
- {
-
- n_g1.resize(100000);
-
- for (int i = 0; i < fNThreads; i++) {
-
- vTracks_local[i].resize(100000);
- vRecoHits_local[i].resize(400000);
-
-
- numberCandidateThread[i] = 0;
- SavedCand[i] = 0;
- SavedHits[i] = 0;
-
- TripForHit[0].resize(TypicalSize);
- TripForHit[1].resize(TypicalSize);
- CandidatesTrack[i].resize(10000);
-
- fT_3[i].reserve(MaxPortionTriplets / fvecLen);
- fhitsl_3[i].reserve(MaxPortionTriplets);
- fhitsm_3[i].reserve(MaxPortionTriplets);
- fhitsr_3[i].reserve(MaxPortionTriplets);
- fu_front3[i].reserve(MaxPortionTriplets / fvecLen);
- fu_back3[i].reserve(MaxPortionTriplets / fvecLen);
- fz_pos3[i].reserve(MaxPortionTriplets / fvecLen);
- fTimeR[i].reserve(MaxPortionTriplets / fvecLen);
- fTimeER[i].reserve(MaxPortionTriplets / fvecLen);
- dx[i].reserve(MaxPortionTriplets / fvecLen);
- dy[i].reserve(MaxPortionTriplets / fvecLen);
- du[i].reserve(MaxPortionTriplets / fvecLen);
- dv[i].reserve(MaxPortionTriplets / fvecLen);
-
- for (int j = 0; j < MaxNStations; j++)
- TripletsLocal1[j][i].resize(400000);
- }
-
- for (int i = 0; i < MaxNStations; i++)
- vGridTime[i].AllocateMemory(fNThreads);
-
-#ifdef _OPENMP
-
- for (unsigned int j = 0; j < hitToBestTrackB.size(); j++) {
- omp_init_lock(&hitToBestTrackB[j]);
- omp_init_lock(&hitToBestTrackF[j]);
- }
-
-#endif
-
-
- for (int i = 0; i < nThreads; i++)
- for (int k = 0; k < MaxNStations; k++)
- nTripletsThread[k][i] = 0;
-
- NTracksIsecAll = 20000;
- NHitsIsecAll = TypicalSize;
-
-
- FirstHit.resize(NTracksIsecAll);
- LastHit.resize(NTracksIsecAll);
- FirstHitIndex.resize(NTracksIsecAll);
- LastHitIndex.resize(NTracksIsecAll);
- // IsUsed.resize(NTracksIsecAll);
- TrackChi2.resize(NTracksIsecAll);
- Neighbour.resize(NTracksIsecAll);
- // IsNext.resize(NTracksIsecAll);
- }
+ L1Algo(int nThreads = 1, int ExpectedHits = 200000);
L1Algo(const L1Algo&) = delete;
L1Algo operator=(const L1Algo&) = delete;
@@ -264,31 +102,30 @@ public:
float fDefaultMass = 0.10565800; // muon mass
- L1Vector<L1Triplet> TripletsLocal1[nSta][nTh];
- L1Vector<L1Branch> CandidatesTrack[nTh];
-
- Tindex portionStopIndex[nSta];
- L1Vector<Tindex> n_g1;
-
-
- int SavedCand[nTh];
- int SavedHits[nTh];
+ L1Vector<L1Triplet> fTriplets[nSta][nTh]; // created triplets at station + thread
- int numberCandidateThread[nTh];
+ // Track candidates created out of adjacent triplets before the final track selection.
+ // The candidates may share any amount of hits.
+ L1Vector<L1Branch> fTrackCandidates[nTh];
- int nTripletsThread[nSta][nTh];
+ Tindex fDupletPortionStopIndex[nSta]; // end of the duplet portions for the station
+ L1Vector<Tindex> fDupletPortionSize; // N duplets in a portion
- //for merger
- L1Vector<unsigned short> FirstHit;
- L1Vector<unsigned short> LastHit;
- L1Vector<THitI> FirstHitIndex;
- L1Vector<THitI> LastHitIndex;
- L1Vector<unsigned short> Neighbour;
- L1Vector<float> TrackChi2;
- // L1Vector<bool> IsNext;
- // L1Vector<bool> IsUsed;
- L1Vector<THitI> vRecoHitsNew;
- L1Vector<L1Track> vTracksNew;
+ //
+ // Temporary vectors used by the clone merger
+ //
+ // vectors that are parallel to fTracks
+ L1Vector<unsigned short> fMergerTrackFirstStation; // first station of a track
+ L1Vector<unsigned short> fMergerTrackLastStation; // last station of a track
+ L1Vector<THitI> fMergerTrackFirstHit; // index of the first tracks hit
+ L1Vector<THitI> fMergerTrackLastHit; // index of the last tracks hit
+ L1Vector<unsigned short> fMergerTrackNeighbour; // track that can be merged with the given track
+ L1Vector<float> fMergerTrackChi2; // chi2 of the merge
+ L1Vector<char> fMergerTrackIsStored; // is the track already stored to the output
+ L1Vector<char> fMergerTrackIsDownstreamNeighbour; // is the track a downstream neighbor of another track
+ // other vectors
+ L1Vector<L1Track> fMergerTracksNew; // vector of tracks after the merge
+ L1Vector<THitI> fMergerRecoHitsNew; // vector of track hits after the merge
#ifdef DRAW
@@ -299,8 +136,8 @@ public:
void Init(const vector<fscal>& geo, const bool UseHitErrors, const bool mCBMmode);
- void SetData(const vector<L1StsHit>& StsHits_, int nStsStrips_, const vector<fscal>& StsZPos_,
- const vector<unsigned char>& SFlag_, const THitI* StsHitsStartIndex_, const THitI* StsHitsStopIndex_,
+ void SetData(vector<L1StsHit>& StsHits_, int nStsStrips_, const vector<fscal>& StsZPos_,
+ L1Vector<unsigned char>& SFlag_, const THitI* StsHitsStartIndex_, const THitI* StsHitsStopIndex_,
const int NhitsGlobal
);
@@ -333,16 +170,16 @@ public:
int NStsStrips; // number of strips
const vector<fscal>* vStsZPos; // all possible z-positions of hits
- const vector<L1StsHit>* vStsHits; // hits as a combination of front-, backstrips and z-position
+ vector<L1StsHit>* vStsHits; // hits as a combination of front-, backstrips and z-position
L1Grid vGrid[MaxNStations]; // hits as a combination of front-, backstrips and z-position
L1Grid vGridTime[MaxNStations];
- const vector<unsigned char>* vSFlag; // information of hits station & using hits in tracks;
+ L1Vector<unsigned char>* fStripFlag; // information of hits station & using hits in tracks;
double CATime; // time of trackfinding
- L1Vector<L1Track> vTracks;
- L1Vector<THitI> vRecoHits;
+ L1Vector<L1Track> fTracks;
+ L1Vector<THitI> fRecoHits;
const THitI *StsHitsStartIndex,
*StsHitsStopIndex; // station-bounders in vStsHits array
@@ -350,7 +187,6 @@ public:
// L1Branch* pointer;
unsigned int NHitsIsecAll;
- unsigned int NTracksIsecAll;
vector<L1StsHit> vStsDontUsedHits_A;
vector<L1StsHit> vStsDontUsedHits_B;
@@ -358,22 +194,20 @@ public:
vector<L1HitPoint> vStsDontUsedHitsxy_A;
vector<L1HitPoint> vStsDontUsedHitsxy_buf;
vector<L1HitPoint> vStsDontUsedHitsxy_B;
- L1Vector<L1Track> vTracks_local[nTh];
- L1Vector<THitI> vRecoHits_local[nTh];
+ L1Vector<L1Track> fTracks_local[nTh];
+ L1Vector<THitI> fRecoHits_local[nTh];
vector<THitI> RealIHit_v;
vector<THitI> RealIHit_v_buf;
vector<THitI> RealIHit_v_buf2;
#ifdef _OPENMP
-
- L1Vector<omp_lock_t> hitToBestTrackF;
- L1Vector<omp_lock_t> hitToBestTrackB;
-
+ L1Vector<omp_lock_t> fHitToBestTrackF;
+ L1Vector<omp_lock_t> fHitToBestTrackB;
#endif
- L1Vector<int> vStripToTrack;
- L1Vector<int> vStripToTrackB;
+ L1Vector<int> fStripToTrack; // front strip to track pointers
+ L1Vector<int> fStripToTrackB; // back strip to track pointers
int fNThreads;
bool fUseHitErrors;
@@ -624,7 +458,7 @@ private:
/// Find doublets on station
void DupletsStaPort( // input
- int istal, int istam, Tindex ip, vector<Tindex>& n_g, Tindex* portionStopIndex_,
+ int istal, int istam, Tindex ip, L1Vector<Tindex>& n_g, Tindex* portionStopIndex_,
// output
L1TrackPar* T_1, L1FieldRegion* fld_1, THitI* hitsl_1,
@@ -755,12 +589,6 @@ private:
L1FieldRegion vtxFieldRegion _fvecalignment; // really doesn't used
L1FieldValue vtxFieldValue _fvecalignment; // field at teh vertex position.
- // vector <L1Triplet> vTriplets; // container for triplets got in finding
- // vector<L1Triplet*> vTripletsP;
- int numPortions[12];
- vector<L1Triplet*>* TripletsLocal[MaxNStations - 2];
-
-
// int TripNumThread;
int fTrackingLevel, fGhostSuppression; // really doesn't used
diff --git a/reco/L1/L1Algo/L1Branch.h b/reco/L1/L1Algo/L1Branch.h
index 10db5245a8f9e0526a6c75eddd530606c4502623..62c05882d4bd6a1670ebc8f52fb35822bf367678 100644
--- a/reco/L1/L1Algo/L1Branch.h
+++ b/reco/L1/L1Algo/L1Branch.h
@@ -18,17 +18,15 @@
*====================================================================
*/
-#ifndef L1Branch_H
-#define L1Branch_H
+#ifndef L1Branch_h
+#define L1Branch_h
-#include "../CbmL1Def.h"
+#include "CbmL1Def.h"
#include <vector>
#include "L1StsHit.h"
-
-class L1Triplet;
-
+#include "L1Vector.h"
struct L1Branch {
L1Branch()
@@ -38,10 +36,10 @@ struct L1Branch {
, NHits(0)
, chi2(0)
, CandIndex(0)
- , StsHits("L1Branch::StsHits")
+ , fStsHits("L1Branch::fStsHits")
{
// L1Branch():Momentum(0),chi2(0),NHits(0),Lengtha(0),ista(0) , StsHits(){
- StsHits.resize(12);
+ fStsHits.reserve(25);
}
// unsigned short int Quality;
@@ -52,9 +50,7 @@ struct L1Branch {
// char Lengtha;
fscal chi2;
int CandIndex;
-
-
- L1Vector<THitI> StsHits;
+ L1Vector<THitI> fStsHits;
// static bool compareCand(const L1Branch *a, const L1Branch *b){
//
diff --git a/reco/L1/L1Algo/L1CAMergeClones.cxx b/reco/L1/L1Algo/L1CAMergeClones.cxx
index e196fdb414d1f4bbc01e6c9eacbb116945bd1ca5..cd548f0c50ac385aedf7f2f41bda1ed65752ac2d 100644
--- a/reco/L1/L1Algo/L1CAMergeClones.cxx
+++ b/reco/L1/L1Algo/L1CAMergeClones.cxx
@@ -202,51 +202,55 @@ void L1Algo::FilterTracks(fvec const r[5], fvec const C[15], fvec const m[5], fv
void L1Algo::CAMergeClones()
{
- // vector<unsigned short> FirstHit;
- // vector<unsigned short> LastHit;
- // vector<THitI> FirstHitIndex;
- // vector<THitI> LastHitIndex;
- // vector<unsigned short> Neighbour;
- // vector<float> TrackChi2;
- vector<bool> IsNext;
- vector<bool> IsUsed;
+
+ L1Vector<unsigned short>& firstStation = fMergerTrackFirstStation;
+ L1Vector<unsigned short>& lastStation = fMergerTrackLastStation;
+ L1Vector<THitI>& firstHit = fMergerTrackFirstHit;
+ L1Vector<THitI>& lastHit = fMergerTrackLastHit;
+ L1Vector<unsigned short>& neighbour = fMergerTrackNeighbour;
+ L1Vector<float>& trackChi2 = fMergerTrackChi2;
+ L1Vector<char>& isStored = fMergerTrackIsStored;
+ L1Vector<char>& isDownstreamNeighbour = fMergerTrackIsDownstreamNeighbour;
+
+ int nTracks = fTracks.size();
+
+ assert(nTracks < std::numeric_limits<unsigned short>::max());
+
+ constexpr unsigned short kNoNeighbour = std::numeric_limits<unsigned short>::max();
// vector< L1Track > vTracksNew;
- vTracksNew.clear();
- vTracksNew.reserve(NTracksIsecAll);
- // vector< THitI > vRecoHitsNew;
- vRecoHitsNew.clear();
- vRecoHitsNew.reserve(vRecoHits.size());
-
- FirstHit.resize(NTracksIsecAll);
- LastHit.resize(NTracksIsecAll);
- FirstHitIndex.resize(NTracksIsecAll);
- LastHitIndex.resize(NTracksIsecAll);
- IsUsed.resize(NTracksIsecAll);
- TrackChi2.resize(NTracksIsecAll);
- Neighbour.resize(NTracksIsecAll);
- IsNext.resize(NTracksIsecAll);
+ fMergerTracksNew.clear();
+ fMergerTracksNew.reserve(nTracks);
+ // vector< THitI > fRecoHitsNew;
+ fMergerRecoHitsNew.clear();
+ fMergerRecoHitsNew.reserve(fRecoHits.size());
+
+ firstStation.resize(nTracks);
+ lastStation.resize(nTracks);
+ firstHit.resize(nTracks);
+ lastHit.resize(nTracks);
+ isStored.resize(nTracks);
+ trackChi2.resize(nTracks);
+ neighbour.resize(nTracks);
+ isDownstreamNeighbour.resize(nTracks);
THitI start_hit = 0;
- unsigned short ista = 0;
#ifdef OMP
#pragma omp parallel for
#endif
- for (unsigned short iTr = 0; iTr < NTracksIsecAll; iTr++) {
- FirstHitIndex[iTr] = start_hit;
- ista = (*vSFlag)[(*vStsHits)[vRecoHits[start_hit]].f] / 4;
- FirstHit[iTr] = ista;
- start_hit += vTracks[iTr].NHits - 1;
- LastHitIndex[iTr] = start_hit;
- ista = (*vSFlag)[(*vStsHits)[vRecoHits[start_hit]].f] / 4;
- LastHit[iTr] = ista;
+ for (int iTr = 0; iTr < nTracks; iTr++) {
+ firstHit[iTr] = start_hit;
+ firstStation[iTr] = (*fStripFlag)[(*vStsHits)[fRecoHits[start_hit]].f] / 4;
+ start_hit += fTracks[iTr].NHits - 1;
+ lastHit[iTr] = start_hit;
+ lastStation[iTr] = (*fStripFlag)[(*vStsHits)[fRecoHits[start_hit]].f] / 4;
start_hit++;
- IsUsed[iTr] = 0;
- Neighbour[iTr] = 50000;
- TrackChi2[iTr] = 100000;
- IsNext[iTr] = 0;
+ isStored[iTr] = false;
+ neighbour[iTr] = kNoNeighbour;
+ trackChi2[iTr] = 100000;
+ isDownstreamNeighbour[iTr] = false;
}
L1KFTrackFitter();
@@ -256,130 +260,74 @@ void L1Algo::CAMergeClones()
L1TrackPar Tf;
L1FieldValue fBm, fBb, fBf _fvecalignment;
L1FieldRegion fld _fvecalignment;
+
+ unsigned char maxLengthForMerge = static_cast<unsigned char>(NStations - 3); // max length for merge
+
#ifdef OMP
#pragma omp parallel for
#endif
- for (int iTr = 0; iTr < static_cast<unsigned short>(NTracksIsecAll); iTr++) {
- if (static_cast<int>(vTracks[iTr].NHits) > (NStations - 3)) continue;
- for (int jTr = 0; jTr < static_cast<unsigned short>(NTracksIsecAll); jTr++) {
- if (static_cast<int>(vTracks[jTr].NHits) > (NStations - 3)) continue;
-
+ for (int iTr = 0; iTr < nTracks; iTr++) {
+ if (fTracks[iTr].NHits > maxLengthForMerge) continue;
+ for (int jTr = 0; jTr < nTracks; jTr++) {
+ if (fTracks[jTr].NHits > maxLengthForMerge) continue;
if (iTr == jTr) continue;
- // if(vTracks[iTr].n != vTracks[jTr].n) continue;
-
+ if (firstStation[iTr] <= lastStation[jTr]) continue;
+ // if(vTracks[iTr].n != vTracks[jTr].n) continue;
// if (fabs(vTracks[iTr].fTrackTime - vTracks[jTr].fTrackTime) > 3) continue;
- unsigned short dist = 0;
- unsigned short stab = 0, staf = 0;
- bool IsNextTemp = 0;
//if((vTracks[iTr].TFirst[4] - vTracks[jTr].TFirst[4])*(vTracks[iTr].TFirst[4] - vTracks[jTr].TFirst[4])
// > 9*(vTracks[iTr].CFirst[14]+vTracks[jTr].CFirst[14]) ) continue;
- if (FirstHit[iTr] > LastHit[jTr]) {
- dist = FirstHit[iTr] - LastHit[jTr];
-
- stab = FirstHit[iTr];
- staf = LastHit[jTr];
- IsNextTemp = 1;
-
- Tb.x = vTracks[iTr].TFirst[0];
- Tb.y = vTracks[iTr].TFirst[1];
- Tb.tx = vTracks[iTr].TFirst[2];
- Tb.ty = vTracks[iTr].TFirst[3];
- Tb.qp = vTracks[iTr].TFirst[4];
- Tb.z = vTracks[iTr].TFirst[5];
- Tb.C00 = vTracks[iTr].CFirst[0];
- Tb.C10 = vTracks[iTr].CFirst[1];
- Tb.C11 = vTracks[iTr].CFirst[2];
- Tb.C20 = vTracks[iTr].CFirst[3];
- Tb.C21 = vTracks[iTr].CFirst[4];
- Tb.C22 = vTracks[iTr].CFirst[5];
- Tb.C30 = vTracks[iTr].CFirst[6];
- Tb.C31 = vTracks[iTr].CFirst[7];
- Tb.C32 = vTracks[iTr].CFirst[8];
- Tb.C33 = vTracks[iTr].CFirst[9];
- Tb.C40 = vTracks[iTr].CFirst[10];
- Tb.C41 = vTracks[iTr].CFirst[11];
- Tb.C42 = vTracks[iTr].CFirst[12];
- Tb.C43 = vTracks[iTr].CFirst[13];
- Tb.C44 = vTracks[iTr].CFirst[14];
-
- Tf.x = vTracks[jTr].TLast[0];
- Tf.y = vTracks[jTr].TLast[1];
- Tf.tx = vTracks[jTr].TLast[2];
- Tf.ty = vTracks[jTr].TLast[3];
- Tf.qp = vTracks[jTr].TLast[4];
- Tf.z = vTracks[jTr].TLast[5];
- Tf.C00 = vTracks[jTr].CLast[0];
- Tf.C10 = vTracks[jTr].CLast[1];
- Tf.C11 = vTracks[jTr].CLast[2];
- Tf.C20 = vTracks[jTr].CLast[3];
- Tf.C21 = vTracks[jTr].CLast[4];
- Tf.C22 = vTracks[jTr].CLast[5];
- Tf.C30 = vTracks[jTr].CLast[6];
- Tf.C31 = vTracks[jTr].CLast[7];
- Tf.C32 = vTracks[jTr].CLast[8];
- Tf.C33 = vTracks[jTr].CLast[9];
- Tf.C40 = vTracks[jTr].CLast[10];
- Tf.C41 = vTracks[jTr].CLast[11];
- Tf.C42 = vTracks[jTr].CLast[12];
- Tf.C43 = vTracks[jTr].CLast[13];
- Tf.C44 = vTracks[jTr].CLast[14];
- //std::cout << "!!!!!!! Chi2 !!!!!! "<<vTracks[iTr].TFirst[0]<<" "<<vTracks[jTr].TLast[0]<<std::endl;
- }
- // if(FirstHit[jTr] > LastHit[iTr])
- // {
- // dist = FirstHit[jTr] - LastHit[iTr];
- //
- // stab = FirstHit[jTr];
- // staf = LastHit[iTr];
- //
- // Tb.x = vTracks[jTr].TFirst[0];
- // Tb.y = vTracks[jTr].TFirst[1];
- // Tb.tx = vTracks[jTr].TFirst[2];
- // Tb.ty = vTracks[jTr].TFirst[3];
- // Tb.qp = vTracks[jTr].TFirst[4];
- // Tb.z = vTracks[jTr].TFirst[5];
- // Tb.C00 = vTracks[jTr].CFirst[0];
- // Tb.C10 = vTracks[jTr].CFirst[1];
- // Tb.C11 = vTracks[jTr].CFirst[2];
- // Tb.C20 = vTracks[jTr].CFirst[3];
- // Tb.C21 = vTracks[jTr].CFirst[4];
- // Tb.C22 = vTracks[jTr].CFirst[5];
- // Tb.C30 = vTracks[jTr].CFirst[6];
- // Tb.C31 = vTracks[jTr].CFirst[7];
- // Tb.C32 = vTracks[jTr].CFirst[8];
- // Tb.C33 = vTracks[jTr].CFirst[9];
- // Tb.C40 = vTracks[jTr].CFirst[10];
- // Tb.C41 = vTracks[jTr].CFirst[11];
- // Tb.C42 = vTracks[jTr].CFirst[12];
- // Tb.C43 = vTracks[jTr].CFirst[13];
- // Tb.C44 = vTracks[jTr].CFirst[14];
- //
- // Tf.x = vTracks[iTr].TLast[0];
- // Tf.y = vTracks[iTr].TLast[1];
- // Tf.tx = vTracks[iTr].TLast[2];
- // Tf.ty = vTracks[iTr].TLast[3];
- // Tf.qp = vTracks[iTr].TLast[4];
- // Tf.z = vTracks[iTr].TLast[5];
- // Tf.C00 = vTracks[iTr].CLast[0];
- // Tf.C10 = vTracks[iTr].CLast[1];
- // Tf.C11 = vTracks[iTr].CLast[2];
- // Tf.C20 = vTracks[iTr].CLast[3];
- // Tf.C21 = vTracks[iTr].CLast[4];
- // Tf.C22 = vTracks[iTr].CLast[5];
- // Tf.C30 = vTracks[iTr].CLast[6];
- // Tf.C31 = vTracks[iTr].CLast[7];
- // Tf.C32 = vTracks[iTr].CLast[8];
- // Tf.C33 = vTracks[iTr].CLast[9];
- // Tf.C40 = vTracks[iTr].CLast[10];
- // Tf.C41 = vTracks[iTr].CLast[11];
- // Tf.C42 = vTracks[iTr].CLast[12];
- // Tf.C43 = vTracks[iTr].CLast[13];
- // Tf.C44 = vTracks[iTr].CLast[14];
- // }
-
- if (dist == 0) continue;
+
+ unsigned short stab = firstStation[iTr];
+
+ Tb.x = fTracks[iTr].TFirst[0];
+ Tb.y = fTracks[iTr].TFirst[1];
+ Tb.tx = fTracks[iTr].TFirst[2];
+ Tb.ty = fTracks[iTr].TFirst[3];
+ Tb.qp = fTracks[iTr].TFirst[4];
+ Tb.z = fTracks[iTr].TFirst[5];
+ Tb.C00 = fTracks[iTr].CFirst[0];
+ Tb.C10 = fTracks[iTr].CFirst[1];
+ Tb.C11 = fTracks[iTr].CFirst[2];
+ Tb.C20 = fTracks[iTr].CFirst[3];
+ Tb.C21 = fTracks[iTr].CFirst[4];
+ Tb.C22 = fTracks[iTr].CFirst[5];
+ Tb.C30 = fTracks[iTr].CFirst[6];
+ Tb.C31 = fTracks[iTr].CFirst[7];
+ Tb.C32 = fTracks[iTr].CFirst[8];
+ Tb.C33 = fTracks[iTr].CFirst[9];
+ Tb.C40 = fTracks[iTr].CFirst[10];
+ Tb.C41 = fTracks[iTr].CFirst[11];
+ Tb.C42 = fTracks[iTr].CFirst[12];
+ Tb.C43 = fTracks[iTr].CFirst[13];
+ Tb.C44 = fTracks[iTr].CFirst[14];
+
+ unsigned short staf = lastStation[jTr];
+
+ Tf.x = fTracks[jTr].TLast[0];
+ Tf.y = fTracks[jTr].TLast[1];
+ Tf.tx = fTracks[jTr].TLast[2];
+ Tf.ty = fTracks[jTr].TLast[3];
+ Tf.qp = fTracks[jTr].TLast[4];
+ Tf.z = fTracks[jTr].TLast[5];
+ Tf.C00 = fTracks[jTr].CLast[0];
+ Tf.C10 = fTracks[jTr].CLast[1];
+ Tf.C11 = fTracks[jTr].CLast[2];
+ Tf.C20 = fTracks[jTr].CLast[3];
+ Tf.C21 = fTracks[jTr].CLast[4];
+ Tf.C22 = fTracks[jTr].CLast[5];
+ Tf.C30 = fTracks[jTr].CLast[6];
+ Tf.C31 = fTracks[jTr].CLast[7];
+ Tf.C32 = fTracks[jTr].CLast[8];
+ Tf.C33 = fTracks[jTr].CLast[9];
+ Tf.C40 = fTracks[jTr].CLast[10];
+ Tf.C41 = fTracks[jTr].CLast[11];
+ Tf.C42 = fTracks[jTr].CLast[12];
+ Tf.C43 = fTracks[jTr].CLast[13];
+ Tf.C44 = fTracks[jTr].CLast[14];
+ //std::cout << "!!!!!!! Chi2 !!!!!! "<<fTracks[iTr].TFirst[0]<<" "<<fTracks[jTr].TLast[0]<<std::endl;
+
//if(((Tf.qp - Tb.qp)*(Tf.qp - Tb.qp)/(Tb.C44+Tf.C44))[0] > 25*10*7) continue;
if (fabs(Tf.t[0] - Tb.t[0]) > 3 * sqrt(Tf.C55[0] + Tb.C55[0])) continue;
// if (fabs (Tf.time[0] - Tb.time[0]) > 500000) continue;
@@ -387,9 +335,13 @@ void L1Algo::CAMergeClones()
vStations[staf].fieldSlice.GetFieldValue(Tf.x, Tf.y, fBf);
vStations[stab].fieldSlice.GetFieldValue(Tb.x, Tb.y, fBb);
+
+ unsigned short dist = firstStation[iTr] - lastStation[jTr];
+
if (dist > 1) stam = staf + 1;
else
stam = staf - 1;
+
fvec zm = vStations[stam].z;
fvec xm = 0.5 * (Tf.x + Tf.tx * (zm - Tf.z) + Tb.x + Tb.tx * (zm - Tb.z));
fvec ym = 0.5 * (Tb.y + Tb.ty * (zm - Tb.z) + Tb.y + Tb.ty * (zm - Tb.z));
@@ -404,56 +356,64 @@ void L1Algo::CAMergeClones()
fvec Chi2Tracks = 0.f;
FilterTracks(&(Tf.x), &(Tf.C00), &(Tb.x), &(Tb.C00), 0, 0, &Chi2Tracks);
if (Chi2Tracks[0] > 50) continue;
- if (Chi2Tracks[0] < TrackChi2[iTr] || Chi2Tracks[0] < TrackChi2[jTr]) {
- if (Neighbour[iTr] < static_cast<unsigned short>(50000)) {
- Neighbour[Neighbour[iTr]] = 50000;
- TrackChi2[Neighbour[iTr]] = 100000;
- IsNext[Neighbour[iTr]] = 0;
+
+ if (Chi2Tracks[0] < trackChi2[iTr] || Chi2Tracks[0] < trackChi2[jTr]) {
+ if (neighbour[iTr] < kNoNeighbour) {
+ neighbour[neighbour[iTr]] = kNoNeighbour;
+ trackChi2[neighbour[iTr]] = 100000;
+ isDownstreamNeighbour[neighbour[iTr]] = false;
}
- if (Neighbour[jTr] < static_cast<unsigned short>(50000)) {
- Neighbour[Neighbour[jTr]] = 50000;
- TrackChi2[Neighbour[jTr]] = 100000;
- IsNext[Neighbour[jTr]] = 0;
+ if (neighbour[jTr] < kNoNeighbour) {
+ neighbour[neighbour[jTr]] = kNoNeighbour;
+ trackChi2[neighbour[jTr]] = 100000;
+ isDownstreamNeighbour[neighbour[jTr]] = false;
}
- Neighbour[iTr] = jTr;
- Neighbour[jTr] = iTr;
- TrackChi2[iTr] = Chi2Tracks[0];
- TrackChi2[jTr] = Chi2Tracks[0];
- IsNext[iTr] = IsNextTemp;
- IsNext[jTr] = (!IsNextTemp);
+ neighbour[iTr] = jTr;
+ neighbour[jTr] = iTr;
+ trackChi2[iTr] = Chi2Tracks[0];
+ trackChi2[jTr] = Chi2Tracks[0];
+ isDownstreamNeighbour[iTr] = true;
+ isDownstreamNeighbour[jTr] = false;
}
}
}
- for (int iTr = 0; iTr < static_cast<unsigned short>(NTracksIsecAll); iTr++) {
- if (IsUsed[iTr]) continue;
-
- vTracksNew.push_back(vTracks[iTr]);
- if (!IsNext[iTr])
- for (THitI HI = FirstHitIndex[iTr]; HI <= LastHitIndex[iTr]; HI++)
- vRecoHitsNew.push_back(vRecoHits[HI]);
-
- if (Neighbour[iTr] < 50000) {
- IsUsed[Neighbour[iTr]] = 1;
- vTracksNew.back().NHits += vTracks[Neighbour[iTr]].NHits;
- for (THitI HI = FirstHitIndex[Neighbour[iTr]]; HI <= LastHitIndex[Neighbour[iTr]]; HI++)
- vRecoHitsNew.push_back(vRecoHits[HI]);
+
+ for (int iTr = 0; iTr < nTracks; iTr++) {
+ if (isStored[iTr]) continue;
+
+ fMergerTracksNew.push_back(fTracks[iTr]);
+ if (!isDownstreamNeighbour[iTr]) {
+ for (THitI HI = firstHit[iTr]; HI <= lastHit[iTr]; HI++) {
+ fMergerRecoHitsNew.push_back(fRecoHits[HI]);
+ }
}
- if (IsNext[iTr])
- for (THitI HI = FirstHitIndex[iTr]; HI <= LastHitIndex[iTr]; HI++)
- vRecoHitsNew.push_back(vRecoHits[HI]);
+ if (neighbour[iTr] < kNoNeighbour) {
+ isStored[neighbour[iTr]] = true;
+ fMergerTracksNew.back().NHits += fTracks[neighbour[iTr]].NHits;
+ for (THitI HI = firstHit[neighbour[iTr]]; HI <= lastHit[neighbour[iTr]]; HI++)
+ fMergerRecoHitsNew.push_back(fRecoHits[HI]);
+ }
+
+ if (isDownstreamNeighbour[iTr]) {
+ for (THitI HI = firstHit[iTr]; HI <= lastHit[iTr]; HI++) {
+ fMergerRecoHitsNew.push_back(fRecoHits[HI]);
+ }
+ }
}
- //vTracks.resize(vTracksNew.size());
- for (unsigned short iTr = 0; iTr < vTracksNew.size(); iTr++)
- vTracks[iTr] = vTracksNew[iTr];
- // vRecoHits.resize(vRecoHitsNew.size());
- for (THitI iHi = 0; iHi < vRecoHitsNew.size(); iHi++)
- vRecoHits[iHi] = vRecoHitsNew[iHi];
+ fTracks.resize(fMergerTracksNew.size());
+ for (unsigned int iTr = 0; iTr < fMergerTracksNew.size(); iTr++) {
+ fTracks[iTr] = fMergerTracksNew[iTr];
+ }
- NHitsIsecAll = vRecoHitsNew.size();
- NTracksIsecAll = vTracksNew.size();
+ assert(fRecoHits.size() == fMergerRecoHitsNew.size());
+ fRecoHits.resize(fMergerRecoHitsNew.size());
+ for (THitI iHi = 0; iHi < fMergerRecoHitsNew.size(); iHi++) {
+ fRecoHits[iHi] = fMergerRecoHitsNew[iHi];
+ }
+ NHitsIsecAll = fMergerRecoHitsNew.size();
//std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!! new "<<vTracksNew.size()<<" old "<< vTracks.size()<<std::endl;
}
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 3ce72a9182496576fe570cede7ce6e3bb21e7a91..4cc2d581bc3b7cd1a9ea7e3d14caff82f7f788bf 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -19,7 +19,6 @@
*/
#include "L1Algo.h"
-#include "L1Branch.h"
#include "L1Extrapolation.h"
#include "L1Filtration.h"
#include "L1Fit.h"
@@ -27,6 +26,7 @@
#include "L1Track.h"
#include "L1TrackPar.h"
//#include "TDHelper.h"
+#include "L1Branch.h"
#include "L1Grid.h"
#include "L1HitArea.h"
#include "L1Portion.h"
@@ -81,7 +81,7 @@ inline void L1Algo::f10( // input
// comment unused parameters, FU, 18.01.21
fvec* /*Event_l*/, fvec* /*d_x*/, fvec* /*d_y*/, fvec* /*d_xy*/, fvec* d_u, fvec* d_v)
{
- const Tindex& end_lh = start_lh + n1_l;
+ const Tindex end_lh = start_lh + n1_l;
for (Tindex ilh = start_lh, i1 = 0; ilh < end_lh; ++ilh, ++i1) {
@@ -403,8 +403,8 @@ inline void L1Algo::f20( // input
n2 = 0; // number of doublet
for (Tindex i1 = 0; i1 < n1; ++i1) // for each singlet
{
- const Tindex& i1_V = i1 / fvecLen;
- const Tindex& i1_4 = i1 % fvecLen;
+ const Tindex i1_V = i1 / fvecLen;
+ const Tindex i1_4 = i1 % fvecLen;
L1TrackPar& T1 = T_1[i1_V];
const int n2Saved = n2;
@@ -455,8 +455,8 @@ inline void L1Algo::f20( // input
// if (fabs(T1_new.t[i1_4]-hitm.time)>sqrt(T1_new.C55[i1_4]+hitm.timeEr*hitm.timeEr)*4) continue;
// if (fabs(T1_new.t[i1_4]-hitm.time)>sqrt(2.9*2.9)*5) continue;
- // const fscal &dt_est2 = Pick_m22[i1_4] * fabs(T1_new.C55[i1_4] + hitm.timeEr*hitm.timeEr);
- // const fscal &dt = hitm.time - T1_new.t[i1_4];
+ // const fscal dt_est2 = Pick_m22[i1_4] * fabs(T1_new.C55[i1_4] + hitm.timeEr*hitm.timeEr);
+ // const fscal dt = hitm.time - T1_new.t[i1_4];
// if ( dt*dt > dt_est2 && dt < 0 ) continue;
@@ -1081,7 +1081,6 @@ inline void L1Algo::f32( // input // TODO not updated after gaps introduction
} // f32
-/// Select triplets. Save them into vTriplets.
inline void L1Algo::f4( // input
Tindex n3, int istal, int istam, int istar, nsL1::vector<L1TrackPar>::TSimd& T_3, vector<THitI>& hitsl_3,
vector<THitI>& hitsm_3, vector<THitI>& hitsr_3,
@@ -1089,144 +1088,106 @@ inline void L1Algo::f4( // input
Tindex& nstaltriplets, vector<THitI>* /*hitsn_3*/, vector<THitI>* /*hitsr_5*/
)
{
- THitI ihitl_priv = 0;
- unsigned int Station = 0;
- unsigned int Thread = 0;
- unsigned int Triplet = 0;
+ /// Select good triplets and save them into fTriplets. Find neighbouring triplets at the next station.
- unsigned int Location = 0;
-
- unsigned char level = 0;
+#ifdef _OPENMP
+ uint Thread = omp_get_thread_num();
+#else
+ uint Thread = 0;
+#endif
+ THitI ihitl_priv = 0;
for (Tindex i3 = 0; i3 < n3; ++i3) {
- const Tindex& i3_V = i3 / fvecLen;
- const Tindex& i3_4 = i3 % fvecLen;
-
+ const Tindex i3_V = i3 / fvecLen;
+ const Tindex i3_4 = i3 % fvecLen;
L1TrackPar& T3 = T_3[i3_V];
// select
fscal& chi2 = T3.chi2[i3_4];
+ const THitI ihitl = hitsl_3[i3] + StsHitsUnusedStartIndex[istal];
+ const THitI ihitm = hitsm_3[i3] + StsHitsUnusedStartIndex[istam];
+ const THitI ihitr = hitsr_3[i3] + StsHitsUnusedStartIndex[istar];
+ L1_ASSERT(ihitl < StsHitsUnusedStopIndex[istal], ihitl << " < " << StsHitsUnusedStopIndex[istal]);
+ L1_ASSERT(ihitm < StsHitsUnusedStopIndex[istam], ihitm << " < " << StsHitsUnusedStopIndex[istam]);
+ L1_ASSERT(ihitr < StsHitsUnusedStopIndex[istar], ihitr << " < " << StsHitsUnusedStopIndex[istar]);
- Station = istal;
-
-#ifdef _OPENMP
- Thread = omp_get_thread_num();
-#else
- Thread = 0;
-#endif
-
- TripletsLocal1[Station][Thread][nTripletsThread[istal][Thread]].SetLevel(0);
-
- TripletsLocal1[Station][Thread][nTripletsThread[istal][Thread]].SetFNeighbour(0);
- TripletsLocal1[Station][Thread][nTripletsThread[istal][Thread]].SetNNeighbours(0);
-
- Triplet = nTripletsThread[istal][Thread];
-
- Location = Triplet + Station * 100000000 + Thread * 1000000;
+ uint Location = L1Triplet::PackTripletID(istal, Thread, fTriplets[istal][Thread].size());
if (ihitl_priv == 0 || ihitl_priv != hitsl_3[i3]) {
- TripForHit[0][hitsl_3[i3] + StsHitsUnusedStartIndex[istal]] = Location;
-
- TripForHit[1][hitsl_3[i3] + StsHitsUnusedStartIndex[istal]] = Location;
+ TripForHit[0][ihitl] = Location;
+ TripForHit[1][ihitl] = Location;
}
ihitl_priv = hitsl_3[i3];
-
#ifdef DO_NOT_SELECT_TRIPLETS
if (isec != TRACKS_FROM_TRIPLETS_ITERATION)
#endif
if (!finite(chi2) || chi2 < 0 || chi2 > TRIPLET_CHI2_CUT) continue;
-
+ // TODO: SG: simplify calculations for qp & Cqp below
// prepare data
fscal MaxInvMomS = MaxInvMom[0];
+ fscal scale = 255 / (MaxInvMomS * 2);
+
fscal qp = MaxInvMomS + T3.qp[i3_4];
if (qp < 0) qp = 0;
if (qp > MaxInvMomS * 2) qp = MaxInvMomS * 2;
- fscal Cqp = 5. * sqrt(fabs(T3.C44[i3_4]));
-
-
- fscal scale = 255 / (MaxInvMom[0] * 2);
-
qp = (static_cast<unsigned int>(qp * scale)) % 256;
- Cqp = (static_cast<unsigned int>(Cqp * scale)) % 256;
- Cqp += 1;
-
- if (Cqp < 0) Cqp = 0;
- if (Cqp > 20) Cqp = 20;
qp = static_cast<unsigned char>(qp);
+ fscal Cqp = 5. * sqrt(fabs(T3.C44[i3_4]));
+ Cqp = (static_cast<unsigned int>(Cqp * scale)) % 256;
+ Cqp += 1.f;
- const THitI& ihitl = hitsl_3[i3] + StsHitsUnusedStartIndex[istal];
- const THitI& ihitm = hitsm_3[i3] + StsHitsUnusedStartIndex[istam];
- const THitI& ihitr = hitsr_3[i3] + StsHitsUnusedStartIndex[istar];
- L1_ASSERT(ihitl < StsHitsUnusedStopIndex[istal], ihitl << " < " << StsHitsUnusedStopIndex[istal]);
- L1_ASSERT(ihitm < StsHitsUnusedStopIndex[istam], ihitm << " < " << StsHitsUnusedStopIndex[istam]);
- L1_ASSERT(ihitr < StsHitsUnusedStopIndex[istar], ihitr << " < " << StsHitsUnusedStopIndex[istar]);
+ if (Cqp < 0.f) Cqp = 0.f;
+ if (Cqp > 20.f) Cqp = 20.f;
- fscal& time = T3.t[i3_4];
+ //fscal& time = T3.t[i3_4];
// int n = T3.n[i3_4];
+ fTriplets[istal][Thread].emplace_back(ihitl, ihitm, ihitr, istal, istam, istar, 0, 0, 0, chi2, qp, Cqp,
+ sqrt(fabs(T3.tx[i3_4])), // TODO: SG: why sqrt(tx)???
+ sqrt(fabs(T3.C22[i3_4])), sqrt(fabs(T3.ty[i3_4])), sqrt(fabs(T3.C33[i3_4])));
- L1Triplet& tr1 = TripletsLocal1[Station][Thread][nTripletsThread[istal][Thread]];
-
-
- tr1.SetLevel(0);
-
-
- tr1.Set(ihitl, ihitm, ihitr, istal, istam, istar, 0, qp, chi2, time, Cqp, 0);
-
- tr1.tx = sqrt(fabs(T3.tx[i3_4]));
- tr1.ty = sqrt(fabs(T3.ty[i3_4]));
- tr1.Ctx = sqrt(fabs(T3.C22[i3_4]));
- tr1.Cty = sqrt(fabs(T3.C33[i3_4]));
-
+ L1Triplet& tr1 = fTriplets[istal][Thread].back();
++nstaltriplets;
-
- nTripletsThread[istal][Thread]++;
-
- Triplet = nTripletsThread[istal][Thread];
-
- Location = Triplet + Station * 100000000 + Thread * 1000000;
-
-
- TripForHit[1][hitsl_3[i3] + StsHitsUnusedStartIndex[istal]] = Location;
-
+ TripForHit[1][ihitl] = Location + 1;
if (istal > (NStations - 4)) continue;
- unsigned int Neighbours = TripForHit[1][ihitm] - TripForHit[0][ihitm];
-
+ unsigned int nNeighbours = TripForHit[1][ihitm] - TripForHit[0][ihitm];
- level = 0;
+ uint neighLocation = TripForHit[0][ihitm];
+ uint neighStation;
+ uint neighThread;
+ uint neighTriplet;
+ L1Triplet::UnpackTripletID(neighLocation, neighStation, neighThread, neighTriplet);
- for (unsigned int iN = 0; iN < Neighbours; ++iN) {
- Location = TripForHit[0][ihitm] + iN;
+ if (nNeighbours > 0) { assert(neighStation == istal + 1 || neighStation == istal + 2); }
+ unsigned char level = 0;
+ for (unsigned int iN = 0; iN < nNeighbours; ++iN, ++neighTriplet, ++neighLocation) {
- Station = Location / 100000000;
- Thread = (Location - Station * 100000000) / 1000000;
- Triplet = (Location - Station * 100000000 - Thread * 1000000);
-
- L1Triplet& curNeighbour = TripletsLocal1[Station][Thread][Triplet];
+ L1Triplet& curNeighbour = fTriplets[neighStation][neighThread][neighTriplet];
if ((curNeighbour.GetMHit() != ihitr)) continue;
- if (tr1.GetFNeighbour() == 0) tr1.SetFNeighbour(Location);
+ if (tr1.GetFNeighbour() == 0) tr1.SetFNeighbour(neighLocation);
- tr1.SetNNeighbours(Location - tr1.GetFNeighbour() + 1);
+ tr1.SetNNeighbours(neighLocation - tr1.GetFNeighbour() + 1);
- const unsigned char& jlevel = curNeighbour.GetLevel();
+ const unsigned char jlevel = curNeighbour.GetLevel();
if (level <= jlevel) level = jlevel + 1;
}
+
tr1.SetLevel(level);
}
}
@@ -1260,19 +1221,18 @@ inline void L1Algo::f5( // input
break;
}
- for (Tindex ip = 0; ip < fNThreads; ++ip) {
- for (Tindex itrip = 0; itrip < nTripletsThread[istal][ip]; ++itrip) {
- L1Triplet* trip = &(TripletsLocal1[istal][ip][itrip]);
- if (istam != trip->GetMSta()) continue;
- if (istar != trip->GetRSta()) continue;
+ for (Tindex iThread = 0; iThread < fNThreads; ++iThread) {
+ for (Tindex itrip = 0; itrip < fTriplets[istal][iThread].size(); ++itrip) {
+ L1Triplet& trip = fTriplets[istal][iThread][itrip];
+ if (istam != trip.GetMSta()) continue;
+ if (istar != trip.GetRSta()) continue;
unsigned char level = 0;
// float chi2 = trip->GetChi2();
- unsigned char qp = trip->GetQp();
- THitI ihitl = trip->GetLHit();
- THitI ihitm = trip->GetMHit();
- THitI ihitr = trip->GetRHit();
+ THitI ihitl = trip.GetLHit();
+ THitI ihitm = trip.GetMHit();
+ THitI ihitr = trip.GetRHit();
L1_ASSERT(ihitl < StsHitsUnusedStopIndex[istal], ihitl << " < " << StsHitsUnusedStopIndex[istal]);
L1_ASSERT(ihitm < StsHitsUnusedStopIndex[istam], ihitm << " < " << StsHitsUnusedStopIndex[istam]);
L1_ASSERT(ihitr < StsHitsUnusedStopIndex[istar], ihitr << " < " << StsHitsUnusedStopIndex[istar]);
@@ -1280,33 +1240,27 @@ inline void L1Algo::f5( // input
vector<unsigned int> neighCands; // save neighbour candidates
neighCands.reserve(8); // ~average is 1-2 for central, up to 5
- unsigned int Neighbours = TripForHit[1][ihitm] - TripForHit[0][ihitm];
-
- for (unsigned int iN = 0; iN < Neighbours; ++iN) {
- // for (iN = first_triplet; iN <= last_triplet; ++iN){
- int Location = TripForHit[0][ihitm] + iN;
-
+ unsigned int nNeighbours = TripForHit[1][ihitm] - TripForHit[0][ihitm];
- int Station = Location / 100000000;
- int Thread = (Location - Station * 100000000) / 1000000;
- int Triplet = (Location - Station * 100000000 - Thread * 1000000);
+ uint neighLocation = TripForHit[0][ihitm];
+ uint neighStation;
+ uint neighThread;
+ uint neighTriplet;
+ L1Triplet::UnpackTripletID(neighLocation, neighStation, neighThread, neighTriplet);
- L1Triplet& triplet = TripletsLocal1[Station][Thread][Triplet];
-
- // if (triplet.GetMSta() != istar) continue; // neighbours should have 2 common hits
- // if (triplet.GetMHit() != ihitr) continue; //!!!
-
- L1Triplet* tripn = &triplet;
+ for (unsigned int iN = 0; iN < nNeighbours; ++iN, ++neighTriplet, ++neighLocation) {
+ // for (iN = first_triplet; iN <= last_triplet; ++iN){
+ L1Triplet& neigh = fTriplets[neighStation][neighThread][neighTriplet];
+ // if (neigh.GetMSta() != istar) continue; // neighbours should have 2 common hits
+ // if (neigh.GetMHit() != ihitr) continue; //!!!
- fscal qp2 = tripn->GetQp();
- fscal Cqp1 = trip->Cqp;
- fscal Cqp2 = tripn->Cqp;
- if (fabs(qp - qp2) > PickNeighbour * (Cqp1 + Cqp2)) continue; // neighbours should have same qp
+ if (fabs(trip.GetQp() - neigh.GetQp()) > PickNeighbour * (trip.GetCqp() + neigh.GetCqp()))
+ continue; // neighbours should have same qp
// calculate level
- unsigned char jlevel = tripn->GetLevel();
+ unsigned char jlevel = neigh.GetLevel();
if (level <= jlevel) level = jlevel + 1;
- if (level == jlevel + 1) neighCands.push_back(Location);
+ if (level == jlevel + 1) neighCands.push_back(neighLocation);
}
// trip->neighbours.resize(0);
@@ -1319,7 +1273,7 @@ inline void L1Algo::f5( // input
// int Thread = (Location -Station*100000000)/1000000;
// int Triplet = (Location- Station*100000000-Thread*1000000);
- // const int nLevel = TripletsLocal1[Station][Thread][Triplet].GetLevel();
+ // const int nLevel = fTriplets[Station][Thread][Triplet].GetLevel();
// if (level == nLevel + 1) trip->neighbours.push_back(Location);
// }
nlevel[level]++;
@@ -1333,7 +1287,7 @@ inline void L1Algo::f5( // input
inline void L1Algo::
DupletsStaPort( /// creates duplets: input: @istal - start station number, @istam - last station number, @ip - index of portion, @&n_g - number of elements in portion, @*portionStopIndex
- int istal, int istam, Tindex ip, vector<Tindex>& n_g, Tindex* portionStopIndex_,
+ int istal, int istam, Tindex ip, L1Vector<Tindex>& n_g, Tindex* portionStopIndex_,
/// output:
L1TrackPar*
T_1, /// @*T_1 - singlets perameters, @*fld_1 - field aproximation, @*hitsl_1- left hits of triplets, @&lmDuplets - existance of a doublet starting from the left hit,
@@ -1670,7 +1624,8 @@ void L1Algo::CATrackFinder()
vector<L1HitPoint>* vStsHitPointsUnused_buf = &vStsDontUsedHitsxy_A;
NHitsIsecAll = 0;
- NTracksIsecAll = 0;
+ fTracks.clear();
+ fRecoHits.clear();
int nDontUsedHits = 0;
@@ -1693,7 +1648,7 @@ void L1Algo::CATrackFinder()
if ((starttime > time) && (time > 0)) starttime = time;
if (ist < NMvdStations) {
- L1StsHit& h = *(const_cast<L1StsHit*>(&((*vStsHits)[ih])));
+ L1StsHit& h = (*vStsHits)[ih];
h.t_reco = 0;
h.t_er = 100;
}
@@ -1738,9 +1693,9 @@ void L1Algo::CATrackFinder()
for (int ist = 0; ist < NStations; ++ist)
for (THitI ih = StsHitsStartIndex[ist]; ih < StsHitsStopIndex[ist]; ++ih) {
- L1StsHit& h = *(const_cast<L1StsHit*>(&((*vStsHits)[ih])));
- SetFUnUsed(const_cast<unsigned char&>((*vSFlag)[h.f]));
- SetFUnUsed(const_cast<unsigned char&>((*vSFlag)[h.b]));
+ L1StsHit& h = (*vStsHits)[ih];
+ SetFUnUsed((*fStripFlag)[h.f]);
+ SetFUnUsed((*fStripFlag)[h.b]);
}
for (int ista = 0; ista < NStations; ++ista) {
@@ -1772,15 +1727,13 @@ void L1Algo::CATrackFinder()
if (fmCBMmode)
if (isec > 1) continue;
- n_g1.assign(n_g1.size(), Portion);
-
- TripForHit[0].assign(nDontUsedHits, 0);
- TripForHit[1].assign(nDontUsedHits, 0);
-
- for (int n = 0; n < nTh; n++)
- for (int j = 0; j < NStations; j++)
- nTripletsThread[j][n] = 0;
+ // n_g1.assign(n_g1.size(), Portion);
+ for (int n = 0; n < nTh; n++) {
+ for (int j = 0; j < NStations; j++) {
+ fTriplets[j][n].clear();
+ }
+ }
/// isec - number of current iterations, fNFindIterations - number of all iterations
#ifdef COUNTERS
Tindex nSinglets = 0;
@@ -1800,6 +1753,17 @@ void L1Algo::CATrackFinder()
vStsHitPointsUnused_buf = vStsHitsUnused_temp2;
}
+ for (int ist = 0; ist < NStations; ++ist) {
+ for (THitI ih = StsHitsUnusedStartIndex[ist]; ih < StsHitsUnusedStopIndex[ist]; ++ih) {
+ //SG!!
+ TripForHit[0][ih] = 0;
+ TripForHit[1][ih] = 0;
+ }
+ }
+ /*
+ TripForHit[0].assign(StsHitsUnusedStopIndex[NStations-1],0);
+ TripForHit[1].assign(StsHitsUnusedStopIndex[NStations-1],0);
+*/
{
// #pragma omp task
{
@@ -1900,26 +1864,17 @@ void L1Algo::CATrackFinder()
{
/// possible left hits of triplets are splited in portions of 16 (4 SIMDs) to use memory faster
- portionStopIndex[NStations - 1] = 0;
- unsigned int ip = 0; //index of curent portion
-
+ fDupletPortionStopIndex[NStations - 1] = 0;
+ fDupletPortionSize.clear();
for (int istal = NStations - 2; istal >= FIRSTCASTATION; istal--) { //start downstream chambers
int NHits_l = StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal];
-
- int NHits_l_P = NHits_l / Portion;
-
- for (int ipp = 0; ipp < NHits_l_P; ipp++) {
- // n_g1[ip++] = Portion;
- n_g1[ip] = (Portion);
-
- ip++;
+ int nPortions = NHits_l / Portion;
+ int rest = NHits_l - nPortions * Portion;
+ for (int ipp = 0; ipp < nPortions; ipp++) {
+ fDupletPortionSize.push_back(Portion);
} // start_lh - portion of left hits
-
- // n_g1[ip++] = NHits_l - NHits_l_P*Portion;
- n_g1[ip] = (NHits_l - NHits_l_P * Portion);
-
- ip++;
- portionStopIndex[istal] = ip;
+ if (rest > 0) fDupletPortionSize.push_back(rest);
+ fDupletPortionStopIndex[istal] = fDupletPortionSize.size();
} // lstations
@@ -1930,7 +1885,7 @@ void L1Algo::CATrackFinder()
/* {
/// possible left hits of triplets are splited in portions of 16 (4 SIMDs) to use memory faster
- portionStopIndex[NStations-1] = 0;
+ fDupletPortionStopIndex[NStations-1] = 0;
unsigned int ip = 0; //index of curent portion
for (int istal = NStations-2; istal >= FIRSTCASTATION; istal--) //start downstream chambers
@@ -1948,7 +1903,7 @@ void L1Algo::CATrackFinder()
n_g1[ip] = nHits - NHits_P*Portion;
ip++;
- portionStopIndex[istal] = ip;
+ fDupletPortionStopIndex[istal] = ip;
}// lstations
// nPortions = ip;
} */
@@ -1990,7 +1945,7 @@ void L1Algo::CATrackFinder()
#pragma omp parallel for firstprivate(T_1, fld_1, hitsl_1, hitsm_2, i1_2, TG_1, fldG_1, hitslG_1, hitsmG_2, \
i1G_2) //schedule(dynamic, 2)
#endif
- for (Tindex ip = portionStopIndex[istal + 1]; ip < portionStopIndex[istal]; ++ip) {
+ for (Tindex ip = fDupletPortionStopIndex[istal + 1]; ip < fDupletPortionStopIndex[istal]; ++ip) {
Tindex n_2; /// number of doublets in portion
int NHitsSta = StsHitsStopIndex[istal] - StsHitsStartIndex[istal];
lmDuplets[istal].resize(NHitsSta);
@@ -2000,7 +1955,7 @@ void L1Algo::CATrackFinder()
i1_2.clear();
- DupletsStaPort(istal, istal + 1, ip, n_g1, portionStopIndex,
+ DupletsStaPort(istal, istal + 1, ip, fDupletPortionSize, fDupletPortionStopIndex,
// output
T_1, fld_1, hitsl_1,
@@ -2029,7 +1984,7 @@ void L1Algo::CATrackFinder()
i1G_2.clear();
DupletsStaPort( // input
- istal, istal + 2, ip, n_g1, portionStopIndex,
+ istal, istal + 2, ip, fDupletPortionSize, fDupletPortionStopIndex,
// output
TG_1, fldG_1, hitslG_1,
@@ -2115,11 +2070,16 @@ void L1Algo::CATrackFinder()
const unsigned char min_best_l = (ilev > min_level) ? ilev + 2 : min_level + 3; // loose maximum
- for (int i = 0; i < fNThreads; ++i)
- numberCandidateThread[i] = 0;
+ for (int i = 0; i < fNThreads; ++i) {
+ fTrackCandidates[i].clear();
+ }
- vStripToTrack.assign(vStripToTrack.size(), -1);
- vStripToTrackB.assign(vStripToTrackB.size(), -1);
+ fStripToTrack.clear();
+ fStripToTrack.resize(NStsStrips, -1);
+ fStripToTrackB.clear();
+ fStripToTrackB.resize(NStsStrips, -1);
+ fStripToTrack.assign(NStsStrips, -1);
+ fStripToTrackB.assign(NStsStrips, -1);
for (int istaF = FIRSTCASTATION; istaF <= NStations - 3 - ilev; ++istaF) {
@@ -2127,18 +2087,17 @@ void L1Algo::CATrackFinder()
#pragma omp parallel for firstprivate(curr_tr, new_tr, best_tr, curr_chi2, best_chi2, best_L, curr_L, \
ndf) // schedule(dynamic, 10)
#endif
- for (Tindex ip = 0; ip < fNThreads; ++ip) {
- for (Tindex itrip = 0; itrip < nTripletsThread[istaF][ip]; ++itrip) {
+ for (Tindex iThread = 0; iThread < fNThreads; ++iThread) {
+ for (Tindex itrip = 0; itrip < fTriplets[istaF][iThread].size(); ++itrip) {
#ifdef _OPENMP
int thread_num = omp_get_thread_num();
#else
int thread_num = 0;
#endif
- L1Triplet& first_trip = (TripletsLocal1[istaF][ip][itrip]);
-
- if (GetFUsed((*vSFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].f]
- | (*vSFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].b]))
+ L1Triplet& first_trip = (fTriplets[istaF][iThread][itrip]);
+ if (GetFUsed((*fStripFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].f]
+ | (*fStripFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].b]))
continue;
@@ -2160,7 +2119,7 @@ void L1Algo::CATrackFinder()
continue; // ghost suppression // find track with 3 hits only if it was created from a chain of triplets, but not from only one triplet
if (!fmCBMmode)
- if ((ilev == 0) && (GetFStation((*vSFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].f]) != 0))
+ if ((ilev == 0) && (GetFStation((*fStripFlag)[(*vStsHitsUnused)[first_trip.GetLHit()].f]) != 0))
continue; // ghost supression // collect only MAPS tracks-triplets CHECK!!!
}
if (first_trip.GetLevel() < ilev)
@@ -2172,10 +2131,9 @@ void L1Algo::CATrackFinder()
curr_tr.chi2 = 0.f;
// curr_tr.Lengtha = 0;
curr_tr.ista = 0;
-
- (curr_tr).StsHits[0] = ((*RealIHitP)[first_trip.GetLHit()]);
-
- (curr_tr).NHits = 1;
+ curr_tr.fStsHits.clear();
+ curr_tr.fStsHits.push_back((*RealIHitP)[first_trip.GetLHit()]);
+ curr_tr.NHits = 1;
curr_L = 1;
curr_chi2 = first_trip.GetChi2();
@@ -2208,33 +2166,33 @@ void L1Algo::CATrackFinder()
}
#endif
best_tr.Set(istaF, best_L, best_chi2, first_trip.GetQpOrig());
- L1Branch& tr = CandidatesTrack[thread_num][numberCandidateThread[thread_num]];
- tr = best_tr;
- tr.CandIndex = numberCandidateThread[thread_num] + thread_num * 10000000;
+ fTrackCandidates[thread_num].push_back(best_tr);
+ L1Branch& tr = fTrackCandidates[thread_num].back();
+ tr.CandIndex = fTrackCandidates[thread_num].size() - 1 + thread_num * 10000000;
bool check = 1;
- for (vector<THitI>::iterator phitIt = tr.StsHits.begin(); /// used strips are marked
- phitIt != tr.StsHits.begin() + tr.NHits; ++phitIt) {
+ for (L1Vector<THitI>::iterator phitIt = tr.fStsHits.begin(); /// used strips are marked
+ phitIt != tr.fStsHits.end(); ++phitIt) {
const L1StsHit& h = (*vStsHits)[*phitIt];
#ifdef _OPENMP
- omp_set_lock(&hitToBestTrackB[h.b]);
+ omp_set_lock(&fHitToBestTrackB[h.b]);
#endif
- int& strip1 = (vStripToTrackB)[h.b];
+ int& strip1 = (fStripToTrackB)[h.b];
if (strip1 != -1) {
int thread = strip1 / 10000000;
int num = (strip1 - thread * 10000000);
- if (L1Branch::compareCand(tr, CandidatesTrack[thread][num])) {
- CandidatesTrack[thread][num].CandIndex = -1;
+ if (L1Branch::compareCand(tr, fTrackCandidates[thread][num])) {
+ fTrackCandidates[thread][num].CandIndex = -1;
strip1 = tr.CandIndex;
}
else {
check = 0;
#ifdef _OPENMP
- omp_unset_lock(&hitToBestTrackB[h.b]);
+ omp_unset_lock(&fHitToBestTrackB[h.b]);
#endif
break;
}
@@ -2242,26 +2200,26 @@ void L1Algo::CATrackFinder()
else
strip1 = tr.CandIndex;
#ifdef _OPENMP
- omp_unset_lock(&hitToBestTrackB[h.b]);
+ omp_unset_lock(&fHitToBestTrackB[h.b]);
#endif
if (check) {
#ifdef _OPENMP
- omp_set_lock(&hitToBestTrackF[h.f]);
+ omp_set_lock(&fHitToBestTrackF[h.f]);
#endif
- int& strip2 = (vStripToTrack)[h.f];
+ int& strip2 = (fStripToTrack)[h.f];
if (strip2 != -1) {
int thread = strip2 / 10000000;
int num = (strip2 - thread * 10000000);
- if (L1Branch::compareCand(tr, CandidatesTrack[thread][num])) {
- CandidatesTrack[thread][num].CandIndex = -1;
+ if (L1Branch::compareCand(tr, fTrackCandidates[thread][num])) {
+ fTrackCandidates[thread][num].CandIndex = -1;
strip2 = tr.CandIndex;
}
else {
check = 0;
#ifdef _OPENMP
- omp_unset_lock(&hitToBestTrackF[h.f]);
+ omp_unset_lock(&fHitToBestTrackF[h.f]);
#endif
break;
}
@@ -2269,20 +2227,20 @@ void L1Algo::CATrackFinder()
else
strip2 = tr.CandIndex;
#ifdef _OPENMP
- omp_unset_lock(&hitToBestTrackF[h.f]);
+ omp_unset_lock(&fHitToBestTrackF[h.f]);
#endif
}
}
- if (check) numberCandidateThread[thread_num]++;
+ if (!check) { fTrackCandidates[thread_num].pop_back(); }
} // itrip
- }
- }
+ } // iThread
+ } // istaF
if (--nlevel == 0) break;
for (int i = 0; i < fNThreads; ++i) {
- SavedCand[i] = 0;
- SavedHits[i] = 0;
+ fTracks_local[i].clear();
+ fRecoHits_local[i].clear();
}
for (int i = 0; i < fNThreads; ++i) {
@@ -2291,16 +2249,16 @@ void L1Algo::CATrackFinder()
#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic, 10) firstprivate(t)
#endif
- for (Tindex iCandidate = 0; iCandidate < numberCandidateThread[i]; ++iCandidate) {
- L1Branch& tr = CandidatesTrack[i][iCandidate];
+ for (Tindex iCandidate = 0; iCandidate < fTrackCandidates[i].size(); ++iCandidate) {
+ L1Branch& tr = fTrackCandidates[i][iCandidate];
bool check = 1;
- if (CandidatesTrack[i][iCandidate].CandIndex != -1) {
- for (vector<THitI>::iterator phIt = tr.StsHits.begin(); /// used strips are marked
- phIt != tr.StsHits.begin() + tr.NHits; ++phIt) {
+ if (fTrackCandidates[i][iCandidate].CandIndex != -1) {
+ for (L1Vector<THitI>::iterator phIt = tr.fStsHits.begin(); /// used strips are marked
+ phIt != tr.fStsHits.end(); ++phIt) {
const L1StsHit& h = (((*vStsHits))[*phIt]);
- if (((vStripToTrackB)[h.b] != tr.CandIndex) || ((vStripToTrack)[h.f] != tr.CandIndex)) {
+ if (((fStripToTrackB)[h.b] != tr.CandIndex) || ((fStripToTrack)[h.f] != tr.CandIndex)) {
check = 0;
break;
}
@@ -2329,18 +2287,15 @@ void L1Algo::CATrackFinder()
#endif
- for (vector<THitI>::iterator phIt = tr.StsHits.begin(); /// used strips are marked
- phIt != tr.StsHits.begin() + tr.NHits; ++phIt) {
- L1StsHit& h = *(const_cast<L1StsHit*>(&(((*vStsHits))[*phIt])));
+ for (L1Vector<THitI>::iterator phIt = tr.fStsHits.begin(); /// used strips are marked
+ phIt != tr.fStsHits.end(); ++phIt) {
+ L1StsHit& h = (*vStsHits)[*phIt];
- SetFUsed(const_cast<unsigned char&>((*vSFlag)[h.f]));
- SetFUsed(const_cast<unsigned char&>((*vSFlag)[h.b]));
+ SetFUsed((*fStripFlag)[h.f]);
+ SetFUsed((*fStripFlag)[h.b]);
-
- vRecoHits_local[num_thread][SavedHits[num_thread]] = (*phIt);
-
- SavedHits[num_thread]++;
+ fRecoHits_local[num_thread].push_back(*phIt);
const L1StsHit& hit = (*vStsHits)[*phIt];
@@ -2359,9 +2314,7 @@ void L1Algo::CATrackFinder()
t.NHits = tr.NHits;
// t.Momentum = tr.Momentum;
t.fTrackTime = sumTime / t.NHits;
-
- vTracks_local[num_thread][SavedCand[num_thread]] = (t);
- SavedCand[num_thread]++;
+ fTracks_local[num_thread].push_back(t);
}
}
}
@@ -2375,30 +2328,33 @@ void L1Algo::CATrackFinder()
Time.Start();
#endif
-
- vector<int> offset_tracks(nTh, NTracksIsecAll);
+ int nTracks = fTracks.size();
+ vector<int> offset_tracks(nTh, nTracks);
vector<int> offset_hits(nTh, NHitsIsecAll);
- NTracksIsecAll += SavedCand[0];
- NHitsIsecAll += SavedHits[0];
+ assert(NHitsIsecAll == fRecoHits.size()); //SG!!
+ nTracks += fTracks_local[0].size();
+ NHitsIsecAll += fRecoHits_local[0].size();
for (int i = 1; i < nTh; ++i) {
- offset_tracks[i] = offset_tracks[i - 1] + SavedCand[i - 1];
- offset_hits[i] = offset_hits[i - 1] + SavedHits[i - 1];
- NTracksIsecAll += SavedCand[i];
- NHitsIsecAll += SavedHits[i];
+ offset_tracks[i] = offset_tracks[i - 1] + fTracks_local[i - 1].size();
+ offset_hits[i] = offset_hits[i - 1] + fRecoHits_local[i - 1].size();
+ nTracks += fTracks_local[i].size();
+ NHitsIsecAll += fRecoHits_local[i].size();
}
-
+ fTracks.resize(nTracks);
+ fRecoHits.resize(NHitsIsecAll);
#ifdef _OPENMP
#pragma omp parallel for
#endif
for (int i = 0; i < nTh; ++i) {
- for (Tindex iC = 0; iC < SavedCand[i]; ++iC)
- vTracks[offset_tracks[i] + iC] = (vTracks_local[i][iC]);
-
- for (Tindex iH = 0; iH < SavedHits[i]; ++iH)
- vRecoHits[offset_hits[i] + iH] = (vRecoHits_local[i][iH]);
+ for (Tindex iC = 0; iC < fTracks_local[i].size(); ++iC) {
+ fTracks[offset_tracks[i] + iC] = fTracks_local[i][iC];
+ }
+ for (Tindex iH = 0; iH < fRecoHits_local[i].size(); ++iH) {
+ fRecoHits[offset_hits[i] + iH] = fRecoHits_local[i][iH];
+ }
}
} //istaf
@@ -2418,7 +2374,7 @@ void L1Algo::CATrackFinder()
vGridTime[ista].UpdateIterGrid(
Nelements, &((*vStsHitsUnused)[StsHitsUnusedStartIndex[ista]]), RealIHitPBuf,
&((*RealIHitP)[StsHitsUnusedStartIndex[ista]]), vStsHitsUnused_buf, vStsHitPointsUnused_buf,
- &((*vStsHitPointsUnused)[StsHitsUnusedStartIndex[ista]]), NHitsOnStation, ista, *this, vSFlag);
+ &((*vStsHitPointsUnused)[StsHitsUnusedStartIndex[ista]]), NHitsOnStation, ista, *this, fStripFlag);
StsHitsUnusedStartIndex[ista] = NHitsOnStationTmp;
StsHitsUnusedStopIndex[ista] = NHitsOnStation;
}
@@ -2508,9 +2464,9 @@ void L1Algo::CATrackFinder()
NHits += vStsHitsUnused->size();
HitSize += vStsHitsUnused->size()*sizeof(L1StsHit);
NStrips+= vStsStrips.size();
- StripSize += vStsStrips.size()*sizeof(fscal) + (*vSFlag).size()*sizeof(unsigned char);
- NStripsB+= (*vSFlagB).size();
- StripSizeB += vStsStripsB.size()*sizeof(fscal) + (*vSFlagB).size()*sizeof(unsigned char);
+ StripSize += vStsStrips.size()*sizeof(fscal) + (*fStripFlag).size()*sizeof(unsigned char);
+ NStripsB+= (*fStripFlagB).size();
+ StripSizeB += vStsStripsB.size()*sizeof(fscal) + (*fStripFlagB).size()*sizeof(unsigned char);
NDup += stat_max_n_dup;
DupSize += stat_max_n_dup*sizeof(/*short*/ int);
NTrip += stat_max_n_trip;
@@ -2518,8 +2474,8 @@ void L1Algo::CATrackFinder()
NBranches += stat_max_n_branches;
BranchSize += stat_max_BranchSize;
- NTracks += vTracks.size();
- TrackSize += sizeof(L1Track)*vTracks.size() + sizeof(THitI)*vRecoHits.size();
+ NTracks += fTracks.size();
+ TrackSize += sizeof(L1Track)*fTracks.size() + sizeof(THitI)*fRecoHits.size();
int k = 1024*NTimes;
cout<<"L1 Event size: \n"
@@ -2587,21 +2543,21 @@ inline void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best
const THitI& ihitr = curr_trip->GetRHit();
- if (!GetFUsed((*vSFlag)[(*vStsHitsUnused)[ihitm].f] | (*vSFlag)[(*vStsHitsUnused)[ihitm].b])) {
+ if (!GetFUsed((*fStripFlag)[(*vStsHitsUnused)[ihitm].f] | (*fStripFlag)[(*vStsHitsUnused)[ihitm].b])) {
// curr_tr.StsHits.push_back((*RealIHitP)[ihitm]);
- curr_tr.StsHits[curr_tr.NHits] = ((*RealIHitP)[ihitm]);
+ curr_tr.fStsHits.push_back((*RealIHitP)[ihitm]);
curr_tr.NHits++;
curr_L++;
}
- if (!GetFUsed((*vSFlag)[(*vStsHitsUnused)[ihitr].f] | (*vSFlag)[(*vStsHitsUnused)[ihitr].b])) {
+ if (!GetFUsed((*fStripFlag)[(*vStsHitsUnused)[ihitr].f] | (*fStripFlag)[(*vStsHitsUnused)[ihitr].b])) {
//curr_tr.StsHits.push_back((*RealIHitP)[ihitr]);
- curr_tr.StsHits[curr_tr.NHits] = ((*RealIHitP)[ihitr]);
+ curr_tr.fStsHits.push_back((*RealIHitP)[ihitr]);
curr_tr.NHits++;
@@ -2654,33 +2610,31 @@ inline void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best
Thread = (Location - Station * 100000000) / 1000000;
Triplet = (Location - Station * 100000000 - Thread * 1000000);
-
- const L1Triplet& new_trip = TripletsLocal1[Station][Thread][Triplet];
-
+ const L1Triplet& new_trip = fTriplets[Station][Thread][Triplet];
if ((new_trip.GetMHit() != curr_trip->GetRHit())) continue;
if ((new_trip.GetLHit() != curr_trip->GetMHit())) continue;
const fscal qp1 = curr_trip->GetQp();
const fscal qp2 = new_trip.GetQp();
fscal dqp = fabs(qp1 - qp2);
- fscal Cqp = curr_trip->Cqp;
- Cqp += new_trip.Cqp;
+ fscal Cqp = curr_trip->GetCqp();
+ Cqp += new_trip.GetCqp();
if (!fmCBMmode)
if (dqp > PickNeighbour * Cqp)
continue; // bad neighbour // CHECKME why do we need recheck it?? (it really change result)
- const fscal& tx1 = curr_trip->tx;
- const fscal& tx2 = new_trip.tx;
+ fscal tx1 = curr_trip->GetTx();
+ fscal tx2 = new_trip.GetTx();
fscal dtx = fabs(tx1 - tx2);
- fscal Ctx = curr_trip->Ctx;
- Ctx += new_trip.Ctx;
+ fscal Ctx = curr_trip->GetCtx();
+ Ctx += new_trip.GetCtx();
- const fscal& ty1 = curr_trip->ty;
- const fscal& ty2 = new_trip.ty;
+ fscal ty1 = curr_trip->GetTy();
+ fscal ty2 = new_trip.GetTy();
fscal dty = fabs(ty1 - ty2);
- fscal Cty = curr_trip->Cty;
- Cty += new_trip.Cty;
+ fscal Cty = curr_trip->GetCty();
+ Cty += new_trip.GetCty();
if (fGlobal || fmCBMmode) {
if (dty > 6 * Cty) continue;
@@ -2688,8 +2642,8 @@ inline void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best
}
- if (GetFUsed((*vSFlag)[(*vStsHitsUnused)[new_trip.GetLHit()].f]
- | (*vSFlag)[(*vStsHitsUnused)[new_trip.GetLHit()].b])) { //hits are used
+ if (GetFUsed((*fStripFlag)[(*vStsHitsUnused)[new_trip.GetLHit()].f]
+ | (*fStripFlag)[(*vStsHitsUnused)[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;
@@ -2708,7 +2662,7 @@ inline void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best
fscal new_chi2 = curr_chi2;
// add new hit
- new_tr[ista].StsHits[new_tr[ista].NHits] = ((*RealIHitP)[new_trip.GetLHit()]);
+ new_tr[ista].fStsHits.push_back((*RealIHitP)[new_trip.GetLHit()]);
new_tr[ista].NHits++;
new_L += 1;
dqp = dqp / Cqp * 5.; // CHECKME: understand 5, why no sqrt(5)?
diff --git a/reco/L1/L1Algo/L1Grid.cxx b/reco/L1/L1Algo/L1Grid.cxx
index 15e839ddc5b4146c3fb99f8dd5b77484aa66540d..4dda96c73ca0cbac9ea6705e6574a06abc1966fc 100644
--- a/reco/L1/L1Algo/L1Grid.cxx
+++ b/reco/L1/L1Algo/L1Grid.cxx
@@ -39,7 +39,7 @@ inline void memset(T* dest, T i, size_t num)
void L1Grid::UpdateIterGrid(unsigned int Nelements, L1StsHit* hits, vector<THitI>* indicesBuf, THitI* indices,
vector<L1StsHit>* hits2, vector<L1HitPoint>* pointsBuf, L1HitPoint* points,
- int& NHitsOnStation, char iS, L1Algo& Algo, const vector<unsigned char>* vSFlag)
+ int& NHitsOnStation, char iS, L1Algo& Algo, const L1Vector<unsigned char>* vSFlag)
{
fFirstHitInBin.assign(fN + 2, 0);
diff --git a/reco/L1/L1Algo/L1Grid.h b/reco/L1/L1Algo/L1Grid.h
index 854de912605bc6ab3dde927fe34dbfa6e64a8965..e5b249e58c0894a45aed81669d515fed6cab2e0b 100644
--- a/reco/L1/L1Algo/L1Grid.h
+++ b/reco/L1/L1Algo/L1Grid.h
@@ -22,6 +22,7 @@
#endif
#include "L1HitPoint.h"
#include "L1StsHit.h"
+#include "L1Vector.h"
using namespace std;
@@ -84,9 +85,6 @@ public:
// }
// }
- void CreatePar(L1HitPoint* points, THitI nhits, L1Vector<L1HitPoint>* pointsBuf, L1Vector<L1StsHit>* hitsBuf,
- const L1StsHit* hits, L1Vector<THitI>* indices, L1Vector<THitI>* indicesBuf, char iS, L1Algo& Algo,
- THitI n);
void StoreHits(THitI nhits, const L1StsHit* hits, char iS, L1Algo& Algo, THitI n, L1StsHit* hitsBuf1,
const L1StsHit* hits1, THitI* indices1);
@@ -146,7 +144,7 @@ public:
void UpdateIterGrid(unsigned int Nelements, L1StsHit* hits, vector<THitI>* indicesBuf, THitI* indices,
vector<L1StsHit>* hits2, vector<L1HitPoint>* pointsBuf, L1HitPoint* points, int& NHitsOnStation,
- char iS, L1Algo& Algo, const vector<unsigned char>* vSFlag);
+ char iS, L1Algo& Algo, const L1Vector<unsigned char>* vSFlag);
private:
diff --git a/reco/L1/L1Algo/L1StsHit.h b/reco/L1/L1Algo/L1StsHit.h
index 7b31bbcf6f583423e8145fb42e3eb6e9b2fc172c..bbe8eae5ef79dca3e6e9597a928cf08a105c3347 100644
--- a/reco/L1/L1Algo/L1StsHit.h
+++ b/reco/L1/L1Algo/L1StsHit.h
@@ -2,8 +2,8 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Valentina Akishina, Igor Kulakov, Sergey Gorbunov [committer], Maksym Zyzak */
-#ifndef _L1StsHit_h_
-#define _L1StsHit_h_
+#ifndef L1StsHit_h
+#define L1StsHit_h
//struct L1Branch;
typedef unsigned /*short*/ int THitI; // strip index type
diff --git a/reco/L1/L1Algo/L1TrackExtender.cxx b/reco/L1/L1Algo/L1TrackExtender.cxx
index 0a660b145fdb5f3f4c41b90964e99853834f425d..2fc66613a1310f37aaf2b60b8524098749ed1871 100644
--- a/reco/L1/L1Algo/L1TrackExtender.cxx
+++ b/reco/L1/L1Algo/L1TrackExtender.cxx
@@ -34,8 +34,8 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
fit.SetParticleMass(GetDefaultParticleMass());
// get hits of current track
- const std::vector<THitI>& hits = t.StsHits; // array of indeses of hits of current track
- const int nHits = t.NHits;
+ const L1Vector<THitI>& hits = t.fStsHits; // array of indeses of hits of current track
+ const int nHits = t.NHits;
const signed short int step = -2 * static_cast<int>(dir) + 1; // increment for station index
const int iFirstHit = (dir) ? nHits - 1 : 0;
@@ -45,9 +45,9 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
const L1StsHit& hit1 = (*vStsHits)[hits[iFirstHit + step]];
const L1StsHit& hit2 = (*vStsHits)[hits[iFirstHit + 2 * step]];
- int ista0 = GetFStation((*vSFlag)[hit0.f]);
- int ista1 = GetFStation((*vSFlag)[hit1.f]);
- int ista2 = GetFStation((*vSFlag)[hit2.f]);
+ int ista0 = GetFStation((*fStripFlag)[hit0.f]);
+ int ista1 = GetFStation((*fStripFlag)[hit1.f]);
+ int ista2 = GetFStation((*fStripFlag)[hit2.f]);
L1Station& sta0 = vStations[ista0];
L1Station& sta1 = vStations[ista1];
@@ -125,7 +125,7 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
for (int i = iFirstHit + step; step * i <= step * iLastHit; i += step) {
const L1StsHit& hit = (*vStsHits)[hits[i]];
ista_prev = ista;
- ista = GetFStation((*vSFlag)[hit.f]);
+ ista = GetFStation((*fStripFlag)[hit.f]);
L1Station& sta = vStations[ista];
@@ -205,15 +205,15 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
const signed short int step = -2 * static_cast<int>(dir) + 1; // increment for station index
const int iFirstHit = (dir) ? 2 : t.NHits - 3;
- // int ista = GetFStation((*vSFlag)[(*vStsHits)[t.StsHits[iFirstHit]].f]) + 2*step; // current station. set to the end of track
+ // int ista = GetFStation((*fStripFlag)[(*vStsHits)[t.StsHits[iFirstHit]].f]) + 2*step; // current station. set to the end of track
- const L1StsHit& hit0 = (*vStsHits)[t.StsHits[iFirstHit]]; // optimize
- const L1StsHit& hit1 = (*vStsHits)[t.StsHits[iFirstHit + step]];
- const L1StsHit& hit2 = (*vStsHits)[t.StsHits[iFirstHit + 2 * step]];
+ const L1StsHit& hit0 = (*vStsHits)[t.fStsHits[iFirstHit]]; // optimize
+ const L1StsHit& hit1 = (*vStsHits)[t.fStsHits[iFirstHit + step]];
+ const L1StsHit& hit2 = (*vStsHits)[t.fStsHits[iFirstHit + 2 * step]];
- const int ista0 = GetFStation((*vSFlag)[hit0.f]);
- const int ista1 = GetFStation((*vSFlag)[hit1.f]);
- const int ista2 = GetFStation((*vSFlag)[hit2.f]);
+ const int ista0 = GetFStation((*fStripFlag)[hit0.f]);
+ const int ista1 = GetFStation((*fStripFlag)[hit1.f]);
+ const int ista2 = GetFStation((*fStripFlag)[hit2.f]);
const L1Station& sta0 = vStations[ista0];
const L1Station& sta1 = vStations[ista1];
@@ -284,7 +284,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
const L1StsHit& hit = (*vStsHitsUnused)[ih];
if (fabs(hit.t_reco - T.t[0]) > sqrt(T.C55[0] + hit.t_er) * 5) continue;
- if (GetFUsed((*vSFlag)[hit.f] | (*vSFlag)[hit.b])) continue; // if used
+ if (GetFUsed((*fStripFlag)[hit.f] | (*fStripFlag)[hit.b])) continue; // if used
fscal xh, yh, zh;
GetHitCoor(hit, xh, yh, zh, sta); // faster
@@ -354,23 +354,22 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
}
// save hits
+ const unsigned int NOldHits = t.NHits;
+ const unsigned int NNewHits = newHits.size();
+ t.fStsHits.resize(NNewHits + NOldHits);
+ t.NHits = (NNewHits + NOldHits);
+
if (dir) { // backward
- const unsigned int NOldHits = t.NHits;
- const unsigned int NNewHits = newHits.size();
- // t.StsHits.resize(NNewHits + NOldHits);
- t.NHits = (NNewHits + NOldHits);
for (int i = NOldHits - 1; i >= 0; i--) {
- t.StsHits[NNewHits + i] = t.StsHits[i];
+ t.fStsHits[NNewHits + i] = t.fStsHits[i];
}
for (unsigned int i = 0, ii = NNewHits - 1; i < NNewHits; i++, ii--) {
- t.StsHits[i] = newHits[ii];
+ t.fStsHits[i] = newHits[ii];
}
}
else { // forward
- const unsigned int NOldHits = t.NHits;
- t.NHits = (newHits.size() + NOldHits);
for (unsigned int i = 0; i < newHits.size(); i++) {
- t.StsHits[NOldHits + i] = (newHits[i]);
+ t.fStsHits[NOldHits + i] = (newHits[i]);
}
}
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 070c3113ca78761270de35b157ecf3822bccff49..385b84c8253d55f3b4c4fd0e7908d167234d27ab 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -30,17 +30,16 @@ using namespace NS_L1TrackFitter;
void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
{
// cout << " Start KF Track Fitter " << endl;
- int start_hit = 0; // for interation in vRecoHits[]
-
- for (unsigned int itrack = 0; itrack < NTracksIsecAll; itrack++) {
- L1Track& t = vTracks[itrack]; // current track
+ int start_hit = 0; // for interation in fRecoHits[]
+ for (unsigned int itrack = 0; itrack < fTracks.size(); itrack++) {
+ L1Track& t = fTracks[itrack]; // current track
// get hits of current track
std::vector<unsigned short int> hits; // array of indeses of hits of current track
hits.clear();
int nHits = t.NHits;
for (int i = 0; i < nHits; i++) {
- hits.push_back(vRecoHits[start_hit++]);
+ hits.push_back(fRecoHits[start_hit++]);
}
L1TrackPar T; // fitting parametr coresponding to current track
@@ -57,9 +56,9 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
const L1StsHit& hit1 = (*vStsHits)[hits[nHits - 2]];
const L1StsHit& hit2 = (*vStsHits)[hits[nHits - 3]];
- int ista0 = (*vSFlag)[hit0.f] / 4;
- int ista1 = (*vSFlag)[hit1.f] / 4;
- int ista2 = (*vSFlag)[hit2.f] / 4;
+ int ista0 = (*fStripFlag)[hit0.f] / 4;
+ int ista1 = (*fStripFlag)[hit1.f] / 4;
+ int ista2 = (*fStripFlag)[hit2.f] / 4;
//cout<<"back: ista012="<<ista0<<" "<<ista1<<" "<<ista2<<endl;
L1Station& sta0 = vStations[ista0];
@@ -128,7 +127,7 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
for (int i = nHits - 2; i >= 0; i--) {
// if( fabs(T.qp[0])>2. ) break; // iklm. Don't know it need for
const L1StsHit& hit = (*vStsHits)[hits[i]];
- ista = (*vSFlag)[hit.f] / 4;
+ ista = (*fStripFlag)[hit.f] / 4;
L1Station& sta = vStations[ista];
@@ -198,9 +197,9 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
const L1StsHit& hit1 = (*vStsHits)[hits[1]];
const L1StsHit& hit2 = (*vStsHits)[hits[2]];
- int ista0 = GetFStation((*vSFlag)[hit0.f]);
- int ista1 = GetFStation((*vSFlag)[hit1.f]);
- int ista2 = GetFStation((*vSFlag)[hit2.f]);
+ int ista0 = GetFStation((*fStripFlag)[hit0.f]);
+ int ista1 = GetFStation((*fStripFlag)[hit1.f]);
+ int ista2 = GetFStation((*fStripFlag)[hit2.f]);
L1Station& sta0 = vStations[ista0];
L1Station& sta1 = vStations[ista1];
@@ -264,7 +263,7 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
for (int i = 1; i < nHits; i++) {
const L1StsHit& hit = (*vStsHits)[hits[i]];
- ista = (*vSFlag)[hit.f] / 4;
+ ista = (*fStripFlag)[hit.f] / 4;
L1Station& sta = vStations[ista];
fvec u = hit.u;
fvec v = hit.v;
@@ -328,7 +327,7 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
void L1Algo::L1KFTrackFitter()
{
// cout << " Start L1 Track Fitter " << endl;
- int start_hit = 0; // for interation in vRecoHits[]
+ int start_hit = 0; // for interation in fRecoHits[]
// static L1FieldValue fB0, fB1, fB2 _fvecalignment;
// static L1FieldRegion fld _fvecalignment;
@@ -369,13 +368,13 @@ void L1Algo::L1KFTrackFitter()
ZSta[iHit] = sta[iHit].z;
}
- unsigned short N_vTracks = NTracksIsecAll;
+ unsigned short N_vTracks = fTracks.size();
for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
if (N_vTracks - itrack < static_cast<unsigned short>(fvecLen)) nTracks_SIMD = N_vTracks - itrack;
for (i = 0; i < nTracks_SIMD; i++)
- t[i] = &vTracks[itrack + i]; // current track
+ t[i] = &fTracks[itrack + i]; // current track
// get hits of current track
for (i = 0; i < nHits; i++) {
@@ -388,8 +387,8 @@ void L1Algo::L1KFTrackFitter()
int nHitsTrack = t[iVec]->NHits;
int iSta[MaxNStations];
for (i = 0; i < nHitsTrack; i++) {
- const L1StsHit& hit = (*vStsHits)[vRecoHits[start_hit++]];
- const int ista = (*vSFlag)[hit.f] / 4;
+ const L1StsHit& hit = (*vStsHits)[fRecoHits[start_hit++]];
+ const int ista = (*fStripFlag)[hit.f] / 4;
iSta[i] = ista;
w[ista][iVec] = 1.;
if (ista > NMvdStations) w_time[ista][iVec] = 1.;
@@ -750,7 +749,7 @@ void L1Algo::L1KFTrackFitter()
void L1Algo::L1KFTrackFitterMuch()
{
// cout << " Start L1 Track Fitter " << endl;
- int start_hit = 0; // for interation in vRecoHits[]
+ int start_hit = 0; // for interation in fRecoHits[]
// static L1FieldValue fB0, fB1, fB2 _fvecalignment;
// static L1FieldRegion fld _fvecalignment;
@@ -792,13 +791,13 @@ void L1Algo::L1KFTrackFitterMuch()
ZSta[iHit] = sta[iHit].z;
}
- unsigned short N_vTracks = NTracksIsecAll;
+ unsigned short N_vTracks = fTracks.size();
for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvecLen) {
if (N_vTracks - itrack < static_cast<unsigned short>(fvecLen)) nTracks_SIMD = N_vTracks - itrack;
for (i = 0; i < nTracks_SIMD; i++)
- t[i] = &vTracks[itrack + i]; // current track
+ t[i] = &fTracks[itrack + i]; // current track
// get hits of current track
for (i = 0; i < nHits; i++) {
@@ -808,7 +807,6 @@ void L1Algo::L1KFTrackFitterMuch()
for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
for (i = 0; i < NStations; i++) {
-
d_x[i][iVec] = 0;
d_y[i][iVec] = 0;
}
@@ -818,8 +816,8 @@ void L1Algo::L1KFTrackFitterMuch()
int nHitsTrack = t[iVec]->NHits;
int nHitsTrackSts = 0;
for (i = 0; i < nHitsTrack; i++) {
- const L1StsHit& hit = (*vStsHits)[vRecoHits[start_hit++]];
- const int ista = (*vSFlag)[hit.f] / 4;
+ const L1StsHit& hit = (*vStsHits)[fRecoHits[start_hit++]];
+ const int ista = (*fStripFlag)[hit.f] / 4;
if (ista < 8) nHitsTrackSts++;
iSta[i] = ista;
w[ista][iVec] = 1.;
diff --git a/reco/L1/L1Algo/L1Triplet.cxx b/reco/L1/L1Algo/L1Triplet.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..bfe9b9ccceec4e95932f8d9e5b9ffe20afaedd2d
--- /dev/null
+++ b/reco/L1/L1Algo/L1Triplet.cxx
@@ -0,0 +1,22 @@
+// @file L1Triplet.cxx
+// @author Sergey Gorbunov
+// @author VAlentina Akishina
+// @since 2021-05-18
+// @copyright Copyright (C) 2021 GSI Darmstadt
+// @licence SPDX-License-Identifier: GPL-3.0-only
+// @date 2021-05-18
+
+#include "L1Triplet.h"
+
+#include <iostream>
+
+void L1Triplet::Print()
+{
+ /// print the tracklet parameters
+ std::cout << " Triplet: station L/M/R " << GetLSta() << "/" << GetMSta() << "/" << GetRSta() << "\n"
+ << " hit L/M/R " << fHitL << "/" << fHitM << "/" << fHitR << "\n"
+ << " level " << fLevel << " first neighbor " << fFirstNeighbour << " Nneighbors " << fNneighbours
+ << "\n"
+ << " qp " << fQp << " Cqp " << fCqp << " chi2 " << fChi2 << "\n"
+ << " tx " << fTx << " Ctx " << fCtx << " ty " << fTy << " Cty " << fCty << std::endl;
+}
diff --git a/reco/L1/L1Algo/L1Triplet.h b/reco/L1/L1Algo/L1Triplet.h
index 50c324ba265c8ceeb5e6aadf5da182ffabd18aa7..bd7f95ca32cec907849b7f08fe3a2a87057039db 100644
--- a/reco/L1/L1Algo/L1Triplet.h
+++ b/reco/L1/L1Algo/L1Triplet.h
@@ -1,158 +1,113 @@
/* Copyright (C) 2019-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: Florian Uhlig [committer], Valentina Akishina */
+ Authors: Florian Uhlig [committer], Valentina Akishina, Sergey Gorbunov */
-#ifndef L1Triplet_H
-#define L1Triplet_H
-#include "../CbmL1Def.h"
-class L1Triplet {
-public:
- // static bool compareChain(const L1Triplet *a, const L1Triplet *b){
- // return a->chi2Chain < b->chi2Chain;
- // }
- // int sta;
- // int por;
- // int tr;
- // int loc;
-
-private:
- unsigned int first_neighbour;
- char last_neighbour;
- unsigned char bl; // level
+#ifndef L1Triplet_h
+#define L1Triplet_h
- unsigned char st; // staL (4b), staM-1-staL (2b), staR-2-staL (2b)
- fscal b1; // qp (8b)
- fscal chi2double;
+// @file L1Triplet.h
+// @author Sergey Gorbunov
+// @author Valentina Akishina
+// @date 2021-05-18
- THitI w0; // left hit (16b) index in vStsHits array
- THitI w1; // middle hit(16b)
- THitI w2; // right hit(16b)
- // unsigned char b0; // chi2(5b) + level(3b)
+#include "CbmL1Def.h"
+/// L1Triplet class represents a short 3-hit track segment called a "triplet".
+///
+class L1Triplet {
public:
- // std::vector<unsigned int> neighbours;
-
- static bool compare(const L1Triplet* a, const L1Triplet* b)
- { // sort in back order
- return ((a->GetLHit() > b->GetLHit())) || ((a->GetLHit() == b->GetLHit()) && (a->GetMHit() < b->GetMHit()))
- || ((a->GetLHit() == b->GetLHit()) && (a->GetMHit() == b->GetMHit()) && (a->GetRHit() < b->GetRHit()));
- // return a.GetLHit() > b.GetLHit();
- }
-
- static bool compare(const L1Triplet& a, const L1Triplet& b)
- { // sort in back order
- return ((a.GetLHit() > b.GetLHit())) || ((a.GetLHit() == b.GetLHit()) && (a.GetMHit() < b.GetMHit()))
- || ((a.GetLHit() == b.GetLHit()) && (a.GetMHit() == b.GetMHit()) && (a.GetRHit() < b.GetRHit()));
- // return a.GetLHit() > b.GetLHit();
- }
-
- // static bool compareLevel(const L1Triplet *a, const L1Triplet *b){ // sort in back order
- // return ( a->GetLevel() > b->GetLevel() );
- //
- // }
-
- static bool compareLevel(const L1Triplet& a, const L1Triplet& b)
- { // sort in back order
- return (a.GetLevel() > b.GetLevel());
- }
-
- static bool compareLevelT(const L1Triplet a, const L1Triplet b)
- { // sort in back order
- return (a.GetLevel() > b.GetLevel());
+ /// default constructor
+ L1Triplet() {}
+
+ /// constructor
+ L1Triplet(unsigned int iHitL, unsigned int iHitM, unsigned int iHitR, unsigned int iStaL, unsigned int iStaM,
+ unsigned int iStaR, unsigned char Level, unsigned int firstNeighbour, char nNeighbours, fscal Chi2,
+ fscal Qp, fscal Cqp, fscal tx, fscal Ctx, fscal ty, fscal Cty)
+ : fChi2(Chi2)
+ , fQp(Qp)
+ , fCqp(Cqp)
+ , fTx(tx)
+ , fCtx(Ctx)
+ , fTy(ty)
+ , fCty(Cty)
+ , fFirstNeighbour(firstNeighbour)
+ , fHitL(iHitL)
+ , fHitM(iHitM)
+ , fHitR(iHitR)
+ , fNneighbours(nNeighbours)
+ , fLevel(Level)
+ , fSta((iStaL << 4) + ((iStaM - iStaL - 1) << 2) + (iStaR - iStaL - 2))
+ {
}
+ /// Setters and getters
- fscal Cqp;
- fscal tx, ty, Ctx, Cty;
-
+ void SetLevel(unsigned char Level) { fLevel = Level; }
+ unsigned char GetLevel() const { return fLevel; }
- L1Triplet()
- : first_neighbour(0)
- , last_neighbour(0)
- , bl()
- , st()
- , b1(0)
- , chi2double(0)
- , w0()
- , w1()
- , w2()
-
- , Cqp(0)
- , tx(0)
- , ty(0)
- , Ctx(0)
- , Cty(0) {};
-
- void Set(unsigned int iHitL, unsigned int iHitM, unsigned int iHitR, unsigned int iStaL, unsigned int iStaM,
- unsigned int iStaR, unsigned char Level, fscal Qp, fscal Chi2, fscal /*time_*/ = 0, fscal _Cqp = 0,
- int /*_n*/ = 0)
- {
- w0 = iHitL;
- w1 = iHitM;
- w2 = iHitR;
-
- // Chi2 = sqrt(fabs(Chi2))*31./3.5 ; // 3.5 == 31
- chi2double = Chi2;
- // if( Chi2>31 || !finite(Chi2) ) Chi2 = 31;
- // b0 = ( (static_cast<unsigned char>( Chi2 ))<<3 ) + (Level%8);
- b1 = Qp;
- bl = Level;
-
- st = (iStaL << 4) + ((iStaM - iStaL - 1) << 2) + (iStaR - iStaL - 2);
- // time = time_;
- Cqp = _Cqp;
- // n = _n;
- }
+ THitI GetLHit() const { return fHitL; }
+ THitI GetMHit() const { return fHitM; }
+ THitI GetRHit() const { return fHitR; }
- void SetLevel(unsigned char Level) { bl = Level; }
+ void SetNNeighbours(char n) { fNneighbours = n; }
+ char GetNNeighbours() const { return fNneighbours; }
+ void SetFNeighbour(unsigned int n) { fFirstNeighbour = n; }
+ unsigned int GetFNeighbour() const { return fFirstNeighbour; }
- THitI GetLHit() const { return w0; }
- THitI GetMHit() const { return w1; }
- THitI GetRHit() const { return w2; }
- char GetNNeighbours() const { return last_neighbour; }
- void SetNNeighbours(unsigned char n) { last_neighbour = n; }
- unsigned int GetFNeighbour() const { return first_neighbour; }
- void SetFNeighbour(unsigned int n) { first_neighbour = n; }
+ fscal GetQp() const { return fQp; }
+ fscal GetChi2() const { return fChi2; }
+ fscal GetTime() const { return -111.; }
- // unsigned int GetFirstNeighbour() const {
- // return w1;
- // }
- // unsigned int GetNNeighbours() const {
- // return w2;
- // }
+ fscal GetQpOrig() const { return fQp; } //TODO: SG: return unscaled qp!!
- unsigned char GetLevel() const
- {
- // return b0%8;
- return bl;
- }
+ int GetLSta() const { return fSta >> 4; }
+ int GetMSta() const { return ((fSta % 16) >> 2) + GetLSta() + 1; }
+ int GetRSta() const { return (fSta % 4) + GetLSta() + 2; }
- fscal GetQp() const { return b1; }
+ fscal GetCqp() const { return fCqp; }
+ fscal GetTx() const { return fTx; }
+ fscal GetCtx() const { return fCtx; }
+ fscal GetTy() const { return fTy; }
+ fscal GetCty() const { return fCty; }
- fscal GetChi2() const
+ /// pack station, thread and triplet indices to an unique triplet ID
+ static uint PackTripletID(uint Station, uint Thread, uint Triplet)
{
- // float x = (b0>>3)*3.5/31.;
- // return x*x;
- return chi2double;
+ return Station * 100000000 + Thread * 1000000 + Triplet;
}
- fscal Time() const
+ /// unpack the triplet ID to its station, thread, triplet index
+ static void UnpackTripletID(uint ID, uint& Station, uint& Thread, uint& Triplet)
{
- // float x = (b0>>3)*3.5/31.;
- // return x*x;
- // return time;
- return -111.;
+ Station = ID / 100000000;
+ Thread = (ID - Station * 100000000) / 1000000;
+ Triplet = (ID - Station * 100000000 - Thread * 1000000);
}
- fscal GetQpOrig() { return b1; }
+ /// print the tracklet parameters
+ void Print();
- int GetLSta() const { return st >> 4; }
-
- int GetMSta() const { return ((st % 16) >> 2) + GetLSta() + 1; }
-
- int GetRSta() const { return (st % 4) + GetLSta() + 2; }
+private:
+ fscal fChi2 = 0.f; // chi^2
+ fscal fQp = 0.f; // q/p packed
+ fscal fCqp = 0.f; // covariance of q/p, packed
+ fscal fTx = 0.f; // tx at the left hit
+ fscal fCtx = 0.f; // covariance of tx
+ fscal fTy = 0.f; // ty at the left hit
+ fscal fCty = 0.f; // covariance of ty
+
+ unsigned int fFirstNeighbour = 0; // ID of the first neighbouring triplet
+ THitI fHitL = 0; // left hit index (16b) in vStsHits array
+ THitI fHitM = 0; // middle hit index (16b)
+ THitI fHitR = 0; // right hit index (16b)
+ char fNneighbours = 0; // n of neighbouring triplets
+ unsigned char fLevel = 0; // triplet level
+ // == its possible position on the longest track candidate it belongs to.
+ // level 0 = rightmost triplet of a track candidate
+ // level k = k-ths triplet along the track counting upstream, from right to left.
+ unsigned char fSta = 0; // packed station numbers: staL (4b), staM-1-staL (2b), staR-2-staL (2b)
};
#endif
diff --git a/reco/L1/L1Algo/L1Vector.h b/reco/L1/L1Algo/L1Vector.h
new file mode 100644
index 0000000000000000000000000000000000000000..2b512df27054cb7959d021cc2577b81c0a4863d6
--- /dev/null
+++ b/reco/L1/L1Algo/L1Vector.h
@@ -0,0 +1,177 @@
+/* Copyright (C) 2021-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov [committer] */
+
+#ifndef L1Vector_H
+#define L1Vector_H
+
+// @file L1Vector.h
+// @author Sergey Gorbunov
+// @date 2021-06-16
+
+
+#include "CbmL1Def.h"
+#ifndef FAST_CODE
+#include <FairLogger.h>
+#endif
+#include <sstream>
+
+using namespace std;
+
+/// L1Vector class is a wrapper around std::vector.
+/// It does the following:
+/// 1. gives names to vectors for better debugging
+/// 2. controls the out-of-range access in debug mode
+/// 3. supresses methods that are currently not controlled
+/// 4. warns when slow memory operations are called,
+/// i.e. when the preallocated capacity is reached and the entire vector should be copied to a new place
+/// 5. blocks usage of boolean vectors, as they have a special
+/// space-optimized but slow implementation in std. (Use L1Vector<char> instead).
+///
+template<class T>
+class L1Vector : private std::vector<T> {
+public:
+ typedef std::vector<T> Tbase;
+
+ L1Vector(const char* name = "no name") : Tbase(), fName(name) {};
+
+ L1Vector(const L1Vector& v) : Tbase(), fName(v.fName)
+ {
+ Tbase::reserve(v.capacity());
+ Tbase::assign(v.begin(), v.end());
+ }
+
+ void SetName(const std::string& s) { fName = s; }
+
+ void SetName(const std::basic_ostream<char>& s)
+ {
+ // helps to set a composed name in a single line via:
+ // SetName(std::stringstream()<<"my name "<<..whatever..);
+ fName = dynamic_cast<const std::stringstream&>(s).str();
+ }
+
+ std::string GetName() const
+ {
+ std::string s = " L1Vector<";
+ s += fName + "> ";
+ return s;
+ }
+
+ template<typename... Tinput>
+ void resize(std::size_t count, Tinput... value)
+ {
+#ifndef FAST_CODE
+ if (!Tbase::empty() && (count > Tbase::capacity())) {
+ LOG(WARNING) << "L1Vector \"" << fName << "\"::resize(" << count << "): allocated capacity of "
+ << Tbase::capacity() << " is reached, re-allocate and copy." << std::endl;
+ }
+#endif
+ Tbase::resize(count, value...);
+ }
+
+ template<typename Tinput>
+ void push_back(Tinput value)
+ {
+#ifndef FAST_CODE
+ if (Tbase::size() >= Tbase::capacity()) {
+ LOG(WARNING) << "L1Vector \"" << fName << "\"::push_back(): allocated capacity of " << Tbase::capacity()
+ << " is reached, re-allocate and copy." << std::endl;
+ }
+#endif
+ Tbase::push_back(value);
+ }
+
+ template<typename... Tinput>
+ void emplace_back(Tinput&&... value)
+ {
+#ifndef FAST_CODE
+ if (Tbase::size() >= Tbase::capacity()) {
+ LOG(WARNING) << "L1Vector \"" << fName << "\"::emplace_back(): allocated capacity of " << Tbase::capacity()
+ << " is reached, re-allocate and copy." << std::endl;
+ }
+#endif
+ Tbase::emplace_back(value...);
+ }
+
+ T& operator[](std::size_t pos)
+ {
+#ifndef FAST_CODE
+ if (pos >= Tbase::size()) {
+ LOG(ERROR) << "L1Vector \"" << fName << "\": trying to access element " << pos
+ << " outside of the vector of the size of " << Tbase::size() << std::endl;
+ exit(0);
+ }
+#endif
+ return Tbase::operator[](pos);
+ }
+
+ const T& operator[](std::size_t pos) const
+ {
+#ifndef FAST_CODE
+ if (pos >= Tbase::size()) {
+ LOG(ERROR) << "L1Vector \"" << fName << "\": trying to access element " << pos
+ << " outside of the vector of the size of " << Tbase::size() << std::endl;
+ exit(0);
+ }
+#endif
+ return Tbase::operator[](pos);
+ }
+
+ T& back()
+ {
+#ifndef FAST_CODE
+ if (Tbase::size() == 0) {
+ LOG(ERROR) << "L1Vector \"" << fName << "\": trying to access element of an empty vector" << std::endl;
+ exit(0);
+ }
+#endif
+ return Tbase::back();
+ }
+
+ const T& back() const
+ {
+#ifndef FAST_CODE
+ if (Tbase::size() == 0) {
+ LOG(ERROR) << "L1Vector \"" << fName << "\": trying to access element of an empty vector" << std::endl;
+ exit(0);
+ }
+#endif
+ return Tbase::back();
+ }
+
+ void assign(std::size_t count, const T& value)
+ {
+#ifndef FAST_CODE
+ if (!Tbase::empty() && (count > Tbase::capacity())) {
+ LOG(WARNING) << "L1Vector \"" << fName << "\"::assign(" << count << "): allocated capacity of "
+ << Tbase::capacity() << " is reached, re-allocate and copy." << std::endl;
+ }
+#endif
+ Tbase::assign(count, value);
+ }
+
+ using Tbase::begin;
+ using Tbase::capacity;
+ using Tbase::clear;
+ using Tbase::end;
+ using Tbase::iterator;
+ using Tbase::pop_back;
+ using Tbase::reserve;
+ using Tbase::size;
+
+private:
+ std::string fName;
+ using Tbase::at;
+};
+
+///
+/// std::vector<bool> has a special implementation that is space-optimized
+/// and therefore slow and not thread-safe.
+/// That is why one should use L1Vector<char> instead.
+///
+template<>
+class L1Vector<bool> { /// Make sure that L1Vector<bool> is not used
+};
+
+
+#endif
diff --git a/reco/L1/L1AlgoInputData.cxx b/reco/L1/L1AlgoInputData.cxx
index 45203ea1074295c262f113a8da161091ed05783e..3902e270e4dc6d65f06972a6c079fe895c2c3f48 100644
--- a/reco/L1/L1AlgoInputData.cxx
+++ b/reco/L1/L1AlgoInputData.cxx
@@ -37,8 +37,8 @@ void L1AlgoInputData::SetData( const vector< L1StsHit > & StsHits_,
vStsStrips.resize(StsStrips_.size());
vStsStripsB.resize(StsStripsB_.size());
vStsZPos.resize(StsZPos_.size());
- vSFlag.resize(SFlag_.size());
- vSFlagB.resize(SFlagB_.size());
+ fStripFlag.resize(SFlag_.size());
+ fStripFlagB.resize(SFlagB_.size());
for(unsigned int i=0; i<StsHits_.size(); ++i ) {vStsHits[i] = StsHits_[i];
@@ -46,8 +46,8 @@ void L1AlgoInputData::SetData( const vector< L1StsHit > & StsHits_,
for(unsigned int i=0; i<StsStrips_.size(); ++i ) vStsStrips[i] = StsStrips_[i];
for(unsigned int i=0; i<StsStripsB_.size(); ++i ) vStsStripsB[i] = StsStripsB_[i];
for(unsigned int i=0; i<StsZPos_.size(); ++i ) vStsZPos[i] = StsZPos_[i];
- for(unsigned int i=0; i<SFlag_.size(); ++i ) vSFlag[i] = SFlag_[i];
- for(unsigned int i=0; i<SFlagB_.size(); ++i ) vSFlagB[i] = SFlagB_[i];
+ for(unsigned int i=0; i<SFlag_.size(); ++i ) fStripFlag[i] = SFlag_[i];
+ for(unsigned int i=0; i<SFlagB_.size(); ++i ) fStripFlagB[i] = SFlagB_[i];
for(unsigned int i=0; i<MaxNStations+1; ++i) StsHitsStartIndex[i] = StsHitsStartIndex_[i];
@@ -87,7 +87,7 @@ bool L1AlgoInputData::ReadHitsFromFile(const char work_dir[100], const int maxNE
NStsStrips = 0;
vStsZPos.clear();
- vSFlag.clear();
+ fStripFlag.clear();
// check correct position in file
char s[] = "Event: ";
@@ -120,16 +120,17 @@ bool L1AlgoInputData::ReadHitsFromFile(const char work_dir[100], const int maxNE
cout << "vStsZPos[" << n << "]"
<< " have been read." << endl;
}
- // read algo->vSFlag
+ // read algo->fStripFlag
fadata >> n;
- // cout << n<< " vSFlagB"<<endl;
+ // cout << n<< " fStripFlagB"<<endl;
+ fStripFlag.reserve(n);
for (int i = 0; i < n; i++) {
int element;
fadata >> element;
- vSFlag.push_back(static_cast<unsigned char>(element));
+ fStripFlag.push_back(static_cast<unsigned char>(element));
}
if (iVerbose >= 4) {
- cout << "vSFlag[" << n << "]"
+ cout << "fStripFlag[" << n << "]"
<< " have been read." << endl;
}
// read algo->vStsHits
@@ -205,16 +206,16 @@ void L1AlgoInputData::PrintHits()
std::cout << vStsZPos[i] << std::endl;
}
- n = vSFlag.size();
+ n = fStripFlag.size();
std::cout << n << std::endl;
for (int i = 0; i < n; i++){
- std::cout << static_cast<int>(vSFlag[i]) << std::endl;
+ std::cout << static_cast<int>(fStripFlag[i]) << std::endl;
}
- n = vSFlagB.size();
+ n = fStripFlagB.size();
std::cout << n << std::endl;
for (int i = 0; i < n; i++){
- std::cout << static_cast<int>(vSFlagB[i]) << std::endl;
+ std::cout << static_cast<int>(fStripFlagB[i]) << std::endl;
}
n = vStsHits.size();
diff --git a/reco/L1/L1AlgoInputData.h b/reco/L1/L1AlgoInputData.h
index 5185cd29eb367a387fc0ea55d74d706f02a4b3ae..13b6e7057a65b7ed027fead54f45a33853f86a4c 100644
--- a/reco/L1/L1AlgoInputData.h
+++ b/reco/L1/L1AlgoInputData.h
@@ -12,6 +12,7 @@
#include <vector>
#include "L1StsHit.h"
+#include "L1Vector.h"
using std::istream;
using std::vector;
@@ -19,7 +20,7 @@ using std::vector;
class L1AlgoInputData {
public:
- L1AlgoInputData() : vStsHits(), NStsStrips(0), vStsZPos(), vSFlag()
+ L1AlgoInputData() : vStsHits(), NStsStrips(0), vStsZPos(), fStripFlag("L1AlgoInputData::fStripFlag")
// MaxNStations(12)
{
@@ -30,10 +31,10 @@ public:
}
~L1AlgoInputData() {};
- const vector<L1StsHit>& GetStsHits() const { return vStsHits; }
+ vector<L1StsHit>& GetStsHits() { return vStsHits; }
int GetNStsStrips() const { return NStsStrips; }
const vector<fscal>& GetStsZPos() const { return vStsZPos; }
- const L1Vector<unsigned char>& GetSFlag() const { return vSFlag; }
+ L1Vector<unsigned char>& GetSFlag() { return fStripFlag; }
const THitI* GetStsHitsStartIndex() const { return StsHitsStartIndex; }
const THitI* GetStsHitsStopIndex() const { return StsHitsStopIndex; }
@@ -59,7 +60,7 @@ public:
vStsHits.clear();
NStsStrips = 0;
vStsZPos.clear();
- //vSFlag.clear();
+ fStripFlag.clear();
{
for (int i = 0; i < MaxNStations + 1; ++i)
@@ -83,7 +84,7 @@ public:
int NStsStrips; // Number of strips in sts
vector<fscal> vStsZPos; // all possible z-positions of hits
- L1Vector<unsigned char> vSFlag; // information of hits station & using hits in tracks;
+ L1Vector<unsigned char> fStripFlag; // information of hits station & using hits in tracks;
THitI StsHitsStartIndex[MaxNStations + 1],
StsHitsStopIndex[MaxNStations + 1]; // station-bounders in vStsHits array