diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index f39bba9a8b75ed8849c55ffe7a14302cad8d3a9f..ec33fcef9da3b10adf655121a030fd5cf4d09f8f 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -1254,7 +1254,7 @@ InitStatus CbmL1::Init()
   algo->Init(geo, fUseHitErrors, fTrackingMode, fMissingHits);
   geo.clear();
 
-  algo->fRadThick.reset(algo->NStations);
+  algo->fRadThick.reset(algo->GetNstations());
 
   // Read STS  MVD TRD MuCh ToF Radiation Thickness table
   TString stationName = "Radiation Thickness [%], Station";
@@ -1283,13 +1283,13 @@ InitStatus CbmL1::Init()
           for (int iB2 = 0; iB2 < NBins; iB2++) {
             algo->fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB, iB2);
             // Correction for holes in material map
-            if (algo->fRadThick[iSta].table[iB][iB2] < algo->vStations[iSta].materialInfo.RadThick[0])
+            if (algo->fRadThick[iSta].table[iB][iB2] < algo->GetStations()[iSta].materialInfo.RadThick[0])
               if (iB2 > 0 && iB2 < NBins - 1)
                 algo->fRadThick[iSta].table[iB][iB2] = TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1),
                                                                   0.01 * hStaRadLen->GetBinContent(iB, iB2 + 1));
             // Correction for the incorrect harcoded value of RadThick of MVD stations
             if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015) algo->fRadThick[iSta].table[iB][iB2] = 0.0015;
-            //              algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0];
+            //              algo->fRadThick[iSta].table[iB][iB2] = algo->GetStations()[iSta].materialInfo.RadThick[0];
             sumRF += algo->fRadThick[iSta].table[iB][iB2];
             nRF++;
           }
@@ -1310,7 +1310,7 @@ InitStatus CbmL1::Init()
                                       100);  // mvd should be in +-100 cm square
         algo->fRadThick[iSta].table.resize(1);
         algo->fRadThick[iSta].table[0].resize(1);
-        algo->fRadThick[iSta].table[0][0] = algo->vStations[iSta].materialInfo.RadThick[0];
+        algo->fRadThick[iSta].table[0][0] = algo->GetStations()[iSta].materialInfo.RadThick[0];
       }
     }
   }
@@ -1337,8 +1337,8 @@ InitStatus CbmL1::Init()
         algo->fRadThick[iSta].table[iB].resize(NBins);
         for (int iB2 = 0; iB2 < NBins; iB2++) {
           algo->fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB, iB2);
-          if (algo->fRadThick[iSta].table[iB][iB2] < algo->vStations[iSta].materialInfo.RadThick[0])
-            algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0];
+          if (algo->fRadThick[iSta].table[iB][iB2] < algo->GetStations()[iSta].materialInfo.RadThick[0])
+            algo->fRadThick[iSta].table[iB][iB2] = algo->GetStations()[iSta].materialInfo.RadThick[0];
           // cout << " iSta " << iSta << " iB " << iB << "iB2 " << iB2 << " RadThick " << algo->fRadThick[iSta].table[iB][iB2] << endl;
           sumRF += algo->fRadThick[iSta].table[iB][iB2];
           nRF++;
@@ -1359,7 +1359,7 @@ InitStatus CbmL1::Init()
       algo->fRadThick[iSta].SetBins(1, 100);
       algo->fRadThick[iSta].table.resize(1);
       algo->fRadThick[iSta].table[0].resize(1);
-      algo->fRadThick[iSta].table[0][0] = algo->vStations[iSta].materialInfo.RadThick[0];
+      algo->fRadThick[iSta].table[0][0] = algo->GetStations()[iSta].materialInfo.RadThick[0];
     }
   }
 
@@ -1393,7 +1393,7 @@ InitStatus CbmL1::Init()
           for (int iB2 = 0; iB2 < NBins; iB2++) {
             algo->fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB, iB2);
             // Correction for holes in material map
-            //if(algo->fRadThick[iSta].table[iB][iB2] < algo->vStations[iSta].materialInfo.RadThick[0])
+            //if(algo->fRadThick[iSta].table[iB][iB2] < algo->GetStations()[iSta].materialInfo.RadThick[0])
 
             if (iB2 > 0 && iB2 < NBins - 1)
               algo->fRadThick[iSta].table[iB][iB2] = TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1),
@@ -1402,7 +1402,7 @@ InitStatus CbmL1::Init()
 
             if (algo->fRadThick[iSta].table[iB][iB2] > 0.0015) hole = algo->fRadThick[iSta].table[iB][iB2];
             if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015) algo->fRadThick[iSta].table[iB][iB2] = hole;
-            //              algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0];
+            //              algo->fRadThick[iSta].table[iB][iB2] = algo->GetStations()[iSta].materialInfo.RadThick[0];
           }
         }
       }
@@ -1419,7 +1419,7 @@ InitStatus CbmL1::Init()
         algo->fRadThick[iSta].SetBins(1, 100);
         algo->fRadThick[iSta].table.resize(1);
         algo->fRadThick[iSta].table[0].resize(1);
-        algo->fRadThick[iSta].table[0][0] = algo->vStations[iSta].materialInfo.RadThick[0];
+        algo->fRadThick[iSta].table[0][0] = algo->GetStations()[iSta].materialInfo.RadThick[0];
       }
     }
 
@@ -1457,7 +1457,7 @@ InitStatus CbmL1::Init()
           for (int iB2 = 0; iB2 < NBins; iB2++) {
             algo->fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB, iB2);
             // Correction for holes in material map
-            //if(algo->fRadThick[iSta].table[iB][iB2] < algo->vStations[iSta].materialInfo.RadThick[0])
+            //if(algo->fRadThick[iSta].table[iB][iB2] < algo->GetStations()[iSta].materialInfo.RadThick[0])
 
             if (iB2 > 0 && iB2 < NBins - 1)
               algo->fRadThick[iSta].table[iB][iB2] = TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1),
@@ -1466,7 +1466,7 @@ InitStatus CbmL1::Init()
 
             if (algo->fRadThick[iSta].table[iB][iB2] > 0.0015) hole = algo->fRadThick[iSta].table[iB][iB2];
             if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015) algo->fRadThick[iSta].table[iB][iB2] = hole;
-            //              algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0];
+            //              algo->fRadThick[iSta].table[iB][iB2] = algo->GetStations()[iSta].materialInfo.RadThick[0];
           }
         }
       }
@@ -1484,8 +1484,8 @@ InitStatus CbmL1::Init()
         algo->fRadThick[iSta].SetBins(1, 100);
         algo->fRadThick[iSta].table.resize(1);
         algo->fRadThick[iSta].table[0].resize(1);
-        algo->fRadThick[iSta].table[0][0] = algo->vStations[iSta].materialInfo.RadThick[0];
-        cout << "TRD material: " << algo->vStations[iSta].materialInfo.RadThick[0] << endl;
+        algo->fRadThick[iSta].table[0][0] = algo->GetStations()[iSta].materialInfo.RadThick[0];
+        cout << "TRD material: " << algo->GetStations()[iSta].materialInfo.RadThick[0] << endl;
       }
     }
 
@@ -1518,7 +1518,7 @@ InitStatus CbmL1::Init()
           for (int iB2 = 0; iB2 < NBins; iB2++) {
             algo->fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB, iB2);
             // Correction for holes in material map
-            //if(algo->fRadThick[iSta].table[iB][iB2] < algo->vStations[iSta].materialInfo.RadThick[0])
+            //if(algo->fRadThick[iSta].table[iB][iB2] < algo->GetStations()[iSta].materialInfo.RadThick[0])
 
             if (iB2 > 0 && iB2 < NBins - 1)
               algo->fRadThick[iSta].table[iB][iB2] = TMath::Min(0.01 * hStaRadLen->GetBinContent(iB, iB2 - 1),
@@ -1527,7 +1527,7 @@ InitStatus CbmL1::Init()
 
             if (algo->fRadThick[iSta].table[iB][iB2] > 0.0015) hole = algo->fRadThick[iSta].table[iB][iB2];
             if (algo->fRadThick[iSta].table[iB][iB2] < 0.0015) algo->fRadThick[iSta].table[iB][iB2] = hole;
-            //              algo->fRadThick[iSta].table[iB][iB2] = algo->vStations[iSta].materialInfo.RadThick[0];
+            //              algo->fRadThick[iSta].table[iB][iB2] = algo->GetStations()[iSta].materialInfo.RadThick[0];
           }
         }
       }
@@ -1545,7 +1545,7 @@ InitStatus CbmL1::Init()
         algo->fRadThick[iSta].SetBins(1, 100);
         algo->fRadThick[iSta].table.resize(1);
         algo->fRadThick[iSta].table[0].resize(1);
-        algo->fRadThick[iSta].table[0][0] = algo->vStations[iSta].materialInfo.RadThick[0];
+        algo->fRadThick[iSta].table[0][0] = algo->GetStations()[iSta].materialInfo.RadThick[0];
       }
     }
   return kSUCCESS;
@@ -1670,7 +1670,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
         h.n = mcp.event;
 #endif
         const int ista       = (*algo->fStripFlag)[h.f] / 4;
