From 2b80e8e070e8a84dd69c1b67ee3d9b1e0450874e Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Tue, 21 Jun 2022 14:28:03 +0200
Subject: [PATCH] L1: introduced electron flag to ca iterations, added
 L1Constants::phys namespace

---
 reco/L1/CbmL1.cxx                  |  16 ++-
 reco/L1/L1Algo/L1Algo.h            |   8 +-
 reco/L1/L1Algo/L1CAIteration.cxx   |  28 ++--
 reco/L1/L1Algo/L1CAIteration.h     |  92 ++++++++++----
 reco/L1/L1Algo/L1CATrackFinder.cxx | 197 +++++++++++------------------
 reco/L1/L1Algo/L1Constants.h       |   8 ++
 reco/L1/L1Algo/L1TrackExtender.cxx |   2 +-
 7 files changed, 177 insertions(+), 174 deletions(-)

diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 530d781b0a..d4b8fccad2 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -706,7 +706,7 @@ InitStatus CbmL1::Init()
   trackingIterFastPrim.SetMaxDZ(0);
   trackingIterFastPrim.SetTargetPosSigmaXY(1, 1);
   trackingIterFastPrim.SetMinLevelTripletStart(0);
-  trackingIterFastPrim.SetPrimary(true);
+  trackingIterFastPrim.SetPrimaryFlag(true);
 
   auto trackingIterAllPrim = L1CAIteration("AllPrimIter");
   trackingIterAllPrim.SetTrackChi2Cut(10.f);
@@ -720,7 +720,7 @@ InitStatus CbmL1::Init()
   trackingIterAllPrim.SetMaxDZ(0.1);
   trackingIterAllPrim.SetTargetPosSigmaXY(1, 1);
   trackingIterAllPrim.SetMinLevelTripletStart(0);
-  trackingIterAllPrim.SetPrimary(true);
+  trackingIterAllPrim.SetPrimaryFlag(true);
 
   auto trackingIterFastPrim2 = L1CAIteration("FastPrim2Iter");
   trackingIterFastPrim2.SetTrackChi2Cut(10.f);
@@ -734,7 +734,7 @@ InitStatus CbmL1::Init()
   trackingIterFastPrim2.SetMaxDZ(0);
   trackingIterFastPrim2.SetTargetPosSigmaXY(5, 5);
   trackingIterFastPrim2.SetMinLevelTripletStart(0);
-  trackingIterFastPrim2.SetPrimary(true);
+  trackingIterFastPrim2.SetPrimaryFlag(true);
 
   auto trackingIterAllSec = L1CAIteration("AllSecIter");
   trackingIterAllSec.SetTrackChi2Cut(10.f);
@@ -748,7 +748,7 @@ InitStatus CbmL1::Init()
   trackingIterAllSec.SetMaxDZ(0.1);
   trackingIterAllSec.SetTargetPosSigmaXY(10, 10);
   trackingIterAllSec.SetMinLevelTripletStart(1);
-  trackingIterAllSec.SetPrimary(false);
+  trackingIterAllSec.SetPrimaryFlag(false);
 
   auto trackingIterFastPrimJump = L1CAIteration("FastPrimJumpIter");
   trackingIterFastPrimJump.SetTrackChi2Cut(10.f);
@@ -762,7 +762,7 @@ InitStatus CbmL1::Init()
   trackingIterFastPrimJump.SetMaxDZ(0);
   trackingIterFastPrimJump.SetTargetPosSigmaXY(5, 5);
   trackingIterFastPrimJump.SetMinLevelTripletStart(0);
-  trackingIterFastPrimJump.SetPrimary(true);
+  trackingIterFastPrimJump.SetPrimaryFlag(true);
 
   auto trackingIterAllPrimJump = L1CAIteration("AllPrimJumpIter");
   trackingIterAllPrimJump.SetTrackChi2Cut(10.f);
@@ -776,7 +776,7 @@ InitStatus CbmL1::Init()
   trackingIterAllPrimJump.SetMaxDZ(0.1);
   trackingIterAllPrimJump.SetTargetPosSigmaXY(5, 5);
   trackingIterAllPrimJump.SetMinLevelTripletStart(0);
-  trackingIterAllPrimJump.SetPrimary(true);
+  trackingIterAllPrimJump.SetPrimaryFlag(true);
 
   auto trackingIterAllSecJump = L1CAIteration("AllSecJumpIter");
   trackingIterAllSecJump.SetTrackChi2Cut(10.f);
@@ -803,7 +803,8 @@ InitStatus CbmL1::Init()
   trackingIterAllPrimE.SetMaxDZ(0.1);
   trackingIterAllPrimE.SetTargetPosSigmaXY(1, 1);
   trackingIterAllPrimE.SetMinLevelTripletStart(0);
-  trackingIterAllPrimE.SetPrimary(true);
+  trackingIterAllPrimE.SetPrimaryFlag(true);
+  trackingIterAllPrimE.SetElectronFlag(true);
 
   auto trackingIterAllSecE = L1CAIteration("AllSecEIter");
   trackingIterAllSecE.SetTrackChi2Cut(10.f);
@@ -817,6 +818,7 @@ InitStatus CbmL1::Init()
   trackingIterAllSecE.SetMaxDZ(0.1);
   trackingIterAllSecE.SetMinLevelTripletStart(1);
   trackingIterAllSecE.SetTargetPosSigmaXY(10, 10);
