CA: de-SIMDIze and clean up the triplet constructor

reco/L1/CMakeLists.txt

@@ -34,6 +34,7 @@ set(SRCS
   L1Algo/utils/CaAlgoRandom.cxx
   L1Algo/L1Algo.cxx
   L1Algo/L1TrackPar.cxx
+  L1Algo/L1TripletConstructor.cxx
   L1Algo/L1CaTrackFinder.cxx
   L1Algo/L1CaTrackFinderSlice.cxx
   L1Algo/L1BranchExtender.cxx
@@ -205,7 +206,6 @@ install(FILES CbmL1Counters.h
   L1Algo/L1Branch.h
   L1Algo/L1HitPoint.h
   L1Algo/L1Hit.h
-  L1Algo/L1Portion.h
   L1Algo/L1Triplet.h
   L1Algo/L1Event.h
   L1Algo/L1EventMatch.h

reco/L1/L1Algo/L1Algo.cxx

@@ -34,31 +34,6 @@ void L1Algo::SetNThreads(unsigned int n)
 
     fTrackCandidates[i].SetName(std::stringstream() << "L1Algo::fTrackCandidates[" << i << "]");
 
-
-    fTripletPar[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-
-    fTripletHitsL[i].SetName(std::stringstream() << "L1Algo::fhitsl_3[" << i << "]");
-    fTripletHitsM[i].SetName(std::stringstream() << "L1Algo::fhitsm_3[" << i << "]");
-    fTripletHitsR[i].SetName(std::stringstream() << "L1Algo::fhitsr_3[" << i << "]");
-
-    fTripletHitsL[i].reserve(L1Constants::size::kMaxPortionTriplets);
-    fTripletHitsM[i].reserve(L1Constants::size::kMaxPortionTriplets);
-    fTripletHitsR[i].reserve(L1Constants::size::kMaxPortionTriplets);
-
-    fTripletHitR_Ufront[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-    fTripletHitR_Uback[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-    fTripletHitR_Z[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-    fTripletHitR_Time[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-    fTripletHitR_TimeErr[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-    fTripletHitR_dUfront[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-    fTripletHitR_dUback[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-
-    fTripletHitR_X[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-    fTripletHitR_Y[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-    fTripletHitR_dX[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-    fTripletHitR_dXY[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-    fTripletHitR_dY[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-
     for (unsigned int j = 0; j < L1Constants::size::kMaxNstations; j++) {
       fTriplets[j][i].SetName(std::stringstream() << "L1Algo::fTriplets[" << i << "][" << j << "]");
     }
@@ -128,12 +103,6 @@ void L1Algo::ResetSliceData()
   fSliceRecoHits.clear();
   fSliceRecoHits.reserve(2 * nHits);
 
-  for (int iS = 0; iS < L1Constants::size::kMaxNstations; iS++) {
-    int nHitsStation = fSliceHitIds[iS].size();
-    fSingletPortionSize[iS].clear();
-    fSingletPortionSize[iS].reserve(2 * nHitsStation);
-  }
-
   for (int iT = 0; iT < fNThreads; iT++) {
     fSliceRecoTracks_thread[iT].clear();
     fSliceRecoTracks_thread[iT].reserve(nHits / 10);

reco/L1/L1Algo/L1Algo.h

@@ -37,7 +37,6 @@ class L1AlgoDraw;
 #include "L1HitPoint.h"
 #include "L1InputData.h"
 #include "L1Parameters.h"
-#include "L1Portion.h"
 #include "L1Station.h"
 #include "L1Track.h"
 #include "L1TrackPar.h"
@@ -88,8 +87,7 @@ public:
   // *************************
 
   friend class CbmL1;  /// TODO: Remove this dependency
-  friend class ParalleledDup;
-  friend class ParalleledTrip;
+  friend class L1TripletConstructor;
 #ifdef DRAW
   friend class L1AlgoDraw;
 #endif
@@ -261,82 +259,6 @@ public:
     flag = iStation * 4 + (flag % 4);
   }
 
-  /// Prepare the portion of left hits data
-  void findSingletsStep0(  // input
-    Tindex start_lh, Tindex n1_l, L1HitPoint* Hits_l,
-    // output
-    L1HitIndex_t* hitsl, fvec* x_l, fvec* y_l, fvec* z_l, fvec* t_l, fvec* dx_l, fvec* dxy_l, fvec* dy_l, fvec* dt2_l);
-
-  /// Get the field approximation. Add the target to parameters estimation. Propagate to middle station.
-  void findSingletsStep1(  // input
-    int istal, int istam, Tindex n1_V, fvec* x_l, fvec* y_l, fvec* zPos_l, fvec* t_l, fvec* dx_l, fvec* dxy_l,
-    fvec* dy_l, fvec* dt2_l,
-    // output
-    L1TrackPar* T_1, L1FieldRegion* fld_1);
-
-  /// Find the doublets. Reformat data in the portion of doublets.
-  void findDoubletsStep0(  // input
-    Tindex n1, const L1Station& stal, L1HitPoint* vHits_l, const L1Station& stam, L1HitPoint* vHits_m, L1TrackPar* T_1,
-    L1HitIndex_t* hitsl_1,
-    // output
-    Tindex& n2, L1Vector<L1HitIndex_t>& i1_2,
-#ifdef DOUB_PERFORMANCE
-    L1Vector<L1HitIndex_t>& hitsl_2,
-#endif  // DOUB_PERFORMANCE
-    L1Vector<L1HitIndex_t>& hitsm_2);
-
-  /// Add the middle hits to parameters estimation. Propagate to right station.
-  /// Find the triplets (right hit). Reformat data in the portion of triplets.
-  void findTripletsStep0(  // input
-    L1HitPoint* vHits_r, int istal, int istam, int istar, L1HitPoint* vHits_m, L1TrackPar* T_1, L1FieldRegion* fld_1,
-    L1HitIndex_t* hitsl_1,
-
-    Tindex n2, L1Vector<L1HitIndex_t>& hitsm_2, L1Vector<L1HitIndex_t>& i1_2,
-
-    // output
-    Tindex& n3, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3, L1Vector<L1HitIndex_t>& hitsm_3,
-    L1Vector<L1HitIndex_t>& hitsr_3,
-
-    L1Vector<fvec>& x_3, L1Vector<fvec>& y_3, L1Vector<fvec>& z_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dx_3,
-    L1Vector<fvec>& dxy_3, L1Vector<fvec>& dy_3, L1Vector<fvec>& dt2_3);
-
-  /// Add the right hits to parameters estimation.
-  void findTripletsStep1(  // input
-    Tindex n3_V, const L1Station& star, L1Vector<fvec>& x_3, L1Vector<fvec>& y_3, L1Vector<fvec>& z_3,
-    L1Vector<fvec>& t_3, L1Vector<fvec>& dx_3, L1Vector<fvec>& dxy_3, L1Vector<fvec>& dy_3, L1Vector<fvec>& dt2_3,
-    // output
-    L1Vector<L1TrackPar>& T_3);
-
-  /// Refit Triplets.
-  void findTripletsStep2(  // input
-    Tindex n3, int istal, int istam, int istar, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3,
-    L1Vector<L1HitIndex_t>& hitsm_3, L1Vector<L1HitIndex_t>& hitsr_3, int nIterations = 0);
-
-  /// Select triplets. Save them into vTriplets.
-  void findTripletsStep3(  // input
-    Tindex n3, int istal, int istam, int istar, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3,
-    L1Vector<L1HitIndex_t>& hitsm_3, L1Vector<L1HitIndex_t>& hitsr_3);
-
-
-  /// Find neighbours of triplets. Calculate level of triplets.
-  void f5(  // input
-    // output
-    int* nlevel);
-
-
-  /// Find doublets on station
-  void CreatePortionOfDoublets(  // input
-    int istal, int istam, Tindex iSingletPortion, Tindex singletPortionSize,
-    // output
-    L1TrackPar* T_1, L1FieldRegion* fld_1, L1HitIndex_t* hitsl_1,
-    //
-    Tindex& n_2, L1Vector<L1HitIndex_t>& i1_2, L1Vector<L1HitIndex_t>& hitsm_2);
-
-  /// Find triplets on station
-  void CreatePortionOfTriplets(  // input
-    int istal, int istam, int istar, L1TrackPar* T_1, L1FieldRegion* fld_1, L1HitIndex_t* hitsl_1, Tindex& n_2,
-    L1Vector<L1HitIndex_t>& i1_2, L1Vector<L1HitIndex_t>& hitsm_2);
-
 #ifdef DRAW
   L1AlgoDraw* draw {nullptr};
   void DrawRecoTracksTime(const L1Vector<CbmL1Track>& tracks);
@@ -432,10 +354,6 @@ public:
   /// The candidates may share any amount of hits.
   L1Vector<L1Branch> fTrackCandidates[L1Constants::size::kMaxNthreads] {"L1Algo::fTrackCandidates"};
 
-  L1Vector<Tindex> fSingletPortionSize[L1Constants::size::kMaxNstations] {
-    "L1Algo::fSingletPortionSize"};  ///< Number of doublets in a portion
-
-
   ///< indices of the sub-slice hits
   L1Vector<L1HitIndex_t> fSliceHitIds[L1Constants::size::kMaxNstations] {"L1Algo::fSliceHitIds"};
 
@@ -476,36 +394,6 @@ public:
   L1Vector<int> fHitNtriplets {"L1Algo::fHitNtriplets"};        /// link hit ->n triplets { hit, *, *}
 
 
-  //  fvec u_front[Portion/fvecLen], u_back[Portion/fvecLen];
-  //  fvec zPos[Portion/fvecLen];
-  //  fvec fHitTime[Portion/fvecLen];
-
-  L1Vector<L1TrackPar> fTripletPar[L1Constants::size::kMaxNthreads] {"fTripletPar"};
-
-  L1Vector<L1HitIndex_t> fTripletHitsL[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitsL"};
-  L1Vector<L1HitIndex_t> fTripletHitsM[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitsM"};
-  L1Vector<L1HitIndex_t> fTripletHitsR[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitsR"};
-
-  L1Vector<fvec> fTripletHitR_Ufront[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_Ufront"};
-  L1Vector<fvec> fTripletHitR_Uback[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_Uback"};
-  L1Vector<fvec> fTripletHitR_Z[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_Z"};
-  L1Vector<fvec> fTripletHitR_Time[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_Time"};
-  L1Vector<fvec> fTripletHitR_TimeErr[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_TimeErr"};
-  L1Vector<fvec> fTripletHitR_dUfront[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dUfront"};
-  L1Vector<fvec> fTripletHitR_dUback[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dUback"};
-
-  L1Vector<fvec> fTripletHitR_X[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_X"};
-  L1Vector<fvec> fTripletHitR_Y[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_Y"};
-  L1Vector<fvec> fTripletHitR_dX[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dX"};
-  L1Vector<fvec> fTripletHitR_dXY[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dXY"};
-  L1Vector<fvec> fTripletHitR_dY[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dY"};
-
-
-  //   Tindex NHits_l[L1Constants::size::kMaxNstations];
-  //   Tindex NHits_l_P[L1Constants::size::kMaxNstations];
-  /// ----- Output data -----
-
-
 private:
   L1CloneMerger fCloneMerger {*this};  ///< Object of L1Algo clones merger algorithm
 

reco/L1/L1Algo/L1BaseStationInfo.cxx

@@ -205,6 +205,7 @@ void L1BaseStationInfo::SetFieldFunction(
     fL1Station.fieldSlice.cy[j] = A[j][N + 1] / A[j][j];
     fL1Station.fieldSlice.cz[j] = A[j][N + 2] / A[j][j];
   }
+  fL1Station.fieldSlice.z = fZref;
 
   fInitController.SetFlag(EInitKey::kFieldSlice);
 }
@@ -285,7 +286,7 @@ void L1BaseStationInfo::SetXmax(double aSize)
 //
 void L1BaseStationInfo::SetYmax(double aSize)
 {
-  fYmax = aSize;
+  fYmax           = aSize;
   fL1Station.Ymax = aSize;
   fInitController.SetFlag(EInitKey::kYmax);
 }

reco/L1/L1Algo/L1BaseStationInfo.h

@@ -99,6 +99,7 @@ public:
   /// Gets array of the coefficients for the field z-component approximation
   const fvec* GetFieldSliceCz() const { return fL1Station.fieldSlice.cz; }
 
+
   /// Gets field status: 0 - station is outside the field, 1 - station is inside the field
   int GetFieldStatus() const { return fL1Station.fieldStatus; }
 

reco/L1/L1Algo/L1Constants.h

@@ -49,8 +49,6 @@ namespace L1Constants
 
     // TODO: Clarify the meaning of these coefficients
     constexpr int kCoeff                 = 64 / 4;                        ///< TODO:
-    constexpr int kSingletPortionSize    = 1024 / kCoeff;                 ///< portion of left hits
-    constexpr int kSingletPortionSizeVec = 1024 / kCoeff / fvec::size();  ///< portion of left hits per one vector word
     constexpr int kMaxPortionDoublets    = 10000 / 5 * 64 / 2 / kCoeff;   ///< Max size of the doublets portion
     constexpr int kMaxPortionTriplets    = 10000 * 5 * 64 / 2 / kCoeff;   ///< Max size of the triplets portion
     constexpr int kMaxPortionTripletsP   = kMaxPortionTriplets / fvec::size();  ///< Max size of the triplets portion

reco/L1/L1Algo/L1Field.cxx

@@ -11,6 +11,7 @@
 #include <iostream>
 #include <sstream>
 
+#include "L1TrackPar.h"
 #include "L1Utils.h"
 
 //
@@ -60,10 +61,11 @@ std::ostream& operator<<(std::ostream& out, const L1FieldValue& B)
 L1FieldSlice::L1FieldSlice()
 {
   for (int i = 0; i < L1Constants::size::kMaxNFieldApproxCoefficients; ++i) {
-    cx[i] = L1Utils::kNaN;
-    cy[i] = L1Utils::kNaN;
-    cz[i] = L1Utils::kNaN;
+    cx[i] = undef::kFvc;
+    cy[i] = undef::kFvc;
+    cz[i] = undef::kFvc;
   }
+  z = undef::kFvc;
 }
 
 //----------------------------------------------------------------------------------------------------------------------
@@ -76,6 +78,7 @@ void L1FieldSlice::CheckConsistency() const
     L1Utils::CheckSimdVectorEquality(cy[i], "L1FieldSlice: cy");
     L1Utils::CheckSimdVectorEquality(cz[i], "L1FieldSlice: cz");
   }
+  L1Utils::CheckSimdVectorEquality(z, "L1FieldSlice: z");
 }
 
 //----------------------------------------------------------------------------------------------------------------------
@@ -117,6 +120,12 @@ void L1FieldSlice::GetFieldValue(const fvec& x, const fvec& y, L1FieldValue& B)
         + cz[16] * x4y + cz[17] * x3y2 + cz[18] * x2y3 + cz[19] * xy4 + cz[20] * y5;
 }
 
+void L1FieldSlice::GetFieldValueForLine(const L1TrackPar& t, L1FieldValue& B) const
+{
+  fvec dz = z - t.z;
+  GetFieldValue(t.x + t.tx * dz, t.y + t.ty * dz, B);
+}
+
 //----------------------------------------------------------------------------------------------------------------------
 //
 std::string L1FieldSlice::ToString(int indentLevel) const

reco/L1/L1Algo/L1Field.h

@@ -10,6 +10,9 @@
 #include "L1Constants.h"
 #include "L1Def.h"
 #include "L1SimdSerializer.h"
+#include "L1Undef.h"
+
+class L1TrackPar;
 
 class L1FieldValue {
 public:
@@ -76,6 +79,11 @@ public:
   /// \param B  the L1FieldValue output
   void GetFieldValue(const fvec& x, const fvec& y, L1FieldValue& B) const;
 
+  /// Gets field value for the intersection with a straight track
+  /// \param t  straight track
+  /// \param B  the L1FieldValue output
+  void GetFieldValueForLine(const L1TrackPar& t, L1FieldValue& B) const;
+
   /// String representation of class contents
   /// \param indentLevel      number of indent characters in the output
   std::string ToString(int indentLevel = 0) const;
@@ -91,6 +99,8 @@ public:
   fvec cz
     [L1Constants::size::kMaxNFieldApproxCoefficients];  ///< Polynomial coefficients for z-component of the field value
 
+  fvec z {undef::kFvc};  ///< z coordinate of the slice
+
   /// Serialization function
   friend class boost::serialization::access;
   template<class Archive>
@@ -99,6 +109,7 @@ public:
     ar& cx;
     ar& cy;
     ar& cz;
+    ar& z;
   }
 } _fvecalignment;
 

reco/L1/L1Algo/L1Portion.h

-/* Copyright (C) 2010-2021 Frankfurt Institute for Advanced Studies, Goethe-Universität Frankfurt, Frankfurt
-   SPDX-License-Identifier: GPL-3.0-only
-   Authors: Igor Kulakov [committer] */
-
-#ifndef L1Portion_H
-#define L1Portion_H
-
-#include <vector>
-
-#include "L1TrackPar.h"
-
-using std::vector;
-// class for portions organization
-
-template<typename T>
-class L1Portion;
-
-template<>
-class L1Portion<L1TrackPar> {
-public:
-  typedef L1TrackPar T;
-  typedef L1Vector<T> vType;
-
-  L1Portion() {};
-  L1Portion(int size) { reserve(size); };
-  L1Portion(int size, int size2) : dataSize(size2)
-  {
-    reserve(size);
-    //     reserve2(size2);
-  };
-  vType& operator[](int i) { return a[i]; }
-  void resize(int size) { a.resize(size); };
-  void reserve(int size) { a.reserve(size); };
-  void reserve2(int size)
-  {
-    for (unsigned int i = 0; i < a.size(); i++)
-      a[i].reserve(size);
-  };
-  void push_back(vType& v) { a.push_back(v); };
-  void add_void()
-  {
-    vType v("L1Portion<L1TrackPar>::add_void");
-    //     v.resize(dataSize);
-    a.push_back(v);
-    a[a.size() - 1].reserve(dataSize);
-  };
-  void add_void(int i)
-  {
-    T t;
-    a[i].push_back(t);
-  };
-
-  int CalcSize()
-  {
-    int size = 0;
-    for (unsigned int i = 0; i < a.size(); i++)
-      size += a[i].size();
-    return size * sizeof(T);
-  };
-
-private:
-  vector<vType> a {"L1Portion<L1TrackPar>::a"};
-  //   int mainSize; // size of a
-  int dataSize {0};  // size of vType
-};
-
-template<>
-class L1Portion<L1FieldRegion> {
-public:
-  typedef L1FieldRegion T;
-  typedef L1Vector<T> vType;
-
-  L1Portion() {};
-  L1Portion(int size) { reserve(size); };
-  L1Portion(int size, int size2) : dataSize(size2)
-  {
-    reserve(size);
-    //     reserve2(size2);
-  };
-
-  vType& operator[](int i) { return a[i]; }
-  void resize(int size) { a.resize(size); };
-  void reserve(int size) { a.reserve(size); };
-  void reserve2(int size)
-  {
-    for (unsigned int i = 0; i < a.size(); i++)
-      a[i].reserve(size);
-  };
-  void push_back(vType& v) { a.push_back(v); };
-  void add_void()
-  {
-    vType v("L1Portion<L1FieldRegion>::add_void");
-    //     v.resize(dataSize);
-    a.push_back(v);
-    a[a.size() - 1].reserve(dataSize);
-  };
-  void add_void(int i)
-  {
-    T t;
-    a[i].push_back(t);
-  };
-
-  int CalcSize()
-  {
-    int size = 0;
-    for (unsigned int i = 0; i < a.size(); i++)
-      size += a[i].size();
-    return size * sizeof(T);
-  };
-
-private:
-  vector<vType> a {"L1Portion<L1FieldRegion>::a"};
-  //   int mainSize; // size of a
-  int dataSize {0};  // size of vType
-};
-
-
-template<typename T>
-class L1Portion {
-public:
-  typedef vector<T> vType;
-
-  L1Portion() {};
-  L1Portion(int size) { reserve(size); };
-  L1Portion(int size, int size2) : dataSize(size2)
-  {
-    reserve(size);
-    //     reserve2(size2);
-  };
-  vType& operator[](int i) { return a[i]; }
-  void resize(int size) { a.resize(size); };
-  void reserve(int size) { a.reserve(size); };
-  void reserve2(int size)
-  {
-    for (int i = 0; i < a.size(); i++)
-      a[i].reserve(size);
-  };
-  void push_back(vType& v) { a.push_back(v); };
-  void add_void()
-  {
-    vType v("L1Portion<T>::add_void");
-    //     v.resize(dataSize);
-    a.push_back(v);
-    a[a.size() - 1].reserve(dataSize);
-  };
-  void add_void(int i)
-  {
-    T t;
-    a[i].push_back(t);
-  };
-
-  int CalcSize()
-  {
-    int size = 0;
-    for (int i = 0; i < a.size(); i++)
-      size += a[i].size();
-    return size * sizeof(T);
-  };
-
-private:
-  vector<vType> a {"L1Portion<T>::a"};
-  //   int mainSize; // size of a
-  int dataSize {0};  // size of vType
-};
-
-#endif  // L1Portion_H

reco/L1/L1Algo/L1Triplet.h

@@ -22,7 +22,7 @@ public:
   /// constructor
   L1Triplet(unsigned int iHitL, unsigned int iHitM, unsigned int iHitR, unsigned int iStaL, unsigned int iStaM,
             unsigned int iStaR, unsigned char Level, unsigned int firstNeighbour, char nNeighbours, fscal Chi2,
-            fscal Qp, fscal Cqp, fscal tx, fscal Ctx, fscal ty, fscal Cty)
+            fscal Qp, fscal Cqp, fscal tx, fscal Ctx, fscal ty, fscal Cty, bool isMomentumFitted)
     : fChi2(Chi2)
     , fQp(Qp)
     , fCqp(Cqp)
@@ -37,6 +37,7 @@ public:
     , fNneighbours(nNeighbours)
     , fLevel(Level)
     , fSta((iStaL << 4) + ((iStaM - iStaL - 1) << 2) + (iStaR - iStaL - 2))
+    , fIsMomentumFitted(isMomentumFitted)
   {
   }
 

reco/L1/L1Algo/L1TripletConstructor.h

+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov [committer] */
+
+#ifndef L1TripletConstructor_h
+#define L1TripletConstructor_h
+
+#include "L1Algo.h"
+#include "L1Field.h"
+#include "L1Fit.h"
+#include "L1HitPoint.h"
+#include "L1Station.h"
+#include "L1TrackPar.h"
+#include "L1Triplet.h"
+#include "L1Vector.h"
+
+
+/// Construction of triplets for the CA tracker
+///
+class L1TripletConstructor {
+public:
+  ///                             ------  Constructors and destructor ------
+
+  /// Constructor
+  /// \param nThreads  Number of threads for multi-threaded mode
+  L1TripletConstructor(L1Algo* algo);
+
+  /// Copy constructor
+  L1TripletConstructor(const L1TripletConstructor&) = delete;
+
+  /// Move constructor
+  L1TripletConstructor(L1TripletConstructor&&) = delete;
+
+  /// Copy assignment operator
+  L1TripletConstructor& operator=(const L1TripletConstructor&) = delete;
+
+  /// Move assignment operator
+  L1TripletConstructor& operator=(L1TripletConstructor&&) = delete;
+
+  /// Destructor
+  ~L1TripletConstructor() = default;
+
+
+  ///                             ------  FUNCTIONAL PART ------
+
+  void InitStations(int istal, int istam, int istar);
+
+  void CreateTripletsForHit(int istal, int istam, int istar, L1HitIndex_t ihl);
+
+  const L1Vector<L1Triplet>& GetTriplets() const { return fTriplets; }
+
+  /// Find the doublets. Reformat data in the portion of doublets.
+  void FindDoublets();
+  void FitDoublets();
+
+  /// Find triplets on station
+  void FindTriplets();
+
+  /// Add the middle hits to parameters estimation. Propagate to right station.
+  /// Find the triplets (right hit). Reformat data in the portion of triplets.
+  void findRightHit();
+
+  /// Fit Triplets
+  void fitTriplets();
+
+  /// Select triplets. Save them into vTriplets.
+  void storeTriplets();
+
+  void CollectHits(const L1TrackPar& Tr, const int iSta, const double chi2Cut, const int iMC,
+                   L1Vector<L1HitIndex_t>& collectedHits, int maxNhits);
+
+private:
+  /// left station
+  const L1Station& staL() { return *fStaL; }
+  /// middle station
+  const L1Station& staM() { return *fStaM; }
+  /// right station
+  const L1Station& staR() { return *fStaR; }
+
+private:
+  bool fIsInitialized {false};  ///< is initialized;
+  L1Algo* fAlgo {nullptr};      ///< pointer to L1Algo object
+
+  int fIstaL {-1};  ///< left station index
+  int fIstaM {-1};  ///< middle station index
+  int fIstaR {-1};  ///< right station index
+
+  const L1Station* fStaL {nullptr};  ///< left station
+  const L1Station* fStaM {nullptr};  ///< mid station
+  const L1Station* fStaR {nullptr};  ///< right station
+
+  L1HitPoint* fHitsL {nullptr};  ///< hits on the left station
+  L1HitPoint* fHitsM {nullptr};  ///< hits on the middle station
+  L1HitPoint* fHitsR {nullptr};  ///< hits on the right station
+
+  const L1Station* fFld0Sta[2];  // two stations for approximating the field between the target and the left hit
+  const L1Station* fFld1Sta[3];  // three stations for approximating the field between the left and the right hit
+
+  L1HitIndex_t fIhitL;
+  L1TrackPar fTrL;
+  L1FieldRegion fFldL;
+
+  L1Vector<L1HitIndex_t> fHitsM_2 {"L1TripletConstructor::fHitsM_2"};
+  L1Vector<L1TrackPar> fTracks_2 {"L1TripletConstructor::fTracks_2"};
+
+  L1Vector<L1TrackPar> fTracks_3 {"L1TripletConstructor::fTracks_3"};
+  L1Vector<L1HitIndex_t> fHitsM_3 {"L1TripletConstructor::fHitsM_3"};
+  L1Vector<L1HitIndex_t> fHitsR_3 {"L1TripletConstructor::fHitsR_3"};
+
+  L1Vector<L1Triplet> fTriplets {"L1TripletConstructor::fTriplets"};
+
+  L1Fit fFit;
+
+} _fvecalignment;
+
+#endif