-        const L1Station& sta = algo->vStations[ista];
+        const L1Station& sta = algo->GetStations()[ista];
         if (std::find(strips.begin(), strips.end(), h.f) != strips.end()) {  // separate strips
 
           (*algo->fStripFlag).push_back((*algo->fStripFlag)[h.f]);
@@ -2008,7 +2008,7 @@ void CbmL1::IdealTrackFinder()
     L1Track algoTr;
     algoTr.NHits = 0;
 
-    L1Vector<int> hitIndices("CbmL1::hitIndices", algo->NStations, -1);
+    L1Vector<int> hitIndices("CbmL1::hitIndices", algo->GetNstations(), -1);
 
     for (unsigned int iH = 0; iH < MC.StsHits.size(); iH++) {
       const int hitI      = MC.StsHits[iH];
@@ -2020,7 +2020,7 @@ void CbmL1::IdealTrackFinder()
     }
 
 
-    for (int iH = 0; iH < algo->NStations; iH++) {
+    for (int iH = 0; iH < algo->GetNstations(); iH++) {
       const int hitI = hitIndices[iH];
       if (hitI < 0) continue;
 
@@ -2031,7 +2031,7 @@ void CbmL1::IdealTrackFinder()
 
     if (algoTr.NHits < 3) continue;
 
-    for (int iH = 0; iH < algo->NStations; iH++) {
+    for (int iH = 0; iH < algo->GetNstations(); iH++) {
       const int hitI = hitIndices[iH];
       if (hitI < 0) continue;
       algo->fRecoHits.push_back(hitI);
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index d5046dc68a04f15cb02e4967c3a2cc0ea6bcc976..6c602bfa69958c1d239156892df6b89f31c5a55c 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -318,9 +318,9 @@ void CbmL1::EfficienciesPerformance()
     //     }
   }
 
-  int sta_nhits[algo->NStations], sta_nfakes[algo->NStations];
+  int sta_nhits[algo->GetNstations()], sta_nfakes[algo->GetNstations()];
 
-  for (int i = 0; i < algo->NStations; i++) {
+  for (int i = 0; i < algo->GetNstations(); i++) {
     sta_nhits[i]  = 0;
     sta_nfakes[i] = 0;
   }
@@ -457,7 +457,7 @@ void CbmL1::EfficienciesPerformance()
     if (fVerbose > 1) {
       ntra.PrintEff();
       cout << "Number of true and fake hits in stations: " << endl;
-      for (int i = 0; i < algo->NStations; i++) {
+      for (int i = 0; i < algo->GetNstations(); i++) {
         cout << sta_nhits[i] - sta_nfakes[i] << "+" << sta_nfakes[i] << "   ";
       }
       cout << endl;
@@ -798,8 +798,8 @@ void CbmL1::HistoPerformance()  // TODO: check if works correctly. Change vHitRe
       CbmL1HitStore& h2 = vHitStore[prtra->StsHits[1]];
       h_ghost_fstation->Fill(h1.iStation);
       h_ghost_r->Fill(sqrt(fabs(h1.x * h1.x + h1.y * h1.y)));
-      double z1 = algo->vStations[h1.iStation].z[0];
-      double z2 = algo->vStations[h2.iStation].z[0];
+      double z1 = algo->GetStations()[h1.iStation].z[0];
+      double z2 = algo->GetStations()[h2.iStation].z[0];
       if (fabs(z2 - z1) > 1.e-4) {
         h_ghost_tx->Fill((h2.x - h1.x) / (z2 - z1));
         h_ghost_ty->Fill((h2.y - h1.y) / (z2 - z1));
@@ -990,8 +990,8 @@ void CbmL1::HistoPerformance()  // TODO: check if works correctly. Change vHitRe
       //      CbmL1HitStore &ph21 = vHitStore[mtra.StsHits[0]];
       //      CbmL1HitStore &ph22 = vHitStore[mtra.StsHits[1]];
 
-      //      double z21 = algo->vStations[ph21.iStation].z[0];
-      //      double z22 = algo->vStations[ph22.iStation].z[0];
+      //      double z21 = algo->GetStations()[ph21.iStation].z[0];
+      //      double z22 = algo->GetStations()[ph22.iStation].z[0];
       //      if( fabs(z22-z21)>1.e-4 ){
       //        h_notfound_tx->Fill((ph22.x-ph21.x)/(z22-z21));
       //        h_notfound_ty->Fill((ph22.y-ph21.y)/(z22-z21));
@@ -1191,7 +1191,7 @@ void CbmL1::TrackFitPerformance()
       for (unsigned int iMCPoint = 0; iMCPoint < mc.Points.size(); iMCPoint++) {
         const int iMCP    = mc.Points[iMCPoint];
         CbmL1MCPoint& mcP = vMCPoints[iMCP];
-        L1Station& st     = algo->vStations[mcP.iStation];
+        const L1Station& st     = algo->GetStations()[mcP.iStation];
         z[ih]             = st.z[0];
         if (ih > 0 && (z[ih] - z[ih - 1]) < 0.1) continue;
         st.fieldSlice.GetFieldValue(mcP.x, mcP.y, B[ih]);
@@ -1302,7 +1302,7 @@ void CbmL1::TrackFitPerformance()
       for (unsigned int iMCPoint = 0; iMCPoint < mc.Points.size(); iMCPoint++) {
         const int iMCP    = mc.Points[iMCPoint];
         CbmL1MCPoint& mcP = vMCPoints[iMCP];
-        L1Station& st     = algo->vStations[mcP.iStation];
+        const L1Station& st     = algo->GetStations()[mcP.iStation];
         z[ih]             = st.z[0];
         if (ih > 0 && (z[ih] - z[ih - 1]) < 0.1) continue;
         st.fieldSlice.GetFieldValue(mcP.x, mcP.y, B[ih]);
@@ -1373,7 +1373,7 @@ void CbmL1::TrackFitPerformance()
             if (ih >= mc.Points.size()) continue;  //If nofMCPoints in track < 3
             const int iMCP    = mc.Points[ih];
             CbmL1MCPoint& mcP = vMCPoints[iMCP];
-            L1Station& st     = algo->vStations[mcP.iStation];
+            const L1Station& st     = algo->GetStations()[mcP.iStation];
             z[ih]             = st.z[0];
             st.fieldSlice.GetFieldValue(mcP.x, mcP.y, B[ih]);
           };
@@ -1382,12 +1382,12 @@ void CbmL1::TrackFitPerformance()
           L1Extrapolate(trPar, mc.z, trPar.qp, fld);
           // add material
           const int fSta = vHitStore[it->StsHits[0]].iStation;
-          const int dir  = int((mc.z - algo->vStations[fSta].z[0]) / fabs(mc.z - algo->vStations[fSta].z[0]));
-          //         if (abs(mc.z - algo->vStations[fSta].z[0]) > 10.) continue; // can't extrapolate on large distance
+          const int dir  = int((mc.z - algo->GetStations()[fSta].z[0]) / fabs(mc.z - algo->GetStations()[fSta].z[0]));
+          //         if (abs(mc.z - algo->GetStations()[fSta].z[0]) > 10.) continue; // can't extrapolate on large distance
           for (int iSta = fSta /*+dir*/;
-               (iSta >= 0) && (iSta < NStation) && (dir * (mc.z - algo->vStations[iSta].z[0]) > 0); iSta += dir) {
+               (iSta >= 0) && (iSta < NStation) && (dir * (mc.z - algo->GetStations()[iSta].z[0]) > 0); iSta += dir) {
             //           cout << iSta << " " << dir << endl;
-            fit.L1AddMaterial(trPar, algo->vStations[iSta].materialInfo, trPar.qp, 1);
+            fit.L1AddMaterial(trPar, algo->GetStations()[iSta].materialInfo, trPar.qp, 1);
             if (iSta + dir == NMvdStations - 1) fit.L1AddPipeMaterial(trPar, trPar.qp, 1);
           }
         }
@@ -1436,7 +1436,7 @@ void CbmL1::TrackFitPerformance()
           int ih = 1;
           for (unsigned int iHit = 0; iHit < it->StsHits.size(); iHit++) {
             const int iStation = vHitStore[it->StsHits[iHit]].iStation;
-            L1Station& st      = algo->vStations[iStation];
+            const L1Station& st      = algo->GetStations()[iStation];
             z[ih]              = st.z[0];
             st.fieldSlice.GetFieldValue(vHitStore[it->StsHits[iHit]].x, vHitStore[it->StsHits[iHit]].y, B[ih]);
             ih++;
@@ -1447,18 +1447,18 @@ void CbmL1::TrackFitPerformance()
           // add material
           const int fSta = vHitStore[it->StsHits[0]].iStation;
 
-          const int dir = (mc.z - algo->vStations[fSta].z[0]) / abs(mc.z - algo->vStations[fSta].z[0]);
-          //         if (abs(mc.z - algo->vStations[fSta].z[0]) > 10.) continue; // can't extrapolate on large distance
+          const int dir = (mc.z - algo->GetStations()[fSta].z[0]) / abs(mc.z - algo->GetStations()[fSta].z[0]);
+          //         if (abs(mc.z - algo->GetStations()[fSta].z[0]) > 10.) continue; // can't extrapolate on large distance
 
           for (int iSta = fSta + dir;
-               (iSta >= 0) && (iSta < NStation) && (dir * (mc.z - algo->vStations[iSta].z[0]) > 0); iSta += dir) {
+               (iSta >= 0) && (iSta < NStation) && (dir * (mc.z - algo->GetStations()[iSta].z[0]) > 0); iSta += dir) {
 
-            z[0]     = algo->vStations[iSta].z[0];
+            z[0]     = algo->GetStations()[iSta].z[0];
             float dz = z[1] - z[0];
-            algo->vStations[iSta].fieldSlice.GetFieldValue(trPar.x - trPar.tx * dz, trPar.y - trPar.ty * dz, B[0]);
+            algo->GetStations()[iSta].fieldSlice.GetFieldValue(trPar.x - trPar.tx * dz, trPar.y - trPar.ty * dz, B[0]);
             fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
 
-            L1Extrapolate(trPar, algo->vStations[iSta].z[0], trPar.qp, fld);
+            L1Extrapolate(trPar, algo->GetStations()[iSta].z[0], trPar.qp, fld);
             fit.L1AddMaterial(trPar, algo->fRadThick[iSta].GetRadThick(trPar.x, trPar.y), trPar.qp, 1);
             fit.EnergyLossCorrection(trPar, algo->fRadThick[iSta].GetRadThick(trPar.x, trPar.y), trPar.qp, fvec(1.f),
                                      fvec(1.f));
@@ -1613,7 +1613,7 @@ void CbmL1::FieldApproxCheck()
 
     const int M   = 5;  // polinom order
     const int N   = (M + 1) * (M + 2) / 2;
-    L1Station& st = algo->vStations[ist];
+    const L1Station& st = algo->GetStations()[ist];
     for (int i = 0; i < N; i++) {
       FSl.cx[i] = st.fieldSlice.cx[i][0];
       FSl.cy[i] = st.fieldSlice.cy[i][0];
@@ -1982,8 +1982,8 @@ void CbmL1::InputPerformance()
           pullTsts->Fill((sh->GetTime() - mcTime) / sh->GetTimeError());
         }
         else {  // errors used in TF
-          pullXsts->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
-          pullYsts->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
+          pullXsts->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->GetStations()[NMvdStations].XYInfo.C00[0]));
+          pullYsts->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->GetStations()[NMvdStations].XYInfo.C11[0]));
         }
 
         resXsts->Fill((hitPos.X() - mcPos.X()) * 10 * 1000);
@@ -2028,8 +2028,8 @@ void CbmL1::InputPerformance()
       //       if (hitErr.X() != 0) pullX->Fill( (hitPos.X() - mcPos.X()) / sh->GetDx() ); // qa errors
       //       if (hitErr.Y() != 0) pullY->Fill( (hitPos.Y() - mcPos.Y()) / sh->GetDy() );
       if (hitErr.X() != 0)
-        pullXmvd->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->vStations[0].XYInfo.C00[0]));  // errors used in TF
-      if (hitErr.Y() != 0) pullYmvd->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[0].XYInfo.C11[0]));
+        pullXmvd->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->GetStations()[0].XYInfo.C00[0]));  // errors used in TF
+      if (hitErr.Y() != 0) pullYmvd->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->GetStations()[0].XYInfo.C11[0]));
 
       resXmvd->Fill((hitPos.X() - mcPos.X()) * 10 * 1000);
       resYmvd->Fill((hitPos.Y() - mcPos.Y()) * 10 * 1000);