+  trackingIterAllSecE.SetElectronFlag(true);
 
   // Select default track finder
   if (L1Algo::TrackingMode::kMcbm == fTrackingMode) {
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index ee87dc9cff..0084cc1387 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -105,8 +105,6 @@ public:
   /// \return particle mass squared
   float GetDefaultParticleMass2() const { return fDefaultMass * fDefaultMass; }
 
-  float fDefaultMass = 0.10565800;  // muon mass
-  // TODO: make fDefaultMass a private member (S.Zharko)
 
   /// pack station, thread and triplet indices to an unique triplet ID
   static unsigned int PackTripletId(unsigned int Station, unsigned int Thread, unsigned int Triplet)
@@ -228,6 +226,7 @@ private:
   int fNfieldStations {0};       ///< number of stations in the field region
   //alignas(16) L1StationsArray_t fStations {};  ///< array of L1Station objects
   //alignas(16) L1MaterialArray_t fRadThick {};  ///< material for each station
+  float fDefaultMass {L1Constants::phys::kMuonMass};  ///< mass of the propagated particle [GeV/c2]
 
 public:
   /// Gets total number of stations used in tracking
@@ -313,8 +312,9 @@ public:
 
 
   /// --- data used during finding iterations
+  int isec {0};  // iteration TODO: to be dispatched (S.Zharko, 21.06.2022)
+  const L1CAIteration* fpCurrentIteration {nullptr};  ///< pointer to the current CA track finder iteration
 
-  int isec {0};  // iteration
   L1Vector<L1Hit>* vHitsUnused {nullptr};
   L1Vector<L1HitIndex_t>* RealIHitP {nullptr};
   L1Vector<L1HitIndex_t>* RealIHitPBuf {nullptr};
@@ -625,7 +625,7 @@ private:
 
   /// ----- Different parameters of CATrackFinder -----
 
-  Tindex FIRSTCASTATION {0};  //first station used in CA
+  Tindex fFirstCAstation {0};  //first station used in CA
 
   // fNFindIterations - set number of interation for trackfinding
   // itetation of finding:
diff --git a/reco/L1/L1Algo/L1CAIteration.cxx b/reco/L1/L1Algo/L1CAIteration.cxx
index da51146069..022507e46d 100644
--- a/reco/L1/L1Algo/L1CAIteration.cxx
+++ b/reco/L1/L1Algo/L1CAIteration.cxx
@@ -1,4 +1,4 @@
-/* Copyright (C) 2016-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
    SPDX-License-Identifier: GPL-3.0-only
    Authors: Sergey Gorbunov, Sergei Zharko [committer] */
 
@@ -24,7 +24,6 @@ L1CAIteration::L1CAIteration() noexcept { LOG(debug) << "L1CAIteration: Default
 L1CAIteration::L1CAIteration(const L1CAIteration& other) noexcept
   // Basic fields
   : fName(other.fName)
-  , fControlFlags(other.fControlFlags)
   // Cuts
   , fTrackChi2Cut(other.fTrackChi2Cut)
   , fTripletChi2Cut(other.fTripletChi2Cut)
@@ -38,8 +37,10 @@ L1CAIteration::L1CAIteration(const L1CAIteration& other) noexcept
   , fTargetPosSigmaX(other.fTargetPosSigmaX)
   , fTargetPosSigmaY(other.fTargetPosSigmaY)
   , fMinLevelTripletStart(other.fMinLevelTripletStart)
+  , fFirstStationIndex(other.fFirstStationIndex)
+  , fIsPrimary(other.fIsPrimary)
+  , fIsElectron(other.fIsElectron)
 {
-  LOG(debug) << "L1CAIteration: Copy constructor called: " << &other << " was copied into " << this;
 }
 //
 //----------------------------------------------------------------------------------------------------------------------
@@ -47,26 +48,21 @@ L1CAIteration::L1CAIteration(const L1CAIteration& other) noexcept
 L1CAIteration::L1CAIteration(L1CAIteration&& other) noexcept
 {
   this->Swap(other);
-  LOG(debug) << "L1CAIteration: Move constructor called: " << &other << " was moved into " << this;
 }
 //
 //----------------------------------------------------------------------------------------------------------------------
 //
-L1CAIteration::L1CAIteration(const std::string& name) noexcept : fName(name)
-{
-  LOG(debug) << "L1CAIteration: Constructor from name called for " << this;
-}
+L1CAIteration::L1CAIteration(const std::string& name) noexcept : fName(name) {}
 //
 //----------------------------------------------------------------------------------------------------------------------
 //
-L1CAIteration::~L1CAIteration() noexcept { LOG(debug) << "L1CAIteration: Destructor called for " << this; }
+L1CAIteration::~L1CAIteration() noexcept {}
 //
 //----------------------------------------------------------------------------------------------------------------------
 //
 L1CAIteration& L1CAIteration::operator=(const L1CAIteration& other) noexcept
 {
   if (this != &other) { L1CAIteration(other).Swap(*this); }
-  LOG(debug) << "L1CAIteration: Copy operator= called: " << &other << " was copied into " << this;
   return *this;
 }
 //
@@ -78,7 +74,6 @@ L1CAIteration& L1CAIteration::operator=(L1CAIteration&& other) noexcept
     L1CAIteration tmp(std::move(other));
     this->Swap(tmp);
   }
-  LOG(debug) << "L1CAIteration: Move operator= called: " << &other << " was moved into " << this;
   return *this;
 }
 //
