From 2cce5b1b415232f9738210c3c4972b3b7130d7a0 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Sat, 5 Feb 2022 19:17:03 +0100
Subject: [PATCH] L1Algo init interface updates
---
reco/L1/CMakeLists.txt | 6 +-
reco/L1/CbmL1.cxx | 93 ++++--
reco/L1/CbmL1.h | 26 +-
reco/L1/L1Algo/L1Algo.cxx | 36 ++
reco/L1/L1Algo/L1Algo.h | 8 +-
reco/L1/L1Algo/L1BaseStationInfo.cxx | 482 +++++++++++++++++++++++++++
reco/L1/L1Algo/L1BaseStationInfo.h | 446 ++++++-------------------
reco/L1/L1Algo/L1CAIteration.cxx | 25 ++
reco/L1/L1Algo/L1CAIteration.h | 121 +++++++
reco/L1/L1Algo/L1Def.h | 5 +-
reco/L1/L1Algo/L1InitManager.cxx | 291 ++++++++++++++++
reco/L1/L1Algo/L1InitManager.h | 167 ++++++++--
reco/L1/L1Algo/L1Parameters.h | 279 ++++++++--------
reco/L1/L1Algo/L1Station.h | 107 ++++--
reco/L1/L1Algo/L1Vector.h | 5 +-
15 files changed, 1510 insertions(+), 587 deletions(-)
create mode 100644 reco/L1/L1Algo/L1BaseStationInfo.cxx
create mode 100644 reco/L1/L1Algo/L1CAIteration.cxx
create mode 100644 reco/L1/L1Algo/L1CAIteration.h
create mode 100644 reco/L1/L1Algo/L1InitManager.cxx
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 6c6b0fefc5..1861b9a600 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -140,7 +140,9 @@ L1Algo/L1EventMatch.cxx
L1Algo/L1MCEvent.cxx
L1Algo/L1Fit.cxx
CbmL1MCTrack.cxx
-
+L1Algo/L1CAIteration.cxx
+L1Algo/L1BaseStationInfo.cxx
+L1Algo/L1InitManager.cxx
L1Algo/utils/L1AlgoDraw.cxx
L1Algo/utils/L1AlgoEfficiencyPerformance.cxx
L1Algo/utils/L1AlgoPulls.cxx
@@ -270,6 +272,8 @@ Install(FILES CbmL1Counters.h
L1Algo/L1UMeasurementInfo.h
L1Algo/L1XYMeasurementInfo.h
L1Algo/L1BaseStationInfo.h
+ L1Algo/L1InitManager.h
+ L1Algo/L1CAIteration.h
L1Algo/L1Parameters.h
L1Algo/utils/L1Functions.h
vectors/vec_arithmetic.h
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 5daeaef357..f340ac7334 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -768,36 +768,80 @@ InitStatus CbmL1::Init()
algo->SetL1Parameters(fL1Parameters);
+
+#ifdef FEATURING_L1ALGO_INIT
/********************************************************************************************************************
- * *
* EXPERIMENTAL FEATURE: usage of L1InitManager for L1Algo initialization *
- * *
- * Stage 1: Repeat initialization steps in parallel to the original code *
- * Stage 2: Compare two initialization procedures *
- * Stage 3: Remove old initialization *
- * *
********************************************************************************************************************/
- // Get reference to the L1Algo initialization manager
+ // NOTE: The
+
+ // Step 0: Get reference to the L1Algo initialization manager
L1InitManager * initMan = algo->GetL1InitManager();
+
+ // Step 1: Initialize magnetic field function
+ // Set magnetic field slices
+ auto fieldGetterFcn =
+ [](const double (&inPos)[3], double (&outB)[3]) {CbmKF::Instance()->GetMagneticField()->GetFieldValue(inPos, outB);};
+ initMan->SetFieldFunction(fieldGetterFcn);
+
+
+ /// TODO: temporary for tests, must be initialized somewhere in run_reco.C or similar
+ fActiveTrackingDetectorIDs = {L1DetectorID::kMvd, L1DetectorID::kSts };
+
+ /// Step 2: initialize IDs of detectors active in tracking
+ initMan->SetActiveDetectorIDs(fActiveTrackingDetectorIDs);
+
+
constexpr double PI = 3.14159265358; // TODO: why cmath is not used?
- // Fill STS info:
- for (int iSt = 0; iSt < NStsStations; ++iSt) {
+ /// Step 2: initialize IDs
+ initMan->SetStationsNumberCrosscheck(L1DetectorID::kMvd, NMvdStations);
+ initMan->SetStationsNumberCrosscheck(L1DetectorID::kSts, NStsStations);
+ initMan->SetStationsNumberCrosscheck(L1DetectorID::kMuch, NMuchStations);
+
+
+ // Setup MVD stations info
+ for (int iSt = 0; iSt < NMvdStations; ++iSt) { // NOTE: example using in-stack defined objects
+ CbmMvdDetector* mvdDetector = CbmMvdDetector::Instance();
+ CbmMvdStationPar* mvdStationPar = mvdDetector->GetParameterFile();
+
+ CbmKFTube& t = CbmKF::Instance()->vMvdMaterial[iSt];
+ auto stationInfo = L1BaseStationInfo(L1DetectorID::kMvd, iSt);
+ stationInfo.SetStationType(1); // MVD // TODO: to be exchanged with specific flags (timeInfo, fieldInfo etc.)
+ stationInfo.SetTimeInfo(0);
+ stationInfo.SetZ(t.z);
+ stationInfo.SetMaterial(t.dz, t.RadLength);
+ stationInfo.SetXmax(t.R);
+ stationInfo.SetYmax(t.R);
+ stationInfo.SetRmin(t.r);
+ stationInfo.SetRmax(t.R);
+ fscal mvdFrontPhi = 0;
+ fscal mvdBackPhi = PI / 2.;
+ fscal mvdFrontSigma = mvdStationPar->GetXRes(iSt) / 10000;
+ fscal mvdBackSigma = mvdStationPar->GetYRes(iSt) / 10000;
+ stationInfo.SetFrontBackStripsGeometry(mvdFrontPhi, mvdFrontSigma, mvdBackPhi, mvdBackSigma);
+ initMan->AddStation(stationInfo);
+ }
+
+ // Setup STS stations info
+ for (int iSt = 0; iSt < NStsStations; ++iSt) { // NOTE: example using smart pointers
auto cbmSts = CbmStsSetup::Instance()->GetStation(iSt);
- auto stsStation = new L1BaseStationInfo();
- stsStation->SetStationID(iSt);
- stsStation->SetStationType(0); // STS
+ std::unique_ptr<L1BaseStationInfo> stationInfo(new L1BaseStationInfo(L1DetectorID::kSts, iSt));
+ // TODO: replace with std::make_unique, when C++14 is avaliable!!!!
+ // auto stsStation = std::make_unique<L1BaseStationInfo>(L1DetectorID::kSts, iSt);
+ stationInfo->SetStationType(0); // STS
+ stationInfo->SetTimeInfo(0);
// Setup station geometry and material
- stsStation->SetZ(cbmSts->GetZ());
+ stationInfo->SetZ(cbmSts->GetZ());
double stsXmax = cbmSts->GetXmax();
double stsYmax = cbmSts->GetYmax();
- stsStation->SetXmax(stsXmax);
- stsStation->SetYmax(stsYmax);
- stsStation->SetRmin(0);
- stsStation->SetRmax(stsXmax > stsYmax ? stsXmax : stsYmax);
- stsStation->SetMaterial(cbmSts->GetSensorD(), cbmSts->GetRadLength());
+ stationInfo->SetXmax(stsXmax);
+ stationInfo->SetYmax(stsYmax);
+ stationInfo->SetRmin(0);
+ stationInfo->SetRmax(stsXmax > stsYmax ? stsXmax : stsYmax);
+ stationInfo->SetMaterial(cbmSts->GetSensorD(), cbmSts->GetRadLength());
// Setup strips geometry
// TODO: why fscal instead of double in initialization?
@@ -805,22 +849,15 @@ InitStatus CbmL1::Init()
fscal stsBackPhi = cbmSts->GetSensorRotation() + cbmSts->GetSensorStereoAngle(1) * PI / 180.;
fscal stsFrontSigma = cbmSts->GetSensorPitch(0) / sqrt(12);
fscal stsBackSigma = stsFrontSigma;
- stsStation->SetFrontBackStripsGeometry(stsFrontPhi, stsFrontSigma, stsBackPhi, stsBackSigma);
-
- // Setup magnetic field
- // NOTE: Such solution is needed to prevent L1Algo from FairRoot dependencies
- auto getFieldValueFcn = [](const double (&inXYZ)[3], double (&outB)[3]) {
- CbmKF::Instance()->GetMagneticField()->GetFieldValue(inXYZ, outB);
- };
- stsStation->SetFieldSlice(getFieldValueFcn);
- initMan->AddStation(stsStation);
- delete stsStation;
+ stationInfo->SetFrontBackStripsGeometry(stsFrontPhi, stsFrontSigma, stsBackPhi, stsBackSigma);
+ initMan->AddStation(stationInfo);
}
initMan->PrintStations(/*vebosity = */ 1);
/********************************************************************************************************************
********************************************************************************************************************/
+#endif // FEATURING_L1ALGO_INIT
algo->Init(geo, fUseHitErrors, fTrackingMode, fMissingHits);
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 4428b32d51..c14f056a6c 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -76,10 +76,22 @@ class CbmL1HitStore {
public:
int ExtIndex;
int iStation;
- double x, y, time, dx, dy, dt, dxy;
+ double x;
+ double y;
+ double time;
+ double dx;
+ double dy;
+ double dt;
+ double dxy;
int Det;
};
+
+/// Enumeration for the detector subsystems used in L1 tracking
+/// Note: L1DetectorID has a forward declaration in L1InitManager.h and L1BaseStationInfo.h
+enum class L1DetectorID { kMvd, kSts, kMuch, kTof, kTrd };
+
+
// TODO: insert documentation!
//
/// L1Algo runtime constants modification can be performed in run_reco.C. Example:
@@ -124,6 +136,13 @@ public:
L1Parameters* GetL1Parameters() { return &fL1Parameters; }
+ /// Gets a set of active detectors used in tracking
+ // TODO: think about return (value, reference or const reference?)
+ std::set<L1DetectorID> GetActiveTrackingDetectorIDs() const { return fActiveTrackingDetectorIDs; }
+ /// Sets a vector of detectors used in tracking
+ void SetActiveTrackingDetectorIDs(const std::set<L1DetectorID>& init) { fActiveTrackingDetectorIDs = init; }
+
+
void SetStsMaterialBudgetFileName(TString fileName) { fStsMatBudgetFileName = fileName; }
void SetMvdMaterialBudgetFileName(TString fileName) { fMvdMatBudgetFileName = fileName; }
void SetMuchMaterialBudgetFileName(TString fileName) { fMuchMatBudgetFileName = fileName; }
@@ -224,6 +243,11 @@ private:
static CbmL1* fInstance;
L1Parameters fL1Parameters;
+ std::set<L1DetectorID> fActiveTrackingDetectorIDs {
+ L1DetectorID::kMvd,
+ L1DetectorID::kSts
+ }; ///< Set of detectors active in tracking
+
L1AlgoInputData* fData {nullptr};
L1Vector<CbmL1MCPoint> vMCPoints {"CbmL1::vMCPoints"};
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 3c2c4f3c15..55ee25a1e2 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -225,6 +225,42 @@ void L1Algo::Init(const L1Vector<fscal>& geo, const bool UseHitErrors, const Tra
#ifndef TBB2
std::cout << "L1Algo initialized" << std::endl;
#endif // TBB2
+
+ //
+ // NEW INITIALIZATION (BETA)
+ //
+
+
+ LOG(debug) << "\\";
+ LOG(debug) << " \\";
+ LOG(debug) << " Original L1Station vector content";
+ LOG(debug) << " /";
+ LOG(debug) << "/";
+ int nStations = fL1InitManager.GetStationsNumber();
+ for (const auto& aStation: vStations) {
+ int idx = &aStation - vStations;
+ if (idx == nStations) {
+ break;
+ }
+ LOG(debug) << "Station Global No: " << idx;
+ aStation.Print();
+ }
+ LOG(debug) << "\\";
+ LOG(debug) << " \\";
+ LOG(debug) << " New L1Station vector content";
+ LOG(debug) << " /";
+ LOG(debug) << "/";
+ fL1InitManager.TransferL1StationArray(fStationsNew);
+ nStations = fL1InitManager.GetStationsNumber();
+ for (const auto& aStation: fStationsNew) {
+ int idx = &aStation - &fStationsNew.front();
+ if (idx == nStations) {
+ break;
+ }
+ LOG(debug) << "Station Global No: " << idx;
+ aStation.Print();
+
+ }
}
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index 9b4685bb30..1bfbf4c89a 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -42,6 +42,7 @@ class L1AlgoDraw;
#include <iostream>
#include <limits>
#include <map>
+#include <array>
#include "L1Branch.h"
#include "L1Field.h"
@@ -125,14 +126,14 @@ public:
static unsigned int TripletId2Thread(unsigned int ID)
{
constexpr unsigned int kMoveThread = L1Parameters::kTripletBits;
- constexpr unsigned int kThreadMask = (1 << L1Parameters::kThreadBits) - 1;
+ constexpr unsigned int kThreadMask = (1u << L1Parameters::kThreadBits) - 1u;
return (ID >> kMoveThread) & kThreadMask;
}
/// unpack the triplet ID to its triplet index
static unsigned int TripletId2Triplet(unsigned int ID)
{
- constexpr unsigned int kTripletMask = (1 << L1Parameters::kTripletBits) - 1;
+ constexpr unsigned int kTripletMask = (1u << L1Parameters::kTripletBits) - 1u;
return ID & kTripletMask;
}
@@ -206,7 +207,10 @@ public:
int NStsStations {0}; ///< number of sts stations
int fNfieldStations {0}; ///< number of stations in the field region
+
+ // TODO: Replace _fvecalignment with C++11 alignas(16) attibute, see vStationsNew
L1Station vStations[L1Parameters::kMaxNstations] _fvecalignment; // station info
+ alignas(16) std::array<L1Station, L1Parameters::kMaxNstations> fStationsNew;
L1Vector<L1Material> fRadThick {"fRadThick"}; // material for each station
int NStsStrips {0}; // number of strips
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.cxx b/reco/L1/L1Algo/L1BaseStationInfo.cxx
new file mode 100644
index 0000000000..5beae9b526
--- /dev/null
+++ b/reco/L1/L1Algo/L1BaseStationInfo.cxx
@@ -0,0 +1,482 @@
+/* Copyright (C) 2016-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+/************************************************************************************************************
+ * @file L1BaseStationInfo.cxx
+ * @bried Realization of L1BaseStationInfo class large methods
+ * @since 18.01.2022
+ *
+ * The class is implemented to connect concrete experimental setup (CBM or BMN in CbmL1
+ * or BmnL1) with general L1Tracking algorithm. Each derived class must contain general
+ * algorithms sutable for the particular station type.
+ *
+ ***********************************************************************************************************/
+
+// FairRoot
+#include <FairLogger.h>
+
+// L1Algo core
+#include "L1BaseStationInfo.h"
+#include "L1Parameters.h"
+#include "L1Def.h"
+
+// C++ STL
+#include <utility>
+#include <iomanip>
+
+//
+// CONSTRUCTORS AND DESTRUCTOR
+//
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+L1BaseStationInfo::L1BaseStationInfo() noexcept
+{
+ LOG(DEBUG) << "L1BaseStationInfo: Default constructor called for " << this << '\n'; // Temporary
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+L1BaseStationInfo::L1BaseStationInfo(L1DetectorID detectorID, int stationID) noexcept
+ : fDetectorID(detectorID)
+ , fStationID(stationID)
+{
+ LOG(DEBUG) << "L1BaseStationInfo: Constructor (detectorID, stationID) called for " << this << '\n'; // Temporary
+ fInitFlags[kEDetectorID] = true;
+ fInitFlags[kEStationID] = true;
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+L1BaseStationInfo::~L1BaseStationInfo() noexcept
+{
+ LOG(DEBUG) << "L1BaseStationInfo: Destructor called for " << this << '\n'; // Temporary
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+L1BaseStationInfo::L1BaseStationInfo(const L1BaseStationInfo& other) noexcept
+ : fDetectorID(other.fDetectorID)
+ , fStationID(other.fStationID)
+ , fXmax(other.fXmax)
+ , fYmax(other.fYmax)
+ , fZPos(other.fZPos)
+ , fL1Station(other.fL1Station)
+ , fInitFlags(other.fInitFlags)
+{
+ LOG(debug) << "L1BaseStationInfo: Copy constructor called: " << &other << " was copied into " << this;
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+L1BaseStationInfo::L1BaseStationInfo(L1BaseStationInfo&& other) noexcept
+{
+ LOG(debug) << "L1BaseStationInfo: Move constructor called: " << &other << " was moved into " << this;
+ this->Swap(other);
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+L1BaseStationInfo& L1BaseStationInfo::operator=(const L1BaseStationInfo& other) noexcept
+{
+ LOG(debug) << "L1BaseStationInfo: Copy operator= called for " << &other << " was copied into" << this;
+ if (this != &other) {
+ L1BaseStationInfo(other).Swap(*this);
+ }
+ return *this;
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+L1BaseStationInfo& L1BaseStationInfo::operator=(L1BaseStationInfo&& other) noexcept
+{
+ LOG(debug) << "L1BaseStationInfo: Move operator= called for " << &other << " was copied into" << this;
+ if (this != &other) {
+ L1BaseStationInfo tmp(std::move(other));
+ this->Swap(tmp);
+ }
+ return *this;
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::Swap(L1BaseStationInfo& other) noexcept
+{
+ std::swap(fDetectorID, other.fDetectorID);
+ std::swap(fStationID, other.fStationID);
+ std::swap(fXmax, other.fXmax);
+ std::swap(fYmax, other.fYmax);
+ std::swap(fZPos, other.fZPos);
+ std::swap(fL1Station, other.fL1Station);
+ std::swap(fInitFlags, other.fInitFlags);
+}
+
+//
+// BASIC METHODS
+//
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::Print(int verbosity) const
+{
+ if (verbosity == 0) {
+ LOG(info)
+ << "L1BaseStationInfo object: {stationID, detectorID, address} = {" << fStationID << ", "
+ << static_cast<int>(fDetectorID) << ", " << this << '}';
+ }
+ else if(verbosity > 0) {
+ LOG(info) << "L1BaseStationInfo object: at " << this;
+ LOG(info) << "- Station ID: " << fStationID;
+ LOG(info) << "- Detector ID: " << static_cast<int>(fDetectorID);
+ LOG(info) << "- L1Station fields:";
+ LOG(info) << "--- Station type ID: " << fL1Station.type;
+ LOG(info) << "--- Time info ID: " << fL1Station.timeInfo;
+ LOG(info) << "--- z position: " << fL1Station.z[0];
+ LOG(info) << "--- Rmin: " << fL1Station.Rmin[0];
+ LOG(info) << "--- Rmax: " << fL1Station.Rmax[0];
+ LOG(info) << "--- Thickness (X), cm: " << fL1Station.materialInfo.thick[0];
+ LOG(info) << "--- Radiational length (X0), cm: " << fL1Station.materialInfo.RL[0];
+ if (verbosity > 1) {
+ 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::kMaxNFieldApproxCoefficients; ++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) << "- Additional fields:";
+ LOG(info) << "--- Xmax: " << fXmax;
+ LOG(info) << "--- Ymax: " << fYmax;
+ }
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::Reset()
+{
+ L1BaseStationInfo other;
+ std::swap(*this, other);
+}
+
+
+//
+// GETTERS
+//
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+const L1Station& L1BaseStationInfo::GetL1Station() const
+{
+ bool isStationInitialized = IsInitialized();
+ if (!isStationInitialized) {
+ LOG(error)
+ << "L1BaseStationInfo::GetL1Station: attempt to get an L1Staion object from uninitialized L1BaseStation with "
+ << "stationID = "<< fStationID << " and detectorID = " << static_cast<int>(fDetectorID);
+ assert((!isStationInitialized));
+ }
+ return fL1Station;
+}
+
+
+//
+// SETTERS
+//
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::SetStationID(int inID)
+{
+ if (! fInitFlags[kEStationID]) {
+ fStationID = inID;
+ fInitFlags[kEStationID] = true;
+ }
+ else {
+ LOG(warn) << "L1BaseStationInfo::SetStationID: Attempt of station ID redifinition";
+ }
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::SetDetectorID(L1DetectorID inID)
+{
+ if (! fInitFlags[kEDetectorID]) {
+ fDetectorID = inID;
+ fInitFlags[kEDetectorID] = true;
+ }
+ else {
+ LOG(warn) << "L1BaseStationInfo::SetDetectorID: Attempt of detector ID redifinition";
+ }
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::SetStationType(int inType)
+{
+ if (! fInitFlags[kEtype]) {
+ fL1Station.type = inType;
+ fInitFlags[kEtype] = true;
+ }
+ else {
+ LOG(warn) << "L1BaseStationInfo::SetStationType: Attempt of station type redifinition";
+ }
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::SetTimeInfo(int inTimeInfo)
+{
+ fL1Station.timeInfo = inTimeInfo;
+ fInitFlags[kEtimeInfo] = true;
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::SetZ(double inZ)
+{
+ fL1Station.z = inZ; // setting simd vector of single-precision floats, which is passed to high performanced L1Algo
+ fZPos = inZ; // setting precised value to use in field approximation etc
+ fInitFlags[kEz] = true;
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::SetRmin(double inRmin)
+{
+ fL1Station.Rmin = inRmin;
+ fInitFlags[kERmin] = true;
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::SetRmax(double inRmax)
+{
+ fL1Station.Rmax = inRmax;
+ fInitFlags[kERmax] = true;
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::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
+ fL1Station.materialInfo.thick = inThickness;
+ fL1Station.materialInfo.RL = inRL;
+ fL1Station.materialInfo.RadThick = fL1Station.materialInfo.thick / fL1Station.materialInfo.RL;
+ fL1Station.materialInfo.logRadThick = log(fL1Station.materialInfo.RadThick);
+ fInitFlags[kEmaterialInfoThick] = true;
+ fInitFlags[kEmaterialInfoRL] = true;
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::SetFieldSlice(const double* Cx, const double* Cy, const double* Cz)
+{
+ if (fInitFlags[kEfieldSlice]) {
+ LOG(warn)
+ << "L1BaseStationInfo::SetFieldSlice: Attempt to redifine field slice for station with detectorID = "
+ << static_cast<int>(fDetectorID) << " and stationID = " << fStationID << ". Redifinition ignored";
+ return;
+ }
+
+ for (int idx = 0; idx < L1Parameters::kMaxNFieldApproxCoefficients; ++idx) {
+ fL1Station.fieldSlice.cx[idx] = Cx[idx];
+ fL1Station.fieldSlice.cy[idx] = Cy[idx];
+ fL1Station.fieldSlice.cz[idx] = Cz[idx];
+ }
+
+ fInitFlags[kEfieldSlice] = true;
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::SetFieldSlice(const std::function<void(const double (&xyz)[3], double (&B)[3])>& getFieldValue)
+{
+ if (fInitFlags[kEfieldSlice]) {
+ LOG(warn)
+ << "L1BaseStationInfo::SetFieldSlice: Attempt to redifine field slice for station with detectorID = "
+ << static_cast<int>(fDetectorID) << " and stationID = " << fStationID << ". Redifinition ignored";
+ return;
+ }
+
+#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;
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::SetFrontBackStripsGeometry(double frontPhi, double frontSigma, double backPhi, double backSigma)
+{
+ // 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 cFront = cos(frontPhi);
+ double sFront = sin(frontPhi);
+ double cBack = cos(backPhi);
+ double sBack = sin(backPhi);
+
+ // NOTE: Here additional double variables are used to save the precission
+
+ fL1Station.frontInfo.cos_phi = cFront;
+ fL1Station.frontInfo.sin_phi = sFront;
+ fL1Station.frontInfo.sigma2 = frontSigma * frontSigma;
+
+ fL1Station.backInfo.cos_phi = cBack;
+ fL1Station.backInfo.sin_phi = sBack;
+ fL1Station.backInfo.sigma2 = backSigma * backSigma;
+
+ //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(backPhi - frontPhi) < 0.1) ? cos(backPhi - frontPhi) : (cFront * sBack - sFront * cBack);
+ det *= det;
+ fL1Station.XYInfo.C00 = (sBack * sBack * frontSigma * frontSigma + sFront * sFront * backSigma * backSigma) / det;
+ fL1Station.XYInfo.C10 = -(sBack * cBack * frontSigma * frontSigma + sFront * cFront * backSigma * backSigma) / det;
+ fL1Station.XYInfo.C11 = (cBack * cBack * frontSigma * frontSigma + cFront * cFront * backSigma * backSigma) / det;
+
+ fL1Station.xInfo.cos_phi = -sFront / (cFront * sBack - cBack * sFront);
+ fL1Station.xInfo.sin_phi = sBack / (cFront * sBack - cBack * sFront);
+ fL1Station.xInfo.sigma2 = fL1Station.XYInfo.C00;
+
+ fL1Station.yInfo.cos_phi = cBack / (cBack * sFront - cFront * sBack);
+ fL1Station.yInfo.sin_phi = -cFront / (cBack * sFront - cFront * sBack);
+ fL1Station.yInfo.sigma2 = fL1Station.XYInfo.C11;
+ //-----------------------------------------------------------------------------------------------------//
+
+ fInitFlags[kEstripsFrontPhi] = true;
+ fInitFlags[kEstripsFrontSigma] = true;
+ fInitFlags[kEstripsBackPhi] = true;
+ fInitFlags[kEstripsBackSigma] = true;
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::SetXmax(double aSize)
+{
+ fXmax = aSize;
+ fInitFlags[kEXmax] = true;
+}
+
+//----------------------------------------------------------------------------------------------------------------------//
+//
+void L1BaseStationInfo::SetYmax(double aSize)
+{
+ fYmax = aSize;
+ fInitFlags[kEYmax] = true;
+}
+
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.h b/reco/L1/L1Algo/L1BaseStationInfo.h
index f298ea9217..f9cc5a4b01 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.h
+++ b/reco/L1/L1Algo/L1BaseStationInfo.h
@@ -8,40 +8,28 @@
* @since 18.12.2021
*
* The class is implemented to connect concrete experimental setup (CBM or BMN in CbmL1
- * or BmnL1) with general L1Tracking algorithm. Each derived class must contain general
- * algorithms sutable for the particular station type.
+ * or BmnL1) with general L1Algo tracking.
*
***********************************************************************************************************/
#ifndef L1BaseStationInfo_h
#define L1BaseStationInfo_h 1
-// L1 core
-#include "L1Def.h"
-#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
+#include "L1Station.h"
+enum class L1DetectorID;
+
+/// A base class which provides interface to L1Algo station geometry
class L1BaseStationInfo {
private:
enum
{ // Here we list all the fields, which must be initialized by user
// Basic fields initialization
+ kEDetectorID,
kEStationID,
kEXmax,
kEYmax,
@@ -63,277 +51,54 @@ private:
};
public:
- L1BaseStationInfo()
- {
- std::cout << ">>>>>> L1BaseStationInfo: Constructor called\n"; // Temporary
- }
- virtual ~L1BaseStationInfo()
- {
- std::cout << ">>>>>> L1BaseStationInfo: Destructor called\n"; // Temporary
- }
-
- //-------------------------------------------------------------------------------------------------------//
- // Basic methods //
- //-------------------------------------------------------------------------------------------------------//
-
- /// Returns the name of the station (TODO: connect official naming)
- std::string GetTypeName() const { return fTypeName; }
- /// Checks if all the necessary fields are initialized by user
- bool IsInitialized() const { return fInitFlags.size() == fInitFlags.count(); }
- /// Transfers all gathered data to L1Algo (TODO)
- void TransferL1Station()
- { /**********/
- }
- void TransferData()
- { /*********/
- }
-
- //-------------------------------------------------------------------------------------------------------//
- // Interface for L1Station object initialization //
- //-------------------------------------------------------------------------------------------------------//
+ /// Default constructor
+ L1BaseStationInfo() noexcept;
+ /// Constructor from stationID and detectorID
+ L1BaseStationInfo(L1DetectorID detetorID, int stationID) noexcept;
+ /// Destructor
+ ~L1BaseStationInfo() noexcept;
+ /// Copy constructor
+ L1BaseStationInfo(const L1BaseStationInfo& other) noexcept;
+ /// Move constructor
+ L1BaseStationInfo(L1BaseStationInfo&& other) noexcept;
+ /// Copy operator=
+ L1BaseStationInfo& operator=(const L1BaseStationInfo& other) noexcept;
+ /// Move operator=
+ L1BaseStationInfo& operator=(L1BaseStationInfo&& other) noexcept;
+ /// Swap method for easier implementation of above ones (NOTE: all the fields must be accounted carefully here)
+ void Swap(L1BaseStationInfo& other) noexcept;
//
- // Setters
+ // BASIC METHODS
//
- /// 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)
- {
- 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)
- {
- fL1Station.timeInfo = inTimeInfo;
- fInitFlags[kEtimeInfo] = true;
- }
- // ???
-
- /// Sets nominal z position of the station
- void SetZ(double 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
- void SetRmin(double inRmin)
- {
- fL1Station.Rmin = inRmin;
- fInitFlags[kERmin] = true;
- }
- /// Sets max transverse size of the station
- void SetRmax(double inRmax)
- {
- fL1Station.Rmax = inRmax;
- fInitFlags[kERmax] = true;
- }
-
- /// Sets station thickness and radiational length
- 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
- fL1Station.materialInfo.thick = inThickness;
- fL1Station.materialInfo.RL = inRL;
- fL1Station.materialInfo.RadThick = fL1Station.materialInfo.thick / fL1Station.materialInfo.RL;
- fL1Station.materialInfo.logRadThick = log(fL1Station.materialInfo.RadThick);
- fInitFlags[kEmaterialInfoThick] = true;
- fInitFlags[kEmaterialInfoRL] = true;
- }
-
- /// Sets arrays of the approximation coefficients for the field x, y and z-components, respectively
- void SetFieldSlice(const double* Cx, const double* Cy, const double* Cz)
- {
- for (int idx = 0; idx < L1Parameters::kMaxNFieldApproxCoefficients; ++idx) {
- fL1Station.fieldSlice.cx[idx] = Cx[idx];
- fL1Station.fieldSlice.cy[idx] = Cy[idx];
- fL1Station.fieldSlice.cz[idx] = Cz[idx];
- }
- 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;
+ /// Checks if all the necessary fields are initialized by user
+ bool IsInitialized() const { return fInitFlags.size() == fInitFlags.count(); }
+ /// Prints registered fields
+ /// verbosity = 0: print only station id, detector id and address in one line
+ /// verbosity = 1: print basic L1Station fields
+ /// verbosity = 2: print all L1Staion fields
+ void Print(int verbosity = 0) const;
+ /// Print initialization flags (debug function)
+ void PrintInit() const;
+ /// Resets fields to the default values
+ void Reset();
+ /// Less operator for L1BaseStationInfo object to perform their sorts. Sorting is carried out first by fDetectorID,
+ /// and then by fStationID
+ bool operator<(const L1BaseStationInfo& right) const
+ {
+ return fDetectorID != right.fDetectorID ? fDetectorID < right.fDetectorID : fStationID < right.fStationID;
}
- /// Sets stereo angles and sigmas for front and back strips, automatically set frontInfo, backInfo,
- /// XYInfo, xInfo and yInfo
- void SetFrontBackStripsGeometry(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?
- // TODO: move into .cxx file
-
- //----- 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) {
- // 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;
- }
- /// 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
+ // GETTERS
//
/// Gets station ID
int GetStationID() const { return fStationID; }
+ /// Gets detector ID
+ L1DetectorID GetDetectorID() const { return fDetectorID; }
/// Gets station type
int GetStationType() const { return fL1Station.type; }
/// Gets SIMD vectorized z position of the station
@@ -346,11 +111,11 @@ public:
fvec GetRmax() const { return fL1Station.Rmax; }
/// Gets material thickness
- fvec GetMaterialThick() const { return fL1Station.materialInfo.thick; }
+ fvec GetMaterialThickness() const { return fL1Station.materialInfo.thick; }
// TODO: investigate if thick, RL and RadThick are useful, may be we should have only
// GetMateralLogRadThick?
/// Gets the radiation length of the station material
- fvec GetMaterialRL() const { return fL1Station.materialInfo.RL; } // TODO: -//-
+ fvec GetMaterialRadLength() const { return fL1Station.materialInfo.RL; } // TODO: -//-
/// Gets the relative material thickness in units of the radiational length
fvec GetMaterialRadThick() const { return fL1Station.materialInfo.RadThick; } // TODO: Rename?
/// Gets log of the relative material thickness in units of the radiational length
@@ -375,81 +140,68 @@ public:
/// Gets maximum distance between station center and its edge in y direction
double GetYmax() const { return fYmax; }
+ /// Gets a reference to L1Station info field of the L1BaseStation info
+ const L1Station& GetL1Station() const;
+ /// Gets a reference to Bitset object of initialization bits
+ const std::bitset<L1BaseStationInfo::kEND>& GetInitFlags() const { return fInitFlags; }
- //-------------------------------------------------------------------------------------------------------//
- // Other methods //
- //-------------------------------------------------------------------------------------------------------//
+
- /// Prints registered fields (TODO: tmp, may be one needs to wrap it into debug directives)
- void Print()
- {
- 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];
- 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::kMaxNFieldApproxCoefficients; ++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) << " Additiona fields:";
- LOG(INFO) << " Xmax: " << fXmax;
- LOG(INFO) << " Ymax: " << fYmax;
- LOG(INFO) << "";
- LOG(INFO) << "===================================================================================";
- }
- void Reset() {
- L1BaseStationInfo other;
- std::swap(*this, other);
- }
+ //
+ // SETTERS
+ //
-private:
-
- 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
+ /// Sets station ID
+ void SetStationID(int inID);
+ /// Sets detector ID
+ void SetDetectorID(L1DetectorID inID );
+
+ /// Sets type of station
+ void SetStationType(int inType); // TODO: this is a temporary solution since
+ /// Sets flag: 0 - time information is not provided by this detector type
+ // 1 - time information is provided by the detector and can be used in tracking
+ void SetTimeInfo(int inTimeInfo);
+ /// Sets nominal z position of the station
+ void SetZ(double inZ);
+ /// Sets min transverse size of the station
+ void SetRmin(double inRmin);
+ /// Sets max transverse size of the station
+ void SetRmax(double inRmax);
+ /// Sets station thickness and radiational length
+ /// \param thickness thickness of station
+ /// \param radiationalLength radiational length of station
+ void SetMaterial(double thickness, double radiationalLength);
+ /// Sets field approximation at the station plain
+ /// \param Cx Array of approximation coefficients for x field component
+ /// \param Cy Array of approximation coefficients for y field component
+ /// \param Cz Array of approximation coefficients for z field component
+ void SetFieldSlice(const double* Cx, const double* Cy, const double* Cz);
+ /// 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);
+ /// Sets stereo angles and sigmas for front and back strips
+ /// \param f_phi Stereoangle of front strips
+ /// \param f_sigma Sigma of front strips
+ /// \param b_phi Stereoangle of back strips
+ /// \param b_sigma Sigma of back strips
+ void SetFrontBackStripsGeometry(double fPhi, double fSigma, double bPhi, double bSigma);
+ /// Sets maximum distance between station center and its edge in x direction
+ void SetXmax(double aSize);
+ /// Sets maximum distance between station center and its edge in y direction
+ void SetYmax(double aSize);
+
+private:
+ L1DetectorID fDetectorID {static_cast<L1DetectorID>(0)}; ///< Detector ID
+ 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
- // be transfered?
+ L1Station fL1Station {}; ///< L1Station structure, describes a station in L1Algo
+ std::bitset<L1BaseStationInfo::kEND> fInitFlags; ///< Class fileds initialization flags
};
+/// swap function for two L1BaseStationInfo objects, expected to be used instead of std::swap
+inline void swap(L1BaseStationInfo& a, L1BaseStationInfo& b) noexcept { a.Swap(b); }
#endif // L1BaseStationInfo_h
diff --git a/reco/L1/L1Algo/L1CAIteration.cxx b/reco/L1/L1Algo/L1CAIteration.cxx
new file mode 100644
index 0000000000..33a5acc536
--- /dev/null
+++ b/reco/L1/L1Algo/L1CAIteration.cxx
@@ -0,0 +1,25 @@
+/* Copyright (C) 2016-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+#include <FairLogger.h>
+#include "L1CAIteration.h"
+
+L1CAIteration::L1CAIteration(const std::string& name, L1CAIterationType type)
+: fName(name),
+ fType(type)
+{
+ this->Print();
+}
+
+void L1CAIteration::Print() const
+{
+ LOG(INFO) << "CA Track Finder Iteration: " << fName;
+ LOG(INFO) << "\t\t-type" << static_cast<int>(fType);
+ LOG(INFO) << "\tTrack cuts:";
+ LOG(INFO) << "\t\t- chi2 cut" << fTrackChi2Cut;
+ LOG(INFO) << "\tTriplet cuts:";
+ LOG(INFO) << "\t\t- chi2 cut" << fTripletChi2Cut;
+ LOG(INFO) << "\tDoublet cuts:";
+ LOG(INFO) << "\t\t- chi2 cut" << fDoubletChi2Cut;
+}
diff --git a/reco/L1/L1Algo/L1CAIteration.h b/reco/L1/L1Algo/L1CAIteration.h
new file mode 100644
index 0000000000..f75ae3bbd2
--- /dev/null
+++ b/reco/L1/L1Algo/L1CAIteration.h
@@ -0,0 +1,121 @@
+/* Copyright (C) 2016-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+
+/************************************************************************************************************
+ * @file L1CAIteration.h
+ * @brief Class for L1 CA Track Finder Iteration
+ * @since 01.15.2022
+ *
+ ***********************************************************************************************************/
+
+#ifndef L1CAIteration_h
+#define L1CAIteration_h 1
+
+/// All possible iteration types of Track Finder algorithm to select
+///
+enum class L1CAIterationType {
+ 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
+};
+
+
+/// Class describes L1 Track finder iteration
+///
+class L1CAIteration {
+public:
+ L1CAIteration() = default;
+ ~L1CAIteration() = default;
+ L1CAIteration(const L1CAIteration& /*other*/) = default;
+ L1CAIteration(L1CAIteration&& /*other*/) = default;
+ L1CAIteration& operator=(const L1CAIteration& /*other*/) = default;
+ L1CAIteration& operator=(L1CAIteration&& /*other*/) = default;
+
+ L1CAIteration(const std::string& /*name*/, L1CAIterationType /*type*/);
+
+ //
+ // BASIC METHODS
+ //
+
+ /// Prints iteration options
+ void Print() const;
+
+ //
+ // SETTERS
+ //
+
+ //-------------------------------------------------------------------------------------------------------------------
+ // Basic fields setters
+ //-------------------------------------------------------------------------------------------------------------------
+
+ /// Sets name for this iteration, the name should be unique
+ // TODO: develope naming system and checkers
+ void SetName(const std::string& name) { fName = name; }
+ /// Sets type of the iteration
+ void SetType(L1CAIterationType type) { fType = type; }
+
+ //-------------------------------------------------------------------------------------------------------------------
+ // Cuts setters
+ //-------------------------------------------------------------------------------------------------------------------
+
+ /// Sets track chi2 upper cut
+ void SetTrackChi2Cut(float input) { fTrackChi2Cut = input; }
+ /// Sets triplet chi2 upper cut
+ void SetTripletChi2Cut(float input) { fTripletChi2Cut = input; }
+ /// Sets doublet chi2 upper cut
+ void SetDoubletChi2Cut(float input) { fDoubletChi2Cut = input; }
+
+
+ //
+ // GETTERS
+ //
+
+ //-------------------------------------------------------------------------------------------------------------------
+ // Basic fields getters
+ //-------------------------------------------------------------------------------------------------------------------
+
+ /// Gets the name of the iteration
+ std::string GetName() const { return fName; }
+ /// Gets the type of the iteration
+ L1CAIterationType GetType() { return fType; }
+
+ //-------------------------------------------------------------------------------------------------------------------
+ // Cuts getters
+ //-------------------------------------------------------------------------------------------------------------------
+
+ /// Gets track chi2 upper cut
+ float GetTrackChi2Cut() const { return fTrackChi2Cut; }
+ /// Gets triplet chi2 upper cut
+ float GetTripletChi2Cut() const { return fTripletChi2Cut; }
+ /// Gets doublet chi2 upper cut
+ float GetDoubletChi2Cut() const { return fDoubletChi2Cut; }
+
+private:
+
+ //-------------------------------------------------------------------------------------------------------------------
+ // Basic fields
+ //-------------------------------------------------------------------------------------------------------------------
+
+ std::string fName; ///< Iteration name
+ L1CAIterationType fType; ///< Iteration type
+
+ //-------------------------------------------------------------------------------------------------------------------
+ // Track finder dependent cuts
+ //-------------------------------------------------------------------------------------------------------------------
+
+ float fTrackChi2Cut; ///> track chi2 upper cut
+ float fTripletChi2Cut; ///> triplet chi2 upper cut
+ float fDoubletChi2Cut; ///> doublet chi2 upper cut
+
+
+
+};
+
+#endif // L1CAIteration_h
diff --git a/reco/L1/L1Algo/L1Def.h b/reco/L1/L1Algo/L1Def.h
index f6974adc5d..904817e62f 100644
--- a/reco/L1/L1Algo/L1Def.h
+++ b/reco/L1/L1Algo/L1Def.h
@@ -65,10 +65,7 @@ typedef int index_type;
/// Hash for unordered_map with enum class keys
struct EnumClassHash {
- template <typename T> size_t operator()(T t) const
- {
- return static_cast<size_t>(t);
- }
+ template <typename T> int operator()(T t) const { return static_cast<int>(t); }
};
#endif
diff --git a/reco/L1/L1Algo/L1InitManager.cxx b/reco/L1/L1Algo/L1InitManager.cxx
new file mode 100644
index 0000000000..590ffdd57a
--- /dev/null
+++ b/reco/L1/L1Algo/L1InitManager.cxx
@@ -0,0 +1,291 @@
+/* Copyright (C) 2016-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+/************************************************************************************************************
+ * @file L1InitManager.cxx
+ * @bried Input data management class for L1Algo
+ * @since 19.01.2022
+ *
+ ***********************************************************************************************************/
+
+#include "L1InitManager.h"
+#include <algorithm>
+
+
+//-----------------------------------------------------------------------------------------------------------------------
+//
+void L1InitManager::AddStation(const L1BaseStationInfo& inStation)
+{
+ // Check if other fields were defined already
+ // Active detector IDs
+ if (!fInitFlags[L1InitManager::kEactiveDetectorIDs]) {
+ LOG(error)
+ << "L1InitManager::AddStation: station initialization called before the active detectors set had been initialized";
+ assert((fInitFlags[L1InitManager::kEactiveDetectorIDs]));
+ }
+
+ // Number of stations check
+ if (!fInitFlags[L1InitManager::kEstationsNumberCrosscheck]) {
+ LOG(error)
+ << "L1InitManager::AddStation: station initialization called before the numbers of stations for each detector "
+ << "had been initialized";
+ assert((fInitFlags[L1InitManager::kEstationsNumberCrosscheck]));
+ }
+
+ // Field function
+ if (!fInitFlags[L1InitManager::kEfieldFunction]) {
+ LOG(error)
+ << "L1InitManager::AddStation: station initialization called before the magnetic field function was intialized";
+ assert((fInitFlags[L1InitManager::kEfieldFunction]));
+ }
+
+ // Check activeness of this station type
+ bool isDetectorActive = fActiveDetectorIDs.find(inStation.GetDetectorID()) != fActiveDetectorIDs.end();
+ if (isDetectorActive) {
+ // initialize magnetic field slice
+ L1BaseStationInfo inStationCopy = L1BaseStationInfo(inStation); // make a copy of station so it can be initialized
+ inStationCopy.SetFieldSlice(fFieldFunction);
+ bool isStationInitialized = inStationCopy.IsInitialized();
+ if (!isStationInitialized) {
+ LOG(debug) << "L1InitManager::AddStation: station init flags (original)" << inStation.GetInitFlags();
+ LOG(debug) << "L1InitManager::AddStation: station init flags (copy) " << inStation.GetInitFlags();
+ LOG(error)
+ << "L1InitManager::AddStation: Trying to add incompletely initialized object with detectorID = "
+ << static_cast<int>(inStationCopy.GetDetectorID()) << " and stationID = " << inStationCopy.GetStationID();
+ assert((isStationInitialized));
+ }
+ // insert the station in a set
+ auto insertionResult = fStationsInfo.insert(std::move(inStationCopy));
+ if (!insertionResult.second) {
+ LOG(error)
+ << "L1InitManager::AddStation: attempt to insert a dublicating L1BaseStationInfo with StationID = "
+ << inStation.GetStationID() << " and DetectorID = " << static_cast<int>(inStation.GetDetectorID()) << ":";
+ LOG(error) << ">>> Already inserted L1BaseStationInfo object:";
+ insertionResult.first->Print();
+ LOG(error) << ">>> A dublicating L1BaseStationInfo object:";
+ inStation.Print();
+ assert((insertionResult.second)); // TODO: rewrite the assertion
+ }
+ }
+ LOG(debug)
+ << "L1InitManager: adding a station with stationID = " << inStation.GetStationID() << " and detectorID = "
+ << static_cast<int>(inStation.GetDetectorID()) << ". Is active: " << isDetectorActive;
+
+
+}
+
+//-----------------------------------------------------------------------------------------------------------------------
+//
+void L1InitManager::Init() const
+{ // To be implemented
+ // Plans:
+ // 1. Must make a final check of the inititalization and turn on a corresponding trigger in L1Algo class to accept
+ // the incoming data
+}
+
+//-----------------------------------------------------------------------------------------------------------------------
+//
+void L1InitManager::PrintStations(int verbosityLevel) const
+{
+ if (verbosityLevel < 1) {
+ for (auto& station : fStationsInfo) {
+ LOG(info) << "----------- station: ";
+ LOG(info) << "\ttype = " << station.GetStationType(); // TMP
+ LOG(info) << "\tz = " << station.GetZdouble();
+ }
+ }
+ else {
+ for (auto& station : fStationsInfo) {
+ station.Print();
+ }
+ }
+}
+
+//-----------------------------------------------------------------------------------------------------------------------
+//
+void L1InitManager::TransferL1StationArray(std::array<L1Station, L1Parameters::kMaxNstations>& destinationArray)
+{
+ /// First of all, we must check if L1Station was properly initialized
+ // TODO: actually, false condition will never reached (must thing about it, may be remove assertions from
+ // CheckStationInfo and leave only warnings and flag)
+ bool ifStationsInitialized = CheckStationsInfo();
+ if (!ifStationsInitialized) {
+ LOG(error) << "L1InitManager::TransferL1StationArray: attempt to pass unitialized L1Station array to L1Algo core";
+ assert((ifStationsInitialized));
+ }
+
+ /// Check if destinationArraySize is enough for the transfer
+ int totalStationsNumber = this->GetStationsNumber();
+ bool ifDestinationArraySizeOk = totalStationsNumber <= static_cast<int>(destinationArray.size());
+ if (!ifDestinationArraySizeOk) {
+ LOG(error)
+ << "L1InitManager::TransferL1StationArray: destination array size (" << destinationArray.size()
+ << ") is smaller then actual number of active tracking stations (" << totalStationsNumber << ")";
+ assert((ifDestinationArraySizeOk));
+ }
+
+ auto destinationArrayIterator = destinationArray.begin();
+ for (const auto& item: fStationsInfo) {
+ *destinationArrayIterator = std::move(item.GetL1Station());
+ ++destinationArrayIterator;
+ }
+ LOG(info) << "L1InitManager: L1Station vector was successfully transfered to L1Algo core :)";
+}
+
+//
+// GETTERS
+//
+
+//-----------------------------------------------------------------------------------------------------------------------
+//
+int L1InitManager::GetStationsNumber(L1DetectorID detectorID) const
+{
+ auto ifDetectorIdDesired = [&detectorID](const L1BaseStationInfo& station) {
+ return station.GetDetectorID() == detectorID;
+ };
+ return std::count_if(fStationsInfo.begin(), fStationsInfo.end(), ifDetectorIdDesired);
+}
+
+//
+// SETTERS
+//
+
+//-----------------------------------------------------------------------------------------------------------------------
+//
+void L1InitManager::SetActiveDetectorIDs(const std::set<L1DetectorID>& detectorIDs)
+{
+ // TODO: To think about redifinition possibilities: should it be allowed or not?
+ fActiveDetectorIDs = detectorIDs;
+ fInitFlags[L1InitManager::kEactiveDetectorIDs] = true;
+}
+
+//-----------------------------------------------------------------------------------------------------------------------
+//
+void L1InitManager::SetFieldFunction(const std::function<void(const double (&xyz)[3], double (&B)[3])>& fieldFunction)
+{
+ if (!fInitFlags[L1InitManager::kEfieldFunction]) {
+ fFieldFunction = fieldFunction;
+ fInitFlags[L1InitManager::kEfieldFunction] = true;
+ }
+ else {
+ LOG(warn) << "L1InitManager::SetFieldFunction: attemt to reinitialize the field function. Ignored";
+ }
+}
+
+//-----------------------------------------------------------------------------------------------------------------------
+//
+void L1InitManager::SetStationsNumberCrosscheck(L1DetectorID detectorID, int nStations)
+{
+ // NOTE: We add and check only those detectors which will be active (?)
+ // For INACTIVE detectors the initialization code for it inside CbmL1/BmnL1 can (and must) be still in,
+ // but it will be ignored inside L1InitManager.
+ if (fActiveDetectorIDs.find(detectorID) != fActiveDetectorIDs.end()) {
+ fStationsNumberCrosscheck[detectorID] = nStations;
+ }
+
+ // Check if all the station numbers for active detectors are initialized now:
+ LOG(debug) << "SetStationsNumberCrosscheck called for detectorID = " << static_cast<int>(detectorID);
+ if (!fInitFlags[L1InitManager::kEstationsNumberCrosscheck]) {
+ bool ifInitialized = true;
+ for (auto item: fActiveDetectorIDs) {
+ if (fStationsNumberCrosscheck.find(item) == fStationsNumberCrosscheck.end()) {
+ LOG(warn)
+ << "L1InitManager::SetStationsNumberCrosscheck: uninitialized number of stations for detectorID = "
+ << static_cast<int>(item);
+ ifInitialized = false;
+ break;
+ }
+ }
+ fInitFlags[L1InitManager::kEstationsNumberCrosscheck] = ifInitialized;
+ }
+ LOG(debug) << "InitResult: " << fInitFlags[L1InitManager::kEstationsNumberCrosscheck];
+}
+
+//-----------------------------------------------------------------------------------------------------------------------
+//
+void L1InitManager::SetReferencePrimaryVertexPoints(double z0, double z1, double z2)
+{
+ if (fInitFlags[L1InitManager::kEprimaryVertexField]) {
+ LOG(warn)
+ << "L1InitManager::SetReferencePrimaryVertexPoints: attempt to redefine reference points for field calculation "
+ << "near primary vertex. Ignore";
+ return;
+ }
+
+ // Check for field function
+ if (!fInitFlags[L1InitManager::kEfieldFunction]) {
+ LOG(error)
+ << "L1InitManager::SetReferencePrimaryVertexPoints: attempt to set reference points for field calculation near "
+ << "primary vertex before the magnetic field function intialization";
+ assert((fInitFlags[L1InitManager::kEfieldFunction]));
+ }
+
+ constexpr int numberOfDimensions {3};
+ constexpr int numberOfReferencePoints {3};
+
+ std::array<double, numberOfReferencePoints> inputZ = { z0, z1, z2 }; // tmp array to store input assigned with index
+ std::array<L1FieldValue, numberOfReferencePoints> B = {};
+ std::array<fvec, numberOfReferencePoints> z = { z0, z1, z2 };
+ for (int idx = 0; idx < numberOfReferencePoints; ++idx) {
+ double point[numberOfDimensions] = {0., 0., inputZ[idx] };
+ double field[numberOfDimensions] = {};
+ fFieldFunction(point, field);
+ B[idx].x = field[0];
+ B[idx].y = field[1];
+ B[idx].z = field[2];
+ }
+ fPrimaryVertexFieldRegion.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
+ fPrimaryVertexFieldValue = B[0];
+
+ fInitFlags[L1InitManager::kEprimaryVertexField] = true;
+}
+
+
+
+//-----------------------------------------------------------------------------------------------------------------------
+//
+bool L1InitManager::CheckStationsInfo()
+{
+ if (!fInitFlags[L1InitManager::kEstationsInfo]) {
+ bool ifInitPassed = true;
+
+ if (!fInitFlags[L1InitManager::kEifStationNumbersChecked]) {
+ for (const auto& itemDetector: fActiveDetectorIDs) {
+ int actualStationsNumber = GetStationsNumber(itemDetector);
+ int expectedStationsNumber = fStationsNumberCrosscheck.at(itemDetector);
+ if (actualStationsNumber != expectedStationsNumber) {
+ LOG(error)
+ << "L1InitManager::CheckStationsInfo: Incorrect number of L1BaseStationInfo objects passed to the L1Manager "
+ << "for L1DetectorID = " << static_cast<int>(itemDetector) << ": " << actualStationsNumber << " of "
+ << expectedStationsNumber << " expected";
+ ifInitPassed = false;
+ }
+ }
+ fInitFlags[L1InitManager::kEifStationNumbersChecked] = ifInitPassed;
+ }
+ if (!ifInitPassed) {
+ LOG(error) << "L1InitManager::CheckStationsInfo: initialization failed";
+ assert((ifInitPassed));
+ }
+
+ // Check for maximum allowed number of stations
+ int totalStationsNumber = GetStationsNumber();
+ if (totalStationsNumber > L1Parameters::kMaxNstations) {
+ LOG(fatal)
+ << "Actual total number of registered stations (" << totalStationsNumber << ") is larger then designed one ("
+ << L1Parameters::kMaxNstations << "). Please, select another set of active tracking detectors";
+ // TODO: We have to provide an instruction of how to increase the kMaxNstations number keeping the code consistent
+ assert((totalStationsNumber <= L1Parameters::kMaxNstations));
+ }
+
+ fInitFlags[L1InitManager::kEstationsInfo] = true;
+ }
+ else {
+ LOG(warn) << "L1InitManager: L1BaseStationInfo set has already been initialized";
+ }
+ // NOTE: we return a flag here to reduce a number of calls outside the funcition. In other hands we keep this flag
+ // to be consistent with other class fields initialization rules
+ return fInitFlags[L1InitManager::kEstationsInfo];
+}
+
diff --git a/reco/L1/L1Algo/L1InitManager.h b/reco/L1/L1Algo/L1InitManager.h
index e234f1dac9..efc1e34833 100644
--- a/reco/L1/L1Algo/L1InitManager.h
+++ b/reco/L1/L1Algo/L1InitManager.h
@@ -8,20 +8,71 @@
* @since 24.12.2021
*
***********************************************************************************************************/
+#ifndef L1InitManager_h
+#define L1InitManager_h 1
#include "L1BaseStationInfo.h"
#include "L1Parameters.h"
-#include "L1Vector.h"
+#include "L1Field.h"
//#include <string>
-#include <unordered_map>
+#include <bitset>
+#include <set>
+#include <memory> //unique_ptr
+
+/// Forward declaration of the tracking detectors scoped enumeration. Concrete realization of this enumeration must be
+/// determined in the concrete setup class (i.e. CbmL1/BmnL1)
+enum class L1DetectorID;
/// Initialization manager for L1Aglo
///
-/// The L1InitManager class provides interface for L1Algo initialization
-
+/// ==== Expected initialization steps ==== (TODO: keep this instruction up-to-date)
+///
+/// 0. Implement your enum class L1DetectorID with detector systems used for tracking:
+/// IMPORTANT!!! Detectors must be sorted in the beam direction
+/// /// Concrete L1DetectorID implementation for CBM
+/// enum class L1DetectorID {
+/// kMvd,
+/// kSts,
+/// kMuch,
+/// kTof,
+/// kTrd
+/// };
+///
+/// 1. Get a pointer to the :L1InitManager field of L1Algo:
+///
+/// L1InitManager * initMan = algo->GetL1InitManager();
+///
+/// 2. Initialize a set of L1DetectorID's for detectors active in tracking:
+///
+/// std::set<L1DetectorID> activeTrackingDetectors { L1Detector::kMvd, L1Detector::kSts };
+/// initMan->SetActiveDetectorIDs(activeTrackingDetectors);
+///
+/// 3. Initialize number of stations for each detector:
+///
+/// initMan->SetStationsNumberCrosscheck(L1DetectorID::kMvd, NMvdStations)
+/// initMan->SetStationsNumberCrosscheck(L1DetectorID::kMvd, NStsStations);
+///
+/// 3. Initialize each station using L1BaseStationInfo:
+///
class L1InitManager {
+private:
+ enum {
+ kEactiveDetectorIDs, ///< If the detector sequence is set
+ kEstationsNumberCrosscheck, ///< If the crosscheck station numbers were setup
+ kEfieldFunction, ///< If magnetic field getter funciton is set
+ kEprimaryVertexField, ///< If magnetic field value and region defined at primary vertex
+ kEifStationNumbersChecked, ///< If the station number was already checked
+ kEstationsInfo, ///< If all the planned stations were added to the manager
+ kEL1StationTransfered, ///< If the L1Station vector was already transfered to destination array
+ kEND
+ };
public:
+
+ //
+ // CONSTRUCTORS AND DESTRUCTOR
+ //
+
/// Default constructor
L1InitManager() = default;
/// Destructor
@@ -35,40 +86,84 @@ public:
/// Move assignment operator is forbidden
L1InitManager& operator=(L1InitManager&& /*other*/) = delete;
+
+ //
+ // BASIC METHODS
+ //
+
+
/// 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
- }
+ void AddStation(const L1BaseStationInfo& inStation);
/// 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();
- }
- }
- }
+ void AddStation(const L1BaseStationInfo* inStationPtr) { AddStation(*inStationPtr); }
+ /// Adds another station of a given type using std::unique_ptr-wraped pointer to L1BaseStationInfo
+ // TODO: Test weaknesses here
+ void AddStation(const std::unique_ptr<L1BaseStationInfo>& inStationUniquePtr) { AddStation(*inStationUniquePtr); }
+ /// Initializes L1Algo: transfers all caputred data to the L1 tracking core
+ void Init() const;
+ /// Prints a list of stations
+ void PrintStations(int verbosityLevel = 0) const;
+
+ /// Transfer an array of L1Stations formed inside a set of L1BaseStationInfo to a destination std::array
+ void TransferL1StationArray(std::array<L1Station, L1Parameters::kMaxNstations>& destinationArray);
+
+ //
+ // GETTERS
+ //
+
+ /// Gets a set of actie detectors for this analysis
+ std::set<L1DetectorID> GetActiveDetectorIDs() const { return fActiveDetectorIDs; }
+ /// Gets a total number of stations (NOTE: this number includes both active and unactive stations!)
+ int GetStationsNumber() const { return static_cast<int>(fStationsInfo.size()); }
+ /// Gets a number of stations for a particualr detector ID
+ int GetStationsNumber(L1DetectorID detectorID) const;
+ /// Gets a L1FieldRegion object at primary vertex
+ const L1FieldRegion& GetPrimaryVertexFieldRegion() const { return fPrimaryVertexFieldRegion; }
+ /// Gets a L1FieldValue object at primary vertex
+ const L1FieldValue& GetPrimaryVertexFieldValue() const { return fPrimaryVertexFieldValue; }
+
+ //
+ // SETTERS
+ //
+
+ /// Sets a set of active tracking detector IDs
+ void SetActiveDetectorIDs(const std::set<L1DetectorID>& detectorIDs);
+ /// Sets a magnetic field function, which will be applied for all the stations
+ void SetFieldFunction(const std::function<void(const double (&xyz)[3], double (&B)[3])>& fieldFcn);
+ /// Sets a number of stations for a particular tracking detector ID to provide initialization cross-check
+ void SetStationsNumberCrosscheck(L1DetectorID detectorID, int nStations);
+ /// Sets z coordinates referecne points in a vicinity of primary vertex to calculate L1FieldValue and L1FieldReference
+ /// values
+ // TODO: Consider posibility for linear approximation
+ void SetReferencePrimaryVertexPoints(double z0, double z1, double z2);
+
+
private:
- // TODO: TEST IT!!!!!!!!!!!!!!!!!!!!!
- 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>
+ /// Checker for L1BaseStationInfo set
+ /// \return true if all L1BaseStationInfo objects were initialized properly. Similar effect can be achieved by
+ /// calling the fInitFlags[L1InitManager::kEstationsInfo] flag
+ bool CheckStationsInfo();
+
+ /* Basic fields */
+
+ std::bitset<L1InitManager::kEND> fInitFlags {}; ///< Initialization flags
+ std::set<L1DetectorID> fActiveDetectorIDs {}; ///< Set of tracking detectors, active during this analysis session
+
+ /* Stations related fields */
+
+ std::set<L1BaseStationInfo> fStationsInfo {}; ///< Set of L1BaseStationInfo objects
+ /// Map of station numbers used for initialization crosscheck
+ std::unordered_map<L1DetectorID, int, EnumClassHash> fStationsNumberCrosscheck {};
+ /// A function which returns magnetic field vector B in a radius-vector xyz
+ std::function<void(const double (&xyz)[3], double (&B)[3])> fFieldFunction {[](const double (&)[3], double (&)[3]){}};
+ // NOTE: Stations of daetectors which will not be assigned as active, will not be included in the tracking!!!!!!!
+ // NOTE: fTotalNumberOfStations is excess field for logic, but it's imortant to track L1Algo initialization
+
+ /* Vertex related fields */
+ L1FieldValue fPrimaryVertexFieldValue {}; ///> L1FieldValue object at primary vertex
+ L1FieldRegion fPrimaryVertexFieldRegion {}; ///> L1FieldRegion object at primary vertex
+
};
+
+#endif // L1InitManager_h
diff --git a/reco/L1/L1Algo/L1Parameters.h b/reco/L1/L1Algo/L1Parameters.h
index 59a3ea5e1e..3bc052549a 100644
--- a/reco/L1/L1Algo/L1Parameters.h
+++ b/reco/L1/L1Algo/L1Parameters.h
@@ -30,44 +30,46 @@
// TODO: Create independent class of L1CATrackFinderIterations ???
-/// All possible iteration species of Track Finder algorithm
-enum class L1CATrackFinderIter {
- 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
-};
-
-// Type for mapping L1CATrackFinder parameters, which are dependent from the current track finder iteration
-template <typename T>
-using L1CATrackFinderIterParameterMap_t = std::unordered_map<L1CATrackFinderIter, T, EnumClassHash>;
-
-const L1CATrackFinderIterParameterMap_t<std::string> kTrackFinderIterNames = {
- { L1CATrackFinderIter::kFastPrim, "Primary tracks, fast" },
- { L1CATrackFinderIter::kAllPrim, "Primary tracks, all" },
- { L1CATrackFinderIter::kAllPrimJump, "Primary tracks with jumped triplets, all" },
- { L1CATrackFinderIter::kAllSec, "Secondary tracks, all" },
- { L1CATrackFinderIter::kAllPrimE, "Primary electron tracks, all" },
- { L1CATrackFinderIter::kAllSecE, "Secondary electron tracks, all" },
- { L1CATrackFinderIter::kFastPrimJump, "Primary tracks with jumped triplets, fast" },
- { L1CATrackFinderIter::kAllSecJump, "Secondary tracks with jumped triplets, all" }
-};
-
-const L1CATrackFinderIterParameterMap_t<std::string> kTrackFinderIterNamesShort = {
- { L1CATrackFinderIter::kFastPrim, "FastPrim" },
- { L1CATrackFinderIter::kAllPrim, "AllPrim" },
- { L1CATrackFinderIter::kAllPrimJump, "AllPrimJump" },
- { L1CATrackFinderIter::kAllSec, "AllSec" },
- { L1CATrackFinderIter::kAllPrimE, "AllPrimE" },
- { L1CATrackFinderIter::kAllSecE, "AllSecE" },
- { L1CATrackFinderIter::kFastPrimJump, "FastPrimJump" },
- { L1CATrackFinderIter::kAllSecJump, "AllSecJump" }
-};
+/////// > Will be moved to a separate class
+/// All possible iteration species of Track Finder algorithm
+// enum class L1CATrackFinderIter {
+// 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
+// };
+//
+// // Type for mapping L1CATrackFinder parameters, which are dependent from the current track finder iteration
+// template <typename T>
+// using L1CATrackFinderIterParameterMap_t = std::unordered_map<L1CATrackFinderIter, T, EnumClassHash>;
+//
+// ////// > Will be moved to separate class
+// const L1CATrackFinderIterParameterMap_t<std::string> kTrackFinderIterNames = {
+// { L1CATrackFinderIter::kFastPrim, "Primary tracks, fast" },
+// { L1CATrackFinderIter::kAllPrim, "Primary tracks, all" },
+// { L1CATrackFinderIter::kAllPrimJump, "Primary tracks with jumped triplets, all" },
+// { L1CATrackFinderIter::kAllSec, "Secondary tracks, all" },
+// { L1CATrackFinderIter::kAllPrimE, "Primary electron tracks, all" },
+// { L1CATrackFinderIter::kAllSecE, "Secondary electron tracks, all" },
+// { L1CATrackFinderIter::kFastPrimJump, "Primary tracks with jumped triplets, fast" },
+// { L1CATrackFinderIter::kAllSecJump, "Secondary tracks with jumped triplets, all" }
+// };
+//
+// const L1CATrackFinderIterParameterMap_t<std::string> kTrackFinderIterNamesShort = {
+// { L1CATrackFinderIter::kFastPrim, "FastPrim" },
+// { L1CATrackFinderIter::kAllPrim, "AllPrim" },
+// { L1CATrackFinderIter::kAllPrimJump, "AllPrimJump" },
+// { L1CATrackFinderIter::kAllSec, "AllSec" },
+// { L1CATrackFinderIter::kAllPrimE, "AllPrimE" },
+// { L1CATrackFinderIter::kAllSecE, "AllSecE" },
+// { L1CATrackFinderIter::kFastPrimJump, "FastPrimJump" },
+// { L1CATrackFinderIter::kAllSecJump, "AllSecJump" }
+// };
class L1Parameters {
public:
@@ -79,15 +81,17 @@ public:
// order instead of this?
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 kStationBits {6u}; ///> Amount of bits to code one station
+ static constexpr unsigned int kThreadBits {6u}; ///> Amount of bits to code one thread
+ static constexpr unsigned int kTripletBits {32u - kStationBits - kThreadBits}; ///> Amount of bits to code one triple
+
+ static constexpr int kMaxNstations { 1u << kStationBits }; ///> Max number of stations, 2^6 = 64
+ static constexpr int kMaxNthreads { 1u << kThreadBits }; ///> Max number of threads, 2^6 = 64
+ static constexpr int kMaxNtriplets { 1u << kTripletBits }; ///> Max number of triplets, 2^20 = 1,048,576
- 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
+ static constexpr int kStandardIOWidth {15}; ///> Width of one output entry, passed to the std::setw()
- static constexpr short kStandardIOWidth = 15; ///> Width of one output entry, passed to the std::setw()
+ static constexpr int kAssertionLevel {0}; ///> Assertion level
public:
//-------------------------------------------------------------------------------------------------------//
@@ -141,73 +145,74 @@ public:
/// 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<L1CATrackFinderIter>& iterations) {
- fTrackFinderIterSequence = iterations;
- }
-
- /// Sets track Chi2 upper cut
- /// \param tfIteration Track Finder iteration of the L1 algorithm run
- /// \param chi2Cut Upper cut on track chi2 during selected iteration
- void SetTrackChi2Cut(L1CATrackFinderIter tfIteration, float chi2Cut) { fTrackChi2Cut[tfIteration] = chi2Cut; }
- void SetTrackChi2Cut(float chi2Cut) { SetSameValueToMap(chi2Cut, fTrackChi2Cut); }
- void SetTrackChi2Cut(const L1CATrackFinderIterParameterMap_t<float>& newValues)
- {
- SetMappedValuesToMap(newValues, fTrackChi2Cut);
- }
- /// Sets triplet Chi2 upper cut
- /// \param tfIteration Track Finder iteration of the L1 algorithm run
- /// \param chi2Cut Upper cut on triplet chi2 during selected iteration
- void SetTripletChi2Cut(L1CATrackFinderIter tfIteration, float chi2Cut) { fTripletChi2Cut[tfIteration] = chi2Cut; }
- void SetTripletChi2Cut(float chi2Cut) { SetSameValueToMap(chi2Cut, fTripletChi2Cut); }
- void SetTripletChi2Cut(const L1CATrackFinderIterParameterMap_t<float>& newValues)
- {
- SetMappedValuesToMap(newValues, fTripletChi2Cut);
- }
- /// Sets triplet Chi2 upper cut
- /// \param tfIteration Track Finder iteration of the L1 algorithm run
- /// \param chi2Cut Upper cut on triplet chi2 during selected iteration
- void SetDoubletChi2Cut(L1CATrackFinderIter tfIteration, float chi2Cut) { fDoubletChi2Cut[tfIteration] = chi2Cut; }
- void SetDoubletChi2Cut(float chi2Cut) { SetSameValueToMap(chi2Cut, fDoubletChi2Cut); }
- void SetDoubletChi2Cut(const L1CATrackFinderIterParameterMap_t<float>& newValues)
- {
- SetMappedValuesToMap(newValues, fDoubletChi2Cut);
- }
+ // /// Sets track finder iteration sequence
+ // void SetTrackFinderIterSequence(const L1Vector<L1CATrackFinderIter>& iterations) {
+ // fTrackFinderIterSequence = iterations;
+ // }
+ //
+ // /// Sets track Chi2 upper cut
+ // /// \param tfIteration Track Finder iteration of the L1 algorithm run
+ // /// \param chi2Cut Upper cut on track chi2 during selected iteration
+ // void SetTrackChi2Cut(L1CATrackFinderIter tfIteration, float chi2Cut) { fTrackChi2Cut[tfIteration] = chi2Cut; }
+ // void SetTrackChi2Cut(float chi2Cut) { SetSameValueToMap(chi2Cut, fTrackChi2Cut); }
+ // void SetTrackChi2Cut(const L1CATrackFinderIterParameterMap_t<float>& newValues)
+ // {
+ // SetMappedValuesToMap(newValues, fTrackChi2Cut);
+ // }
+ // /// Sets triplet Chi2 upper cut
+ // /// \param tfIteration Track Finder iteration of the L1 algorithm run
+ // /// \param chi2Cut Upper cut on triplet chi2 during selected iteration
+ // void SetTripletChi2Cut(L1CATrackFinderIter tfIteration, float chi2Cut) { fTripletChi2Cut[tfIteration] = chi2Cut; }
+ // void SetTripletChi2Cut(float chi2Cut) { SetSameValueToMap(chi2Cut, fTripletChi2Cut); }
+ // void SetTripletChi2Cut(const L1CATrackFinderIterParameterMap_t<float>& newValues)
+ // {
+ // SetMappedValuesToMap(newValues, fTripletChi2Cut);
+ // }
+ // /// Sets triplet Chi2 upper cut
+ // /// \param tfIteration Track Finder iteration of the L1 algorithm run
+ // /// \param chi2Cut Upper cut on triplet chi2 during selected iteration
+ // void SetDoubletChi2Cut(L1CATrackFinderIter tfIteration, float chi2Cut) { fDoubletChi2Cut[tfIteration] = chi2Cut; }
+ // void SetDoubletChi2Cut(float chi2Cut) { SetSameValueToMap(chi2Cut, fDoubletChi2Cut); }
+ // void SetDoubletChi2Cut(const L1CATrackFinderIterParameterMap_t<float>& newValues)
+ // {
+ // SetMappedValuesToMap(newValues, fDoubletChi2Cut);
+ // }
+ //
// ***** Getters ***** //
unsigned int GetMaxDoubletsPerSinglet() const { return fMaxDoubletsPerSinglet; }
unsigned int GetMaxTripletPerDoublets() const { return fMaxTripletPerDoublets; }
- /// Gets track Chi2 upper cut
- /// \param tfIteration Track Finder iteration of the L1 algorithm run
- /// \return Upper cut on track chi2 during selected iteration
- float GetTrackChi2Cut(L1CATrackFinderIter tfIteration) const { return fTrackChi2Cut.at(tfIteration); }
- /// Gets triplet Chi2 upper cut
- /// \param tfIteration Track Finder iteration of the L1 algorithm run
- /// \return Upper cut on triplet chi2 during selected iteration
- float GetTripletChi2Cut(L1CATrackFinderIter tfIteration) const { return fTripletChi2Cut.at(tfIteration); }
- /// Gets triplet Chi2 upper cut
- /// \param tfIteration Track Finder iteration of the L1 algorithm run
- /// \return Upper cut on triplet chi2 during selected iteration
- float GetDoubletChi2Cut(L1CATrackFinderIter tfIteration) const { return fDoubletChi2Cut.at(tfIteration); }
+ // /// Gets track Chi2 upper cut
+ // /// \param tfIteration Track Finder iteration of the L1 algorithm run
+ // /// \return Upper cut on track chi2 during selected iteration
+ // float GetTrackChi2Cut(L1CATrackFinderIter tfIteration) const { return fTrackChi2Cut.at(tfIteration); }
+ // /// Gets triplet Chi2 upper cut
+ // /// \param tfIteration Track Finder iteration of the L1 algorithm run
+ // /// \return Upper cut on triplet chi2 during selected iteration
+ // float GetTripletChi2Cut(L1CATrackFinderIter tfIteration) const { return fTripletChi2Cut.at(tfIteration); }
+ // /// Gets triplet Chi2 upper cut
+ // /// \param tfIteration Track Finder iteration of the L1 algorithm run
+ // /// \return Upper cut on triplet chi2 during selected iteration
+ // float GetDoubletChi2Cut(L1CATrackFinderIter tfIteration) const { return fDoubletChi2Cut.at(tfIteration); }
//-------------------------------------------------------------------------------------------------------//
// Auxilary methods //
//-------------------------------------------------------------------------------------------------------//
public:
- void PrintTrackFinderIterSequence() const
- {
- LOG(INFO) << "== L1Algo track finder iterations sequence ==========";
- 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) << "=====================================================";
- }
+ //void PrintTrackFinderIterSequence() const
+ //{
+ // LOG(INFO) << "== L1Algo track finder iterations sequence ==========";
+ // 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:
/// Aux method for printing mapped Track Finder parameters
@@ -224,45 +229,45 @@ private:
unsigned int fMaxDoubletsPerSinglet {150}; ///<
unsigned int fMaxTripletPerDoublets {15}; ///<
- L1Vector<L1CATrackFinderIter> fTrackFinderIterSequence {}; ///< Sequence of iterations in the track finder algorithm
+ // L1Vector<L1CATrackFinderIter> fTrackFinderIterSequence {}; ///< Sequence of iterations in the track finder algorithm
- // TODO: Develope naming system for these constants
- L1CATrackFinderIterParameterMap_t<float> fTrackChi2Cut {
- { L1CATrackFinderIter::kFastPrim, 10.f },
- { L1CATrackFinderIter::kAllPrim, 10.f },
- { L1CATrackFinderIter::kAllPrimJump, 10.f },
- { L1CATrackFinderIter::kAllSec, 10.f },
- { L1CATrackFinderIter::kAllPrimE, 10.f },
- { L1CATrackFinderIter::kAllSecE, 10.f },
- { L1CATrackFinderIter::kFastPrimJump, 10.f },
- { L1CATrackFinderIter::kAllSecJump, 10.f }
- };
- // TODO: What is the meaning of the following numbers?
- // 23.4450 <- 7.815 * 3
- // 18.7560 <- 6.252 * 3
- // 11.3449 * 2. / 3.
- //
- L1CATrackFinderIterParameterMap_t<float> fTripletChi2Cut {
- { L1CATrackFinderIter::kFastPrim, 23.4450f }, // kFastPrimIter
- { L1CATrackFinderIter::kAllPrim, 23.4450f },
- { L1CATrackFinderIter::kAllPrimE, 23.4450f },
- { L1CATrackFinderIter::kAllPrimJump, 18.7560f },
- { L1CATrackFinderIter::kAllSec, 18.7560f },
- { L1CATrackFinderIter::kAllSecE, 18.7560f },
- { L1CATrackFinderIter::kFastPrimJump, 21.1075f },
- { L1CATrackFinderIter::kAllSecJump, 21.1075f }
- };
+ // // TODO: Develope naming system for these constants
+ // L1CATrackFinderIterParameterMap_t<float> fTrackChi2Cut {
+ // { L1CATrackFinderIter::kFastPrim, 10.f },
+ // { L1CATrackFinderIter::kAllPrim, 10.f },
+ // { L1CATrackFinderIter::kAllPrimJump, 10.f },
+ // { L1CATrackFinderIter::kAllSec, 10.f },
+ // { L1CATrackFinderIter::kAllPrimE, 10.f },
+ // { L1CATrackFinderIter::kAllSecE, 10.f },
+ // { L1CATrackFinderIter::kFastPrimJump, 10.f },
+ // { L1CATrackFinderIter::kAllSecJump, 10.f }
+ // };
+ // // TODO: What is the meaning of the following numbers?
+ // // 23.4450 <- 7.815 * 3
+ // // 18.7560 <- 6.252 * 3
+ // // 11.3449 * 2. / 3.
+ // //
+ // L1CATrackFinderIterParameterMap_t<float> fTripletChi2Cut {
+ // { L1CATrackFinderIter::kFastPrim, 23.4450f }, // kFastPrimIter
+ // { L1CATrackFinderIter::kAllPrim, 23.4450f },
+ // { L1CATrackFinderIter::kAllPrimE, 23.4450f },
+ // { L1CATrackFinderIter::kAllPrimJump, 18.7560f },
+ // { L1CATrackFinderIter::kAllSec, 18.7560f },
+ // { L1CATrackFinderIter::kAllSecE, 18.7560f },
+ // { L1CATrackFinderIter::kFastPrimJump, 21.1075f },
+ // { L1CATrackFinderIter::kAllSecJump, 21.1075f }
+ // };
- L1CATrackFinderIterParameterMap_t<float> fDoubletChi2Cut {
- { L1CATrackFinderIter::kFastPrim, 11.3449 * 2.f / 3.f },
- { L1CATrackFinderIter::kAllPrim, 11.3449 * 2.f / 3.f },
- { L1CATrackFinderIter::kAllPrimJump, 11.3449 * 2.f / 3.f },
- { L1CATrackFinderIter::kAllSec, 11.3449 * 2.f / 3.f },
- { L1CATrackFinderIter::kAllPrimE, 11.3449 * 2.f / 3.f },
- { L1CATrackFinderIter::kAllSecE, 11.3449 * 2.f / 3.f },
- { L1CATrackFinderIter::kFastPrimJump, 11.3449 * 2.f / 3.f },
- { L1CATrackFinderIter::kAllSecJump, 11.3449 * 2.f / 3.f }
- };
+ // L1CATrackFinderIterParameterMap_t<float> fDoubletChi2Cut {
+ // { L1CATrackFinderIter::kFastPrim, 11.3449 * 2.f / 3.f },
+ // { L1CATrackFinderIter::kAllPrim, 11.3449 * 2.f / 3.f },
+ // { L1CATrackFinderIter::kAllPrimJump, 11.3449 * 2.f / 3.f },
+ // { L1CATrackFinderIter::kAllSec, 11.3449 * 2.f / 3.f },
+ // { L1CATrackFinderIter::kAllPrimE, 11.3449 * 2.f / 3.f },
+ // { L1CATrackFinderIter::kAllSecE, 11.3449 * 2.f / 3.f },
+ // { L1CATrackFinderIter::kFastPrimJump, 11.3449 * 2.f / 3.f },
+ // { L1CATrackFinderIter::kAllSecJump, 11.3449 * 2.f / 3.f }
+ // };
};
#endif // L1Parameters_h
diff --git a/reco/L1/L1Algo/L1Station.h b/reco/L1/L1Algo/L1Station.h
index fd3770e428..b2305966b3 100644
--- a/reco/L1/L1Algo/L1Station.h
+++ b/reco/L1/L1Algo/L1Station.h
@@ -13,23 +13,24 @@
class L1Station {
public:
- L1Station()
- : type(0)
- , timeInfo(0)
- , z(0)
- , Rmin(0)
- , Rmax(0)
- , Sy(0)
- , /// z-position, min & max station radius, field integral
- materialInfo()
- , fieldSlice()
- , frontInfo()
- , backInfo()
- , xInfo()
- , yInfo()
- , XYInfo()
- {
- }
+ // TODO: rewrite initialization (default values to fields, remove constructor (check struct))
+ //L1Station()
+ // : type(0)
+ // , timeInfo(0)
+ // , z(0)
+ // , Rmin(0)
+ // , Rmax(0)
+ // , Sy(0)
+ // , /// z-position, min & max station radius, field integral
+ // materialInfo()
+ // , fieldSlice()
+ // , frontInfo()
+ // , backInfo()
+ // , xInfo()
+ // , yInfo()
+ // , XYInfo()
+ //{
+ //}
// TODO: test this constructors
//L1Station(const L1Station &) = default;
@@ -37,19 +38,65 @@ public:
//L1Station(L1Station &&) = default;
//L1Station & operator=(L1Station &&) = default;
- int type;
- int timeInfo;
- fvec z;
- fvec Rmin;
- fvec Rmax;
- fvec Sy;
- L1MaterialInfo materialInfo;
- L1FieldSlice fieldSlice;
- L1UMeasurementInfo frontInfo;
- L1UMeasurementInfo backInfo;
- L1UMeasurementInfo xInfo;
- L1UMeasurementInfo yInfo;
- L1XYMeasurementInfo XYInfo;
+ int type {0};
+ int timeInfo {0};
+ fvec z {0};
+ fvec Rmin {0};
+ fvec Rmax {0};
+ fvec Sy {0};
+ L1MaterialInfo materialInfo {};
+ L1FieldSlice fieldSlice {};
+ L1UMeasurementInfo frontInfo {};
+ L1UMeasurementInfo backInfo {};
+ L1UMeasurementInfo xInfo {};
+ L1UMeasurementInfo yInfo {};
+ L1XYMeasurementInfo XYInfo {};
+
+
+ void Print() const
+ {
+ LOG(info) << "==== L1Station object at " << this;
+ LOG(info) << "- L1Station fields:";
+ LOG(info) << "--- Station type ID: " << type;
+ LOG(info) << "--- z position: " << z[0];
+ LOG(info) << "--- Rmin: " << Rmin[0];
+ LOG(info) << "--- Rmax: " << Rmax[0];
+ LOG(info) << "--- Thickness (X), cm: " << materialInfo.thick[0];
+ LOG(info) << "--- Radiational length (X0), cm: " << materialInfo.RL[0];
+ LOG(info) << "--- X / X0: " << materialInfo.RadThick[0];
+ LOG(info) << "--- log(X / X0): " << materialInfo.logRadThick[0];
+ LOG(info) << "--- Field approximation coefficients:";
+ LOG(info) << " idx CX CY CZ";
+ for (int idx = 0; idx < 21; ++idx) {
+ LOG(info)
+ << std::setw(9) << std::setfill(' ') << idx << ' ' << std::setw(10) << std::setfill(' ')
+ << fieldSlice.cx[idx][0] << ' ' << std::setw(10) << std::setfill(' ')
+ << fieldSlice.cy[idx][0] << ' ' << std::setw(10) << std::setfill(' ')
+ << fieldSlice.cz[idx][0];
+ }
+ LOG(info) << "--- Strips geometry:";
+ LOG(info) << "----- Front:";
+ LOG(info) << "------- cos(phi): " << frontInfo.cos_phi[0];
+ LOG(info) << "------- sin(phi): " << frontInfo.sin_phi[0];
+ LOG(info) << "------- sigma2: " << frontInfo.sigma2[0];
+ LOG(info) << "----- Back:";
+ LOG(info) << "------- cos(phi): " << backInfo.cos_phi[0];
+ LOG(info) << "------- sin(phi): " << backInfo.sin_phi[0];
+ LOG(info) << "------- sigma2: " << backInfo.sigma2[0];
+ LOG(info) << "----- XY cov matrix:";
+ LOG(info) << "------- C00: " << XYInfo.C00[0];
+ LOG(info) << "------- C10: " << XYInfo.C10[0];
+ LOG(info) << "------- C11: " << XYInfo.C11[0];
+ LOG(info) << "----- X layer:";
+ LOG(info) << "------- cos(phi): " << xInfo.cos_phi[0];
+ LOG(info) << "------- sin(phi): " << xInfo.sin_phi[0];
+ LOG(info) << "------- sigma2: " << xInfo.sigma2[0];
+ LOG(info) << "----- Y layer:";
+ LOG(info) << "------- cos(phi): " << yInfo.cos_phi[0];
+ LOG(info) << "------- sin(phi): " << yInfo.sin_phi[0];
+ LOG(info) << "------- sigma2: " << yInfo.sigma2[0];
+ LOG(info) << "";
+ }
} _fvecalignment;
diff --git a/reco/L1/L1Algo/L1Vector.h b/reco/L1/L1Algo/L1Vector.h
index 81b88d25c1..51ecf90d7a 100644
--- a/reco/L1/L1Algo/L1Vector.h
+++ b/reco/L1/L1Algo/L1Vector.h
@@ -45,7 +45,10 @@ public:
{
}
-
+ /// Constructor to make initializations like L1Vector<int> myVector {"MyVector", {1, 2, 3}}
+ L1Vector(const std::string& name, std::initializer_list<T> init): Tbase(init), fName(name)
+ {
+ }
L1Vector(const L1Vector& v) : Tbase() { *this = v; }
--
GitLab