@@ -2095,8 +2095,8 @@ void CbmL1::InputPerformance()
         pullTmuch->Fill((h.t - mcTime) / sh->GetTimeError());
       }
       else {  // errors used in TF
-        pullXmuch->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
-        pullYmuch->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
+        pullXmuch->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->GetStations()[NMvdStations].XYInfo.C00[0]));
+        pullYmuch->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->GetStations()[NMvdStations].XYInfo.C11[0]));
       }
 
       resXmuch->Fill((h.x - mcPos.X()) * 10 * 1000);
@@ -2162,8 +2162,8 @@ void CbmL1::InputPerformance()
         pullTtrd->Fill((h.t - mcTime) / sh->GetTimeError());
       }
       else {  // errors used in TF
-        pullXtrd->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
-        pullYtrd->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
+        pullXtrd->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->GetStations()[NMvdStations].XYInfo.C00[0]));
+        pullYtrd->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->GetStations()[NMvdStations].XYInfo.C11[0]));
       }
 
       resXtrd->Fill((h.x - mcPos.X()) * 10 * 1000);
@@ -2230,8 +2230,8 @@ void CbmL1::InputPerformance()
         pullTtof->Fill((sh->GetTime() - mcTime) / sh->GetTimeError());
       }
       else {  // errors used in TF
-        pullXtof->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->vStations[NMvdStations].XYInfo.C00[0]));
-        pullYtof->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->vStations[NMvdStations].XYInfo.C11[0]));
+        pullXtof->Fill((hitPos.X() - mcPos.X()) / sqrt(algo->GetStations()[NMvdStations].XYInfo.C00[0]));
+        pullYtof->Fill((hitPos.Y() - mcPos.Y()) / sqrt(algo->GetStations()[NMvdStations].XYInfo.C11[0]));
       }
 
       resXtof->Fill((h.x - mcPos.X()) * 10 * 1000);
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 176daae07aef3172afcd125ad66ec54a66b27d85..34a9a262d91322d9caaf712914c796811394f6d7 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -110,7 +110,7 @@ struct TmpHit {
   /// \param st       reference to the station info object
   // TODO: Probably, L1Station& st parameter should be constant. Do we really want to modify a station here? (S.Zharko)
   void CreateHitFromPoint(const CbmL1MCPoint& point, int det, int nTmpHits, int nStripF, int ip, int& NStrips,
-                          L1Station& st)
+                          const L1Station& st)
   {
     ExtIndex = 0;
     Det      = det;
@@ -156,7 +156,7 @@ struct TmpHit {
   /// \param point  constant reference to the input MC-point
   /// \param st     reference to the station info object
   // TODO: Probably, L1Station& st parameter should be constant. Do we really want to modify a station here? (S.Zharko)
-  void SetHitFromPoint(const CbmL1MCPoint& point, L1Station& st)
+  void SetHitFromPoint(const CbmL1MCPoint& point, const L1Station& st)
   {
     x       = 0.5 * (point.xIn + point.xOut) + gRandom->Gaus(0, dx);
     y       = 0.5 * (point.yIn + point.yOut) + gRandom->Gaus(0, dy);
@@ -266,7 +266,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
           CbmL1MCPoint MC;
           if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 1)) {
             MC.iStation     = -1;
-            L1Station* sta  = algo->vStations; // TODO: Wrap it into interface algo->GetStations() (S.Zharko)
+            const L1Station* sta  = algo->GetStations().begin(); 
             double bestDist = 1.e20;
             for (Int_t iSt = 0; iSt < NMvdStations; iSt++) {
               // use z_in since z_out is sometimes very wrong
@@ -303,7 +303,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
           CbmL1MCPoint MC;
           if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 0)) {
             MC.iStation     = -1;
-            L1Station* sta  = algo->vStations + NMvdStations;
+            const L1Station* sta  = algo->GetStations().begin() + NMvdStations;
             double bestDist = 1.e20;
             for (Int_t iSt = 0; iSt < NStsStations; iSt++) {
               // use z_in since z_out is sometimes very wrong
@@ -339,7 +339,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
           CbmL1MCPoint MC;
           if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 2)) {
             MC.iStation    = -1;
-            L1Station* sta = algo->vStations + NMvdStations + NStsStations;
+            const L1Station* sta = algo->GetStations().begin() + NMvdStations + NStsStations;
             for (Int_t iSt = 0; iSt < NMuchStations; iSt++) {
               if (MC.z > sta[iSt].z[0] - 2.5) { MC.iStation = NMvdStations + NStsStations + iSt; }
             }
@@ -364,7 +364,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
           CbmL1MCPoint MC;
           if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 3)) {
             MC.iStation    = -1;
-            L1Station* sta = algo->vStations + NMvdStations + NStsStations + NMuchStations;
+            const L1Station* sta = algo->GetStations().begin() + NMvdStations + NStsStations + NMuchStations;
             for (Int_t iSt = 0; iSt < NTrdStations; iSt++) {
               if (MC.z > sta[iSt].z[0] - 4.0) { MC.iStation = NMvdStations + NStsStations + NMuchStations + iSt; }
             }
@@ -416,7 +416,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
             Int_t IND_Track = trk_it->second;
 
             MC.iStation    = -1;
-            L1Station* sta = algo->vStations + NMvdStations + NStsStations + NMuchStations + NTrdStations;
+            const L1Station* sta = algo->GetStations().begin() + NMvdStations + NStsStations + NMuchStations + NTrdStations;
             for (Int_t iSt = 0; iSt < NTOFStation; iSt++)
               MC.iStation = (MC.z > sta[iSt].z[0] - 15)
                               ? (NMvdStations + NStsStations + NMuchStations + NTrdStations + iSt)
@@ -443,7 +443,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
             if (!ReadMCPoint(&MC, TofPointToTrack[iTofSta][iMC], iFile, iEvent, 4)) {
 
               MC.iStation    = -1;
-              L1Station* sta = algo->vStations + NMvdStations + NStsStations + NMuchStations + NTrdStations;
+              const L1Station* sta = algo->GetStations().begin() + NMvdStations + NStsStations + NMuchStations + NTrdStations;
               for (Int_t iSt = 0; iSt < NTOFStation; iSt++)
                 MC.iStation = (MC.z > sta[iSt].z[0] - 15)
                                 ? (NMvdStations + NStsStations + NMuchStations + NTrdStations + iSt)
@@ -534,7 +534,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
         th.y = pos.Y();
         th.z = pos.Z();
 
-        L1Station& st = algo->vStations[th.iStation];
+        const L1Station& st = algo->GetStations()[th.iStation];
         th.u_front    = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
         th.u_back     = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
       }
@@ -575,7 +575,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
       //            << p.time << " mc " << p.ID << " p " << p.p << endl;
       TmpHit th;
       int DetId = 1;
-      th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->vStations[p.iStation]);
+      th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->GetStations()[p.iStation]);
       tmpHits.push_back(th);
       nStsHits++;
     }
@@ -652,7 +652,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
         th.du = mh->GetDu();
         th.dv = mh->GetDv();
 
-        L1Station& st = algo->vStations[th.iStation];
+        const L1Station& st = algo->GetStations()[th.iStation];
         th.u_front    = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
         th.u_back     = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
       }
@@ -736,7 +736,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
 
       TmpHit th;
       int DetId = 2;
-      th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->vStations[p.iStation]);
+      th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->GetStations()[p.iStation]);
 
       tmpHits.push_back(th);
       nMuchHits++;
@@ -791,7 +791,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
         th.dv = mh->GetDy();
 
 
-        L1Station& st = algo->vStations[th.iStation];
+        const L1Station& st = algo->GetStations()[th.iStation];
         th.u_front    = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
         th.u_back     = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
       }
@@ -817,7 +817,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
               if (trk_it == dFEI2vMCPoints.end()) continue;
               th.iMC = trk_it->second;
               if ((1 == fMuchUseMcHit) && (th.iMC > -1))
-                th.SetHitFromPoint(vMCPoints[th.iMC], algo->vStations[th.iStation]);
+                th.SetHitFromPoint(vMCPoints[th.iMC], algo->GetStations()[th.iStation]);
             }
           }
         }
@@ -841,7 +841,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
 
       TmpHit th;
       int DetId = 3;
-      th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->vStations[p.iStation]);
+      th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->GetStations()[p.iStation]);
       tmpHits.push_back(th);
       nTrdHits++;
     }