@@ -104,7 +99,6 @@ void L1CAIteration::Swap(L1CAIteration& other) noexcept
 {
   // Basic fields
   std::swap(fName, other.fName);
-  std::swap(fControlFlags, other.fControlFlags);
   // Cuts
   std::swap(fTrackChi2Cut, other.fTrackChi2Cut);
   std::swap(fTripletChi2Cut, other.fTripletChi2Cut);
@@ -118,6 +112,9 @@ void L1CAIteration::Swap(L1CAIteration& other) noexcept
   std::swap(fTargetPosSigmaX, other.fTargetPosSigmaX);
   std::swap(fTargetPosSigmaY, other.fTargetPosSigmaY);
   std::swap(fMinLevelTripletStart, other.fMinLevelTripletStart);
+  std::swap(fFirstStationIndex, other.fFirstStationIndex);
+  std::swap(fIsPrimary, other.fIsPrimary);
+  std::swap(fIsElectron, other.fIsElectron);
 }
 //
 //----------------------------------------------------------------------------------------------------------------------
@@ -128,8 +125,8 @@ std::string L1CAIteration::ToString(int indentLevel) const
   constexpr char indentChar = '\t';
   std::string indent(indentLevel, indentChar);
   aStream << indent << "L1CAIteration: " << fName << '\n';
-  aStream << indent << indentChar
-          << "Is primary:                   " << fControlFlags[static_cast<int>(ControlFlag::kePrimary)] << '\n';
+  aStream << indent << indentChar << "Is primary:   " << fIsPrimary << '\n';
+  aStream << indent << indentChar << "Is electron:  " << fIsElectron << '\n';
   aStream << indent << indentChar << "Track chi2 cut:               " << fTrackChi2Cut << '\n';
   aStream << indent << indentChar << "Triplet chi2 cut:             " << fTripletChi2Cut << '\n';
   aStream << indent << indentChar << "Doublet chi2 cut:             " << fDoubletChi2Cut << '\n';
@@ -141,7 +138,8 @@ std::string L1CAIteration::ToString(int indentLevel) const
   aStream << indent << indentChar << "Max DZ:                       " << fMaxDZ << '\n';
   aStream << indent << indentChar << "Target position sigma X [cm]: " << fTargetPosSigmaX << '\n';
   aStream << indent << indentChar << "Target position sigma Y [cm]: " << fTargetPosSigmaY << '\n';
-  aStream << indent << indentChar << "Min level for triplet start:  " << fMinLevelTripletStart;
+  aStream << indent << indentChar << "Min level for triplet start:  " << fMinLevelTripletStart << '\n';
+  aStream << indent << indentChar << "First tracking station index: " << fFirstStationIndex;
 
   return aStream.str();
 }
diff --git a/reco/L1/L1Algo/L1CAIteration.h b/reco/L1/L1Algo/L1CAIteration.h
index 538a916690..455d124eba 100644
--- a/reco/L1/L1Algo/L1CAIteration.h
+++ b/reco/L1/L1Algo/L1CAIteration.h
@@ -1,4 +1,4 @@
-/* Copyright (C) 2016-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
    SPDX-License-Identifier: GPL-3.0-only
    Authors: Sergey Gorbunov, Sergei Zharko [committer] */
 
