From ddec1f6e67b9483a5db1d46c140df7e82446fda1 Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Wed, 15 Mar 2023 19:36:07 +0000
Subject: [PATCH] L1: improve the read of unsorted data

---
 reco/L1/CMakeLists.txt                  |  1 +
 reco/L1/CbmL1ReadEvent.cxx              | 82 ++++++++++++-------------
 reco/L1/L1Algo/L1Algo.cxx               | 35 ++++++-----
 reco/L1/L1Algo/L1Algo.h                 |  5 +-
 reco/L1/L1Algo/L1CaTrackFinder.cxx      | 47 +++++++-------
 reco/L1/L1Algo/L1CaTrackFinderSlice.cxx | 11 ++--
 reco/L1/L1Algo/L1Grid.cxx               |  8 +--
 reco/L1/L1Algo/L1IODataManager.cxx      | 36 +++++------
 reco/L1/L1Algo/L1IODataManager.h        | 20 ++++--
 reco/L1/L1Algo/L1InputData.cxx          |  5 +-
 reco/L1/L1Algo/L1InputData.h            | 45 +++++++-------
 reco/L1/L1Algo/utils/L1AlgoDraw.cxx     |  9 +--
 12 files changed, 156 insertions(+), 148 deletions(-)

diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 6396001755..571ed8bbb5 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -94,6 +94,7 @@ set(NO_DICT_SRCS
 )
 
 set(HEADERS
+  CbmL1Constants.h
   CbmL1CATrdTrackFinderSA.h
   CbmL1MCPoint.h
   CbmL1Hit.h
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index a6e37619fc..fdd5f17244 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -64,23 +64,24 @@ static bool compareMcPointZ(const int& a, const int& b)
 /// The structure is used to sort hits before writing them into L1 input arrays
 ///
 struct TmpHit {
-  int iStripF;         ///< index of front strip
-  int iStripB;         ///< index of back strip
-  int iStation;        ///< index of station
-  int ExtIndex;        ///< index of hit in the external TClonesArray array (NOTE: negative for MVD)
-  double u;            ///< position of hit along axis orthogonal to front strips [cm]
-  double v;            ///< position of hit along axis orthogonal to back strips [cm]
-  double x;            ///< position of hit along x axis [cm]
-  double y;            ///< position of hit along y axis [cm]
-  double z;            ///< position of hit along z axis [cm]
-  double dx;           ///< hit position uncertainty along x axis [cm]
-  double dy;           ///< hit position uncertainty along y axis [cm]
-  double dxy;          ///< hit position covariance along x and y axes [cm2]
-  double du;           ///< hit position uncertainty along axis orthogonal to front strips [cm]
-  double dv;           ///< hit position uncertainty along axis orthogonal to back strips [cm]
-  int iMC;             ///< index of MCPoint in the fvMCPoints array
-  double time = 0.;    ///< time of the hit [ns]
-  double dt   = 1.e4;  ///< time error of the hit [ns]
+  int iStripF;          ///< index of front strip
+  int iStripB;          ///< index of back strip
+  int iStation;         ///< index of station
+  int64_t fDataStream;  ///< global index of detector module
+  int ExtIndex;         ///< index of hit in the external TClonesArray array (NOTE: negative for MVD)
+  double u;             ///< position of hit along axis orthogonal to front strips [cm]
+  double v;             ///< position of hit along axis orthogonal to back strips [cm]
+  double x;             ///< position of hit along x axis [cm]
+  double y;             ///< position of hit along y axis [cm]
+  double z;             ///< position of hit along z axis [cm]
+  double dx;            ///< hit position uncertainty along x axis [cm]
+  double dy;            ///< hit position uncertainty along y axis [cm]
+  double dxy;           ///< hit position covariance along x and y axes [cm2]
+  double du;            ///< hit position uncertainty along axis orthogonal to front strips [cm]
+  double dv;            ///< hit position uncertainty along axis orthogonal to back strips [cm]
+  int iMC;              ///< index of MCPoint in the fvMCPoints array
+  double time = 0.;     ///< time of the hit [ns]
+  double dt   = 1.e4;   ///< time error of the hit [ns]
   int Det;
 
   /// Creates a hit from the CbmL1MCPoint object
@@ -93,12 +94,12 @@ struct TmpHit {
   void CreateHitFromPoint(const CbmL1MCPoint& point, int ip, int det, int& NStrips, const L1Station& st, double du_,
                           double dv_, double dt_, bool doSmear)
   {
-    ExtIndex = 0;
-    Det      = det;
-    iStation = point.iStation;
-
-    dt   = dt_;
-    time = point.time;
+    ExtIndex    = 0;
+    Det         = det;
+    iStation    = point.iStation;
+    fDataStream = (static_cast<int64_t>(Det) << 60);
+    dt          = dt_;
+    time        = point.time;
 
     iStripF = NStrips;
     iStripB = iStripF;
@@ -663,10 +664,10 @@ void CbmL1::ReadEvent(CbmEvent* event)
       Int_t hitIndex = (readFromEvent ? event->GetIndex(ECbmDataType::kMvdHit, j) : j);
 
       TmpHit th;
+      CbmMvdHit* h = L1_DYNAMIC_CAST<CbmMvdHit*>(fpMvdHits->At(hitIndex));
       {
-        CbmMvdHit* h = L1_DYNAMIC_CAST<CbmMvdHit*>(fpMvdHits->At(hitIndex));
-        th.ExtIndex  = hitIndex;
-        th.iStation  = h->GetStationNr();
+        th.ExtIndex = hitIndex;
+        th.iStation = h->GetStationNr();
 
         int stIdx = fpAlgo->GetParameters()->GetStationIndexActive(h->GetStationNr(), L1DetectorID::kMvd);
         if (stIdx == -1) continue;
@@ -699,6 +700,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
       if (1 == fMvdUseMcHit && th.iMC >= 0) {
         th.SetHitFromPoint(fvMCPoints[th.iMC], th.iMC, fpAlgo->GetParameters()->GetStation(th.iStation));
       }
+      th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
       //if(  th.iMC >=0 ) // DEBUG !!!!
       {
         tmpHits.push_back(th);
@@ -739,13 +741,13 @@ void CbmL1::ReadEvent(CbmEvent* event)
       // ***********************************
 
       TmpHit th;
+      CbmStsHit* h = L1_DYNAMIC_CAST<CbmStsHit*>(fpStsHits->At(hitIndex));
 
       // Fill reconstructed information
       {
-        CbmStsHit* h = L1_DYNAMIC_CAST<CbmStsHit*>(fpStsHits->At(hitIndex));
-        th.ExtIndex  = hitIndex;
-        th.Det       = 1;
-        int stIdx    = fpAlgo->GetParameters()->GetStationIndexActive(
+        th.ExtIndex = hitIndex;
+        th.Det      = 1;
+        int stIdx   = fpAlgo->GetParameters()->GetStationIndexActive(
           CbmStsSetup::Instance()->GetStationNumber(h->GetAddress()), L1DetectorID::kSts);
 
         if (stIdx == -1) continue;
@@ -786,7 +788,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
       if (1 == fStsUseMcHit && th.iMC >= 0) {
         th.SetHitFromPoint(fvMCPoints[th.iMC], th.iMC, fpAlgo->GetParameters()->GetStation(th.iStation));
       }
-
+      th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
       tmpHits.push_back(th);
       nStsHits++;
 
@@ -819,11 +821,9 @@ void CbmL1::ReadEvent(CbmEvent* event)
 
     for (int j = 0; j < nEnt; j++) {
       TmpHit th;
+      Int_t hitIndex     = (event ? event->GetIndex(ECbmDataType::kMuchPixelHit, j) : j);
+      CbmMuchPixelHit* h = static_cast<CbmMuchPixelHit*>(fpMuchPixelHits->At(hitIndex));
       {
-        Int_t hitIndex = (event ? event->GetIndex(ECbmDataType::kMuchPixelHit, j) : j);
-
-        CbmMuchPixelHit* h = static_cast<CbmMuchPixelHit*>(fpMuchPixelHits->At(hitIndex));
-
         th.ExtIndex = hitIndex;
         th.Det      = 2;
 
@@ -889,7 +889,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
       //    }
       //  }
       //}
-
+      th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
       tmpHits.push_back(th);
       nMuchHits++;
 
@@ -912,7 +912,6 @@ void CbmL1::ReadEvent(CbmEvent* event)
       double du = 100.e-4;
       double dt = 3.9;
       th.CreateHitFromPoint(p, ip, DetId, NStrips, fpAlgo->GetParameters()->GetStation(p.iStation), du, du, dt, true);
-
       tmpHits.push_back(th);
       nMuchHits++;
     }
@@ -935,8 +934,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
       TmpHit th;
 
       Int_t hitIndex = (event ? event->GetIndex(ECbmDataType::kTrdHit, iHit) : iHit);
-
-      CbmTrdHit* h = L1_DYNAMIC_CAST<CbmTrdHit*>(fpTrdHits->At(hitIndex));
+      CbmTrdHit* h   = L1_DYNAMIC_CAST<CbmTrdHit*>(fpTrdHits->At(hitIndex));
 
       if ((L1Algo::TrackingMode::kGlobal == fTrackingMode) && (int) h->GetClassType() != 1) {
         // SGtrd2d!! skip TRD-1D hit
@@ -1014,7 +1012,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
         }
         if (mcP < 1.) continue;
       }
-
+      th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
       tmpHits.push_back(th);
       nTrdHits++;
     }  // for fpTrdHits
@@ -1104,7 +1102,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
       if (1 == fTofUseMcHit && th.iMC > -1) {
         th.SetHitFromPoint(fvMCPoints[th.iMC], th.iMC, fpAlgo->GetParameters()->GetStation(th.iStation));
       }
-
+      th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
       tmpHits.push_back(th);
       nTofHits++;
 
@@ -1195,7 +1193,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
 
     // TODO: Here one should fill in the fvExternalHits[iHit].mcPointIds
 
-    fIODataManager.PushBackHit(h);
+    fIODataManager.PushBackHit(h, th.fDataStream);
 
     fvHitDebugInfo.push_back(s);
     fvHitPointIndexes.push_back(th.iMC);
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 4980c5ef9a..4477cf7128 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -91,7 +91,10 @@ void L1Algo::ReceiveInputData(L1InputData&& inputData)
 //
 void L1Algo::ResetSliceData()
 {
-  int nHits = fSliceHitIds.size();
+  int nHits = 0;
+  for (int iS = 0; iS < fParameters.GetNstationsActive(); iS++) {
+    nHits += fSliceHitIds[iS].size();
+  }
 
   fGridHits.reset(nHits);
   fGridHitsBuf.reset(nHits);
@@ -120,23 +123,23 @@ void L1Algo::ResetSliceData()
   fSliceRecoHits.clear();
   fSliceRecoHits.reserve(2 * nHits);
 
-  for (int i = 0; i < L1Constants::size::kMaxNstations; i++) {
-    int nHitsStation = fSliceHitIdsStopIndex[i] - fSliceHitIdsStartIndex[i] + 1;
-    fSingletPortionSize[i].clear();
-    fSingletPortionSize[i].reserve(2 * nHitsStation);
+  for (int iS = 0; iS < L1Constants::size::kMaxNstations; iS++) {
+    int nHitsStation = fSliceHitIds[iS].size();
+    fSingletPortionSize[iS].clear();
+    fSingletPortionSize[iS].reserve(2 * nHitsStation);
   }
 
-  for (int i = 0; i < fNThreads; i++) {
-    fSliceRecoTracks_thread[i].clear();
-    fSliceRecoTracks_thread[i].reserve(nHits / 10);
-    fSliceRecoHits_thread[i].clear();
-    fSliceRecoHits_thread[i].reserve(nHits);
-    fTrackCandidates[i].clear();
-    fTrackCandidates[i].reserve(nHits / 10);
-    for (unsigned int j = 0; j < L1Constants::size::kMaxNstations; j++) {
-      int nHitsStation = fSliceHitIdsStopIndex[i] - fSliceHitIdsStartIndex[i] + 1;
-      fTriplets[j][i].clear();
-      fTriplets[j][i].reserve(2 * nHitsStation);
+  for (int iT = 0; iT < fNThreads; iT++) {
+    fSliceRecoTracks_thread[iT].clear();
+    fSliceRecoTracks_thread[iT].reserve(nHits / 10);
+    fSliceRecoHits_thread[iT].clear();
+    fSliceRecoHits_thread[iT].reserve(nHits);
+    fTrackCandidates[iT].clear();
+    fTrackCandidates[iT].reserve(nHits / 10);
+    for (unsigned int iS = 0; iS < L1Constants::size::kMaxNstations; iS++) {
+      int nHitsStation = fSliceHitIds[iS].size();
+      fTriplets[iS][iT].clear();
+      fTriplets[iS][iT].reserve(2 * nHitsStation);
     }
   }
 }
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index 4e8f2dcfbb..3c89d2cf89 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -441,9 +441,8 @@ public:
     "L1Algo::fSingletPortionSize"};  ///< Number of doublets in a portion
 
 
-  L1Vector<L1HitIndex_t> fSliceHitIds {"L1Algo::fSliceHitIds"};                   ///< indices of the sub-slice hits
-  L1HitIndex_t fSliceHitIdsStartIndex[L1Constants::size::kMaxNstations + 1] {0};  ///< start of station hit inices
-  L1HitIndex_t fSliceHitIdsStopIndex[L1Constants::size::kMaxNstations + 1] {0};   ///< stop of station hit inices
+  ///< indices of the sub-slice hits
+  L1Vector<L1HitIndex_t> fSliceHitIds[L1Constants::size::kMaxNstations] {"L1Algo::fSliceHitIds"};
 
   L1Vector<L1Hit> fGridHits {"L1Algo::fGridHits"};        ///< hits, ordered with respect to their grid bins
   L1Vector<L1Hit> fGridHitsBuf {"L1Algo::fGridHitsBuf"};  ///< hits, ordered with respect to their grid bins
diff --git a/reco/L1/L1Algo/L1CaTrackFinder.cxx b/reco/L1/L1Algo/L1CaTrackFinder.cxx
index 846ec6d4e0..eba4866229 100644
--- a/reco/L1/L1Algo/L1CaTrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CaTrackFinder.cxx
@@ -41,8 +41,10 @@ void L1Algo::CaTrackFinder()
   fRecoHits.reserve(2 * fInputData.GetNhits());
   fRecoTracks.reserve(2 * fInputData.GetNhits() / fParameters.GetNstationsActive());
 
-  fSliceHitIds.clear();
-  fSliceHitIds.reserve(fInputData.GetNhits());
+  for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
+    fSliceHitIds[iS].clear();
+    fSliceHitIds[iS].reserve(fInputData.GetNhits());
+  }
 
   fvHitKeyFlags.reset(fInputData.GetNhitKeys(), 0);
   fHitTimeInfo.reset(fInputData.GetNhits());
@@ -60,18 +62,19 @@ void L1Algo::CaTrackFinder()
 
   // calculate possible event time for the hits (fHitTimeInfo array)
 
-  for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
+  const int nDataStreams = fInputData.GetNdataStreams();
 
-    const L1Station& st = fParameters.GetStation(iS);
-    int nStationHits    = fInputData.GetNhits(iS);
+  for (int iStream = 0; iStream < nDataStreams; ++iStream) {
 
     fscal maxTimeBeforeHit = std::numeric_limits<fscal>::min();
 
-    // loop in the reverse order to fill L1HitTimeInfo::fMinTimeAfterHit fields
+    int nStreamHits = fInputData.GetStreamNhits(iStream);
+
+    for (int ih = 0; ih < nStreamHits; ++ih) {
 
-    for (int ih = 0; ih < nStationHits; ++ih) {
-      L1HitIndex_t caHitId = fInputData.GetStartHitIndex(iS) + ih;
+      L1HitIndex_t caHitId = fInputData.GetStreamStartIndex(iStream) + ih;
       const L1Hit& h       = fInputData.GetHit(caHitId);
+      const L1Station& st  = fParameters.GetStation(h.iSt);
 
       auto [x, y] = st.ConvUVtoXY(h.u, h.v);
       fscal dx    = x - targX;
@@ -98,9 +101,10 @@ void L1Algo::CaTrackFinder()
     }
 
     fscal minTimeAfterHit = std::numeric_limits<fscal>::max();
+    // loop in the reverse order to fill L1HitTimeInfo::fMinTimeAfterHit fields
 
-    for (int ih = nStationHits - 1; ih >= 0; --ih) {
-      L1HitIndex_t caHitId = fInputData.GetStartHitIndex(iS) + ih;
+    for (int ih = nStreamHits - 1; ih >= 0; --ih) {
+      L1HitIndex_t caHitId = fInputData.GetStreamStartIndex(iStream) + ih;
       L1HitTimeInfo& info  = fHitTimeInfo[caHitId];
       if (minTimeAfterHit > info.fEventTimeMin) { minTimeAfterHit = info.fEventTimeMin; }
       info.fMinTimeAfterHit = minTimeAfterHit;
@@ -112,10 +116,10 @@ void L1Algo::CaTrackFinder()
   bool areDataLeft = true;  // is the whole TS processed
   int nSubSlices   = 0;
 
-  L1HitIndex_t sliceFirstHit[L1Constants::size::kMaxNstations] {0};
+  L1HitIndex_t sliceFirstHit[nDataStreams];
 
-  for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
-    sliceFirstHit[iS] = fInputData.GetStartHitIndex(iS);
+  for (int iStream = 0; iStream < nDataStreams; ++iStream) {
+    sliceFirstHit[iStream] = fInputData.GetStreamStartIndex(iStream);
   }
 
   while (areDataLeft) {
@@ -123,16 +127,17 @@ void L1Algo::CaTrackFinder()
     nSubSlices++;
 
     // select the sub-slice hits
-    fSliceHitIds.clear();
+    for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
+      fSliceHitIds[iS].clear();
+    }
 
     areDataLeft = false;
 
-    // TODO: skip empty regions and start the subslice with the earliest hit
+    // TODO: SG: skip empty regions and start the subslice with the earliest hit
 
-    for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
-      fSliceHitIdsStartIndex[iS] = fSliceHitIds.size();
+    for (int iStream = 0; iStream < nDataStreams; ++iStream) {
 
-      for (L1HitIndex_t caHitId = sliceFirstHit[iS]; caHitId < fInputData.GetStopHitIndex(iS); ++caHitId) {
+      for (L1HitIndex_t caHitId = sliceFirstHit[iStream]; caHitId < fInputData.GetStreamStopIndex(iStream); ++caHitId) {
         L1HitTimeInfo& info = fHitTimeInfo[caHitId];
         const L1Hit& h      = fInputData.GetHit(caHitId);
         if (fvHitKeyFlags[h.f] || fvHitKeyFlags[h.b]) {  // the hit is already reconstructed
@@ -147,16 +152,14 @@ void L1Algo::CaTrackFinder()
         }
         else {
           if (tsStart <= info.fEventTimeMax) {  // the hit belongs to the sub-slice
-            fSliceHitIds.push_back(caHitId);
+            fSliceHitIds[h.iSt].push_back(caHitId);
             if (info.fMaxTimeBeforeHit < tsStart + tsLength - tsOverlap) {
               // this hit and all hits before are before the overlap
-              sliceFirstHit[iS] = caHitId + 1;
+              sliceFirstHit[iStream] = caHitId + 1;
             }
           }
         }  // else the hit has been alread processed in previous sub-slices
       }
-
-      fSliceHitIdsStopIndex[iS] = fSliceHitIds.size();
     }
 
     CaTrackFinderSlice();
diff --git a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
index f57dc23beb..360cb31b1f 100644
--- a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
+++ b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
@@ -1552,7 +1552,9 @@ void L1Algo::CaTrackFinderSlice()
   static Tindex stat_nCalls[fNFindIterations]                                        = {0};  // n calls of CAFindTrack
   static Tindex stat_nTrCandidates[fNFindIterations]                                 = {0};
 
-  stat_nStartHits += fSliceHitIds.size();
+  for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
+    stat_nStartHits += fSliceHitIds[iS].size();
+  }
 #endif
 
   /********************************/ /**
@@ -1567,7 +1569,6 @@ void L1Algo::CaTrackFinderSlice()
 
 
   L1HitIndex_t nGridHitsFilled = 0;
-
   for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
     const L1Station& st    = fParameters.GetStation(iS);
     fMaxDx[iS]             = 0.;
@@ -1583,8 +1584,8 @@ void L1Algo::CaTrackFinderSlice()
     fscal gridMinY = -0.1;
     fscal gridMaxY = 0.1;
 
-    for (L1HitIndex_t ih = fSliceHitIdsStartIndex[iS]; ih < fSliceHitIdsStopIndex[iS]; ++ih) {
-      const L1Hit& h = fInputData.GetHit(fSliceHitIds[ih]);
+    for (L1HitIndex_t ih = 0; ih < fSliceHitIds[iS].size(); ++ih) {
+      const L1Hit& h = fInputData.GetHit(fSliceHitIds[iS][ih]);
 
       auto [dxx, dxy, dyy] = st.FormXYCovarianceMatrix(h.du2, h.dv2);
       fscal dx             = sqrt(dxx);
@@ -1621,7 +1622,7 @@ void L1Algo::CaTrackFinderSlice()
       vGridTime[iS].BuildBins(-1, 1, -0.6, 0.6, starttime, lasttime, xStep, yStep, (lasttime - starttime + 1));
     */
 
-    int nSliceHits = fSliceHitIdsStopIndex[iS] - fSliceHitIdsStartIndex[iS];
+    int nSliceHits = fSliceHitIds[iS].size();
     fscal sizeY    = gridMaxY - gridMinY;
     fscal sizeX    = gridMaxX - gridMinX;
     int nBins2D    = 1 + nSliceHits;
diff --git a/reco/L1/L1Algo/L1Grid.cxx b/reco/L1/L1Algo/L1Grid.cxx
index 835ed39f97..67b2208f98 100644
--- a/reco/L1/L1Algo/L1Grid.cxx
+++ b/reco/L1/L1Algo/L1Grid.cxx
@@ -167,9 +167,7 @@ void L1Grid::BuildBins(float xMin, float xMax, float yMin, float yMax, float tMi
 
 void L1Grid::StoreHits(L1Algo& algo, int iS, L1HitIndex_t& nGridHitsFilled)
 {
-
-  L1HitIndex_t firstHitIndex = algo.fSliceHitIdsStartIndex[iS];
-  L1HitIndex_t nHits         = algo.fSliceHitIdsStopIndex[iS] - firstHitIndex;
+  L1HitIndex_t nHits = algo.fSliceHitIds[iS].size();
 
   algo.fGridHitStartIndex[iS] = nGridHitsFilled;
 
@@ -179,7 +177,7 @@ void L1Grid::StoreHits(L1Algo& algo, int iS, L1HitIndex_t& nGridHitsFilled)
 #pragma omp parallel for
 #endif
   for (L1HitIndex_t ih = 0; ih < nHits; ih++) {
-    L1HitIndex_t caHitId = algo.fSliceHitIds[firstHitIndex + ih];
+    L1HitIndex_t caHitId = algo.fSliceHitIds[iS][ih];
     const L1Hit& h       = algo.GetInputData().GetHit(caHitId);
     fscal x, y;
     std::tie(x, y) = algo.GetHitCoorOnGrid(h, iS);
@@ -218,7 +216,7 @@ void L1Grid::StoreHits(L1Algo& algo, int iS, L1HitIndex_t& nGridHitsFilled)
 
 #pragma omp parallel for
   for (L1HitIndex_t ih = 0; ih < nHits; ih++) {
-    L1HitIndex_t caHitId = algo.fSliceHitIds[firstHitIndex + ih];
+    L1HitIndex_t caHitId = algo.fSliceHitIds[iS][ih];
     const L1Hit& h       = algo.GetInputData().GetHit(caHitId);
     fscal x, y;
     std::tie(x, y) = algo.GetHitCoorOnGrid(h, iS);
diff --git a/reco/L1/L1Algo/L1IODataManager.cxx b/reco/L1/L1Algo/L1IODataManager.cxx
index 1c389b282b..dbe5a6a95d 100644
--- a/reco/L1/L1Algo/L1IODataManager.cxx
+++ b/reco/L1/L1Algo/L1IODataManager.cxx
@@ -42,7 +42,7 @@ bool L1IODataManager::SendInputData(L1Algo* pAlgo)
 void L1IODataManager::ReadInputData(const std::string& fileName)
 {
   // Reset input data object
-  this->ResetInputData();
+  ResetInputData();
   LOG(info) << "L1: Input data will be read from file \"" << fileName << "\"";
 
   // Open input binary file
@@ -65,37 +65,29 @@ void L1IODataManager::ResetInputData() noexcept
 {
   L1InputData tmp;
   fInputData.Swap(tmp);
+  fLastStreamId = -1;
+  fInputData.fStreamStartIndices.reserve(2000);
+  fInputData.fStreamStopIndices.reserve(2000);
 }
 
 // ---------------------------------------------------------------------------------------------------------------------
 //
 void L1IODataManager::InitData()
 {
-  // Sort hits within stations
-  L1Vector<L1Hit> fvHitsNew;
-  fvHitsNew.reset(fInputData.fvHits.size());
+  // set the end pointers to data streams
 
-  std::fill(fvHitIndex.begin(), fvHitIndex.end(), 0);
-  std::fill(fInputData.fvStartHitIndexes.begin(), fInputData.fvStartHitIndexes.end(), 0);
+  //std::cout << "N data streams: " << fInputData.fStreamStartIndices.size() << std::endl;
 
-  // Count number of hits in each station
-  for (const auto& hit : fInputData.fvHits) {
-    ++fInputData.fvStartHitIndexes[hit.iSt + 1];
+  if (fInputData.fStreamStartIndices.size() > 3000) {
+    LOG(warning) << "L1: unexpected order of input data: too many data streams!!! ";
+    fInputData.fStreamStartIndices.reduce(3000);
   }
-
-  // Fill bordering numbers of hits for each station
-  for (int iSt = 0; iSt < fpParameters->GetNstationsActive(); ++iSt) {
-    fInputData.fvStartHitIndexes[iSt + 1] += fInputData.fvStartHitIndexes[iSt];
-  }
-
-  // Save ordered hits to new vector
-  for (const auto& hit : fInputData.fvHits) {
-    int iSt                                                            = hit.iSt;
-    fvHitsNew[fInputData.fvStartHitIndexes[iSt] + (fvHitIndex[iSt]++)] = hit;
+  int nStreams = fInputData.fStreamStartIndices.size();
+  fInputData.fStreamStopIndices.reset(nStreams);
+  for (int i = 0; i < nStreams - 1; i++) {
+    fInputData.fStreamStopIndices[i] = fInputData.fStreamStartIndices[i + 1];
   }
-
-  // Swap contents of old and new hits vector
-  std::swap(fvHitsNew, fInputData.fvHits);
+  if (nStreams > 0) { fInputData.fStreamStopIndices[nStreams - 1] = fInputData.fHits.size(); }
 }
 
 // ---------------------------------------------------------------------------------------------------------------------
diff --git a/reco/L1/L1Algo/L1IODataManager.h b/reco/L1/L1Algo/L1IODataManager.h
index c744493b46..05c2c3f131 100644
--- a/reco/L1/L1Algo/L1IODataManager.h
+++ b/reco/L1/L1Algo/L1IODataManager.h
@@ -56,14 +56,23 @@ public:
   /// Reserve number of hits
   /// \param  nHits  Number of hits to be stored
   /// \note   If one does not call this method, the underlying vector of hits will be filled with the time penalty
-  void ReserveNhits(L1HitIndex_t nHits) { fInputData.fvHits.reserve(nHits); }
+  void ReserveNhits(L1HitIndex_t nHits) { fInputData.fHits.reserve(nHits); }
 
   /// Resets the input data block
   void ResetInputData() noexcept;
 
   /// Pushes back a hit
   /// \param  hit  An L1Hit object
-  void PushBackHit(const L1Hit& hit) { fInputData.fvHits.push_back(hit); }
+  void PushBackHit(const L1Hit& hit, int64_t streamId)
+  {
+    if (fInputData.fStreamStartIndices.size() == 0 || fLastStreamId != streamId) {  // new data stream
+      fLastStreamId = streamId;
+      fInputData.fStreamStartIndices.push_back(fInputData.fHits.size());
+      // for a case.. it is fixed later in InitData()
+      fInputData.fStreamStopIndices.push_back(fInputData.fHits.size());
+    }
+    fInputData.fHits.push_back(hit);
+  }
 
   /// Sets the number of hit keys
   /// \param  nKeys  Number of hit keys
@@ -101,10 +110,9 @@ private:
 
   L1InputData fInputData {};  ///< Object of input data
 
-  const L1Parameters* fpParameters = nullptr;  ///< Pointer to the tracking parameters object
+  const L1Parameters* fpParameters {nullptr};  ///< Pointer to the tracking parameters object
 
-  /// @brief Utility array for sorting hits by stations
-  std::array<L1HitIndex_t, L1Constants::size::kMaxNstations + 1> fvHitIndex = {0};
+  int64_t fLastStreamId {-1};  ///< data stream Id of the last hit added
 };
 
 
@@ -122,7 +130,7 @@ inline bool L1IODataManager::CheckInputData() const
   if constexpr (Level == 0) { return true; }  // Level = 0 -> do nothing
   else if constexpr (Level > 0) {             // Level = 1 and higher
     // ----- Check if the hits container is not empty ------------------------------------------------------------------
-    if (fInputData.fvHits.size() == 0) {
+    if (fInputData.fHits.size() == 0) {
       LOG(warn) << "L1IODataManager [check input]: Sample contains empty hits, tracking will not be executed";
       return false;
     }
diff --git a/reco/L1/L1Algo/L1InputData.cxx b/reco/L1/L1Algo/L1InputData.cxx
index e301b6a86e..d9087691f0 100644
--- a/reco/L1/L1Algo/L1InputData.cxx
+++ b/reco/L1/L1Algo/L1InputData.cxx
@@ -17,8 +17,9 @@ L1InputData::L1InputData() {}
 // ---------------------------------------------------------------------------------------------------------------------
 //
 L1InputData::L1InputData(const L1InputData& other)
-  : fvHits(other.fvHits)
-  , fvStartHitIndexes(other.fvStartHitIndexes)
+  : fHits(other.fHits)
+  , fStreamStartIndices(other.fStreamStartIndices)
+  , fStreamStopIndices(other.fStreamStopIndices)
   , fNhitKeys(other.fNhitKeys)
 {
 }
diff --git a/reco/L1/L1Algo/L1InputData.h b/reco/L1/L1Algo/L1InputData.h
index be6756666b..ae5c8a930b 100644
--- a/reco/L1/L1Algo/L1InputData.h
+++ b/reco/L1/L1Algo/L1InputData.h
@@ -13,9 +13,6 @@
 #include <boost/serialization/access.hpp>
 #include <boost/serialization/array.hpp>
 
-#include <array>
-
-#include "L1Constants.h"
 #include "L1Hit.h"
 #include "L1Vector.h"
 
@@ -56,35 +53,38 @@ public:
   L1InputData& operator=(L1InputData&& other) noexcept;
 
   /// Gets hits sample size
-  L1HitIndex_t GetSampleSize() const { return fvHits.size(); }
+  L1HitIndex_t GetSampleSize() const { return fHits.size(); }
 
 
   // ** Accessors **
 
+  /// Gets number of data streams
+  int GetNdataStreams() const { return fStreamStartIndices.size(); }
+
   /// Gets reference to hit by its index
   /// \param  index  Index of hit in the hits sample
-  const L1Hit& GetHit(L1HitIndex_t index) const { return fvHits[index]; }
+  const L1Hit& GetHit(L1HitIndex_t index) const { return fHits[index]; }
 
   /// Gets reference to hits vector
-  const L1Vector<L1Hit>& GetHits() const { return fvHits; }
+  const L1Vector<L1Hit>& GetHits() const { return fHits; }
 
   /// Gets number of hits in the hits vector
-  L1HitIndex_t GetNhits() const { return fvHits.size(); }
+  L1HitIndex_t GetNhits() const { return fHits.size(); }
 
   /// Gets total number of stored keys
   int GetNhitKeys() const { return fNhitKeys; }
 
   /// Gets index of the first hit in the sorted hits vector
-  /// \param iStation  Index of the tracking station in the active stations array
-  L1HitIndex_t GetStartHitIndex(int iStation) const { return fvStartHitIndexes[iStation]; }
+  /// \param iStream  Index of the data stream
+  L1HitIndex_t GetStreamStartIndex(int iStream) const { return fStreamStartIndices[iStream]; }
 
   /// Gets index of (the last + 1) hit in the sorted hits vector
-  /// \param iStation  Index of the tracking station in the active stations array
-  L1HitIndex_t GetStopHitIndex(int iStation) const { return fvStartHitIndexes[iStation + 1]; }
+  /// \param iStream  Index of the data stream
+  L1HitIndex_t GetStreamStopIndex(int iStream) const { return fStreamStopIndices[iStream]; }
 
-  /// Gets n hits for the station
-  /// \param iStation  Index of the tracking station in the active stations array
-  L1HitIndex_t GetNhits(int iStation) const { return fvStartHitIndexes[iStation + 1] - fvStartHitIndexes[iStation]; }
+  /// Gets n hits for the data stream
+  /// \param iStream  Index of the data stream
+  L1HitIndex_t GetStreamNhits(int iStream) const { return fStreamStopIndices[iStream] - fStreamStartIndices[iStream]; }
 
 
 private:
@@ -95,8 +95,9 @@ private:
   template<class Archive>
   void serialize(Archive& ar, const unsigned int /*versino*/)
   {
-    ar& fvHits;
-    ar& fvStartHitIndexes;
+    ar& fHits;
+    ar& fStreamStartIndices;
+    ar& fStreamStopIndices;
     ar& fNhitKeys;
   }
 
@@ -105,10 +106,11 @@ private:
   // ***************************
 
   /// @brief Sample of input hits
-  L1Vector<L1Hit> fvHits = {"L1InputData::fvHits"};
+  L1Vector<L1Hit> fHits {"L1InputData::fHits"};
 
-  /// @brief Index of the first hit in the sorted hits vector for a given station
-  std::array<L1HitIndex_t, L1Constants::size::kMaxNstations + 1> fvStartHitIndexes = {0};
+  /// @brief Index of the first hit in the sorted hits vector for a given data stream
+  L1Vector<L1HitIndex_t> fStreamStartIndices {"L1InputData::fStreamStartIndices"};
+  L1Vector<L1HitIndex_t> fStreamStopIndices {"L1InputData::fStreamStopIndices"};
 
   /// @brief Number of hit keys used for rejecting fake STS hits
   int fNhitKeys = -1;
@@ -123,8 +125,9 @@ private:
 //
 [[gnu::always_inline]] inline void L1InputData::Swap(L1InputData& other) noexcept
 {
-  std::swap(fvHits, other.fvHits);
-  std::swap(fvStartHitIndexes, other.fvStartHitIndexes);
+  std::swap(fHits, other.fHits);
+  std::swap(fStreamStartIndices, other.fStreamStartIndices);
+  std::swap(fStreamStopIndices, other.fStreamStopIndices);
   std::swap(fNhitKeys, other.fNhitKeys);
 }
 
diff --git a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
index c11bfb6a5e..67e72ae7c9 100644
--- a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
+++ b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
@@ -74,14 +74,15 @@ void L1AlgoDraw::InitL1Draw(L1Algo* algo_)
   //   algo = CbmL1::Instance()->algo;
   algo = algo_;
 
-  vHits.reserve(algo->GetInputData().GetNhits());
+  vHits.clear();
+  vHits.reserve(algo->fGridHits.size());
   for (unsigned int i = 0; i < algo->GetInputData().GetNhits(); i++) {
-    vHits.push_back(algo->GetInputData().GetHit(i));
+    vHits.push_back(algo->fGridHits[i]);
   }
   NStations = algo->GetParameters()->GetNstationsActive();
   for (int i = 0; i < NStations; i++) {
-    HitsStartIndex[i] = algo->GetInputData().GetStartHitIndex(i);
-    HitsStopIndex[i]  = algo->GetInputData().GetStopHitIndex(i);
+    HitsStartIndex[i] = algo->fGridHitStartIndex[i];
+    HitsStopIndex[i]  = algo->fGridHitStopIndex[i];
     vStations[i]      = algo->GetParameters()->GetStation(i);
   }
 }
-- 
GitLab