@@ -902,7 +902,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
       th.du = fabs(mh->GetDx());
       th.dv = fabs(mh->GetDy());
 
-      L1Station& st = algo->vStations[th.iStation];
+      const L1Station& st = algo->GetStations()[th.iStation];
       th.u_front    = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
       th.u_back     = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
 
@@ -921,7 +921,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
 
             if ((1 == fTrdUseMcHit) && (th.iMC > -1))  //SG!!! replace hits by MC points
 
-              th.SetHitFromPoint(vMCPoints[th.iMC], algo->vStations[th.iStation]);
+              th.SetHitFromPoint(vMCPoints[th.iMC], algo->GetStations()[th.iStation]);
 
             if (L1Algo::TrackingMode::kGlobal == fTrackingMode) {  //SG!!! replace hits by MC points
 
@@ -987,7 +987,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
       TmpHit th;
 
       int DetId = 4;
-      th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->vStations[p.iStation]);
+      th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips, algo->GetStations()[p.iStation]);
       tmpHits.push_back(th);
 
       nTofHits++;
@@ -1058,7 +1058,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
       if (stIdx == -1) continue;
       th.iStation = stIdx;
 
-      L1Station& st = algo->vStations[th.iStation];
+      const L1Station& st = algo->GetStations()[th.iStation];
       th.u_front    = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
       th.u_back     = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
 
@@ -1082,7 +1082,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
           Double_t dtrck          = dFEI(iFile, iEvent, pt->GetTrackID());
           DFEI2I::iterator trk_it = dFEI2vMCPoints.find(dtrck);
           if (trk_it != dFEI2vMCPoints.end()) th.iMC = TofPointToTrack[sttof][trk_it->second];
-          if ((1 == fTofUseMcHit) && (th.iMC > -1)) th.SetHitFromPoint(vMCPoints[th.iMC], algo->vStations[th.iStation]);
+          if ((1 == fTofUseMcHit) && (th.iMC > -1)) th.SetHitFromPoint(vMCPoints[th.iMC], algo->GetStations()[th.iStation]);
         }
       }
 
diff --git a/reco/L1/CbmL1TrackFitter.cxx b/reco/L1/CbmL1TrackFitter.cxx
index acc6bd4c663c3851180e8a0ea23a909002971916..375150ebfb422b41341947f2bb763c29fb772d4d 100644
--- a/reco/L1/CbmL1TrackFitter.cxx
+++ b/reco/L1/CbmL1TrackFitter.cxx
@@ -35,7 +35,7 @@ void CbmL1::TrackFitter(vector<CbmL1Track>& Tracks, CbmL1Vtx* V)
 {
   TStopwatch timer;
   timer.Start();
-  /*
+  /* NOTE: algo->vStations is deprecated; one have to change it with algo->GetStations() (see L1Algo.h for documentation) (S.Zharko)
   //int nt=0;
   for( vector<CbmL1Track>::iterator i = Tracks.begin(); i!=Tracks.end(); ++i)
     {
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 2d9196640b488619db661f6907b51d0c1ce130e1..f0251aaeee8a8bb79c70139aff3e1585471d4cff 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -100,7 +100,8 @@ void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const Tra
     vtxFieldValue = B[0];
   }
   //vStations.clear();
-  NStations    = static_cast<int>(geo[ind++]);
+  L1Station vStations[L1Parameters::kMaxNstations] _fvecalignment;
+  int NStations    = static_cast<int>(geo[ind++]);
   NMvdStations = static_cast<int>(geo[ind++]);  // TODO: get rid of NMbdStations (S. Zh.)
   NStsStations = static_cast<int>(geo[ind++]);  // TODO: get rid of NStsStations (S. Zh.)
 
@@ -241,12 +242,14 @@ void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const Tra
   int nStsStationsNew   = fInitManager.GetStationsNumber(static_cast<L1DetectorID>(1));
   int nFieldStationsNew = nMvdStationsNew + nStsStationsNew;
 
+  fNstations = fInitManager.GetStationsNumber();
+
   // Get field near target
   L1FieldValue vtxFieldValueNew   = fInitManager.GetTargetFieldValue();
   L1FieldRegion vtxFieldRegionNew = fInitManager.GetTargetFieldRegion();
 
   // Fill L1Station array
-  fInitManager.TransferL1StationArray(fStationsNew);
+  fInitManager.TransferL1StationArray(fStations);
 
   LOG(info) << "**********************************************************************";
   LOG(info) << "*  New L1Algo initialization cross check  (tmp log, to be removed!)  *";
@@ -257,7 +260,7 @@ void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const Tra
   LOG(info) << "\tSTS:   " << NStsStations;
   LOG(info) << "\tField: " << fNfieldStations;
   LOG(info) << "** Number of stations (new) **";
-  LOG(info) << "\tTotal: " << nStationsNew;
+  LOG(info) << "\tTotal: " << fNstations;
   LOG(info) << "\tMVD:   " << nMvdStationsNew;
   LOG(info) << "\tSTS:   " << nStsStationsNew;
   LOG(info) << "\tField: " << nFieldStationsNew;
@@ -272,15 +275,14 @@ void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const Tra
 
   LOG(info) << "** Original L1Station array content **";
   int nStations = fInitManager.GetStationsNumber();
-  for (int iSt = 0; iSt < nStations; ++iSt) {
+  for (int iSt = 0; iSt < NStations; ++iSt) {
     LOG(info) << "Station Global No: " << iSt << ' ' << vStations[iSt].ToString(/*verbosity = */ 0);
   }
-  //LOG(info) << "** New L1Station array content **";
-  //nStations = fInitManager.GetStationsNumber();
-  //for (int iSt = 0; iSt < nStations; ++iSt) {
-  //  LOG(info) << "Station Global No: " << iSt;
-  //  LOG(info) << '\n' << fStationsNew[iSt].ToString(/*verbosity = */ 0);
-  //}
+  LOG(info) << "** New L1Station array content **";
+  for (int iSt = 0; iSt < fNstations; ++iSt) {
+    LOG(info) << "Station Global No: " << iSt;
+    LOG(info) << '\n' << fStations[iSt].ToString(/*verbosity = */ 0);
+  }
 
   // Print L1Parameters
   fParameters.Print(/*verbosity=*/0);