@@ -23,96 +23,129 @@
 /// reconstructed during current iteration are removed from the further iterations.
 ///
 class L1CAIteration {
-  /// Enumeration ControlFlag is used to keep flags, which controls behaviour of the track finder
-  /// iterations loop
-  ///
-  enum class ControlFlag
-  {
-    kePrimary,  ///> true - track is primary, false - track is secondary (not primary)
-    keEnd
-  };
-  using ControlFlags_t = std::bitset<static_cast<int>(ControlFlag::keEnd)>;
-
 public:
   /// Default constructor
   L1CAIteration() noexcept;
+  
   /// Copy constructor
   L1CAIteration(const L1CAIteration& other) noexcept;
+  
   /// Move constructor
   L1CAIteration(L1CAIteration&& other) noexcept;
+  
   /// Constructor from L1CAIteration type
   L1CAIteration(const std::string& name) noexcept;
+  
   /// Destructor
   ~L1CAIteration() noexcept;
+  
   /// Copy assignment operator
   L1CAIteration& operator=(const L1CAIteration& other) noexcept;
+  
   /// Move assignment operator
   L1CAIteration& operator=(L1CAIteration&& other) noexcept;
 
   /// Gets doublet chi2 upper cut
   float GetDoubletChi2Cut() const { return fDoubletChi2Cut; }
+  
+  /// Gets station index of the first station used in tracking
+  int GetFirstStationIndex() const { return fFirstStationIndex; }
+
   /// Gets correction for accounting overlaping and iff z
   float GetMaxDZ() const { return fMaxDZ; }
+  
   /// Gets max considered q/p for tracks
   float GetMaxInvMom() const { return fMaxInvMom; }
+  
   /// Gets max slope (tx\ty) in 3D hit position of a triplet
   float GetMaxSlope() const { return fMaxSlope; }
+  
   /// Gets max slope (tx\ty) in primary vertex
   float GetMaxSlopePV() const { return fMaxSlopePV; }
+  
   /// Gets min level of the triplet start
   int GetMinLevelTripletStart() const { return fMinLevelTripletStart; }
+  
   /// Gets the name of the iteration
   const std::string& GetName() const { return fName; }
+  
   /// Gets size of region [TODO: units??] to attach new hits to the created track
   float GetPickGather() const { return fPickGather; }
+  
   /// Gets min value of dp/dp_error, for which two tiplets are neighbours
   float GetPickNeighbour() const { return fPickNeighbour; }
+  
   /// Gets sigma target position in X direction [cm]
   float GetTargetPosSigmaX() const { return fTargetPosSigmaX; }
+  
   /// Gets sigma target position in Y direction [cm]
   float GetTargetPosSigmaY() const { return fTargetPosSigmaY; }
+  
   /// Gets track chi2 upper cut
   float GetTrackChi2Cut() const { return fTrackChi2Cut; }
+  
   /// Gets triplet chi2 upper cut
   float GetTripletChi2Cut() const { return fTripletChi2Cut; }
 
-  /// flag check: primary tracks - true, secondary tracks - false
-  bool IsPrimary() const { return fControlFlags[static_cast<int>(ControlFlag::kePrimary)]; }
+  /// flag check: electrons/positrons - true, heavy charged - false
+  bool GetElectronFlag() const { return fIsElectron; }
+
+  /// Checks flag: true - only primary tracks are searched, false - [all or only secondary?] 
+  bool GetPrimaryFlag() const { return fIsPrimary; }
 
   /// Prints iteration options
   void Print(int verbosityLevel = 0) const;
 
   /// Sets doublet chi2 upper cut
   void SetDoubletChi2Cut(float input) { fDoubletChi2Cut = input; }
+  
+  /// Sets flag: electron tracks - true, heavy ion tracks - false
+  void SetElectronFlag(bool flag) { fIsElectron = flag; }
+
+  /// Sets index of first station used in tracking
+  void SetFirstStationIndex(int index) { fFirstStationIndex = index; } 
+
   /// Sets correction for accounting overlaping and iff z
   void SetMaxDZ(float input) { fMaxDZ = input; }
+  
   /// Sets max considered q/p for tracks
   void SetMaxInvMom(float input) { fMaxInvMom = input; }
+  
   /// Sets max slope (tx\ty) in 3D hit position of a triplet
   void SetMaxSlope(float input) { fMaxSlope = input; }
+  
   /// Sets max slope (tx\ty) in primary vertex
   void SetMaxSlopePV(float input) { fMaxSlopePV = input; }
+  
   /// Sets min level of the triplet start
   void SetMinLevelTripletStart(int input) { fMinLevelTripletStart = input; }
+  
   /// Sets name of the iteration
   void SetName(const std::string& name) { fName = name; }
+  
   /// Sets size of region [TODO: units??] to attach new hits to the created track
   void SetPickGather(float input) { fPickGather = input; }
+  
   /// Sets min value of dp/dp_error, for which two tiplets are neighbours
   void SetPickNeighbour(float input) { fPickNeighbour = input; }
+  
   /// Sets flag: primary tracks - true, secondary tracks - false
-  void SetPrimary(bool flag) { fControlFlags[static_cast<int>(ControlFlag::kePrimary)] = flag; }
+  void SetPrimaryFlag(bool flag) { fIsPrimary = flag; }
+  
   /// Sets sigma of target positions in XY plane
   /// \param  sigmaX  Sigma value in X direction [cm]
   /// \param  sigmaX  Sigma value in Y direction [cm]
   void SetTargetPosSigmaXY(float sigmaX, float sigmaY);
+  
   /// Sets track chi2 upper cut
   void SetTrackChi2Cut(float input) { fTrackChi2Cut = input; }
+  
   /// Sets triplet chi2 upper cut
   void SetTripletChi2Cut(float input) { fTripletChi2Cut = input; }
 
   /// Swap method
   void Swap(L1CAIteration& other) noexcept;
+  
   /// String representation of the class contents
   /// \param indentLevel  Level of indentation for the text (in terms of \t symbols)
   std::string ToString(int indentLevel = 0) const;
@@ -120,25 +153,30 @@ public:
 private:
   /** Basic fields **/
   std::string fName {""};           ///< Iteration name
-  ControlFlags_t fControlFlags {};  ///< bitset flags to control iteration behaviour
 
   /** Track finder dependent cuts **/
   // TODO: Iteratively change the literals to floats (S.Zharko)
   // NOTE: For each new cut one should not forget to create a setter and a getter, insert the value
   //       initialization in the copy constructor and the Swap operator as well as a string repre-
   //       sentation to the ToString method (S.Zharko)
-  float fTrackChi2Cut {10.f};                   ///> Track chi2 upper cut
-  float fTripletChi2Cut {21.1075f};             ///> Triplet chi2 upper cut
-  float fDoubletChi2Cut {11.3449 * 2.f / 3.f};  ///> Doublet chi2 upper cut
-  float fPickGather {3.0};       ///> Size of region to attach new hits to the created track [TODO: units??]
-  float fPickNeighbour {5.0};    ///> Min value of dp/dp_error, for which two tiplets are neighbours
-  float fMaxInvMom {1.0 / 0.5};  ///> Max considered q/p for tracks
-  float fMaxSlopePV {1.1};       ///> Max slope (tx\ty) in primary vertex
-  float fMaxSlope {2.748};       ///> Max slope (tx\ty) in 3D hit position of a triplet
-  float fMaxDZ {0.f};            ///> Correction for accounting overlaping and iff z [TODO: units??]
-  float fTargetPosSigmaX {0};                   ///> Constraint on target position in X direction [cm]
-  float fTargetPosSigmaY {0};                   ///> Constraint on target position in Y direction [cm]
-  int fMinLevelTripletStart {0};  ///> Min level for starting a triplet. Track length = fMinLevelTripletStart + 3
+  float fTrackChi2Cut {10.f};                   ///< Track chi2 upper cut
+  float fTripletChi2Cut {21.1075f};             ///< Triplet chi2 upper cut
+  float fDoubletChi2Cut {11.3449 * 2.f / 3.f};  ///< Doublet chi2 upper cut
+  float fPickGather {3.0};                      ///< Size of region to attach new hits to the created track [TODO: units??]
+  float fPickNeighbour {5.0};                   ///< Min value of dp/dp_error, for which two tiplets are neighbours
+  float fMaxInvMom {1.0 / 0.5};                 ///< Max considered q/p for tracks
+  float fMaxSlopePV {1.1};                      ///< Max slope (tx\ty) in primary vertex
+  float fMaxSlope {2.748};                      ///< Max slope (tx\ty) in 3D hit position of a triplet
+  float fMaxDZ {0.f};                           ///< Correction for accounting overlaping and iff z [TODO: units??]
+  float fTargetPosSigmaX {0};                   ///< Constraint on target position in X direction [cm]
+  float fTargetPosSigmaY {0};                   ///< Constraint on target position in Y direction [cm]
+  int   fMinLevelTripletStart {0};              ///< Min level for starting a triplet. Track length = fMinLevelTripletStart + 3
+  int   fFirstStationIndex {0};                 ///< First station, used for tracking
+
+  bool  fIsPrimary  {false};                    ///< Flag: true - only primary tracks are searched for
+  bool  fIsElectron {false};                    ///< Flag: true - only electrons are searched for
+
+
   // ^ TODO: invent more proper name
 };
 
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 0799c88db2..5b461674dc 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -122,7 +122,6 @@ inline void L1Algo::findSingletsStep1(  /// input 1st stage of singlet search
   int istal,
   int istam,    /// indexes of left and middle stations of a triplet
   Tindex n1_V,  ///
-
   fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* HitTime_l, fvec* HitTimeEr,
   // output
   //                 L1TrackPar *T_1,
@@ -155,11 +154,11 @@ inline void L1Algo::findSingletsStep1(  /// input 1st stage of singlet search
 
   L1Fit fit;
 
-  if ((isec != kAllPrimEIter) && (isec != kAllSecEIter)) {
-    fit.SetParticleMass(fDefaultMass);  // muon
+  if (fpCurrentIteration->GetElectronFlag()) {
+    fit.SetParticleMass(L1Constants::phys::kElectronMass);
   }
   else {
-    fit.SetParticleMass(0.000511f);  // electron
+    fit.SetParticleMass(fDefaultMass);
   }
 
   for (int i1_V = 0; i1_V < n1_V; i1_V++) {
@@ -1199,7 +1198,7 @@ inline void L1Algo::f5(  // input
   if (isec != TRACKS_FROM_TRIPLETS_ITERATION)
 #endif
 
-    for (int istal = fNstations - 4; istal >= FIRSTCASTATION; istal--) {
+    for (int istal = fNstations - 4; istal >= fFirstCAstation; istal--) {
       for (int tripType = 0; tripType < 3; tripType++)  // tT = 0 - 123triplet, tT = 1 - 124triplet, tT = 2 - 134triplet
       {
         if ((((isec != kFastPrimJumpIter) && (isec != kAllPrimJumpIter) && (isec != kAllSecJumpIter))
@@ -1279,10 +1278,19 @@ inline void L1Algo::f5(  // input
 
 inline void L1Algo::DupletsStaPort(
   /// input:
-  int istal, int istam, Tindex ip, L1Vector<Tindex>& n_g, Tindex* portionStopIndex_,
+  int istal, 
+  int istam, 
+  Tindex ip, 
+  L1Vector<Tindex>& n_g, 
+  Tindex* portionStopIndex_,
   /// output:
-  L1TrackPar* T_1, L1FieldRegion* fld_1, L1HitIndex_t* hitsl_1, L1Vector<char>& lmDuplets, Tindex& n_2,
-  L1Vector<L1HitIndex_t>& i1_2, L1Vector<L1HitIndex_t>& hitsm_2
+  L1TrackPar* T_1, 
+  L1FieldRegion* fld_1, 
+  L1HitIndex_t* hitsl_1, 
+  L1Vector<char>& lmDuplets, 
+  Tindex& n_2,
+  L1Vector<L1HitIndex_t>& i1_2, 
+  L1Vector<L1HitIndex_t>& hitsm_2
   ///
 )
 {
@@ -1377,20 +1385,43 @@ inline void L1Algo::DupletsStaPort(
 /// ------------------- Triplets on station ----------------------
 
 inline void
-L1Algo::TripletsStaPort(  /// creates triplets: input: @istal - start station number, @istam - middle station number, @istar - last station number, @ip - index of portion, @&n_g - numer of elements in portion, @*portionStopIndex
-  int istal, int istam, int istar,
-
-  ///@nstaltriplets - , @*portionStopIndex, @*T_1 - track parameters for singlets, @*fld_1 - field approximation for singlets, @&n_2 - number of doublets in portion
-  ///  @&n_2 - number of douplets,@&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
-
-
-  Tindex& nstaltriplets, L1TrackPar* T_1, L1FieldRegion* fld_1, L1HitIndex_t* hitsl_1,
-
-  Tindex& n_2, L1Vector<L1HitIndex_t>& i1_2, L1Vector<L1HitIndex_t>& hitsm_2,
+L1Algo::TripletsStaPort(  /// creates triplets: 
+  /// input: 
+  ///  @istal - start station number, 
+  ///  @istam - middle station number, 
+  ///  @istar - last station number, 
+  ///  @ip - index of portion, 
+  ///  @&n_g - numer of elements in portion, 
+  ///  @*portionStopIndex
+  int istal, 
+  int istam, 
+  int istar,
+
+  /// @nstaltriplets - , 
+  /// @*portionStopIndex, 
+  /// @*T_1 - track parameters for singlets, 
+  /// @*fld_1 - field approximation for singlets, 
+  /// @&n_2 - number of doublets in portion
+  /// @&n_2 - number of douplets,
+  /// @&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
+
+
+  Tindex& nstaltriplets, 
+  L1TrackPar* T_1, 
+  L1FieldRegion* fld_1, 
+  L1HitIndex_t* hitsl_1,
+
+  Tindex& n_2, 
+  L1Vector<L1HitIndex_t>& i1_2, 
+  L1Vector<L1HitIndex_t>& hitsm_2,
 
   const L1Vector<char>& mrDuplets
 
-  /// output: @*vTriplets_part - array of triplets, @*TripStartIndexH, @*TripStopIndexH - start/stop index of a triplet in the array
+  /// output: 
+  // @*vTriplets_part - array of triplets, 
+  // @*TripStartIndexH, 
+  // @*TripStopIndexH - start/stop index of a triplet in the array
 
 )
 {
@@ -1734,12 +1765,10 @@ void L1Algo::CATrackFinder()
   // ---- Loop over Track Finder iterations ----------------------------------------------------------------//
   L1ASSERT(0, fNFindIterations == fParameters.GetCAIterations().size());
   isec = 0;                                                      // TODO: temporary! (S.Zharko)
+  
   for (const auto& caIteration : fParameters.GetCAIterations())  // all finder
   {
-    //if (fTrackingMode == kMcbm) {
-    //  if (isec > 3) { continue; }
-    //}
-    // n_g1.assign(n_g1.size(), Portion);
+    fpCurrentIteration = &caIteration;  // select current iteration
 
     for (int n = 0; n < fNThreads; n++) {
       for (int j = 0; j < fNstations; j++) {
@@ -1781,66 +1810,14 @@ void L1Algo::CATrackFinder()
       {
         // --- SET PARAMETERS FOR THE ITERATION ---
 
-        FIRSTCASTATION = 0;
-
-        // if ( (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter) )
-        //   FIRSTCASTATION = 2;
-
+        fFirstCAstation = caIteration.GetFirstStationIndex();
         fDoubletChi2Cut = caIteration.GetDoubletChi2Cut();  //11.3449 * 2.f / 3.f;  // prob = 0.1
-        //fDoubletChi2Cut = 11.3449 * 2.f / 3.f;  // prob = 0.1
-
         fTripletChi2Cut = caIteration.GetTripletChi2Cut();  //21.1075;  // prob = 0.01%
-        //fTripletChi2Cut = 21.1075;  // prob = 0.01%
-
-        //switch (isec) {
-        //  case kFastPrimIter:
-        //    fTripletChi2Cut = 7.815 * 3;  // prob = 0.05
-        //    break;
-        //  case kAllPrimIter:
-        //  case kAllPrimEIter:
-        //    fTripletChi2Cut = 7.815 * 3;  // prob = 0.05
-        //    break;
-        //  case kAllPrimJumpIter:
-        //    fTripletChi2Cut = 6.252 * 3;  // prob = 0.1
-        //    break;
-        //  case kAllSecIter:
-        //  case kAllSecEIter:
-        //    fTripletChi2Cut = 6.252 * 3;  //2.706; // prob = 0.1
-        //    break;
-        //}
-
-        /// coefficient for size of region for attach new hits to the created track
         fPickGather = caIteration.GetPickGather();  //3.0;
-        //fPickGather = 3.0;
-        //if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecIter)
-        //    || (isec == kAllSecEIter) || (isec == kAllSecJumpIter))
-        //  fPickGather = 4.0;
-
-        // (fPickNeighbour < dp/dp_error)  =>  triplets are neighbours
         fPickNeighbour = caIteration.GetPickNeighbour();  //5.0;
-        //fPickNeighbour = 5.0;
-        //if ( (isec == kFastPrimIter) )
-        //  fPickNeighbour = 5.0*0.5; // TODO understand why works with 0.2
-
         fMaxInvMom = caIteration.GetMaxInvMom();  //1.0 / 0.5;  // max considered q/p
-        //fMaxInvMom = 1.0 / 0.5;  // max considered q/p
-
-        //if (fTrackingMode == kMcbm) fMaxInvMom = 1 / 0.3;  // max considered q/p
-        //if ((isec == kAllPrimJumpIter) || (isec == kAllSecIter) || (isec == kAllSecJumpIter)) fMaxInvMom = 1.0 / 0.1;
-        //if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllSecEIter)) fMaxInvMom = 1. / 0.05;
-
-        //if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllSecEIter))
-        //  if (fTrackingMode == kMcbm) fMaxInvMom = 1 / 0.1;  // max considered q/p
-
-
         fMaxSlopePV = caIteration.GetMaxSlopePV();  //1.1;
-        //fMaxSlopePV = 1.1;
-        //if (  // (isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) ||
-        //  (isec == kAllSecIter) || (isec == kAllSecEIter) || (isec == kAllSecJumpIter))
-        //  fMaxSlopePV = 1.5;
-
         fMaxSlope = caIteration.GetMaxSlope();  //2.748;  // corresponds to 70 grad
-        //fMaxSlope = 2.748;  // corresponds to 70 grad
 
         // define the target
         fTargX = fParameters.GetTargetPositionX();
@@ -1851,28 +1828,11 @@ void L1Algo::CATrackFinder()
         float SigmaTargetY = caIteration.GetTargetPosSigmaY();  // target constraint [cm]
 
         // Select magnetic field. For primary tracks - fVtxFieldValue, for secondary tracks - st.fieldSlice
-        if (caIteration.IsPrimary()) { fTargB = fParameters.GetVertexFieldValue(); }
+        if (caIteration.GetPrimaryFlag()) { fTargB = fParameters.GetVertexFieldValue(); }
         else {
           fParameters.GetStation(0).fieldSlice.GetFieldValue(0, 0, fTargB);
         }  // NOTE: calculates field fTargB in the center of 0th station
 
-
-        //if ((isec == kFastPrimIter) || (isec == kFastPrimIter2) || (isec == kFastPrimJumpIter) || (isec == kAllPrimIter)
-        //    || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter)) {  // target
-        //  fTargB = fVtxFieldValue;
-        //  if ((isec == kFastPrimIter) || (isec == kAllPrimIter) || (isec == kAllPrimEIter))
-        //    SigmaTargetX = SigmaTargetY = 1;  // target
-        //  else
-        //    SigmaTargetX = SigmaTargetY = 5;
-        //}
-        //if ((isec == kAllSecIter) || (isec == kAllSecEIter)
-        //    || (isec == kAllSecJumpIter)) {  //use outer radius of the 1st station as a constraint // ?
-        //  L1Station& st = fParameters.GetStation(0);
-        //  SigmaTargetX = SigmaTargetY = 10;  //st.Rmax[0];
-        //  fTargZ                      = fParameters.GetTargetPositionZ();  // fParameters.GetTargetPositionZ()-1.;
-        //  st.fieldSlice.GetFieldValue(0, 0, fTargB);
-        //}
-
         TargetXYInfo.C00 = SigmaTargetX * SigmaTargetX;
         TargetXYInfo.C10 = 0;
         TargetXYInfo.C11 = SigmaTargetY * SigmaTargetY;
@@ -1881,14 +1841,6 @@ void L1Algo::CATrackFinder()
         /// The reason is that low momentum tracks are too curved and goes not from target direction. That's why sort by hit_y/hit_z is not work idealy
         /// If sort by y then it is max diff between same station's modules (~0.4cm)
         fMaxDZ = caIteration.GetMaxDZ();  //0;
-        //fMaxDZ = 0;
-        //if ((isec == kAllPrimIter) || (isec == kAllPrimEIter) || (isec == kAllPrimJumpIter) || (isec == kAllSecIter)
-        //    || (isec == kAllSecEIter) || (isec == kAllSecJumpIter))
-        //  fMaxDZ = 0.1;
-
-        // TODO: to be removed, because this condition is checked in L1InitManager (S.Zharko)
-        if (fNstations > (int) L1Constants::size::kMaxNstations)
-          cout << " CATrackFinder: Error: Too many Stations" << endl;
       }
 
 #ifndef L1_NO_ASSERT
@@ -1906,7 +1858,7 @@ void L1Algo::CATrackFinder()
       /// possible left hits of triplets are splited in portions of 16 (4 SIMDs) to use memory faster
       fDupletPortionStopIndex[fNstations - 1] = 0;
       fDupletPortionSize.clear();
-      for (int istal = fNstations - 2; istal >= FIRSTCASTATION; istal--) {  //start downstream chambers
+      for (int istal = fNstations - 2; istal >= fFirstCAstation; istal--) {  //start downstream chambers
         int NHits_l   = HitsUnusedStopIndex[istal] - HitsUnusedStartIndex[istal];
         int nPortions = NHits_l / Portion;
         int rest      = NHits_l - nPortions * Portion;
@@ -1928,7 +1880,7 @@ void L1Algo::CATrackFinder()
          fDupletPortionStopIndex[fNstations-1] = 0;
          unsigned int ip = 0;  //index of curent portion
          
-         for (int istal = fNstations-2; istal >= FIRSTCASTATION; istal--)  //start downstream chambers
+         for (int istal = fNstations-2; istal >= fFirstCAstation; istal--)  //start downstream chambers
          {
          int nHits = HitsUnusedStopIndex[istal] - HitsUnusedStartIndex[istal];
          
@@ -1964,20 +1916,23 @@ void L1Algo::CATrackFinder()
     L1FieldRegion fldG_1[vPortion];
     L1HitIndex_t hitslG_1[Portion];
 
-    L1Vector<L1HitIndex_t> hitsm_2(
-      "L1CATrackFinder::hitsm_2");  /// middle hits indexed by number of doublets in portion(i2)
-    L1Vector<L1HitIndex_t> i1_2(
-      "L1CATrackFinder::i1_2");  /// index in portion of singlets(i1) indexed by index in portion of doublets(i2)
-
-
-    L1Vector<L1HitIndex_t> hitsmG_2(
-      "L1CATrackFinder::hitsmG_2");  /// middle hits indexed by number of doublets in portion(i2)
-    L1Vector<L1HitIndex_t> i1G_2(
-      "L1CATrackFinder::i1G_2");  /// index in portion of singlets(i1) indexed by index in portion of doublets(i2)
-    L1Vector<char> lmDuplets[L1Constants::size::kMaxNstations] {
-      "L1CATrackFinder::lmDuplets"};  // is exist a doublet started from indexed by left hit
-    L1Vector<char> lmDupletsG[L1Constants::size::kMaxNstations] {
-      "L1CATrackFinder::lmDupletsG"};  // is exist a doublet started from indexed by left hit
+    /// middle hits indexed by number of doublets in portion(i2)
+    L1Vector<L1HitIndex_t> hitsm_2("L1CATrackFinder::hitsm_2");  
+    
+    /// index in portion of singlets(i1) indexed by index in portion of doublets(i2)
+    L1Vector<L1HitIndex_t> i1_2("L1CATrackFinder::i1_2");
+
+    /// middle hits indexed by number of doublets in portion(i2)
+    L1Vector<L1HitIndex_t> hitsmG_2("L1CATrackFinder::hitsmG_2");
+ 
+    /// index in portion of singlets(i1) indexed by index in portion of doublets(i2)
+    L1Vector<L1HitIndex_t> i1G_2("L1CATrackFinder::i1G_2");
+    
+    /// is exist a doublet started from indexed by left hit
+    L1Vector<char> lmDuplets[L1Constants::size::kMaxNstations] {"L1CATrackFinder::lmDuplets"};
+    
+    /// is exist a doublet started from indexed by left hit
+    L1Vector<char> lmDupletsG[L1Constants::size::kMaxNstations] {"L1CATrackFinder::lmDupletsG"};
 
     for (int i = 0; i < fNstations; i++) {
       lmDuplets[i].SetName(std::stringstream() << "L1CATrackFinder::lmDuplets[" << i << "]");
@@ -1989,7 +1944,7 @@ void L1Algo::CATrackFinder()
     hitsmG_2.reserve(9000);
     i1G_2.reserve(9000);
 
-    for (int istal = fNstations - 2; istal >= FIRSTCASTATION; istal--)  //  //start downstream chambers
+    for (int istal = fNstations - 2; istal >= fFirstCAstation; istal--)  //  //start downstream chambers
     {
 
 #ifdef _OPENMP
@@ -2131,7 +2086,7 @@ void L1Algo::CATrackFinder()
       fStripToTrack.reset(fNstrips, -1);
       fStripToTrackB.reset(fNstrips, -1);
 
-      for (int istaF = FIRSTCASTATION; istaF <= fNstations - 3 - ilev; ++istaF) {
+      for (int istaF = fFirstCAstation; 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,                  \
@@ -2164,10 +2119,12 @@ void L1Algo::CATrackFinder()
               if (isec != TRACKS_FROM_TRIPLETS_ITERATION)
 #endif
               {  // ghost supression !!!
+                // TODO: Primary => 3 hits tracks are saved, other whise 3 hit tracks are thrown away
                 if (isec != kFastPrimIter && isec != kAllPrimIter && isec != kAllPrimEIter && isec != kAllSecEIter)
                   if (first_trip.GetLevel() == 0)
                     continue;  // ghost suppression // find track with 3 hits only if it was created from a chain of triplets, but not from only one triplet
 
+
                 if (kGlobal != fTrackingMode && kMcbm != fTrackingMode) {
                   if ((ilev == 0) && (GetFStation((*fStripFlag)[(*vHitsUnused)[first_trip.GetLHit()].f]) != 0))
                     continue;  // ghost supression // collect only MAPS tracks-triplets  CHECK!!!
diff --git a/reco/L1/L1Algo/L1Constants.h b/reco/L1/L1Algo/L1Constants.h
index 3137ead7ca..fcf4d51a74 100644
--- a/reco/L1/L1Algo/L1Constants.h
+++ b/reco/L1/L1Algo/L1Constants.h
@@ -58,6 +58,14 @@ namespace L1Constants
     //constexpr bool kIfCreateTripletPulls {false};
   }  // end namespace control
 
+  /// Physics constants
+  namespace phys
+  {
+    /* Particle masses used for track fit */
+    constexpr float kMuonMass     {0.10565800f};  ///< Muon mass     [GeV/c2]
+    constexpr float kElectronMass {0.000511f};    ///< Electron mass [GeV/c2]
+  }
+
   /// Miscellaneous constants
   namespace misc
   {
diff --git a/reco/L1/L1Algo/L1TrackExtender.cxx b/reco/L1/L1Algo/L1TrackExtender.cxx
index c865df9494..ece5b290d7 100644
--- a/reco/L1/L1Algo/L1TrackExtender.cxx
+++ b/reco/L1/L1Algo/L1TrackExtender.cxx
@@ -246,7 +246,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
 
   int ista = ista2 + 2 * step;  // skip one station. if there would be hit it has to be found on previous stap
 
-  if (ista2 == FIRSTCASTATION) ista = ista2 + step;
+  if (ista2 == fFirstCAstation) ista = ista2 + step;
 
   const fvec pickGather2 = fPickGather * fPickGather;
 
-- 
GitLab