From b4e190e69282cdb8039b5af736725a4d7699a5f4 Mon Sep 17 00:00:00 2001
From: sergeizharko <zharkosergey94@yandex.ru>
Date: Mon, 20 Dec 2021 19:58:30 +0100
Subject: [PATCH] L1Parameters.h: new file ; L1BaseStationInfo.h: getters and
 setters added; L1Station.h: minor modifications; L1Field.h: magic constant
 '21' changed to the constexpr parameter L1Parameters::kN_FS_COEFFS

---
 reco/L1/CMakeLists.txt             |   1 +
 reco/L1/L1Algo/L1BaseStationInfo.h | 163 +++++++++++++++++++++++++----
 reco/L1/L1Algo/L1Field.h           |   8 +-
 reco/L1/L1Algo/L1Parameters.h      |  25 +++++
 reco/L1/L1Algo/L1Station.h         |  16 ++-
 5 files changed, 187 insertions(+), 26 deletions(-)
 create mode 100644 reco/L1/L1Algo/L1Parameters.h

diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 68e6fe0746..7e4abe3ced 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -270,6 +270,7 @@ Install(FILES CbmL1Counters.h
               L1Algo/L1UMeasurementInfo.h
               L1Algo/L1XYMeasurementInfo.h
               L1Algo/L1BaseStationInfo.h
+              L1Algo/L1Parameters.h
               vectors/vec_arithmetic.h
               vectors/std_alloc.h
         DESTINATION include
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.h b/reco/L1/L1Algo/L1BaseStationInfo.h
index d20e560c16..4b5a05701d 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.h
+++ b/reco/L1/L1Algo/L1BaseStationInfo.h
@@ -2,7 +2,7 @@
    SPDX-License-Identifier: GPL-3.0-only
    Authors: Sergey Gorbunov, Sergei Zharko [committer] */
 
-/**************************************************************************************************************
+/************************************************************************************************************
  * @file L1BaseStationInfo.h
  * @bried A base class for a L1 station interface
  * @since 18.12.2021
@@ -11,12 +11,13 @@
  * or BmnL1) with general L1Tracking algorithm. Each derived class must contain general
  * algorithms sutable for the particular station type.
  * 
- *************************************************************************************************************/
+ ***********************************************************************************************************/
 
 #ifndef L1BaseStationInfo_h
 #define L1BaseStationInfo_h 1
 
 // L1 core
+#include "L1Parameters.h"
 #include "L1Station.h"
 #include "L1Def.h"
 
@@ -24,18 +25,25 @@
 #include <string>
 #include <bitset>
 #include <iomanip>
+//#include <cmath>
 
 class L1BaseStationInfo {
 private:
   enum { // Here we list all the fields, which must be initialized by user
     /// v Basic fields initialization
     /// v L1Station initialization
-    kEz, 
+    kEtype,
+    kEtimeInfo,
+    kEz,
     kERmin, 
     kERmax, 
     kEmaterialInfoThick, 
     kEmaterialInfoRL,
     kEfieldSlice,
+    kEstripsFrontPhi,
+    kEstripsFrontSigma,
+    kEstripsBackPhi,
+    kEstripsBackSigma,
     /// v The last item is equal to the number of bits in fInitFlags
     kEND
   };
@@ -55,26 +63,36 @@ public:
     std::cout << ">>>>>> L1BaseStationInfo: Destructor called\n"; // Temporary
   }
 
-  ///-------------------------------------------------------------------------------------------------------///
-  ///   Basic methods                                                                                       ///
-  ///-------------------------------------------------------------------------------------------------------///
+  ///-----------------------------------------------------------------------------------------------------///
+  ///   Basic methods                                                                                     ///
+  ///-----------------------------------------------------------------------------------------------------///
   
   std::string GetTypeName()   const { return fTypeName; }
-  bool        IsInitialized() const { return fInitFlags.size() == fInitFlags.count(); };
+  // Checks if all the necessary fields are initialized by user
+  bool IsInitialized() const { return fInitFlags.size() == fInitFlags.count(); };
+  // Transfers all gathered data to L1Algo
+  void TransferL1Station() { 
+    /**********/
+  }
+  void TransferData() { /*********/ }    
 
-  ///-------------------------------------------------------------------------------------------------------///
-  ///   Interface for L1Station object initialization                                                       ///
-  ///-------------------------------------------------------------------------------------------------------///
+  ///-----------------------------------------------------------------------------------------------------///
+  ///   Interface for L1Station object initialization                                                     ///
+  ///-----------------------------------------------------------------------------------------------------///
 
   /// Setters
 
-  void SetZ(fscal inZ)        { fL1Station.z = inZ;        fInitFlags[kEz] = true; }
-  void SetRmin(fscal inRmin)  { fL1Station.Rmin = inRmin;  fInitFlags[kERmin] = true; }
-  void SetRmax(fscal inRmax)  { fL1Station.Rmax = inRmax;  fInitFlags[kERmax] = true; }
+  // Sets type of station // TODO: Temporary method, must be selected automatically in the derived class
+  void SetStationType(int inType)  { fL1Station.type = inType;         fInitFlags[kEtype] = true;}
+  void SetTimeInfo(int inTimeInfo) { fL1Station.timeInfo = inTimeInfo; fInitFlags[kEtimeInfo] = true;}
+
+  void SetZ(double inZ)        { fL1Station.z = inZ;        fInitFlags[kEz] = true; }
+  void SetRmin(double inRmin)  { fL1Station.Rmin = inRmin;  fInitFlags[kERmin] = true; }
+  void SetRmax(double inRmax)  { fL1Station.Rmax = inRmax;  fInitFlags[kERmax] = true; }
 
   // Sets material thickness and radiational length, calculates their ratio ("RadThick")
   // and its logarithm ("logRadThick")
-  void SetMaterial     (fscal inThickness, fscal inRL) {
+  void SetMaterial(double inThickness, double inRL) {
 #ifndef L1_NO_ASSERT // check for zero denominator
     L1_ASSERT(inRL, "Attempt of entering zero inRL (radiational length) value");
 #endif 
@@ -87,9 +105,8 @@ public:
   }
 
   // Sets coefficients of the magnetic field approximation
-  // TODO: Do something with N parameter
-  void SetFieldSlice (int N, const fscal * Cx, const fscal * Cy, const fscal * Cz) {
-    for (int idx = 0; idx < N; ++idx) {
+  void SetFieldSlice(const double * Cx, const double * Cy, const double * Cz) {
+    for (int idx = 0; idx < L1Parameters::kN_FS_COEFFS; ++idx) {
       fL1Station.fieldSlice.cx[idx] = Cx[idx];
       fL1Station.fieldSlice.cy[idx] = Cy[idx];
       fL1Station.fieldSlice.cz[idx] = Cz[idx];
@@ -97,6 +114,68 @@ public:
     }
   }
 
+  // Sets stereo angles and sigmas for front and back strips, automatically set frontInfo, backInfo, 
+  // XYInfo, xInfo and yInfo
+  void SetFBStripsGeometry(double f_phi, double f_sigma, double b_phi, double b_sigma) {
+    // TODO: may be is it better to define separate setters for these values and then calculate the 
+    //       rest somewhere below?
+
+    //----- Original code from L1Algo ---------------------------------------------------------------------//
+    double c_f     = cos(f_phi);
+    double s_f     = sin(f_phi);
+    double c_b     = cos(b_phi);
+    double s_b     = sin(b_phi);
+
+    // NOTE: Here additional double variables are used to save the precission
+
+    fL1Station.frontInfo.cos_phi = c_f;
+    fL1Station.frontInfo.sin_phi = s_f;
+    fL1Station.frontInfo.sigma2  = f_sigma * f_sigma;
+
+    fL1Station.backInfo.cos_phi = c_b;
+    fL1Station.backInfo.sin_phi = s_b;
+    fL1Station.backInfo.sigma2  = b_sigma * b_sigma;
+  
+    if (fabs(b_phi - f_phi) < .1) {
+      std::cout << "test";
+    }
+
+    //if (fabs(b_phi - f_phi) < .1) {
+    //  double th  = b_phi - f_phi;
+    //  double det = cos(th);
+    //  det *= det;
+    //  fL1Station.XYInfo.C00 = (s_b * s_b * f_sigma * f_sigma + s_f * s_f * b_sigma * b_sigma) / det;
+    //  fL1Station.XYInfo.C10 = -(s_b * c_b * f_sigma * f_sigma + s_f * c_f * b_sigma * b_sigma) / det;
+    //  fL1Station.XYInfo.C11 = (c_b * c_b * f_sigma * f_sigma + c_f * c_f * b_sigma * b_sigma) / det;
+    //} else {
+    //  double det = c_f * s_b - s_f * c_b;
+    //  det *= det;
+    //  fL1Station.XYInfo.C00 = (s_b * s_b * f_sigma * f_sigma + s_f * s_f * b_sigma * b_sigma) / det;
+    //  fL1Station.XYInfo.C10 = -(s_b * c_b * f_sigma * f_sigma + s_f * c_f * b_sigma * b_sigma) / det;
+    //  fL1Station.XYInfo.C11 = (c_b * c_b * f_sigma * f_sigma + c_f * c_f * b_sigma * b_sigma) / det;
+    //}
+
+    double det = (fabs(b_phi - f_phi) < 0.1) ? cos(b_phi - f_phi) : (c_f * s_b - s_f * c_b);
+    det *= det;
+    fL1Station.XYInfo.C00 = (s_b * s_b * f_sigma * f_sigma + s_f * s_f * b_sigma * b_sigma) / det;
+    fL1Station.XYInfo.C10 = -(s_b * c_b * f_sigma * f_sigma + s_f * c_f * b_sigma * b_sigma) / det;
+    fL1Station.XYInfo.C11 = (c_b * c_b * f_sigma * f_sigma + c_f * c_f * b_sigma * b_sigma) / det;
+
+    fL1Station.xInfo.cos_phi = -s_f / (c_f * s_b - c_b * s_f);
+    fL1Station.xInfo.sin_phi = s_b / (c_f * s_b - c_b * s_f);
+    fL1Station.xInfo.sigma2  = fL1Station.XYInfo.C00;
+
+    fL1Station.yInfo.cos_phi = c_b / (c_b * s_f - c_f * s_b);
+    fL1Station.yInfo.sin_phi = -c_f / (c_b * s_f - c_f * s_b);
+    fL1Station.yInfo.sigma2  = fL1Station.XYInfo.C11;
+    //-----------------------------------------------------------------------------------------------------//
+ 
+    fInitFlags[kEstripsFrontPhi]   = true;
+    fInitFlags[kEstripsFrontSigma] = true;
+    fInitFlags[kEstripsBackPhi]    = true;
+    fInitFlags[kEstripsBackSigma]  = true;
+  }
+
   /// Getters
   fvec GetZ()                    const { return fL1Station.z; }
   fvec GetRmin()                 const { return fL1Station.Rmin; }
@@ -119,17 +198,57 @@ public:
   const fvec * GetFieldSliceCy() const { return fL1Station.fieldSlice.cy; }
   const fvec * GetFieldSliceCz() const { return fL1Station.fieldSlice.cz; }
 
-  ///-------------------------------------------------------------------------------------------------------///
-  ///   Other methods                                                                                       ///
-  ///-------------------------------------------------------------------------------------------------------///
+  // Methods to get strip layers geometry
+
+  ///-----------------------------------------------------------------------------------------------------///
+  ///   Other methods                                                                                     ///
+  ///-----------------------------------------------------------------------------------------------------///
 
-  void Print() { // TODO: use here LOG() macro instead of std::cout
+  void Print() { 
     LOG(INFO) << "== L1Algo: station info ===========================================================";
     LOG(INFO) << "";
     LOG(INFO) << "  L1Station fields:";
-    LOG(INFO) << "    z                    " << fL1Station.z ;
+    LOG(INFO) << "    z position:              " << fL1Station.z[0] ;
+    LOG(INFO) << "    Rmin:                    " << fL1Station.Rmin[0] ;
+    LOG(INFO) << "    Rmax:                    " << fL1Station.Rmax[0] ;
+    LOG(INFO) << "    Thickness (X):           " << fL1Station.materialInfo.thick[0] ;
+    LOG(INFO) << "    Radiational length (X0): " << fL1Station.materialInfo.RL[0] ;
+    LOG(INFO) << "    X / X0:                  " << fL1Station.materialInfo.RadThick[0] ;
+    LOG(INFO) << "    log(X / X0):             " << fL1Station.materialInfo.logRadThick[0] ;
+    LOG(INFO) << "    Field approximation coefficients:";
+    LOG(INFO) << "      idx         CX         CY         CZ";
+    for (int idx = 0; idx < L1Parameters::kN_FS_COEFFS; ++idx) {
+      LOG(INFO) << std::setw(9)  << std::setfill(' ') << idx << ' '
+                << std::setw(10) << std::setfill(' ') << fL1Station.fieldSlice.cx[idx][0] << ' '
+                << std::setw(10) << std::setfill(' ') << fL1Station.fieldSlice.cy[idx][0] << ' '
+                << std::setw(10) << std::setfill(' ') << fL1Station.fieldSlice.cz[idx][0];
+    }
+    LOG(INFO) << "    Strips geometry:";
+    LOG(INFO) << "      Front:";
+    LOG(INFO) << "        cos(phi):            " << fL1Station.frontInfo.cos_phi[0];
+    LOG(INFO) << "        sin(phi):            " << fL1Station.frontInfo.sin_phi[0];
+    LOG(INFO) << "        sigma2:              " << fL1Station.frontInfo.sigma2[0];
+    LOG(INFO) << "      Back:";
+    LOG(INFO) << "        cos(phi):            " << fL1Station.backInfo.cos_phi[0];
+    LOG(INFO) << "        sin(phi):            " << fL1Station.backInfo.sin_phi[0];
+    LOG(INFO) << "        sigma2:              " << fL1Station.backInfo.sigma2[0];
+    LOG(INFO) << "      XY cov matrix:";
+    LOG(INFO) << "        C00:                 " << fL1Station.XYInfo.C00[0];
+    LOG(INFO) << "        C10:                 " << fL1Station.XYInfo.C10[0];
+    LOG(INFO) << "        C11:                 " << fL1Station.XYInfo.C11[0];
+    LOG(INFO) << "      X layer:";
+    LOG(INFO) << "        cos(phi):            " << fL1Station.xInfo.cos_phi[0];
+    LOG(INFO) << "        sin(phi):            " << fL1Station.xInfo.sin_phi[0];
+    LOG(INFO) << "        sigma2:              " << fL1Station.xInfo.sigma2[0];
+    LOG(INFO) << "      Y layer:";
+    LOG(INFO) << "        cos(phi):            " << fL1Station.yInfo.cos_phi[0];
+    LOG(INFO) << "        sin(phi):            " << fL1Station.yInfo.sin_phi[0];
+    LOG(INFO) << "        sigma2:              " << fL1Station.yInfo.sigma2[0];
+    LOG(INFO) << "";
     LOG(INFO) << "===================================================================================";
   }
+  
+
 };
 
 #endif // L1BaseStationInfo_h
diff --git a/reco/L1/L1Algo/L1Field.h b/reco/L1/L1Algo/L1Field.h
index 3e0fa2433d..114313ab58 100644
--- a/reco/L1/L1Algo/L1Field.h
+++ b/reco/L1/L1Algo/L1Field.h
@@ -8,6 +8,8 @@
 #include <iostream>
 
 #include "L1Def.h"
+#include "L1Parameters.h"
+
 using std::cout;
 using std::endl;
 using std::ostream;
@@ -38,11 +40,13 @@ public:
 class L1FieldSlice {
 
 public:
-  fvec cx[21], cy[21], cz[21];  // polinom coeff.
+  fvec cx[L1Parameters::kN_FS_COEFFS];
+  fvec cy[L1Parameters::kN_FS_COEFFS];
+  fvec cz[L1Parameters::kN_FS_COEFFS];  // polinom coeff.
 
   L1FieldSlice()
   {
-    for (int i = 0; i < 21; i++)
+    for (int i = 0; i < L1Parameters::kN_FS_COEFFS; ++i)
       cx[i] = cy[i] = cz[i] = 0;
   }
 
diff --git a/reco/L1/L1Algo/L1Parameters.h b/reco/L1/L1Algo/L1Parameters.h
new file mode 100644
index 0000000000..72b49c2c46
--- /dev/null
+++ b/reco/L1/L1Algo/L1Parameters.h
@@ -0,0 +1,25 @@
+/* Copyright (C) 2016-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+/************************************************************************************************************
+ * @file L1Parameters.h
+ * @brief Parameter container for the L1Algo library
+ * @since 19.12.2021
+ *
+ ***********************************************************************************************************/
+
+#ifndef L1Parameters_h
+#define L1Parameters_h 1
+
+class L1Parameters {
+public:
+  ///-----------------------------------------------------------------------------------------------------///
+  ///   Array sizes                                                                                       ///
+  ///-----------------------------------------------------------------------------------------------------///
+  // Amount of elements in cx, cy and cz coefficient arrays of L1FieldSlice
+  static constexpr int kN_FS_COEFFS {21}; // TODO: May be we should return the order of polynomial?
+
+};
+
+#endif // L1Parameters_h
diff --git a/reco/L1/L1Algo/L1Station.h b/reco/L1/L1Algo/L1Station.h
index c409f632ec..fd3770e428 100644
--- a/reco/L1/L1Algo/L1Station.h
+++ b/reco/L1/L1Algo/L1Station.h
@@ -31,12 +31,24 @@ public:
   {
   }
 
+  // TODO: test this constructors
+  //L1Station(const L1Station &) = default;
+  //L1Station & operator=(const L1Station &) = default;
+  //L1Station(L1Station &&) = default;
+  //L1Station & operator=(L1Station &&) = default;
+
   int type;
   int timeInfo;
-  fvec z, Rmin, Rmax, Sy;
+  fvec z;
+  fvec Rmin;
+  fvec Rmax;
+  fvec Sy;
   L1MaterialInfo materialInfo;
   L1FieldSlice fieldSlice;
-  L1UMeasurementInfo frontInfo, backInfo, xInfo, yInfo;
+  L1UMeasurementInfo frontInfo;
+  L1UMeasurementInfo backInfo;
+  L1UMeasurementInfo xInfo;
+  L1UMeasurementInfo yInfo;
   L1XYMeasurementInfo XYInfo;
 
 } _fvecalignment;
-- 
GitLab