@@ -352,7 +354,7 @@ void L1Algo::SetData(L1Vector<L1Hit>& StsHits_, int nStsStrips_, L1Vector<unsign
 
 void L1Algo::GetHitCoor(const L1Hit& _h, fscal& _x, fscal& _y, char iS)
 {
-  L1Station& sta = vStations[int(iS)];
+  L1Station& sta = fStations[int(iS)];
   fscal u        = _h.u;
   fscal v        = _h.v;
   _x             = (sta.xInfo.sin_phi[0] * u + sta.xInfo.cos_phi[0] * v) / (_h.z - fCbmTargetZ[0]);
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index d5037da6145f97b86e2e1403f39e39e4c43e4e38..5560a8695aa0cd6f18f4289bfdd2dbce7efebbbc 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -219,15 +219,27 @@ public:
 
   void SetNThreads(unsigned int n);
 
-  int NStations {0};        ///< number of all detector stations
+private:
+  int fNstations {0};        ///< number of all detector stations
+public:
+  
+  int GetNstations() const { return fNstations; } 
+
   int NMvdStations {0};     ///< number of mvd stations
   int NStsStations {0};     ///< number of sts stations
   int fNfieldStations {0};  ///< number of stations in the field region
 
 
+
   // TODO: Replace _fvecalignment with C++11 alignas(16) attibute, see vStationsNew (S.Zh.)
-  L1Station vStations[L1Parameters::kMaxNstations] _fvecalignment;  // station info
-  alignas(16) std::array<L1Station, L1Parameters::kMaxNstations> fStationsNew;
+  //L1Station vStations[L1Parameters::kMaxNstations] _fvecalignment;  // station info
+
+private:
+  alignas(16) std::array<L1Station, L1Parameters::kMaxNstations> fStations;
+public:
+  const std::array<L1Station, L1Parameters::kMaxNstations>& GetStations() const { return fStations; }
+  
+  
   L1Vector<L1Material> fRadThick {"fRadThick"};        // material for each station
 
   int NStsStrips {0};                             ///> number of strips
diff --git a/reco/L1/L1Algo/L1CAMergeClones.cxx b/reco/L1/L1Algo/L1CAMergeClones.cxx
index b9f14fadcdafb75744926abc6f97d254d79c9ae6..86da1df512de2e5c75107779c50c7421c0c682b2 100644
--- a/reco/L1/L1Algo/L1CAMergeClones.cxx
+++ b/reco/L1/L1Algo/L1CAMergeClones.cxx
@@ -260,7 +260,7 @@ void L1Algo::CAMergeClones()
   L1FieldValue fBm, fBb, fBf _fvecalignment;
   L1FieldRegion fld _fvecalignment;
 
-  unsigned char maxLengthForMerge = static_cast<unsigned char>(NStations - 3);  // max length for merge
+  unsigned char maxLengthForMerge = static_cast<unsigned char>(fNstations - 3);  // max length for merge
 
 #ifdef OMP
 #pragma omp parallel for
@@ -347,8 +347,8 @@ void L1Algo::CAMergeClones()
       // if (fabs (Tf.time[0] - Tb.time[0]) > 500000) continue;
       unsigned short stam;
 
-      vStations[staf].fieldSlice.GetFieldValue(Tf.x, Tf.y, fBf);
-      vStations[stab].fieldSlice.GetFieldValue(Tb.x, Tb.y, fBb);
+      fStations[staf].fieldSlice.GetFieldValue(Tf.x, Tf.y, fBf);
+      fStations[stab].fieldSlice.GetFieldValue(Tb.x, Tb.y, fBb);
 
       unsigned short dist = firstStation[iTr] - lastStation[jTr];
 
@@ -356,10 +356,10 @@ void L1Algo::CAMergeClones()
       else
         stam = staf - 1;
 
-      fvec zm = vStations[stam].z;
+      fvec zm = fStations[stam].z;
       fvec xm = 0.5 * (Tf.x + Tf.tx * (zm - Tf.z) + Tb.x + Tb.tx * (zm - Tb.z));
       fvec ym = 0.5 * (Tf.y + Tf.ty * (zm - Tf.z) + Tb.y + Tb.ty * (zm - Tb.z));
-      vStations[stam].fieldSlice.GetFieldValue(xm, ym, fBm);
+      fStations[stam].fieldSlice.GetFieldValue(xm, ym, fBm);
       fld.Set(fBb, Tb.z, fBm, zm, fBf, Tf.z);
 
       fvec zMiddle = 0.5 * (Tb.z + Tf.z);
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index cac3f853b518f03219ce810179e10190f2a797c5..556ff6d87d2bc531562261118c5153b03b822d27 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -128,15 +128,19 @@ inline void L1Algo::f11(  /// input 1st stage of singlet search
   // comment unused parameters, FU, 18.01.21
   fvec* /*d_x*/, fvec* /*d_y*/, fvec* /*d_xy*/, fvec* d_u, fvec* d_v)
 {
-  L1Station& stal = vStations[istal];
-  L1Station& stam = vStations[istam];
+  L1Station& stal = fStations[istal];
+  L1Station& stam = fStations[istam];
   fvec zstal      = stal.z;
   fvec zstam      = stam.z;
 
+  //L1Station& stal_CHECK = fStations[istal];
+  //std::cout << "\033[1;34m>>>>>>>>>>>>\033[0m " << istal << ' ' << &stal_CHECK << ' ' << fStations.begin() << ' ' << &stal_CHECK - fStations.begin() << '\n';
+
+
   int istal_global       = 5;  //TODO: SG: what are these stations? Should the numbers depend on the presence of mvd?
   int istam_global       = 6;
-  L1Station& stal_global = vStations[istal_global];
-  L1Station& stam_global = vStations[istam_global];
+  L1Station& stal_global = fStations[istal_global];
+  L1Station& stam_global = fStations[istam_global];
   fvec zstal_global      = stal_global.z;
   fvec zstam_global      = stam_global.z;
 
@@ -185,12 +189,12 @@ inline void L1Algo::f11(  /// input 1st stage of singlet search
     int istac = istal / 2;  // "center" station
     //     if (istal >= NMvdStations) // to make it in the same way for both with and w\o mvd
     //       istac = (istal-NMvdStations)/2+NMvdStations;
-    L1Station& stac = vStations[istac];
+    L1Station& stac = fStations[istac];
     fvec zstac      = stac.z;
 
     int istac_global = istal_global / 2;  // "center" station
 
-    L1Station& stac_global = vStations[istac_global];
+    L1Station& stac_global = fStations[istac_global];
     fvec zstac_global      = stac.z;
 
     stac.fieldSlice.GetFieldValue(fTargX + tx * (zstac - fTargZ), fTargY + ty * (zstac - fTargZ), centB);
@@ -215,7 +219,7 @@ inline void L1Algo::f11(  /// input 1st stage of singlet search
       fld1.Set(m_B, zstam, l_B, zstal, fTargB, fTargZ);
 #else  // if USE_3HITS  // the best now
     L1FieldValue r_B _fvecalignment;
-    L1Station& star = vStations[istam + 1];
+    L1Station& star = fStations[istam + 1];
     fvec zstar      = star.z;
     star.fieldSlice.GetFieldValue(fTargX + tx * (zstar - fTargZ), fTargY + ty * (zstar - fTargZ), r_B);
     fld1.Set(r_B, zstar, m_B, zstam, l_B, zstal);
@@ -327,7 +331,7 @@ inline void L1Algo::f11(  /// input 1st stage of singlet search
         fit.L1AddMaterial(T, fRadThick[ista].GetRadThick(T.x, T.y), fMaxInvMom, 1);
       }
       else {
-        fit.L1AddMaterial(T, vStations[ista].materialInfo, fMaxInvMom, 1);
+        fit.L1AddMaterial(T, fStations[ista].materialInfo, fMaxInvMom, 1);
       }
       if (ista == NMvdStations - 1) fit.L1AddPipeMaterial(T, fMaxInvMom, 1);
     }
@@ -418,7 +422,7 @@ inline void L1Algo::f20(
     const float& time      = T1.t[i1_4];
 
 
-    L1HitAreaTime areaTime(vGridTime[&stam - vStations], T1.x[i1_4] * iz, T1.y[i1_4] * iz,
+    L1HitAreaTime areaTime(vGridTime[&stam - fStations.begin()], T1.x[i1_4] * iz, T1.y[i1_4] * iz,
                            (sqrt(Pick_m22 * (T1.C00 + stam.XYInfo.C00)) + fMaxDZ * fabs(T1.tx))[i1_4] * iz,
                            (sqrt(Pick_m22 * (T1.C11 + stam.XYInfo.C11)) + fMaxDZ * fabs(T1.ty))[i1_4] * iz, time,
                            sqrt(timeError) * 5);
@@ -429,7 +433,7 @@ inline void L1Algo::f20(
     while (true) {
       if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) {
         irm1++;
-        if ((L1HitIndex_t) irm1 >= (StsHitsUnusedStopIndex[&stam - vStations] - StsHitsUnusedStartIndex[&stam - vStations]))
+        if ((L1HitIndex_t) irm1 >= (StsHitsUnusedStopIndex[&stam - fStations.begin()] - StsHitsUnusedStartIndex[&stam - fStations.begin()]))
           break;
         imh = irm1;
       }
@@ -610,7 +614,7 @@ inline void L1Algo::f30(  // input
 
 
   // ---- Add the middle hits to parameters estimation. Propagate to right station. ----
-  if (istar < NStations) {
+  if (istar < fNstations) {
     for (Tindex i2 = 0; i2 < n2;) {
       Tindex n2_4 = 0;
       for (; n2_4 < fvecLen && i2 < n2; n2_4++, i2++) {
@@ -724,7 +728,7 @@ inline void L1Algo::f30(  // input
         if (isec == TRACKS_FROM_TRIPLETS_ITERATION) Pick_r22 = Pick_r2 + 1;
 #endif  // DO_NOT_SELECT_TRIPLETS
         const fscal iz = 1.f / (T2.z[i2_4] - fCbmTargetZ[0]);
-        L1HitAreaTime area(vGridTime[&star - vStations], T2.x[i2_4] * iz, T2.y[i2_4] * iz,
+        L1HitAreaTime area(vGridTime[&star - fStations.begin()], T2.x[i2_4] * iz, T2.y[i2_4] * iz,
                            (sqrt(Pick_r22 * (T2.C00 + stam.XYInfo.C00)) + fMaxDZ * fabs(T2.tx))[i2_4] * iz,
                            (sqrt(Pick_r22 * (T2.C11 + stam.XYInfo.C11)) + fMaxDZ * fabs(T2.ty))[i2_4] * iz, time,
                            sqrt(timeError) * 5);
@@ -922,7 +926,7 @@ inline void L1Algo::f31(  // input
     if (fUseHitErrors) info.sigma2 = du_[i3_V] * du_[i3_V];
 
     if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) {
-      if ((&star - vStations) < fNfieldStations) L1Filter(T_3[i3_V], info, u_front_[i3_V]);
+      if ((&star - fStations.begin()) < fNfieldStations) L1Filter(T_3[i3_V], info, u_front_[i3_V]);
       else {
         L1FilterNoField(T_3[i3_V], info, u_front_[i3_V]);
       }
@@ -936,7 +940,7 @@ inline void L1Algo::f31(  // input
     if (fUseHitErrors) info.sigma2 = dv_[i3_V] * dv_[i3_V];
 
     if (fTrackingMode == kGlobal || fTrackingMode == kMcbm) {
-      if ((&star - vStations) < fNfieldStations) L1Filter(T_3[i3_V], info, u_back_[i3_V]);
+      if ((&star - fStations.begin()) < fNfieldStations) L1Filter(T_3[i3_V], info, u_back_[i3_V]);
       else
         L1FilterNoField(T_3[i3_V], info, u_back_[i3_V]);
     }
@@ -965,7 +969,7 @@ inline void L1Algo::f32(  // input // TODO not updated after gaps introduction
 
   L1Station sta[3];
   for (int is = 0; is < NHits; ++is) {
-    sta[is] = vStations[ista[is]];
+    sta[is] = fStations[ista[is]];
   };
 
   for (int i3 = 0; i3 < n3; ++i3) {
@@ -1195,7 +1199,7 @@ inline void L1Algo::f4(  // input
 
     TripForHit[1][ihitl] = Location + 1;
 
-    if (istal > (NStations - 4)) continue;
+    if (istal > (fNstations - 4)) continue;
 
     unsigned int nNeighbours = TripForHit[1][ihitm] - TripForHit[0][ihitm];
 
@@ -1236,13 +1240,13 @@ inline void L1Algo::f5(  // input
   if (isec != TRACKS_FROM_TRIPLETS_ITERATION)
 #endif
 
-    for (int istal = NStations - 4; istal >= FIRSTCASTATION; istal--) {
+    for (int istal = fNstations - 4; istal >= FIRSTCASTATION; istal--) {
       for (int tripType = 0; tripType < 3; tripType++)  // tT = 0 - 123triplet, tT = 1 - 124triplet, tT = 2 - 134triplet
       {
         if ((((isec != kFastPrimJumpIter) && (isec != kAllPrimJumpIter) && (isec != kAllSecJumpIter))
              && (tripType != 0))
             || (((isec == kFastPrimJumpIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecJumpIter))
-                && (tripType != 0) && (istal == NStations - 4)))
+                && (tripType != 0) && (istal == fNstations - 4)))
           continue;
 
         int istam = istal + 1;
@@ -1339,9 +1343,9 @@ inline void L1Algo::DupletsStaPort(
   ///   @&i1_2 - index of 1st hit in portion indexed by doublet index
   ///   @&hitsm_2 - index of middle hit in hits array indexed by doublet index
 
-  if (istam < NStations) {
-    L1Station& stal = vStations[istal];
-    L1Station& stam = vStations[istam];
+  if (istam < fNstations) {
+    L1Station& stal = fStations[istal];
+    L1Station& stam = fStations[istam];
 
     // prepare data
     L1HitPoint* vStsHits_l = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istal];
@@ -1431,10 +1435,10 @@ L1Algo::TripletsStaPort(  /// creates triplets: input: @istal - start station nu
 
 )
 {
-  if (istar < NStations) {
+  if (istar < fNstations) {
     // prepare data
-    L1Station& stam = vStations[istam];
-    L1Station& star = vStations[istar];
+    L1Station& stam = fStations[istam];
+    L1Station& star = fStations[istar];
 
     L1HitPoint* vStsHits_m = &((*vStsHitPointsUnused)[0]) + StsHitsUnusedStartIndex[istam];
 
@@ -1672,12 +1676,12 @@ void L1Algo::CATrackFinder()
   fRecoHits.clear();
 
   fRecoHits.reserve(2 * vStsHits->size());
-  fTracks.reserve(2 * vStsHits->size() / NStations);
+  fTracks.reserve(2 * vStsHits->size() / fNstations);
 
   int nDontUsedHits = 0;
 
   // #pragma omp parallel for  reduction(+:nDontUsedHits)
-  for (int ista = 0; ista < NStations; ++ista) {
+  for (int ista = 0; ista < fNstations; ++ista) {
     nDontUsedHits += (StsHitsStopIndex[ista] - StsHitsStartIndex[ista]);
     StsHitsUnusedStartIndex[ista] = StsHitsStartIndex[ista];
     StsHitsUnusedStopIndex[ista]  = StsHitsStopIndex[ista];
@@ -1686,7 +1690,7 @@ void L1Algo::CATrackFinder()
   float lasttime  = 0;
   float starttime = std::numeric_limits<float>::max();
 
-  for (int ist = 0; ist < NStations; ++ist)
+  for (int ist = 0; ist < fNstations; ++ist)
     for (L1HitIndex_t ih = StsHitsStartIndex[ist]; ih < StsHitsStopIndex[ist]; ++ih) {
 
       const float& time = (*vStsHits)[ih].t;
@@ -1719,7 +1723,7 @@ void L1Algo::CATrackFinder()
   vStsHitsUnused = &vStsDontUsedHits_Buf;
 
 
-  for (int iS = 0; iS < NStations; ++iS) {
+  for (int iS = 0; iS < fNstations; ++iS) {
 
     vGridTime[iS].BuildBins(-1, 1, -0.6, 0.6, starttime, lasttime, xStep, yStep, (lasttime - starttime + 1));
     int start = StsHitsUnusedStartIndex[iS];
@@ -1734,14 +1738,14 @@ void L1Algo::CATrackFinder()
   }
 
 
-  for (int ist = 0; ist < NStations; ++ist)
+  for (int ist = 0; ist < fNstations; ++ist)
     for (L1HitIndex_t ih = StsHitsStartIndex[ist]; ih < StsHitsStopIndex[ist]; ++ih) {
       L1Hit& h = (*vStsHits)[ih];
       SetFUnUsed((*fStripFlag)[h.f]);
       SetFUnUsed((*fStripFlag)[h.b]);
     }
 
-  for (int ista = 0; ista < NStations; ++ista) {
+  for (int ista = 0; ista < fNstations; ++ista) {
 #ifdef _OPENMP
 #pragma omp parallel for schedule(dynamic, 5)
 #endif
@@ -1781,7 +1785,7 @@ void L1Algo::CATrackFinder()
     // n_g1.assign(n_g1.size(), Portion);
 
     for (int n = 0; n < fNThreads; n++) {
-      for (int j = 0; j < NStations; j++) {
+      for (int j = 0; j < fNstations; j++) {
         fTriplets[j][n].clear();
       }
     }
@@ -1804,7 +1808,7 @@ void L1Algo::CATrackFinder()
       vStsHitPointsUnused_buf                    = vStsHitsUnused_temp2;
     }
     // TODO: Replace NStations with fInitManager.GetStationsNumber() (S.Zharko)
-    for (int ist = 0; ist < NStations; ++ist) {
+    for (int ist = 0; ist < fNstations; ++ist) {
       for (L1HitIndex_t ih = StsHitsUnusedStartIndex[ist]; ih < StsHitsUnusedStopIndex[ist]; ++ih) {
         //SG!!
         TripForHit[0][ih] = 0;
@@ -1812,8 +1816,8 @@ void L1Algo::CATrackFinder()
       }
     }
     /*
-   TripForHit[0].assign(StsHitsUnusedStopIndex[NStations-1],0);
-   TripForHit[1].assign(StsHitsUnusedStopIndex[NStations-1],0);
+   TripForHit[0].assign(StsHitsUnusedStopIndex[fNstations-1],0);
+   TripForHit[1].assign(StsHitsUnusedStopIndex[fNstations-1],0);
 */
     {
       // #pragma omp  task
@@ -1893,7 +1897,7 @@ void L1Algo::CATrackFinder()
         // Select magnetic field. For primary tracks - vtxFieldValue, for secondary tracks - st.fieldSlice
         if (caIteration.IsPrimary()) { fTargB = vtxFieldValue; }
         else {
-          vStations[0].fieldSlice.GetFieldValue(0, 0, fTargB);
+          fStations[0].fieldSlice.GetFieldValue(0, 0, fTargB);
         }  // NOTE: calculates field fTargB in the center of 0th station
 
 
@@ -1907,7 +1911,7 @@ void L1Algo::CATrackFinder()
         //}
         //if ((isec == kAllSecIter) || (isec == kAllSecEIter)
         //    || (isec == kAllSecJumpIter)) {  //use outer radius of the 1st station as a constraint // ?
-        //  L1Station& st = vStations[0];
+        //  L1Station& st = fStations[0];
         //  SigmaTargetX = SigmaTargetY = 10;  //st.Rmax[0];
         //  fTargZ                      = fCbmTargetZ;  // fCbmTargetZ-1.;
         //  st.fieldSlice.GetFieldValue(0, 0, fTargB);
@@ -1927,11 +1931,11 @@ void L1Algo::CATrackFinder()
         //  fMaxDZ = 0.1;
 
         // TODO: to be removed, because this condition is checked in L1InitManager (S.Zharko)
-        if (NStations > (int) L1Parameters::kMaxNstations) cout << " CATrackFinder: Error: Too many Stations" << endl;
+        if (fNstations > (int) L1Parameters::kMaxNstations) cout << " CATrackFinder: Error: Too many Stations" << endl;
       }
 
 #ifndef L1_NO_ASSERT
-      for (int istal = NStations - 1; istal >= 0; istal--) {
+      for (int istal = fNstations - 1; istal >= 0; istal--) {
         L1_ASSERT(StsHitsUnusedStopIndex[istal] >= StsHitsUnusedStartIndex[istal],
                   StsHitsUnusedStopIndex[istal] << " >= " << StsHitsUnusedStartIndex[istal]);
         L1_ASSERT(StsHitsUnusedStopIndex[istal] <= (*vStsHitsUnused).size(),
@@ -1943,9 +1947,9 @@ void L1Algo::CATrackFinder()
 
     {
       /// possible left hits of triplets are splited in portions of 16 (4 SIMDs) to use memory faster
-      fDupletPortionStopIndex[NStations - 1] = 0;
+      fDupletPortionStopIndex[fNstations - 1] = 0;
       fDupletPortionSize.clear();
-      for (int istal = NStations - 2; istal >= FIRSTCASTATION; istal--) {  //start downstream chambers
+      for (int istal = fNstations - 2; istal >= FIRSTCASTATION; istal--) {  //start downstream chambers
         int NHits_l   = StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal];
         int nPortions = NHits_l / Portion;
         int rest      = NHits_l - nPortions * Portion;
@@ -1964,10 +1968,10 @@ void L1Algo::CATrackFinder()
 
     /*    {
          /// possible left hits of triplets are splited in portions of 16 (4 SIMDs) to use memory faster
-         fDupletPortionStopIndex[NStations-1] = 0;
+         fDupletPortionStopIndex[fNstations-1] = 0;
          unsigned int ip = 0;  //index of curent portion
          
-         for (int istal = NStations-2; istal >= FIRSTCASTATION; istal--)  //start downstream chambers
+         for (int istal = fNstations-2; istal >= FIRSTCASTATION; istal--)  //start downstream chambers
          {
          int nHits = StsHitsUnusedStopIndex[istal] - StsHitsUnusedStartIndex[istal];
          
@@ -2016,7 +2020,7 @@ void L1Algo::CATrackFinder()
     L1Vector<char> lmDupletsG[L1Parameters::kMaxNstations] {
       "L1CATrackFinder::lmDupletsG"};  // is exist a doublet started from indexed by left hit
 
-    for (int i = 0; i < NStations; i++) {
+    for (int i = 0; i < fNstations; i++) {
       lmDuplets[i].SetName(std::stringstream() << "L1CATrackFinder::lmDuplets[" << i << "]");
       lmDupletsG[i].SetName(std::stringstream() << "L1CATrackFinder::lmDupletsG[" << i << "]");
     }
@@ -2026,7 +2030,7 @@ void L1Algo::CATrackFinder()
     hitsmG_2.reserve(9000);
     i1G_2.reserve(9000);
 
-    for (int istal = NStations - 2; istal >= FIRSTCASTATION; istal--)  //  //start downstream chambers
+    for (int istal = fNstations - 2; istal >= FIRSTCASTATION; istal--)  //  //start downstream chambers
     {
 
 #ifdef _OPENMP
@@ -2097,7 +2101,7 @@ void L1Algo::CATrackFinder()
     }
 
     //     int nlevels[L1Parameters::kMaxNstations];  // number of triplets with some number of neighbours.
-    //     for (int il = 0; il < NStations; ++il) nlevels[il] = 0;
+    //     for (int il = 0; il < fNstations; ++il) nlevels[il] = 0;
     //
     //      f5(   // input
     //           // output
@@ -2133,7 +2137,7 @@ void L1Algo::CATrackFinder()
     //     int min_level = 1; // min level for start triplet. So min track length = min_level+3.
     //     if (isec == kAllPrimJumpIter) min_level = 1;
     //     if ( (isec == kAllSecIter) || (isec == kAllSecJumpIter) ) min_level = 2; // only the long low momentum tracks
-    // //    if (isec == -1) min_level = NStations-3 - 3; //only the longest tracks
+    // //    if (isec == -1) min_level = fNstations-3 - 3; //only the longest tracks
 
 
     L1Branch curr_tr;
@@ -2149,12 +2153,12 @@ void L1Algo::CATrackFinder()
 
 
     // collect consequtive: the longest tracks, shorter, more shorter and so on
-    for (int ilev = NStations - 3; ilev >= min_level; ilev--) {
+    for (int ilev = fNstations - 3; ilev >= min_level; ilev--) {
       // choose length in triplets number - ilev - the maximum possible triplet level among all triplets in the searched candidate
       TStopwatch Time;
 
       //  how many levels to check
-      int nlevel = (NStations - 2) - ilev + 1;
+      int nlevel = (fNstations - 2) - ilev + 1;
 
       const unsigned char min_best_l = (ilev > min_level) ? ilev + 2 : min_level + 3;  // loose maximum
 
@@ -2166,7 +2170,7 @@ void L1Algo::CATrackFinder()
       fStripToTrack.reset(NStsStrips, -1);
       fStripToTrackB.reset(NStsStrips, -1);
 
-      for (int istaF = FIRSTCASTATION; istaF <= NStations - 3 - ilev; ++istaF) {
+      for (int istaF = FIRSTCASTATION; istaF <= fNstations - 3 - ilev; ++istaF) {
 
 #ifdef _OPENMP
 #pragma omp parallel for firstprivate(curr_tr, new_tr, best_tr, curr_chi2, best_chi2, best_L, curr_L,                  \
@@ -2359,7 +2363,7 @@ void L1Algo::CATrackFinder()
             if (check) {
 #ifdef EXTEND_TRACKS
               if (fTrackingMode != kMcbm)
-                if (tr.NHits != NStations) BranchExtender(tr);
+                if (tr.NHits != fNstations) BranchExtender(tr);
 #endif
               float sumTime = 0;
 
@@ -2390,7 +2394,7 @@ void L1Algo::CATrackFinder()
                 L1HitPoint tempPoint = CreateHitPoint(hit);  //TODO take number of station from hit
 
                 float xcoor, ycoor = 0;
-                L1Station stah = vStations[0];
+                L1Station stah = fStations[0];
                 StripsToCoor(tempPoint.U(), tempPoint.V(), xcoor, ycoor, stah);
                 float zcoor = tempPoint.Z() - fCbmTargetZ[0];
 
@@ -2455,7 +2459,7 @@ void L1Algo::CATrackFinder()
     if (isec < (fNFindIterations - 1)) {
       int NHitsOnStation = 0;
 
-      for (int ista = 0; ista < NStations; ++ista) {
+      for (int ista = 0; ista < fNstations; ++ista) {
         int start                = StsHitsUnusedStartIndex[ista];
         int Nelements            = StsHitsUnusedStopIndex[ista] - start;
         L1Hit* staHits           = nullptr;  // to avoid out-of-range error in ..[start]
diff --git a/reco/L1/L1Algo/L1Filtration.h b/reco/L1/L1Algo/L1Filtration.h
index 0f2e448cef8f54231fd5abacea56446fd0237a17..df438092cb41e48c633c77d91c35e1b45e3054ca 100644
--- a/reco/L1/L1Algo/L1Filtration.h
+++ b/reco/L1/L1Algo/L1Filtration.h
@@ -81,7 +81,7 @@ inline void FilterTime(L1TrackPar& T, fvec t0, fvec dt0, fvec timeInfo = 1., fve
   T.C55 -= K5 * F5;
 }
 
-inline void L1Filter(L1TrackPar& T, L1UMeasurementInfo& info, fvec u, fvec w = 1.)
+inline void L1Filter(L1TrackPar& T, const L1UMeasurementInfo& info, fvec u, fvec w = 1.)
 {
   fvec wi, zeta, zetawi, HCH;
   fvec F0, F1, F2, F3, F4, F5;
@@ -148,7 +148,7 @@ inline void L1Filter(L1TrackPar& T, L1UMeasurementInfo& info, fvec u, fvec w = 1
   T.C55 -= K5 * F5;
 }
 
-inline void L1FilterNoField(L1TrackPar& T, L1UMeasurementInfo& info, fvec u, fvec w = 1.)
+inline void L1FilterNoField(L1TrackPar& T, const L1UMeasurementInfo& info, fvec u, fvec w = 1.)
 {
   fvec wi, zeta, zetawi, HCH;
   fvec F0, F1, F2, F3, F4, F5;
@@ -265,7 +265,7 @@ inline void L1FilterChi2XYC00C10C11(const L1UMeasurementInfo& info, fvec& x, fve
   C11 -= K1 * F1;
 }
 
-inline void L1FilterVtx(L1TrackPar& T, fvec x, fvec y, L1XYMeasurementInfo& info, fvec extrDx, fvec extrDy, fvec J[])
+inline void L1FilterVtx(L1TrackPar& T, fvec x, fvec y, const L1XYMeasurementInfo& info, fvec extrDx, fvec extrDy, fvec J[])
 {
   cnst TWO = 2.;
 
@@ -343,7 +343,7 @@ inline void L1FilterVtx(L1TrackPar& T, fvec x, fvec y, L1XYMeasurementInfo& info
 }
 
 
-inline void L1FilterXY(L1TrackPar& T, fvec x, fvec y, L1XYMeasurementInfo& info)
+inline void L1FilterXY(L1TrackPar& T, fvec x, fvec y, const L1XYMeasurementInfo& info)
 {
   cnst TWO = 2.;
 
diff --git a/reco/L1/L1Algo/L1Fit.h b/reco/L1/L1Algo/L1Fit.h
index a40bbc2096edc2da3007e443dd89c54e206d5ecd..c56928f3433a16e982f8200ac9d3af11ea2ca5fa 100644
--- a/reco/L1/L1Algo/L1Fit.h
+++ b/reco/L1/L1Algo/L1Fit.h
@@ -64,12 +64,12 @@ public:
 
   void L1AddMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w);
 
-  void L1AddMaterial(L1TrackPar& T, L1MaterialInfo& info, fvec qp0, fvec w);
+  void L1AddMaterial(L1TrackPar& T, const L1MaterialInfo& info, fvec qp0, fvec w);
 
 
   void L1AddThickMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream);
 
-  void L1AddHalfMaterial(L1TrackPar& T, L1MaterialInfo& info, fvec qp0);
+  void L1AddHalfMaterial(L1TrackPar& T, const L1MaterialInfo& info, fvec qp0);
 
   void L1AddPipeMaterial(L1TrackPar& T, fvec qp0, fvec w);
 
diff --git a/reco/L1/L1Algo/L1FitMaterial.h b/reco/L1/L1Algo/L1FitMaterial.h
index 24a769613105b32d651d1b75f0d929c3c279a8b0..e5a2410c2da1aad2f77acf06cbac4d4829bfcb5a 100644
--- a/reco/L1/L1Algo/L1FitMaterial.h
+++ b/reco/L1/L1Algo/L1FitMaterial.h
@@ -533,7 +533,7 @@ inline void L1Fit::L1AddMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fvec w)
   T.C33 += w * (ONE + tyty) * a;
 }
 
-inline void L1Fit::L1AddMaterial(L1TrackPar& T, L1MaterialInfo& info, fvec qp0, fvec w)
+inline void L1Fit::L1AddMaterial(L1TrackPar& T, const L1MaterialInfo& info, fvec qp0, fvec w)
 {
   cnst ONE = 1.f;
 
@@ -640,7 +640,7 @@ inline void L1Fit::L1AddThickMaterial(L1TrackPar& T, fvec radThick, fvec qp0, fv
 }
 
 
-inline void L1Fit::L1AddHalfMaterial(L1TrackPar& T, L1MaterialInfo& info, fvec qp0)
+inline void L1Fit::L1AddHalfMaterial(L1TrackPar& T, const L1MaterialInfo& info, fvec qp0)
 {
   cnst ONE   = 1.f;
   fvec tx    = T.tx;
diff --git a/reco/L1/L1Algo/L1TrackExtender.cxx b/reco/L1/L1Algo/L1TrackExtender.cxx
index 8d3bde8615c31f6e0caf8344379849a10a7399e4..8b370d9955107a49267e127c956008a845c3a0be 100644
--- a/reco/L1/L1Algo/L1TrackExtender.cxx
+++ b/reco/L1/L1Algo/L1TrackExtender.cxx
@@ -48,9 +48,9 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
   int ista1 = GetFStation((*fStripFlag)[hit1.f]);
   int ista2 = GetFStation((*fStripFlag)[hit2.f]);
 
-  L1Station& sta0 = vStations[ista0];
-  L1Station& sta1 = vStations[ista1];
-  L1Station& sta2 = vStations[ista2];
+  L1Station& sta0 = fStations[ista0];
+  L1Station& sta1 = fStations[ista1];
+  L1Station& sta2 = fStations[ista2];
 
   fvec u0 = hit0.u;
   fvec v0 = hit0.v;
@@ -126,7 +126,7 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
     ista_prev        = ista;
     ista             = GetFStation((*fStripFlag)[hit.f]);
 
-    L1Station& sta = vStations[ista];
+    L1Station& sta = fStations[ista];
 
     float z_sta = hit.z;
 
@@ -197,7 +197,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
                           const fvec qp0)  // TODO take into account pipe
 {
   L1Vector<L1HitIndex_t> newHits {"L1TrackExtender::newHits"};
-  newHits.reserve(NStations);
+  newHits.reserve(fNstations);
   L1Fit fit;
   fit.SetParticleMass(GetDefaultParticleMass());
 
@@ -213,9 +213,9 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
   const int ista1 = GetFStation((*fStripFlag)[hit1.f]);
   const int ista2 = GetFStation((*fStripFlag)[hit2.f]);
 
-  const L1Station& sta0 = vStations[ista0];
-  const L1Station& sta1 = vStations[ista1];
-  const L1Station& sta2 = vStations[ista2];
+  const L1Station& sta0 = fStations[ista0];
+  const L1Station& sta1 = fStations[ista1];
+  const L1Station& sta2 = fStations[ista2];
 
   fvec u0 = hit0.u;
   fvec v0 = hit0.v;
@@ -250,9 +250,9 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
 
   const fvec pickGather2 = fPickGather * fPickGather;
 
-  for (; (ista < NStations) && (ista >= 0); ista += step) {  // CHECKME why ista2?
+  for (; (ista < fNstations) && (ista >= 0); ista += step) {  // CHECKME why ista2?
 
-    L1Station& sta = vStations[ista];
+    L1Station& sta = fStations[ista];
 
     fvec dz = sta.z - T.z;
 
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index a8354b19327f940aae43f83116e192ba215a7d10..5304cc071ce7053ed9e37ca8bb9057313fb7cbe0 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -61,9 +61,9 @@ void L1Algo::KFTrackFitter_simple()  // TODO: Add pipe.
         int ista2 = (*fStripFlag)[hit2.f] / 4;
 
         //cout<<"back: ista012="<<ista0<<" "<<ista1<<" "<<ista2<<endl;
-        L1Station& sta0 = vStations[ista0];
-        L1Station& sta1 = vStations[ista1];
-        L1Station& sta2 = vStations[ista2];
+        L1Station& sta0 = fStations[ista0];
+        L1Station& sta1 = fStations[ista1];
+        L1Station& sta2 = fStations[ista2];
 
         fvec u0 = hit0.u;
         fvec v0 = hit0.v;
@@ -129,7 +129,7 @@ void L1Algo::KFTrackFitter_simple()  // TODO: Add pipe.
           const L1Hit& hit = (*vStsHits)[hits[i]];
           ista             = (*fStripFlag)[hit.f] / 4;
 
-          L1Station& sta = vStations[ista];
+          L1Station& sta = fStations[ista];
 
           //    L1Extrapolate( T, hit.z, qp0, fld );
           L1ExtrapolateLine(T, hit.z);
@@ -202,9 +202,9 @@ void L1Algo::KFTrackFitter_simple()  // TODO: Add pipe.
         int ista1 = GetFStation((*fStripFlag)[hit1.f]);
         int ista2 = GetFStation((*fStripFlag)[hit2.f]);
 
-        L1Station& sta0 = vStations[ista0];
-        L1Station& sta1 = vStations[ista1];
-        L1Station& sta2 = vStations[ista2];
+        L1Station& sta0 = fStations[ista0];
+        L1Station& sta1 = fStations[ista1];
+        L1Station& sta2 = fStations[ista2];
 
         fvec u0 = hit0.u;
         fvec v0 = hit0.v;
@@ -264,7 +264,7 @@ void L1Algo::KFTrackFitter_simple()  // TODO: Add pipe.
         for (int i = 1; i < nHits; i++) {
           const L1Hit& hit = (*vStsHits)[hits[i]];
           ista             = (*fStripFlag)[hit.f] / 4;
-          L1Station& sta   = vStations[ista];
+          L1Station& sta   = fStations[ista];
           fvec u           = hit.u;
           fvec v           = hit.v;
           fvec x, y;
@@ -339,7 +339,7 @@ void L1Algo::L1KFTrackFitter()
   L1FieldValue fldB01, fldB11, fldB21 _fvecalignment;
   L1FieldRegion fld1 _fvecalignment;
 
-  const int nHits = NStations;
+  const int nHits = fNstations;
   int iVec = 0, i = 0;
   int nTracks_SIMD = fvecLen;
   L1TrackPar T;  // fitting parametr coresponding to current track
@@ -352,7 +352,7 @@ void L1Algo::L1KFTrackFitter()
 
   L1Track* t[fvecLen];
 
-  L1Station* sta = vStations;
+  L1Station* sta = fStations.begin();
   L1Station staFirst, staLast;
   fvec x[L1Parameters::kMaxNstations], u[L1Parameters::kMaxNstations], v[L1Parameters::kMaxNstations],
     y[L1Parameters::kMaxNstations], time[L1Parameters::kMaxNstations], timeEr[L1Parameters::kMaxNstations],
@@ -445,7 +445,7 @@ void L1Algo::L1KFTrackFitter()
       fscal cursy = 0., curSy = 0.;
       for (i = nHitsTrack - 1; i >= 0; i--) {
         const int ista   = iSta[i];
-        L1Station& st    = vStations[ista];
+        L1Station& st    = fStations[ista];
         const fscal curZ = z[ista][iVec];
         fscal dZ         = curZ - prevZ;
         fscal By         = st.fieldSlice.cy[0][0];
@@ -765,7 +765,7 @@ void L1Algo::L1KFTrackFitterMuch()
   L1FieldValue fldB01, fldB11, fldB21 _fvecalignment;
   L1FieldRegion fld1 _fvecalignment;
 
-  const int nHits = NStations;
+  const int nHits = fNstations;
   int iVec = 0, i = 0;
   int nTracks_SIMD = fvecLen;
   L1TrackPar T;  // fitting parametr coresponding to current track
@@ -778,7 +778,7 @@ void L1Algo::L1KFTrackFitterMuch()
 
   L1Track* t[fvecLen];
 
-  L1Station* sta = vStations;
+  L1Station* sta = fStations.begin();
   L1Station staFirst, staLast;
   fvec x[L1Parameters::kMaxNstations], u[L1Parameters::kMaxNstations], v[L1Parameters::kMaxNstations],
     y[L1Parameters::kMaxNstations], time[L1Parameters::kMaxNstations], timeEr[L1Parameters::kMaxNstations],
@@ -813,7 +813,7 @@ void L1Algo::L1KFTrackFitterMuch()
     }
 
     for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
-      for (i = 0; i < NStations; i++) {
+      for (i = 0; i < fNstations; i++) {
         d_x[i][iVec] = 0;
         d_y[i][iVec] = 0;
       }
@@ -884,7 +884,7 @@ void L1Algo::L1KFTrackFitterMuch()
       //   for(i = nHitsTrack - 1; i >= 0; i-- ){
       for (i = 0; i < nHitsTrackSts; i++) {
         const int ista   = iSta[i];
-        L1Station& st    = vStations[ista];
+        L1Station& st    = fStations[ista];
         const fscal curZ = z[ista][iVec];
         fscal dZ         = curZ - prevZ;
         fscal By         = st.fieldSlice.cy[0][0];
diff --git a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
index b83b8eeeacd2116a9b9eec700a0357884db94966..be7519a930fb0294127fd1fd0d3ae35ee1145929 100644
--- a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
+++ b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
@@ -78,11 +78,11 @@ void L1AlgoDraw::InitL1Draw(L1Algo* algo_)
   for (unsigned int i = 0; i < algo->vStsHits->size(); i++) {
     vStsHits.push_back((*algo->vStsHits)[i]);
   }
-  NStations = algo->NStations;
+  NStations = algo->GetNstations();
   for (int i = 0; i < NStations; i++) {
     StsHitsStartIndex[i] = algo->StsHitsStartIndex[i];
     StsHitsStopIndex[i]  = algo->StsHitsStopIndex[i];
-    vStations[i]         = algo->vStations[i];
+    vStations[i]         = algo->GetStations()[i];
   }
 }
 
diff --git a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
index 436f402dd556efdea0ce3a014ae41686fa0bc326..f0927c0cab9cc68e384a963584b2e6a53ebbfa24 100644
--- a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
+++ b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
@@ -97,7 +97,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
   TClonesArray* listMvdHits = 0;
   if (NMvdStations > 0.) listMvdHits = (TClonesArray*) fManger->GetObject("MvdHit");
 
-  static int nHits = CbmL1::Instance()->algo->NStations;
+  static int nHits = CbmL1::Instance()->algo->GetNstations();
   int iVec = 0, i = 0;
   int nTracks_SIMD = fvecLen;
   L1TrackPar T;  // fitting parametr coresponding to current track
@@ -108,7 +108,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
   CbmStsTrack* t[fvecLen];
 
   int ista;
-  L1Station* sta = CbmL1::Instance()->algo->vStations;
+  const L1Station* sta = CbmL1::Instance()->algo->GetStations().begin();
   L1Station staFirst, staLast;
   fvec* x = new fvec[nHits];
   fvec* u = new fvec[nHits];
@@ -405,9 +405,9 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<L1FieldRe
 
   CbmStsTrack* t[fvecLen];
 
-  int nStations    = CbmL1::Instance()->algo->NStations;
+  int nStations    = CbmL1::Instance()->algo->GetNstations();
   int NMvdStations = CbmL1::Instance()->algo->NMvdStations;
-  L1Station* sta   = CbmL1::Instance()->algo->vStations;
+  const L1Station* sta   = CbmL1::Instance()->algo->GetStations().begin();
   fvec* zSta       = new fvec[nStations];
   for (int iSta = 0; iSta < nStations; iSta++)
     zSta[iSta] = sta[iSta].z;
@@ -581,7 +581,7 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<L1F
   CbmStsTrack* t[fvecLen];
 
   int ista;
-  L1Station* sta = CbmL1::Instance()->algo->vStations;
+  const L1Station* sta = CbmL1::Instance()->algo->GetStations().begin();
   L1FieldValue fB[3], fB_temp _fvecalignment;
   fvec zField[3];
 
@@ -653,7 +653,7 @@ void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks,
   CbmStsTrack* t[fvecLen];
 
   int ista;
-  L1Station* sta = CbmL1::Instance()->algo->vStations;
+  const L1Station* sta = CbmL1::Instance()->algo->GetStations().begin();
   L1FieldValue fB[3], fB_temp _fvecalignment;
   fvec zField[3];