diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index 7b4f244a78554358684b28aeed5e9c1a560545b9..9b4685bb30be32ae297ab86886336a72eac754d5 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -36,6 +36,7 @@ class L1AlgoDraw;
//#endif
//#define MERGE_CLONES
+#define FEATURING_L1ALGO_INIT
#include <iomanip>
#include <iostream>
@@ -48,6 +49,8 @@ class L1AlgoDraw;
#include "L1Hit.h"
#include "L1HitPoint.h"
#include "L1HitsSortHelper.h"
+#include "L1InitManager.h"
+#include "L1Parameters.h"
#include "L1Portion.h"
#include "L1Station.h"
#include "L1Track.h"
@@ -135,7 +138,6 @@ public:
L1Vector<L1Triplet> fTriplets[L1Parameters::kMaxNstations][L1Parameters::kMaxNthreads] {
- //<--------
{"L1Algo::fTriplets"}}; // created triplets at station + thread
// Track candidates created out of adjacent triplets before the final track selection.
@@ -359,10 +361,12 @@ public:
// TODO: We should think about, where non-constexpr L1Alog parameters can be modified. At the moment we can create a
// L1Parameters object somewhere outside the L1Algo, fill its fields there and then pass it directly to
// the L1Algo instance.
+ L1InitManager * GetL1InitManager() { return &fL1InitManager; }
private:
/// Object containing L1Parameters. Default consturctor is used
- L1Parameters fL1Parameters;
+ L1Parameters fL1Parameters; ///< Object of L1Algo parameters class
+ L1InitManager fL1InitManager; ///< Object of L1Algo initialization manager class
/// ================================= FUNCTIONAL PART =================================
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.h b/reco/L1/L1Algo/L1BaseStationInfo.h
index ae09356723001e8e95bb65354892edc40f6784f9..f298ea92178bc7df0971190f2e8ef5afd1b0a4bf 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.h
+++ b/reco/L1/L1Algo/L1BaseStationInfo.h
@@ -21,10 +21,18 @@
#include "L1Parameters.h"
#include "L1Station.h"
+// FairRoot
+#include "FairField.h"
+
+// Root
+#include "TMatrixD.h"
+#include "TVectorD.h"
+
// C++ std
#include <bitset>
#include <iomanip>
#include <string>
+#include <functional>
//#include <cmath>
/// A base class which provides interface to L1Algo station geometry
@@ -34,6 +42,9 @@ private:
enum
{ // Here we list all the fields, which must be initialized by user
// Basic fields initialization
+ kEStationID,
+ kEXmax,
+ kEYmax,
// L1Station initialization
kEtype,
kEtimeInfo,
@@ -84,11 +95,31 @@ public:
//
// Setters
//
+ /// Sets station ID
+ void SetStationID(int inID)
+ {
+ if (! fInitFlags[kEStationID]) {
+ fStationID = inID;
+ fInitFlags[kEStationID] = true;
+ }
+#ifndef FAST_CODE
+ else {
+ LOG(WARNING) << __FILE__ << ":" << __LINE__ << " Attempt of station ID redifinition\n";
+ }
+#endif // ! FAST_CODE
+ }
/// Sets type of station
void SetStationType(int inType)
{
- fL1Station.type = inType;
- fInitFlags[kEtype] = true;
+ if (! fInitFlags[kEtype]) {
+ fL1Station.type = inType;
+ fInitFlags[kEtype] = true;
+ }
+#ifndef FAST_CODE
+ else {
+ LOG(WARNING) << __FILE__ << ":" << __LINE__ << " Attempt of station type redifinition\n";
+ }
+#endif // ! FAST_CODE
}
// TODO: Temporary method, a type must be selected automatically in the derived classes
void SetTimeInfo(int inTimeInfo)
@@ -101,7 +132,8 @@ public:
/// Sets nominal z position of the station
void SetZ(double inZ)
{
- fL1Station.z = inZ;
+ fL1Station.z = inZ; // setting simd vector of single-precision floats, which is passed to high performanced L1Algo
+ fZPos = inZ; // setting precise value to use in field approximation etc
fInitFlags[kEz] = true;
}
/// Sets min transverse size of the station
@@ -138,8 +170,90 @@ public:
fL1Station.fieldSlice.cx[idx] = Cx[idx];
fL1Station.fieldSlice.cy[idx] = Cy[idx];
fL1Station.fieldSlice.cz[idx] = Cz[idx];
- fInitFlags[kEfieldSlice] = true;
}
+ fInitFlags[kEfieldSlice] = true;
+ }
+ /// Sets arrays of the approcimation
+ /// \param getField A user function, which gets a xyz array of position coordinates and fills B array
+ /// of magnetic field components in position
+ void SetFieldSlice(const std::function<void(const double (&xyz)[3], double (&B)[3])>& getFieldValue)
+ {
+ std::cout << ">>>>>>> CALL: SetFieldSlice\n";
+#ifndef L1_NO_ASSERT // check for zero denominator
+ L1_ASSERT(fInitFlags[kEz], "Attempt to set magnetic field slice before setting z position of the station");
+ L1_ASSERT(fInitFlags[kEXmax], "Attempt to set magnetic field slice before setting Xmax size of the station");
+ L1_ASSERT(fInitFlags[kEYmax], "Attempt to set magnetic field slice before setting Ymax size of the station");
+#endif
+ constexpr int M = L1Parameters::kMaxFieldApproxPolynomialOrder;
+ constexpr int N = L1Parameters::kMaxNFieldApproxCoefficients;
+ constexpr int D = 3; ///> number of dimensions
+
+ // SLE initialization
+ double A[N][N + D] = {}; // augmented matrix
+ double dx = (fXmax / N / 2 < 1.) ? fXmax / N / 4. : 1.;
+ double dy = (fYmax / N / 2 < 1.) ? fYmax / N / 4. : 1.;
+
+ for (double x = -fXmax; x <= fXmax; x += dx) {
+ for (double y = -fYmax; y <= fYmax; y += dy) {
+ double r = sqrt(fabs(x * x / fXmax / fXmax + y / fYmax * y / fYmax));
+ if (r > 1.) {
+ continue;
+ }
+ double p[D] = {x, y, fZPos};
+ double B[D] = {};
+ getFieldValue(p, B);
+
+ double m[N] = {1};
+ m[0] = 1;
+ for (int i = 1; i <= M; ++i) {
+ int k = (i - 1) * (i) / 2;
+ int l = i * (i + 1) / 2;
+ for (int j = 0; j < i; ++j) {
+ m[l + j] = x * m[k + j];
+ }
+ m[l + i] = y * m[k + i - 1];
+ }
+
+ double w = 1. / (r * r + 1);
+ for (int i = 0; i < N; ++i) {
+ // fill the left part of the matrix
+ for (int j = 0; j < N; ++j) {
+ A[i][j] += w * m[i] * m[j];
+ }
+ // fill the right part of the matrix
+ for (int j = 0; j < D; ++j) {
+ A[i][N + j] += w * B[j] * m[i];
+ }
+ }
+ }
+ }
+
+ // SLE solution (Gaussian elimination)
+ //
+ for (int kCol = 0; kCol < N - 1; ++kCol) {
+ for (int jRow = kCol + 1; jRow < N; ++jRow) {
+ double factor = A[jRow][kCol] / A[kCol][kCol];
+ for (int iCol = kCol; iCol < N + D; ++iCol) {
+ A[jRow][iCol] -= factor * A[kCol][iCol];
+ }
+ }
+ }
+ for (int kCol = N - 1; kCol > 0; --kCol) {
+ for (int jRow = kCol - 1; jRow >= 0; --jRow) {
+ double factor = A[jRow][kCol] / A[kCol][kCol];
+ for (int iCol = kCol; iCol < N + D; ++iCol) {
+ A[jRow][iCol] -= factor * A[kCol][iCol];
+ }
+ }
+ }
+
+ for (int j = 0; j < N; ++j) {
+ fL1Station.fieldSlice.cx[j] = A[j][N + 0] / A[j][j];
+ fL1Station.fieldSlice.cy[j] = A[j][N + 1] / A[j][j];
+ fL1Station.fieldSlice.cz[j] = A[j][N + 2] / A[j][j];
+ }
+
+ fInitFlags[kEfieldSlice] = true;
}
/// Sets stereo angles and sigmas for front and back strips, automatically set frontInfo, backInfo,
@@ -201,13 +315,31 @@ public:
fInitFlags[kEstripsBackPhi] = true;
fInitFlags[kEstripsBackSigma] = true;
}
+ /// Sets maximum distance between station center and its edge in x direction
+ void SetXmax(double aSize)
+ {
+ fXmax = aSize;
+ fInitFlags[kEXmax] = true;
+ }
+ /// Sets maximum distance between station center and its edge in y direction
+ void SetYmax(double aSize)
+ {
+ fYmax = aSize;
+ fInitFlags[kEYmax] = true;
+ }
//
// Getters
//
- /// Gets nominal z position of the station
- fvec GetZ() const { return fL1Station.z; }
+ /// Gets station ID
+ int GetStationID() const { return fStationID; }
+ /// Gets station type
+ int GetStationType() const { return fL1Station.type; }
+ /// Gets SIMD vectorized z position of the station
+ fvec GetZsimdVec() const { return fL1Station.z; }
+ /// Gets double precised z position of the station
+ double GetZdouble() const { return fZPos; }
/// Gets min transverse size of the station
fvec GetRmin() const { return fL1Station.Rmin; }
/// Gets max transverse size of the station
@@ -238,6 +370,11 @@ public:
/// Gets array of the coefficients for the field z-component approximation
const fvec* GetFieldSliceCz() const { return fL1Station.fieldSlice.cz; }
+ /// Gets maximum distance between station center and its edge in x direction
+ double GetXmax() const { return fXmax; }
+ /// Gets maximum distance between station center and its edge in y direction
+ double GetYmax() const { return fYmax; }
+
//-------------------------------------------------------------------------------------------------------//
// Other methods //
@@ -248,7 +385,9 @@ public:
{
LOG(INFO) << "== L1Algo: station info ===========================================================";
LOG(INFO) << "";
+ LOG(INFO) << " Station ID: " << fStationID;
LOG(INFO) << " L1Station fields:";
+ LOG(INFO) << " Station type ID: " << fL1Station.type;
LOG(INFO) << " z position: " << fL1Station.z[0];
LOG(INFO) << " Rmin: " << fL1Station.Rmin[0];
LOG(INFO) << " Rmax: " << fL1Station.Rmax[0];
@@ -286,14 +425,26 @@ public:
LOG(INFO) << " sin(phi): " << fL1Station.yInfo.sin_phi[0];
LOG(INFO) << " sigma2: " << fL1Station.yInfo.sigma2[0];
LOG(INFO) << "";
+ LOG(INFO) << " Additiona fields:";
+ LOG(INFO) << " Xmax: " << fXmax;
+ LOG(INFO) << " Ymax: " << fYmax;
+ LOG(INFO) << "";
LOG(INFO) << "===================================================================================";
}
+ void Reset() {
+ L1BaseStationInfo other;
+ std::swap(*this, other);
+ }
private:
- /// Flags set. A particular flag is up if the corresponding setter was called (experimental)
- std::bitset<L1BaseStationInfo::kEND> fInitFlags;
- std::string fTypeName {"BASE"};
- L1Station fL1Station {};
+
+ std::bitset<L1BaseStationInfo::kEND> fInitFlags; ///< Class fileds initialization flags
+ std::string fTypeName {"BASE"}; ///< Station type name (TODO: probably should be replaced with function of type)
+ L1Station fL1Station {}; ///< L1Station structure, describes a station in L1Algo
+ int fStationID {-1}; ///< Station ID
+ double fXmax {0}; ///< Maximum distance between station center and its edge in x direction
+ double fYmax {0}; ///< Maximum distance between station center and its edge in y direction
+ double fZPos {0}; ///< z position of the station in double precision, used in field approximation
// TODO (!!!!) MUST THINK ABOUT THIS OBJECT LIFE TIME: will it be better to copy
// or to move it to the core? If the second, what should we do with
// getters? Do we need to lock them, when the fL1Station object will
diff --git a/reco/L1/L1Algo/L1InitManager.h b/reco/L1/L1Algo/L1InitManager.h
new file mode 100644
index 0000000000000000000000000000000000000000..d46295478cbbb1b2300e14ee694853cfc51255cd
--- /dev/null
+++ b/reco/L1/L1Algo/L1InitManager.h
@@ -0,0 +1,73 @@
+/* Copyright (C) 2016-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+/************************************************************************************************************
+ * @file L1InitManager.h
+ * @bried Input data management class for L1Algo
+ * @since 24.12.2021
+ *
+ ***********************************************************************************************************/
+
+#include "L1BaseStationInfo.h"
+#include "L1Parameters.h"
+#include "L1Vector.h"
+
+//#include <string>
+#include <unordered_map>
+
+/// Initialization manager for L1Aglo
+///
+/// The L1InitManager class provides interface for L1Algo initialization
+
+class L1InitManager {
+public:
+ /// Default constructor
+ L1InitManager() = default;
+ /// Destructor
+ ~L1InitManager() = default;
+ /// Copy constructor is forbidden
+ L1InitManager(const L1InitManager& /*other*/) = delete;
+ /// Move constructor is forbidden
+ L1InitManager(L1InitManager&& /*other*/) = delete;
+ /// Copy assignment operator is forbidden
+ L1InitManager& operator=(const L1InitManager& /*other*/) = delete;
+ /// Move assignment operator is forbidden
+ L1InitManager& operator=(L1InitManager&& /*other*/) = delete;
+
+ /// Adds another station of a given type using reference to a L1BaseStationInfo object
+ void AddStation(const L1BaseStationInfo& inStation)
+ {
+ fvStationsInfo.push_back(inStation);
+ ++fmStationsNumber[inStation.GetStationType()];
+ // TODO: Think, how we would like to initialize station ID
+ }
+ /// Adds another station of a given type using pointer to a L1BaseStationInfo object
+ void AddStation(const L1BaseStationInfo* inStationPtr)
+ {
+ fvStationsInfo.push_back(*inStationPtr); // Copy is occured
+ ++fmStationsNumber[inStationPtr->GetStationType()];
+ }
+
+ /// Prints out a list of stations
+ void PrintStations(int verbosityLevel = 0)
+ {
+ if (verbosityLevel < 1) {
+ for (auto& station : fvStationsInfo) {
+ LOG(INFO) << "----------- station: ";
+ LOG(INFO) << "\ttype = " << station.GetStationType(); // TMP
+ LOG(INFO) << "\tz = " << station.GetZdouble();
+ }
+ }
+ else {
+ for (auto& station : fvStationsInfo) {
+ station.Print();
+ }
+ }
+ }
+
+private:
+ L1Vector<L1BaseStationInfo> fvStationsInfo {"L1InitManager::fvStationsInfo"}; ///< Vector containing all the stations
+ // TODO: Max stations is probably not optimal for vector initialization. Must think, what to do here.
+ std::unordered_map<int, int> fmStationsNumber; ///< Number of each station type <station type, counter>
+};
diff --git a/reco/L1/L1Algo/L1Parameters.h b/reco/L1/L1Algo/L1Parameters.h
index 718d3bea164b134385f40cc9cd33f437b94bdf39..94de70cf100beb82e5accaaee22901779a438169 100644
--- a/reco/L1/L1Algo/L1Parameters.h
+++ b/reco/L1/L1Algo/L1Parameters.h
@@ -13,6 +13,44 @@
#define L1Parameters_h 1
#include <FairLogger.h>
+#include <string>
+#include <unordered_map>
+#include <iomanip>
+
+#include "L1Vector.h"
+
+//---------------------------------------------------------------------------------------------------------//
+// Track finder iterations definition
+//---------------------------------------------------------------------------------------------------------//
+
+// TODO: May be it is reasonable to create TrackFinderIterations class
+
+/// Different iterations for Track Finder running
+enum class TrackFinderIter {
+ kFastPrim, ///< primary fast tracks
+ kAllPrim, ///< primary all tracks
+ kAllPrimJump, ///< primary all tracks with jumped triplets
+ kAllSec, ///< secondary all tracks
+ kAllPrimE, ///< primary all electron tracks
+ kAllSecE, ///< secondary all electron tracks
+ kFastPrimJump, ///< primary fast tracks with jumped triplets
+ kAllSecJump ///< secondary all tracks with jumped triplets
+};
+
+
+const std::unordered_map<TrackFinderIter, std::string> kTrackFinderIterNames {
+ { TrackFinderIter::kFastPrim, "Primary tracks, fast" },
+ { TrackFinderIter::kAllPrim, "Primary tracks, all" },
+ { TrackFinderIter::kAllPrimJump, "Primary tracks with jumped triplets, all" },
+ { TrackFinderIter::kAllSec, "Secondary tracks, all" },
+ { TrackFinderIter::kAllPrimE, "Primary electron tracks, all" },
+ { TrackFinderIter::kAllSecE, "Secondary electron tracks, all" },
+ { TrackFinderIter::kFastPrimJump, "Primary tracks with jumped triplets, fast" },
+ { TrackFinderIter::kAllSecJump, "Secondary tracks with jumped triplets, all" }
+};
+
+template <typename T>
+using TrackFinderIterParameter_t = std::unordered_map<TrackFinderIter, T>;
class L1Parameters {
public:
@@ -22,19 +60,15 @@ public:
/// Amount of coefficients in field approximations
static constexpr int kMaxNFieldApproxCoefficients {21}; // TODO: May be it is better to use the polynomial
// order instead of this?
- /// Amount of bits to code one station
- static constexpr unsigned int kStationBits {6};
- /// Amount of bits to code one thread
- static constexpr unsigned int kThreadBits {6};
- /// Amount of bits to code one triple
- static constexpr unsigned int kTripletBits {32 - kStationBits - kThreadBits};
-
- /// Max number of stations
- static constexpr unsigned int kMaxNstations {1 << kStationBits}; // 2^6 = 64
- /// Max number of threads
- static constexpr unsigned int kMaxNthreads {1 << kThreadBits}; // 2^6 = 64
- /// Max number of triplets
- static constexpr unsigned int kMaxNtriplets {1 << kTripletBits}; // 2^20 = 1,048,576
+ static constexpr int kMaxFieldApproxPolynomialOrder {5}; ///> Order of polynomial to approximate field in the vicinity of station plane
+
+ static constexpr unsigned int kStationBits {6}; ///> Amount of bits to code one station
+ static constexpr unsigned int kThreadBits {6}; ///> Amount of bits to code one thread
+ static constexpr unsigned int kTripletBits {32 - kStationBits - kThreadBits}; ///> Amount of bits to code one triple
+
+ static constexpr unsigned int kMaxNstations {1 << kStationBits}; ///> Max number of stations, 2^6 = 64
+ static constexpr unsigned int kMaxNthreads {1 << kThreadBits}; ///> Max number of threads, 2^6 = 64
+ static constexpr unsigned int kMaxNtriplets {1 << kTripletBits}; ///> Max number of triplets, 2^20 = 1,048,576
public:
//-------------------------------------------------------------------------------------------------------//
@@ -71,26 +105,90 @@ public:
LOG(INFO) << "===================================================================================";
}
// TODO: change constant names with actual (human) names
- // TODO: create a map with parameters, their description and
+ // TODO: create a map with parameters, their description etc.
public:
//-------------------------------------------------------------------------------------------------------//
- // Runtime constants //
+ // Runtime constants //
//-------------------------------------------------------------------------------------------------------//
/// Sets upper-bound cut on max number of doublets per one singlet
void SetMaxDoubletsPerSinglet(unsigned int value) { fMaxDoubletsPerSinglet = value; }
/// Sets upper-bound cut on max number of triplets per one doublet
void SetMaxTripletPerDoublets(unsigned int value) { fMaxTripletPerDoublets = value; }
+ /// Sets track finder iteration sequence
+ void SetTrackFinderIterSequence(const L1Vector<TrackFinderIter>& iterations) {
+ fTrackFinderIterSequence = iterations;
+ }
+
+
unsigned int GetMaxDoubletsPerSinglet() const { return fMaxDoubletsPerSinglet; }
unsigned int GetMaxTripletPerDoublets() const { return fMaxTripletPerDoublets; }
+ //-------------------------------------------------------------------------------------------------------//
+ // Auxilary methods //
+ //-------------------------------------------------------------------------------------------------------//
+
+ void PrintTrackFinderIterationsSequence() const
+ {
+ LOG(INFO) << "== L1Algo track finder iterations =================================================";
+ LOG(INFO) << "";
+ int idx = 0;
+ for (auto& iteration: fTrackFinderIterSequence) {
+ LOG(INFO) << " " << std::setw(2) << std::setfill(' ') << idx << ' ' << kTrackFinderIterNames.at(iteration);
+ ++idx;
+ }
+ LOG(INFO) << "";
+ LOG(INFO) << "===================================================================================";
+
+ }
private:
- ///
- unsigned int fMaxDoubletsPerSinglet {150};
- ///
- unsigned int fMaxTripletPerDoublets {15};
+ unsigned int fMaxDoubletsPerSinglet {150}; ///<
+ unsigned int fMaxTripletPerDoublets {15}; ///<
+
+ L1Vector<TrackFinderIter> fTrackFinderIterSequence {}; ///< Sequence of iterations in the track finder algorithm
+
+ // TODO: Develope naming system for these constants
+
+ std::unordered_map<TrackFinderIter, float> fTrackChi2Cut {
+ { TrackFinderIter::kFastPrim, 10.f },
+ { TrackFinderIter::kAllPrim, 10.f },
+ { TrackFinderIter::kAllPrimJump, 10.f },
+ { TrackFinderIter::kAllSec, 10.f },
+ { TrackFinderIter::kAllPrimE, 10.f },
+ { TrackFinderIter::kAllSecE, 10.f },
+ { TrackFinderIter::kFastPrimJump, 10.f },
+ { TrackFinderIter::kAllSecJump, 10.f }
+ };
+
+ std::unordered_map<TrackFinderIter, float> fTripletChi2Cut {
+ { TrackFinderIter::kFastPrim, 5.f },
+ { TrackFinderIter::kAllPrim, 5.f },
+ { TrackFinderIter::kAllPrimJump, 5.f },
+ { TrackFinderIter::kAllSec, 5.f },
+ { TrackFinderIter::kAllPrimE, 5.f },
+ { TrackFinderIter::kAllSecE, 5.f },
+ { TrackFinderIter::kFastPrimJump, 5.f },
+ { TrackFinderIter::kAllSecJump, 5.f }
+ };
+
+ std::unordered_map<TrackFinderIter, float> fDoubletChi2Cut {
+ { TrackFinderIter::kFastPrim, 5.f },
+ { TrackFinderIter::kAllPrim, 5.f },
+ { TrackFinderIter::kAllPrimJump, 5.f },
+ { TrackFinderIter::kAllSec, 5.f },
+ { TrackFinderIter::kAllPrimE, 5.f },
+ { TrackFinderIter::kAllSecE, 5.f },
+ { TrackFinderIter::kFastPrimJump, 5.f },
+ { TrackFinderIter::kAllSecJump, 5.f }
+ };
+
+
+
+ //float fTripletChi2Cut {5.f};
+ //float fTripletChi2Cut {5.f};
+
};
#endif // L1Parameters_h