From 14f058f8c222106faf1dafb7fc1c1bf53a1743c1 Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Thu, 19 Oct 2023 21:37:12 +0000
Subject: [PATCH] Ca: move Ca track finder to /algo
---
algo/ca/core/CMakeLists.txt | 16 +
.../ca/core/tracking/CaCloneMerger.cxx | 40 +-
algo/ca/core/tracking/CaCloneMerger.h | 127 +++
.../ca/core/tracking/CaFramework.cxx | 59 +-
algo/ca/core/tracking/CaFramework.h | 358 +++++++
.../ca/core/tracking/CaTrackExtender.cxx | 136 +--
algo/ca/core/tracking/CaTrackExtender.h | 93 ++
.../ca/core/tracking/CaTrackFinder.cxx | 150 +--
algo/ca/core/tracking/CaTrackFinder.h | 66 ++
algo/ca/core/tracking/CaTrackFinderWindow.cxx | 771 ++++++++++++++++
algo/ca/core/tracking/CaTrackFinderWindow.h | 73 ++
.../ca/core/tracking/CaTrackFitter.cxx | 68 +-
algo/ca/core/tracking/CaTrackFitter.h | 66 ++
.../ca/core/tracking/CaTripletConstructor.cxx | 55 +-
algo/ca/core/tracking/CaTripletConstructor.h | 122 +++
reco/KF/ParticleFitter/CbmL1PFFitter.cxx | 43 +-
reco/L1/CMakeLists.txt | 16 +-
reco/L1/CbmL1.cxx | 72 +-
reco/L1/CbmL1.h | 21 +-
reco/L1/CbmL1MCTrack.cxx | 4 +-
reco/L1/CbmL1Performance.cxx | 61 +-
reco/L1/CbmL1ReadEvent.cxx | 20 +-
reco/L1/L1Algo/L1Algo.h | 435 ---------
reco/L1/L1Algo/L1CaTrackFinderSlice.cxx | 871 ------------------
reco/L1/L1Algo/L1CloneMerger.h | 133 ---
reco/L1/L1Algo/L1TripletConstructor.h | 117 ---
reco/L1/L1Algo/utils/L1AlgoDraw.cxx | 11 +-
reco/L1/L1Algo/utils/L1AlgoDraw.h | 10 +-
.../utils/L1AlgoEfficiencyPerformance.h | 3 +-
reco/L1/L1Algo/utils/L1AlgoPulls.cxx | 2 +
reco/L1/L1Algo/utils/L1AlgoPulls.h | 5 +-
.../CbmL1GlobalTrackFinder.cxx | 1 -
.../OffLineInterface/CbmL1StsTrackFinder.cxx | 1 -
reco/L1/catools/CaToolsDebugger.cxx | 8 +-
reco/L1/catools/CaToolsDebugger.h | 1 +
35 files changed, 2110 insertions(+), 1925 deletions(-)
rename reco/L1/L1Algo/L1CloneMerger.cxx => algo/ca/core/tracking/CaCloneMerger.cxx (91%)
create mode 100644 algo/ca/core/tracking/CaCloneMerger.h
rename reco/L1/L1Algo/L1Algo.cxx => algo/ca/core/tracking/CaFramework.cxx (76%)
create mode 100644 algo/ca/core/tracking/CaFramework.h
rename reco/L1/L1Algo/L1BranchExtender.cxx => algo/ca/core/tracking/CaTrackExtender.cxx (60%)
create mode 100644 algo/ca/core/tracking/CaTrackExtender.h
rename reco/L1/L1Algo/L1CaTrackFinder.cxx => algo/ca/core/tracking/CaTrackFinder.cxx (53%)
create mode 100644 algo/ca/core/tracking/CaTrackFinder.h
create mode 100644 algo/ca/core/tracking/CaTrackFinderWindow.cxx
create mode 100644 algo/ca/core/tracking/CaTrackFinderWindow.h
rename reco/L1/L1Algo/L1TrackFitter.cxx => algo/ca/core/tracking/CaTrackFitter.cxx (81%)
create mode 100644 algo/ca/core/tracking/CaTrackFitter.h
rename reco/L1/L1Algo/L1TripletConstructor.cxx => algo/ca/core/tracking/CaTripletConstructor.cxx (93%)
create mode 100644 algo/ca/core/tracking/CaTripletConstructor.h
delete mode 100644 reco/L1/L1Algo/L1Algo.h
delete mode 100644 reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
delete mode 100644 reco/L1/L1Algo/L1CloneMerger.h
delete mode 100644 reco/L1/L1Algo/L1TripletConstructor.h
diff --git a/algo/ca/core/CMakeLists.txt b/algo/ca/core/CMakeLists.txt
index da259ab55b..0adcb851a5 100644
--- a/algo/ca/core/CMakeLists.txt
+++ b/algo/ca/core/CMakeLists.txt
@@ -26,9 +26,18 @@ set(SRCS
${CMAKE_CURRENT_SOURCE_DIR}/pars/CaSearchWindow.cxx
${CMAKE_CURRENT_SOURCE_DIR}/pars/CaStation.cxx
${CMAKE_CURRENT_SOURCE_DIR}/pars/CaStationInitializer.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/tracking/CaCloneMerger.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/tracking/CaFramework.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/tracking/CaTrackExtender.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/tracking/CaTrackFinder.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/tracking/CaTrackFinderWindow.cxx
${CMAKE_CURRENT_SOURCE_DIR}/tracking/CaTrackFit.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/tracking/CaTrackFitter.cxx
+ ${CMAKE_CURRENT_SOURCE_DIR}/tracking/CaTripletConstructor.cxx
)
+SET_SOURCE_FILES_PROPERTIES(${SRCS} PROPERTIES COMPILE_FLAGS "-O3")
+
add_library(CaCore SHARED ${SRCS})
target_include_directories(CaCore
@@ -93,7 +102,14 @@ install(
utils/CaVector.h
utils/CaUtils.h
utils/CaDefines.h
+ tracking/CaCloneMerger.h
+ tracking/CaFramework.h
+ tracking/CaTrackExtender.h
+ tracking/CaTrackFinder.h
+ tracking/CaTrackFinderWindow.h
tracking/CaTrackFit.h
+ tracking/CaTrackFitter.h
+ tracking/CaTripletConstructor.h
DESTINATION
include/
)
diff --git a/reco/L1/L1Algo/L1CloneMerger.cxx b/algo/ca/core/tracking/CaCloneMerger.cxx
similarity index 91%
rename from reco/L1/L1Algo/L1CloneMerger.cxx
rename to algo/ca/core/tracking/CaCloneMerger.cxx
index 707d12ba6c..37074f5a88 100644
--- a/reco/L1/L1Algo/L1CloneMerger.cxx
+++ b/algo/ca/core/tracking/CaCloneMerger.cxx
@@ -2,35 +2,35 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Sergei Zharko, Maksym Zyzak [committer]*/
-/// \file L1CloneMerger.h
+/// \file CloneMerger.h
/// \authors S.Zharko <s.zharko@gsi.de> (interface), M.Zyzak (original algorithm)
-/// \brief A class wrapper over clones merger algorithm for the L1 track finder (implementation)
+/// \brief A class wrapper over clones merger algorithm for the Ca track finder (implementation)
/// \since 22.07.2022 (second version)
-#include "L1CloneMerger.h"
+#include "CaCloneMerger.h"
#include <iostream>
+#include "CaFramework.h"
#include "CaParameters.h"
#include "CaTrack.h"
#include "CaTrackFit.h"
#include "CaVector.h"
-#include "L1Algo.h"
-using cbm::algo::ca::Track; // TMP
-using cbm::algo::ca::Vector; // TMP
+using namespace cbm::algo::ca;
+using namespace cbm::algo;
// ---------------------------------------------------------------------------------------------------------------------
//
-L1CloneMerger::L1CloneMerger(const L1Algo& algo) : frAlgo(algo) {}
+CloneMerger::CloneMerger(const ca::Framework& algo) : frAlgo(algo) {}
// ---------------------------------------------------------------------------------------------------------------------
//
-L1CloneMerger::~L1CloneMerger() {}
+CloneMerger::~CloneMerger() {}
// ---------------------------------------------------------------------------------------------------------------------
//
-void L1CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extRecoHits)
+void CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extRecoHits)
{
Vector<unsigned short>& firstStation = fTrackFirstStation;
Vector<unsigned short>& lastStation = fTrackLastStation;
@@ -93,7 +93,7 @@ void L1CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extRe
ca::FieldRegion fld _fvecalignment;
// Max length for merging
- unsigned char maxLengthForMerge = static_cast<unsigned char>(frAlgo.GetParameters()->GetNstationsActive() - 3);
+ unsigned char maxLengthForMerge = static_cast<unsigned char>(frAlgo.GetParameters().GetNstationsActive() - 3);
for (int iTr = 0; iTr < nTracks; iTr++) {
if (extTracks[iTr].fNofHits > maxLengthForMerge) continue;
@@ -119,8 +119,8 @@ void L1CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extRe
unsigned short stam;
- fBf = frAlgo.GetParameters()->GetStation(staf).fieldSlice.GetFieldValue(Tf.X(), Tf.Y());
- fBb = frAlgo.GetParameters()->GetStation(stab).fieldSlice.GetFieldValue(Tb.X(), Tb.Y());
+ fBf = frAlgo.GetParameters().GetStation(staf).fieldSlice.GetFieldValue(Tf.X(), Tf.Y());
+ fBb = frAlgo.GetParameters().GetStation(stab).fieldSlice.GetFieldValue(Tb.X(), Tb.Y());
unsigned short dist = firstStation[iTr] - lastStation[jTr];
@@ -128,10 +128,10 @@ void L1CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extRe
else
stam = staf - 1;
- fvec zm = frAlgo.GetParameters()->GetStation(stam).fZ;
+ fvec zm = frAlgo.GetParameters().GetStation(stam).fZ;
fvec xm = fvec(0.5) * (Tf.GetX() + Tf.Tx() * (zm - Tf.Z()) + Tb.GetX() + Tb.Tx() * (zm - Tb.Z()));
fvec ym = fvec(0.5) * (Tf.Y() + Tf.Ty() * (zm - Tf.Z()) + Tb.Y() + Tb.Ty() * (zm - Tb.Z()));
- fBm = frAlgo.GetParameters()->GetStation(stam).fieldSlice.GetFieldValue(xm, ym);
+ fBm = frAlgo.GetParameters().GetStation(stam).fieldSlice.GetFieldValue(xm, ym);
fld.Set(fBb, Tb.Z(), fBm, zm, fBf, Tf.Z());
fvec zMiddle = fvec(0.5) * (Tb.Z() + Tf.Z());
@@ -195,8 +195,8 @@ void L1CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extRe
// ---------------------------------------------------------------------------------------------------------------------
//
-void L1CloneMerger::FilterTracks(fvec const r[5], fvec const C[15], fvec const m[5], fvec const V[15], fvec R[5],
- fvec W[15], fvec* chi2)
+void CloneMerger::FilterTracks(fvec const r[5], fvec const C[15], fvec const m[5], fvec const V[15], fvec R[5],
+ fvec W[15], fvec* chi2)
{
fvec S[15];
for (int i = 0; i < 15; i++) {
@@ -241,7 +241,7 @@ void L1CloneMerger::FilterTracks(fvec const r[5], fvec const C[15], fvec const m
// ---------------------------------------------------------------------------------------------------------------------
//
-void L1CloneMerger::InvertCholesky(fvec a[15])
+void CloneMerger::InvertCholesky(fvec a[15])
{
fvec d[5], uud, u[5][5];
for (int i = 0; i < 5; i++) {
@@ -306,7 +306,7 @@ void L1CloneMerger::InvertCholesky(fvec a[15])
// ---------------------------------------------------------------------------------------------------------------------
//
-void L1CloneMerger::MultiplyMS(fvec const C[5][5], fvec const V[15], fvec K[15])
+void CloneMerger::MultiplyMS(fvec const C[5][5], fvec const V[15], fvec K[15])
{
K[0] = C[0][0] * V[0] + C[0][1] * V[1] + C[0][2] * V[3] + C[0][3] * V[6] + C[0][4] * V[10];
@@ -331,7 +331,7 @@ void L1CloneMerger::MultiplyMS(fvec const C[5][5], fvec const V[15], fvec K[15])
// ---------------------------------------------------------------------------------------------------------------------
//
-void L1CloneMerger::MultiplySR(fvec const C[15], fvec const r_in[5], fvec r_out[5])
+void CloneMerger::MultiplySR(fvec const C[15], fvec const r_in[5], fvec r_out[5])
{
r_out[0] = r_in[0] * C[0] + r_in[1] * C[1] + r_in[2] * C[3] + r_in[3] * C[6] + r_in[4] * C[10];
r_out[1] = r_in[0] * C[1] + r_in[1] * C[2] + r_in[2] * C[4] + r_in[3] * C[7] + r_in[4] * C[11];
@@ -342,7 +342,7 @@ void L1CloneMerger::MultiplySR(fvec const C[15], fvec const r_in[5], fvec r_out[
// ---------------------------------------------------------------------------------------------------------------------
//
-void L1CloneMerger::MultiplySS(fvec const C[15], fvec const V[15], fvec K[5][5])
+void CloneMerger::MultiplySS(fvec const C[15], fvec const V[15], fvec K[5][5])
{
K[0][0] = C[0] * V[0] + C[1] * V[1] + C[3] * V[3] + C[6] * V[6] + C[10] * V[10];
K[0][1] = C[0] * V[1] + C[1] * V[2] + C[3] * V[4] + C[6] * V[7] + C[10] * V[11];
diff --git a/algo/ca/core/tracking/CaCloneMerger.h b/algo/ca/core/tracking/CaCloneMerger.h
new file mode 100644
index 0000000000..8acc9fb88d
--- /dev/null
+++ b/algo/ca/core/tracking/CaCloneMerger.h
@@ -0,0 +1,127 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer], Maksym Zyzak */
+
+/// \file CloneMerger.h
+/// \authors S.Zharko <s.zharko@gsi.de> (interface), M.Zyzak (original algorithm)
+/// \brief A class wrapper over clones merger algorithm for the Ca track finder (declaration)
+/// \since 22.07.2022
+
+#pragma once // include this header only once per compilation unit
+
+#include "CaHit.h" // For ca::HitIndex_t
+#include "CaSimd.h" // TEMPORARY FOR fvec, fscal
+#include "CaVector.h"
+
+
+namespace cbm::algo::ca
+{
+ class Track;
+ class Framework;
+
+ namespace
+ {
+ using namespace cbm::algo; // to keep ca:: prefices in the code
+ }
+
+
+ /// Class implements a clones merger algorithm for the CA track finder
+ ///
+ class CloneMerger {
+ public:
+ /// Default constructor
+ CloneMerger(const ca::Framework& algo);
+
+ /// Destructor
+ ~CloneMerger();
+
+ /// Copy constructor
+ CloneMerger(const CloneMerger&) = delete;
+
+ /// Move constructor
+ CloneMerger(CloneMerger&&) = delete;
+
+ /// Copy assignment operator
+ CloneMerger& operator=(const CloneMerger&) = delete;
+
+ /// Move assignment operator
+ CloneMerger& operator=(CloneMerger&&) = delete;
+
+ /// Registers
+
+ /// Executes track clones merging algorithm and updates input containers
+ /// \param extTracks Reference to the external container of reconstructed tracks
+ /// \param extRecoHits Reference to the external container of reconstructed hit indexes
+ void Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>&);
+
+ private:
+ // ***************
+ // ** Functions **
+ // ***************
+
+ ///
+ static void InvertCholesky(fvec a[15]);
+
+ /// Multiplication of two symmetric matryces 5x5
+ /// \param C Left matrix:
+ /// C[0] C[1] C[3] C[6] C[10]
+ /// C[1] C[2] C[4] C[7] C[11]
+ /// C = C[3] C[4] C[5] C[8] C[12]
+ /// C[6] C[7] C[8] C[9] C[13]
+ /// C[10] C[11] C[12] C[13] C[14]
+ /// \param V Right matrix:
+ /// V[0] V[1] V[3] V[6] V[10]
+ /// V[1] V[2] V[4] V[7] V[11]
+ /// V = V[3] V[4] V[5] V[8] V[12]
+ /// V[6] V[7] V[8] V[9] V[13]
+ /// V[10] V[11] V[12] V[13] V[14]
+ /// \param K Output: K = C * V
+ static void MultiplySS(fvec const C[15], fvec const V[15], fvec K[5][5]);
+
+ ///
+ static void MultiplyMS(fvec const C[5][5], fvec const V[15], fvec K[15]);
+
+ ///
+ static void MultiplySR(fvec const C[15], fvec const r_in[5], fvec r_out[5]);
+
+ ///
+ static void FilterTracks(fvec const r[5], fvec const C[15], fvec const m[5], fvec const V[15], fvec R[5],
+ fvec W[15], fvec* chi2);
+
+
+ // ***************
+ // ** Variables **
+ // ***************
+
+ /// First station of a track
+ Vector<unsigned short> fTrackFirstStation {"CloneMerger::fTrackFirstStation"};
+
+ /// Last station of a track
+ Vector<unsigned short> fTrackLastStation {"CloneMerger::fTrackLastStation"};
+
+ /// Index of the first hit of a track
+ Vector<ca::HitIndex_t> fTrackFirstHit {"CloneMerger::fTrackFirstHit"};
+
+ /// Index of the last hit of a track
+ Vector<ca::HitIndex_t> fTrackLastHit {"CloneMerger::fTrackLastHit"};
+
+ /// Index (TODO:??) of a track that can be merge with the given track
+ Vector<unsigned short> fTrackNeighbour {"CloneMerger::fTrackNeighbour"};
+
+ /// Chi2 value of the track merging procedure
+ Vector<float> fTrackChi2 {"CloneMerger::fTrackChi2"};
+
+ /// Flag: is the given track already stored to the output
+ Vector<char> fTrackIsStored {"CloneMerger::fTrackIsStored"};
+
+ /// Flag: is the track a downstream neighbour of another track
+ Vector<char> fTrackIsDownstreamNeighbour {"CloneMerger::fTrackIsDownstreamNeighbour"};
+
+ Vector<Track> fTracksNew {"CaCloneMerger::fTracksNew"}; ///< vector of tracks after the merge
+
+ Vector<ca::HitIndex_t> fRecoHitsNew {"CaCloneMerger::fRecoHitsNew"}; ///< vector of track hits after the merge
+
+ const ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class
+ };
+
+} // namespace cbm::algo::ca
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/algo/ca/core/tracking/CaFramework.cxx
similarity index 76%
rename from reco/L1/L1Algo/L1Algo.cxx
rename to algo/ca/core/tracking/CaFramework.cxx
index 9b3a9f2e27..539483fe87 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/algo/ca/core/tracking/CaFramework.cxx
@@ -2,31 +2,33 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Igor Kulakov [committer], Maksym Zyzak, Valentina Akishina */
-#include "L1Algo.h"
-
-#include "CbmL1.h"
+#include "CaFramework.h"
#include "CaGridEntry.h"
-#include "CaToolsDebugger.h"
+// #include "CaToolsDebugger.h"
+
+using namespace cbm::algo::ca;
-L1Algo::L1Algo()
+namespace
+{
+ using namespace cbm::algo; // to keep ca:: prefices in the code
+}
+
+Framework::Framework()
{
for (unsigned int i = 0; i < constants::size::MaxNstations; i++) {
- fTriplets[i].SetName(std::stringstream() << "L1Algo::fTriplets[" << i << "]");
+ fTriplets[i].SetName(std::stringstream() << "Framework::fTriplets[" << i << "]");
}
}
using cbm::algo::ca::ECounter; // monitor counter key type
using cbm::algo::ca::EDetectorID;
using cbm::algo::ca::ETimer; // monitor timer key type
-using cbm::ca::tools::Debugger;
+//using cbm::ca::tools::Debugger;
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void L1Algo::Init(const TrackingMode mode)
+void Framework::Init(const TrackingMode mode)
{
fTrackingMode = mode;
-
// Monitor settings
fMonitor.SetCounterName(ECounter::RecoTrack, "reco tracks");
fMonitor.SetCounterName(ECounter::RecoHit, "reco hits");
@@ -42,21 +44,23 @@ void L1Algo::Init(const TrackingMode mode)
fMonitor.Reset();
}
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void L1Algo::Finish() { Debugger::Instance().Write(); }
+void Framework::Finish()
+{
+ //Debugger::Instance().Write();
+}
// ---------------------------------------------------------------------------------------------------------------------
//
-void L1Algo::ReceiveInputData(InputData&& inputData)
+void Framework::ReceiveInputData(InputData&& inputData)
{
// ----- Get input data ----------------------------------------------------------------------------------------------
fInputData = std::move(inputData);
}
+
// ---------------------------------------------------------------------------------------------------------------------
//
-void L1Algo::ReceiveParameters(Parameters&& parameters)
+void Framework::ReceiveParameters(Parameters&& parameters)
{
fParameters = std::move(parameters);
@@ -76,34 +80,43 @@ void L1Algo::ReceiveParameters(Parameters&& parameters)
ca::FieldRegion::ForceUseOfOriginalField(fParameters.DevIsUseOfOriginalField());
}
-int L1Algo::GetMcTrackIdForCaHit(int iHit) const
+int Framework::GetMcTrackIdForCaHit(int /*iHit*/) const
{
+ return -1;
+ /*
int hitId = iHit;
int iMcPoint = CbmL1::Instance()->GetHitBestMcRefs()[hitId];
if (iMcPoint < 0) return -1;
return CbmL1::Instance()->GetMcPoints()[iMcPoint].ID;
+ */
}
-int L1Algo::GetMcTrackIdForWindowHit(int iHit) const
+int Framework::GetMcTrackIdForWindowHit(int /*iHit*/) const
{
+ return -1;
+ /*
int hitId = fWindowHits[iHit].Id();
int iMcPoint = CbmL1::Instance()->GetHitBestMcRefs()[hitId];
if (iMcPoint < 0) return -1;
return CbmL1::Instance()->GetMcPoints()[iMcPoint].ID;
+ */
}
-const CbmL1MCTrack* L1Algo::GetMcTrackForWindowHit(int iHit) const
+/*
+const CbmL1MCTrack* Framework::GetMcTrackForWindowHit(int iHit) const
{
+ return nullptr;
int id = GetMcTrackIdForWindowHit(iHit);
if (id < 0) return nullptr;
return &CbmL1::Instance()->GetMcTracks()[id];
}
+*/
-// bool L1Algo::SortTrip(TripSort const& a, TripSort const& b) {
+// bool Framework::SortTrip(TripSort const& a, TripSort const& b) {
// return ( a.trip.GetLevel() > b.trip.GetLevel() );
// }
//
-// bool L1Algo::SortCand(CandSort const& a, CandSort const& b) {
+// bool Framework::SortCand(CandSort const& a, CandSort const& b) {
// if (a.cand.Lengtha != b.cand.Lengtha) return (a.cand.Lengtha > b.cand.Lengtha);
//
// if (a.cand.ista != b.cand.ista ) return (a.cand.ista < b.cand.ista );
@@ -114,8 +127,8 @@ const CbmL1MCTrack* L1Algo::GetMcTrackForWindowHit(int iHit) const
// // return (a.cand.CandIndex > b.cand.CandIndex );
// }
-// inline int L1Algo::PackIndex(const int& a, const int& b, const int& c) {
+// inline int Framework::PackIndex(const int& a, const int& b, const int& c) {
// return (a) + ((b)*10000) + (c*100000000);
// }
//
-// inline int L1Algo::UnPackIndex(const int& i, int& a, int& b, int& c) {
+// inline int Framework::UnPackIndex(const int& i, int& a, int& b, int& c) {
diff --git a/algo/ca/core/tracking/CaFramework.h b/algo/ca/core/tracking/CaFramework.h
new file mode 100644
index 0000000000..a279435e0d
--- /dev/null
+++ b/algo/ca/core/tracking/CaFramework.h
@@ -0,0 +1,358 @@
+/* Copyright (C) 2007-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Maksym Zyzak, Valentina Akishina, Igor Kulakov [committer], Sergei Zharko */
+
+#pragma once // include this header only once per compilation unit
+
+#include <array>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <map>
+
+#include "CaBranch.h"
+#include "CaCloneMerger.h"
+#include "CaConstants.h"
+#include "CaField.h"
+#include "CaGrid.h"
+#include "CaGridEntry.h"
+#include "CaHit.h"
+#include "CaInputData.h"
+#include "CaMeasurementXy.h"
+#include "CaParameters.h"
+#include "CaStation.h"
+#include "CaTimer.h"
+#include "CaTrack.h"
+#include "CaTrackExtender.h"
+#include "CaTrackFinder.h"
+#include "CaTrackFinderWindow.h"
+#include "CaTrackFitter.h"
+#include "CaTrackParam.h"
+#include "CaTrackingMonitor.h"
+#include "CaTriplet.h"
+#include "CaVector.h"
+
+namespace cbm::algo::ca
+{
+ class TripletConstructor;
+
+
+ namespace
+ {
+ using namespace cbm::algo; // to keep ca:: prefices in the code
+ }
+
+ //namespace
+ //{
+ // using cbm::algo::ca::Track; // TMP
+ // using cbm::algo::ca::Vector; // TMP
+ // using cbm::algo::ca::Iteration; // TMP
+ //} // namespace
+
+ // *******************************
+ // ** Types definition (global) **
+ // *******************************
+
+ using CaStationsArray_t = std::array<ca::Station, constants::size::MaxNstations>;
+ using CaMaterialArray_t = std::array<ca::MaterialMap, constants::size::MaxNstations>;
+ using Tindex = int; // TODO: Replace with ca::HitIndex_t, if suitable
+
+ struct CaHitTimeInfo {
+ fscal fEventTimeMin {-std::numeric_limits<fscal>::max() / 2.};
+ fscal fEventTimeMax {std::numeric_limits<fscal>::max() / 2.};
+ fscal fMaxTimeBeforeHit {0.}; //< max event time for hits [0 .. hit] in the station hit array
+ fscal fMinTimeAfterHit {0.}; //< min event time for hits [hit ... ] in the station hit array
+ };
+
+ /// Main class of CA track finder algorithm
+ ///
+ class Framework {
+ public:
+ // *************************
+ // ** Friend classes list **
+ // *************************
+
+ friend class ca::TripletConstructor;
+ friend class ca::TrackExtender;
+ friend class ca::TrackFitter;
+ friend class ca::TrackFinderWindow;
+ friend class ca::TrackFinder;
+
+ ///---------------------------
+ /// Internal enumerations
+
+ enum TrackingMode
+ {
+ kSts,
+ kGlobal,
+ kMcbm
+ };
+
+ // **********************************
+ // ** Member functions declaration **
+ // **********************************
+
+ // ** Constructors and destructor
+
+ /// Constructor
+ Framework();
+
+ /// Copy constructor
+ Framework(const Framework&) = delete;
+
+ /// Move constructor
+ Framework(Framework&&) = delete;
+
+ /// Copy assignment operator
+ Framework& operator=(const Framework&) = delete;
+
+ /// Move assignment operator
+ Framework& operator=(Framework&&) = delete;
+
+ /// Destructor
+ ~Framework() = default;
+
+
+ // ** Functions, which pack and unpack indexes of station and triplet **
+
+ // TODO: move to ca::Triplet class (S.Zharko)
+ /// Packs station and triplet indices to an unique triplet ID
+ /// \param iStation Index of station in the active stations array
+ /// \param iTriplet Index of triplet
+ static unsigned int PackTripletId(unsigned int iStation, unsigned int iTriplet);
+
+ /// Unpacks the triplet ID to its station index
+ /// \param id Unique triplet ID
+ static unsigned int TripletId2Station(unsigned int id);
+
+ /// Unpacks the triplet ID to its triplet index
+ /// \param id Unique triplet ID
+ static unsigned int TripletId2Triplet(unsigned int id);
+
+
+ /// Sets Framework parameters object
+ /// \param other - reference to the Parameters object
+ void SetParameters(const Parameters& other) { fParameters = other; }
+ // TODO: remove it (S.Zharko)
+
+ /// Gets a pointer to the Framework parameters object
+ const Parameters& GetParameters() const { return fParameters; }
+
+ /// Receives input data
+ void ReceiveInputData(InputData&& inputData);
+
+ /// Receives tracking parameters
+ void ReceiveParameters(Parameters&& parameters);
+
+ /// Gets pointer to input data object for external access
+ const InputData& GetInputData() const { return fInputData; }
+
+ int PackIndex(const int& a, const int& b, const int& c);
+
+ int UnPackIndex(const int& i, int& a, int& b, int& c);
+
+ /// \brief Sets a default particle mass for the track fit
+ ///
+ /// The function is used during the reconstruction in order to estimate the multiple scattering and energy loss
+ /// \param mass Default particle mass
+ void SetDefaultParticleMass(float mass) { fDefaultMass = mass; }
+
+ /// Gets default particle mass
+ /// \return particle mass
+ float GetDefaultParticleMass() const { return fDefaultMass; }
+
+ /// Gets default particle mass squared
+ /// \return particle mass squared
+ float GetDefaultParticleMass2() const { return fDefaultMass * fDefaultMass; }
+
+
+ /*********************************************************************************************/ /**
+ * ------ FUNCTIONAL PART ------
+ ************************************************************************************************/
+
+
+ void Init(const TrackingMode mode);
+
+ void Finish();
+
+ void PrintHits();
+
+ float GetMaxInvMom() const { return fMaxInvMom[0]; }
+
+ /// Gets reference to the monitor
+ const TrackingMonitor& GetMonitor() const { return fMonitor; }
+
+ public:
+ /// Gets number of stations before the pipe (MVD stations in CBM)
+ int GetNstationsBeforePipe() const { return fNstationsBeforePipe; }
+
+ /// Gets number of stations situated in field region (MVD + STS in CBM)
+ int GetNfieldStations() const { return fNfieldStations; }
+
+ /// Get mc track ID for a hit (debug tool)
+ int GetMcTrackIdForCaHit(int iHit) const;
+
+ /// Get mc track ID for a hit (debug tool)
+ int GetMcTrackIdForWindowHit(int iGridHit) const;
+
+ // const CbmL1MCTrack* GetMcTrackForWindowHit(int iHit) const;
+
+ private:
+ int fNstationsBeforePipe {0}; ///< number of stations before pipe (MVD stations in CBM)
+ int fNfieldStations {0}; ///< number of stations in the field region
+ float fDefaultMass {constants::phys::MuonMass}; ///< mass of the propagated particle [GeV/c2]
+
+
+ // ***************************
+ // ** Member variables list **
+ // ***************************
+
+ Parameters fParameters; ///< Object of Framework parameters class
+ InputData fInputData; ///< Tracking input data
+
+ Vector<unsigned char> fvHitKeyFlags {
+ "Framework::fvHitKeyFlags"}; ///< List of key flags: has been this hit or cluster already used
+
+ ca::TrackingMonitor fMonitor {"CA Algo"}; ///< Tracking monitor
+
+ public:
+ Vector<CaHitTimeInfo> fHitTimeInfo;
+
+ ca::Grid vGrid[constants::size::MaxNstations]; ///<
+
+ /// ----- Data used during track finding -----
+ ///
+
+ ///\brief hits of the current time window
+ /// it is a portion of fInputData to process in the current time window
+ /// hit.Id is replaced by the hit index in fInputData
+ Vector<ca::Hit> fWindowHits {"Framework::fWindowHits"};
+
+ Vector<unsigned char> fIsWindowHitSuppressed {"Framework::fIsWindowHitSuppressed"};
+
+ ///\brief first index of the station hits in the time window
+ ca::HitIndex_t fStationHitsStartIndex[constants::size::MaxNstations + 1] {0};
+
+ ///\brief number of station hits in the time window
+ ca::HitIndex_t fStationNhits[constants::size::MaxNstations + 1] {0};
+
+ /// ----------
+
+ double fCaRecoTime {0.}; // time of the track finder + fitter
+
+ Vector<Track> fRecoTracks {"Framework::fRecoTracks"}; ///< reconstructed tracks
+ Vector<ca::HitIndex_t> fRecoHits {"Framework::fRecoHits"}; ///< packed hits of reconstructed tracks
+
+ Vector<Track> fSliceRecoTracks {"Framework::fSliceRecoTracks"}; ///< reconstructed tracks in sub-timeslice
+ Vector<ca::HitIndex_t> fSliceRecoHits {"Framework::fSliceRecoHits"}; ///< packed hits of reconstructed tracks
+
+ /// Created triplets vs station index
+ Vector<ca::Triplet> fTriplets[constants::size::MaxNstations] {{"Framework::fTriplets"}};
+
+ /// Track candidates created out of adjacent triplets before the final track selection.
+ /// The candidates may share any amount of hits.
+ Vector<ca::Branch> fTrackCandidates {"Framework::fTrackCandidates"};
+
+ ///< indices of the sub-slice hits
+ Vector<ca::HitIndex_t> fSliceHitIds[constants::size::MaxNstations] {"Framework::fSliceHitIds"};
+
+ Vector<int> fStripToTrack {"Framework::fStripToTrack"}; // strip to track pointers
+
+ TrackingMode fTrackingMode {kSts};
+
+ fvec EventTime {0.f};
+ fvec Err {0.f};
+
+ /// --- data used during finding iterations
+ unsigned int fCurrentIterationIndex {0}; ///< index of the corrent iteration, needed for debug output
+ const Iteration* fpCurrentIteration = nullptr; ///< pointer to the current CA track finder iteration
+
+ Vector<int> fHitFirstTriplet {"Framework::fHitFirstTriplet"}; /// link hit -> first triplet { hit, *, *}
+ Vector<int> fHitNtriplets {"Framework::fHitNtriplets"}; /// link hit ->n triplets { hit, *, *}
+
+
+ public:
+ ca::TrackExtender fTrackExtender {*this}; ///< Object of the track extender algorithm
+ ca::CloneMerger fCloneMerger {*this}; ///< Object of the clone merger algorithm
+ ca::TrackFitter fTrackFitter {*this}; ///< Object of the track extender algorithm
+ ca::TrackFinderWindow fTrackFinderWindow {*this}; ///< Object of the track finder algorithm for the time window
+ ca::TrackFinder fTrackFinder {*this}; ///< Object of the track finder algorithm for the entire time slice
+
+ private:
+ /// ================================= DATA PART =================================
+
+ /// ----- Different parameters of CATrackFinder -----
+
+ Tindex fFirstCAstation {0}; //first station used in CA
+
+ // fNFindIterations - set number of interation for trackfinding
+ // itetation of finding:
+
+ int fNFindIterations = -1; // TODO investigate kAllPrimJumpIter & kAllSecJumpIter
+
+ float fTrackChi2Cut {10.f};
+ float fTripletFinalChi2Cut {10.f};
+ float fTripletChi2Cut {5.f}; // cut for selecting triplets before collecting tracks.per one DoF
+ float fDoubletChi2Cut {5.f};
+ float fTimeCut1 {0.f}; // TODO: please, specify "1" and "2" (S.Zharko)
+ float fTimeCut2 {0.f};
+
+ /// correction in order to take into account overlaping and iff z. if sort by y then it is max diff between same station's modules (~0.4cm)
+ fvec fMaxDZ {constants::Undef<fvec>};
+
+ /// parameters which are different for different iterations. Set in the begin of CA Track Finder
+
+ float fPickGather {constants::Undef<fvec>}; ///< same for attaching additional hits to track
+ float fTripletLinkChi2 {
+ constants::Undef<fvec>}; ///< (dp2/dp_error2 < fTripletLinkChi2) => triplets are neighbours
+ fvec fMaxInvMom {constants::Undef<fvec>}; ///< max considered q/p for tracks
+ fvec fMaxSlopePV {constants::Undef<fvec>}; ///< max slope (tx\ty) in prim vertex
+ float fMaxSlope {constants::Undef<fvec>}; ///< max slope (tx\ty) in 3d hit position of a triplet
+
+ fvec fTargX {constants::Undef<fvec>}; ///< target position x coordinate for the current iteration (modifiable)
+ fvec fTargY {constants::Undef<fvec>}; ///< target position y coordinate for the current iteration (modifiable)
+ fvec fTargZ {constants::Undef<fvec>}; ///< target position z coordinate for the current iteration (modifiable)
+
+ bool fIsTargetField {false}; ///< is the magnetic field present at the target
+
+ ca::FieldValue fTargB _fvecalignment {}; // field in the target point (modifiable, do not touch!!)
+ ca::MeasurementXy<fvec> TargetMeasurement _fvecalignment {}; // target constraint [cm]
+
+ int fGhostSuppression {1}; // NOTE: Should be equal to 0 in TRACKS_FROM_TRIPLETS mode!
+
+ } _fvecalignment;
+
+
+ // ********************************************
+ // ** Inline member functions implementation **
+ // ********************************************
+
+ // ---------------------------------------------------------------------------------------------------------------------
+ //
+ [[gnu::always_inline]] inline unsigned int Framework::PackTripletId(unsigned int iStation, unsigned int iTriplet)
+ {
+#ifndef FAST_CODE
+ assert(iStation < constants::size::MaxNstations);
+ assert(iTriplet < constants::size::MaxNtriplets);
+#endif
+ constexpr unsigned int kMoveStation = constants::size::TripletBits;
+ return (iStation << kMoveStation) + iTriplet;
+ }
+
+ // ---------------------------------------------------------------------------------------------------------------------
+ //
+ [[gnu::always_inline]] inline unsigned int Framework::TripletId2Station(unsigned int id)
+ {
+ constexpr unsigned int kMoveStation = constants::size::TripletBits;
+ return id >> kMoveStation;
+ }
+
+ // ---------------------------------------------------------------------------------------------------------------------
+ //
+ [[gnu::always_inline]] inline unsigned int Framework::TripletId2Triplet(unsigned int id)
+ {
+ constexpr unsigned int kTripletMask = (1u << constants::size::TripletBits) - 1u;
+ return id & kTripletMask;
+ }
+
+} // namespace cbm::algo::ca
diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/algo/ca/core/tracking/CaTrackExtender.cxx
similarity index 60%
rename from reco/L1/L1Algo/L1BranchExtender.cxx
rename to algo/ca/core/tracking/CaTrackExtender.cxx
index fb0490dd0d..0bfbb98c4f 100644
--- a/reco/L1/L1Algo/L1BranchExtender.cxx
+++ b/algo/ca/core/tracking/CaTrackExtender.cxx
@@ -2,34 +2,45 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Igor Kulakov [committer], Maksym Zyzak, Valentina Akishina */
+#include "CaTrackExtender.h"
+
#include <iostream>
#include "CaBranch.h"
#include "CaDefines.h"
+#include "CaFramework.h"
#include "CaGridArea.h"
#include "CaTrack.h"
#include "CaTrackFit.h"
#include "CaTrackParam.h"
#include "CaVector.h"
-#include "L1Algo.h"
// using namespace std;
using cbm::algo::ca::Vector; // TMP!!
using std::cout;
using std::endl;
-/// Fit track
-/// t - track with hits
-/// T - track params
-/// qp0 - momentum for extrapolation
-/// initialize - should be params ititialized. 1 - yes.
-void L1Algo::BranchFitterFast(const ca::Branch& t, TrackParamV& Tout, const bool upstream, const fvec qp0,
- const bool initParams)
+using namespace cbm::algo::ca;
+
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+TrackExtender::TrackExtender(const ca::Framework& algo) : frAlgo(algo) {}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+TrackExtender::~TrackExtender() {}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+
+void TrackExtender::FitBranchFast(const ca::Branch& t, TrackParamV& Tout, const bool upstream, const fvec qp0,
+ const bool initParams)
{
CBMCA_DEBUG_ASSERT(t.NHits >= 3);
ca::TrackFit fit;
- fit.SetParticleMass(GetDefaultParticleMass());
+ fit.SetParticleMass(frAlgo.GetDefaultParticleMass());
fit.SetMask(fmask::One());
fit.SetTrack(Tout);
TrackParamV& T = fit.Tr();
@@ -42,17 +53,17 @@ void L1Algo::BranchFitterFast(const ca::Branch& t, TrackParamV& Tout, const bool
const int iFirstHit = (upstream) ? nHits - 1 : 0;
const int iLastHit = (upstream) ? 0 : nHits - 1;
- const ca::Hit& hit0 = fInputData.GetHit(hits[iFirstHit]);
- const ca::Hit& hit1 = fInputData.GetHit(hits[iFirstHit + step]);
- const ca::Hit& hit2 = fInputData.GetHit(hits[iFirstHit + 2 * step]);
+ const ca::Hit& hit0 = frAlgo.GetInputData().GetHit(hits[iFirstHit]);
+ const ca::Hit& hit1 = frAlgo.GetInputData().GetHit(hits[iFirstHit + step]);
+ const ca::Hit& hit2 = frAlgo.GetInputData().GetHit(hits[iFirstHit + 2 * step]);
int ista0 = hit0.Station();
int ista1 = hit1.Station();
int ista2 = hit2.Station();
- const ca::Station& sta0 = fParameters.GetStation(ista0);
- const ca::Station& sta1 = fParameters.GetStation(ista1);
- const ca::Station& sta2 = fParameters.GetStation(ista2);
+ const ca::Station& sta0 = frAlgo.GetParameters().GetStation(ista0);
+ const ca::Station& sta1 = frAlgo.GetParameters().GetStation(ista1);
+ const ca::Station& sta2 = frAlgo.GetParameters().GetStation(ista2);
fvec x0 = hit0.X();
fvec y0 = hit0.Y();
@@ -100,13 +111,13 @@ void L1Algo::BranchFitterFast(const ca::Branch& t, TrackParamV& Tout, const bool
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
for (int i = iFirstHit + step; step * i <= step * iLastHit; i += step) {
- const ca::Hit& hit = fInputData.GetHit(hits[i]);
+ const ca::Hit& hit = frAlgo.GetInputData().GetHit(hits[i]);
int ista = hit.Station();
- const ca::Station& sta = fParameters.GetStation(ista);
+ const ca::Station& sta = frAlgo.GetParameters().GetStation(ista);
fit.Extrapolate(hit.Z(), fld);
fit.FilterHit(sta, hit);
- auto radThick = fParameters.GetMaterialThickness(ista, T.X(), T.Y());
+ auto radThick = frAlgo.GetParameters().GetMaterialThickness(ista, T.X(), T.Y());
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, fvec(upstream ? 1. : -1.));
@@ -120,50 +131,46 @@ void L1Algo::BranchFitterFast(const ca::Branch& t, TrackParamV& Tout, const bool
} // i
Tout = T;
-} // void L1Algo::BranchFitterFast
+} // void ca::Framework::BranchFitterFast
/// like BranchFitterFast but more precise
-void L1Algo::BranchFitter(const ca::Branch& t, TrackParamV& T, const bool upstream, const fvec qp0,
- const bool initParams)
+void TrackExtender::FitBranch(const ca::Branch& t, TrackParamV& T, const bool upstream, const fvec qp0,
+ const bool initParams)
{
- BranchFitterFast(t, T, upstream, qp0, initParams);
+ FitBranchFast(t, T, upstream, qp0, initParams);
for (int i = 0; i < 1; i++) {
- BranchFitterFast(t, T, !upstream, T.Qp(), false);
- BranchFitterFast(t, T, upstream, T.Qp(), false);
+ FitBranchFast(t, T, !upstream, T.Qp(), false);
+ FitBranchFast(t, T, upstream, T.Qp(), false);
}
-} // void L1Algo::BranchFitter
-
-/// Find additional hits for existing track
-/// t - track with hits
-/// T - track params
-/// qp0 - momentum for extrapolation
-/// initialize - should be params ititialized. 1 - yes.
-void L1Algo::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool upstream, const fvec qp0)
+} // void ca::Framework::BranchFitter
+
+
+void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool upstream, const fvec qp0)
{
- Vector<ca::HitIndex_t> newHits {"L1TrackExtender::newHits"};
- newHits.reserve(fParameters.GetNstationsActive());
+ Vector<ca::HitIndex_t> newHits {"ca::TrackExtender::newHits"};
+ newHits.reserve(frAlgo.GetParameters().GetNstationsActive());
ca::TrackFit fit;
- fit.SetParticleMass(GetDefaultParticleMass());
+ fit.SetParticleMass(frAlgo.GetDefaultParticleMass());
fit.SetMask(fmask::One());
fit.SetTrack(Tout);
fit.SetQp0(qp0);
const signed short int step = -2 * static_cast<int>(upstream) + 1; // increment for station index
const int iFirstHit = (upstream) ? 2 : t.NofHits() - 3;
- // int ista = fInputData.GetHit(t.Hits[iFirstHit]).iSt + 2 * step; // current station. set to the end of track
+ // int ista = frAlgo.GetInputData().GetHit(t.Hits[iFirstHit]).iSt + 2 * step; // current station. set to the end of track
- const ca::Hit& hit0 = fInputData.GetHit(t.Hits()[iFirstHit]); // optimize
- const ca::Hit& hit1 = fInputData.GetHit(t.Hits()[iFirstHit + step]);
- const ca::Hit& hit2 = fInputData.GetHit(t.Hits()[iFirstHit + 2 * step]);
+ const ca::Hit& hit0 = frAlgo.GetInputData().GetHit(t.Hits()[iFirstHit]); // optimize
+ const ca::Hit& hit1 = frAlgo.GetInputData().GetHit(t.Hits()[iFirstHit + step]);
+ const ca::Hit& hit2 = frAlgo.GetInputData().GetHit(t.Hits()[iFirstHit + 2 * step]);
const int ista0 = hit0.Station();
const int ista1 = hit1.Station();
const int ista2 = hit2.Station();
- const ca::Station& sta0 = fParameters.GetStation(ista0);
- const ca::Station& sta1 = fParameters.GetStation(ista1);
- const ca::Station& sta2 = fParameters.GetStation(ista2);
+ const ca::Station& sta0 = frAlgo.GetParameters().GetStation(ista0);
+ const ca::Station& sta1 = frAlgo.GetParameters().GetStation(ista1);
+ const ca::Station& sta2 = frAlgo.GetParameters().GetStation(ista2);
fvec x0 = hit0.X();
fvec y0 = hit0.Y();
@@ -187,13 +194,13 @@ void L1Algo::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool upstream,
int ista = ista2 + 2 * step; // skip one station. if there would be hit it has to be found on previous stap
- if (ista2 == fFirstCAstation) ista = ista2 + step;
+ if (ista2 == frAlgo.fFirstCAstation) ista = ista2 + step;
- const fvec pickGather2 = fPickGather * fPickGather;
+ const fvec pickGather2 = frAlgo.fPickGather * frAlgo.fPickGather;
- for (; (ista < fParameters.GetNstationsActive()) && (ista >= 0); ista += step) { // CHECKME why ista2?
+ for (; (ista < frAlgo.GetParameters().GetNstationsActive()) && (ista >= 0); ista += step) { // CHECKME why ista2?
- const ca::Station& sta = fParameters.GetStation(ista);
+ const ca::Station& sta = frAlgo.GetParameters().GetStation(ista);
fit.Extrapolate(sta.fZ, fld);
@@ -202,18 +209,19 @@ void L1Algo::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool upstream,
TrackParamV& tr = fit.Tr();
- ca::GridArea area(vGrid[ista], tr.X()[0], tr.Y()[0],
- (sqrt(fPickGather * tr.C00()) + vGrid[ista].GetMaxRangeX() + fMaxDZ * abs(tr.Tx()))[0],
- (sqrt(fPickGather * tr.C11()) + vGrid[ista].GetMaxRangeY() + fMaxDZ * abs(tr.Ty()))[0]);
+ ca::GridArea area(
+ frAlgo.vGrid[ista], tr.X()[0], tr.Y()[0],
+ (sqrt(frAlgo.fPickGather * tr.C00()) + frAlgo.vGrid[ista].GetMaxRangeX() + frAlgo.fMaxDZ * abs(tr.Tx()))[0],
+ (sqrt(frAlgo.fPickGather * tr.C11()) + frAlgo.vGrid[ista].GetMaxRangeY() + frAlgo.fMaxDZ * abs(tr.Ty()))[0]);
- if (fParameters.DevIsIgnoreHitSearchAreas()) { area.DoLoopOverEntireGrid(); }
+ if (frAlgo.GetParameters().DevIsIgnoreHitSearchAreas()) { area.DoLoopOverEntireGrid(); }
ca::HitIndex_t ih = 0;
while (area.GetNextObjectId(ih)) { // loop over the hits in the area
- if (fIsWindowHitSuppressed[ih]) { continue; }
- const ca::Hit& hit = fWindowHits[ih];
+ if (frAlgo.fIsWindowHitSuppressed[ih]) { continue; }
+ const ca::Hit& hit = frAlgo.fWindowHits[ih];
if (sta.timeInfo && tr.NdfTime()[0] > -2.) {
fscal dt = hit.T() - tr.Time()[0];
@@ -221,10 +229,8 @@ void L1Algo::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool upstream,
}
//if (GetFUsed((*fStripFlag)[hit.f] | (*fStripFlag)[hit.b])) continue; // if used
- //L1_SHOW(fvHitKeyFlags.size());
- //L1_SHOW(hit.f);
- //L1_SHOW(hit.b);
- if (fvHitKeyFlags[hit.FrontKey()] || fvHitKeyFlags[hit.BackKey()]) { continue; }
+
+ if (frAlgo.fvHitKeyFlags[hit.FrontKey()] || frAlgo.fvHitKeyFlags[hit.BackKey()]) { continue; }
auto [y, C11] = fit.ExtrapolateLineYdY2(hit.Z());
@@ -239,7 +245,7 @@ void L1Algo::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool upstream,
fscal d2 = d_x * d_x + d_y * d_y;
if (d2 > r2_best) continue;
- fscal dxm_est = sqrt(pickGather2 * C00)[0] + vGrid[ista].GetMaxRangeX();
+ fscal dxm_est = sqrt(pickGather2 * C00)[0] + frAlgo.vGrid[ista].GetMaxRangeX();
if (fabs(d_x) > dxm_est) continue;
r2_best = d2;
@@ -248,13 +254,13 @@ void L1Algo::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool upstream,
if (iHit_best < 0) break;
- newHits.push_back(fWindowHits[iHit_best].Id());
+ newHits.push_back(frAlgo.fWindowHits[iHit_best].Id());
- const ca::Hit& hit = fWindowHits[iHit_best];
+ const ca::Hit& hit = frAlgo.fWindowHits[iHit_best];
fit.Extrapolate(hit.Z(), fld);
fit.FilterHit(sta, hit);
- auto radThick = fParameters.GetMaterialThickness(ista, tr.X(), tr.Y());
+ auto radThick = frAlgo.GetParameters().GetMaterialThickness(ista, tr.X(), tr.Y());
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, upstream ? fvec(1.f) : fvec(-1.f));
@@ -288,10 +294,10 @@ void L1Algo::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool upstream,
Tout = fit.Tr();
-} // void L1Algo::FindMoreHits
+} // void ca::Framework::FindMoreHits
/// Try to extrapolate and find additional hits on other stations
-fscal L1Algo::BranchExtender(ca::Branch& t) // TODO Simdize
+fscal TrackExtender::ExtendBranch(ca::Branch& t)
{
// const unsigned int minNHits = 3;
@@ -300,15 +306,13 @@ fscal L1Algo::BranchExtender(ca::Branch& t) // TODO Simdize
// downstream
bool upstream = 0;
- BranchFitter(t, T, upstream);
- // BranchFitterFast (t, T, upstream, 0, true);
+ FitBranch(t, T, upstream, 0.0);
- // if (t.NHits < minNHits) return T.chi2[0];
FindMoreHits(t, T, upstream, T.Qp());
// upstream
upstream = 1;
- BranchFitterFast(t, T, upstream, T.Qp(), false);
+ FitBranchFast(t, T, upstream, T.Qp(), false);
FindMoreHits(t, T, upstream, T.Qp());
diff --git a/algo/ca/core/tracking/CaTrackExtender.h b/algo/ca/core/tracking/CaTrackExtender.h
new file mode 100644
index 0000000000..3e9ee941d6
--- /dev/null
+++ b/algo/ca/core/tracking/CaTrackExtender.h
@@ -0,0 +1,93 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer], Maksym Zyzak */
+
+/// \file TrackExtender.h
+/// \authors S.Zharko <s.zharko@gsi.de> (interface), M.Zyzak (original algorithm)
+/// \brief A class wrapper over clones merger algorithm for the CA track finder (declaration)
+/// \since 22.07.2022
+
+#pragma once // include this header only once per compilation unit
+
+#include "CaBranch.h"
+#include "CaHit.h"
+#include "CaSimd.h"
+#include "CaTrackParam.h"
+#include "CaVector.h"
+
+
+namespace cbm::algo::ca
+{
+ class Track;
+ class Framework;
+
+ namespace
+ {
+ using namespace cbm::algo; // to keep ca:: prefices in the code
+ }
+
+
+ /// Class implements a clones merger algorithm for the CA track finder
+ ///
+ class TrackExtender {
+ public:
+ /// Default constructor
+ TrackExtender(const ca::Framework& algo);
+
+ /// Destructor
+ ~TrackExtender();
+
+ /// Copy constructor
+ TrackExtender(const TrackExtender&) = delete;
+
+ /// Move constructor
+ TrackExtender(TrackExtender&&) = delete;
+
+ /// Copy assignment operator
+ TrackExtender& operator=(const TrackExtender&) = delete;
+
+ /// Move assignment operator
+ TrackExtender& operator=(TrackExtender&&) = delete;
+
+
+ /// Find additional hits for existing track both downstream and upstream
+ /// \return chi2
+ fscal ExtendBranch(ca::Branch& t);
+
+ private:
+ ///-------------------------------
+ /// Private methods
+
+ /// Finds additional hits for already found track in the given direction
+ /// \param t - track with hits
+ /// \param T - track params
+ /// \param dir - 0 - forward, 1 - backward
+ /// \param qp0 - momentum value to linearize the extrapolation
+ void FindMoreHits(ca::Branch& t, TrackParamV& T, const bool dir, const fvec qp0);
+
+ /// Fits the branch. Does few passes over the hits.
+ /// \param t - track branch with hits
+ /// \param T - track parameters
+ /// \param dir - false - forward, true - backward
+ /// \param qp0 - momentum value to linearize the extrapolation
+ /// \param initParams - should params be ititialized. 1 - yes.
+ void FitBranch(const ca::Branch& t, TrackParamV& T, const bool dir, const fvec qp0, const bool initParams = true);
+
+
+ /// Fits the branch. Does only one pass over the hits.
+ /// \param t - track branch with hits
+ /// \param T - track parameters
+ /// \param dir - false - forward, true - backward
+ /// \param qp0 - momentum value to linearize the extrapolation
+ /// \param initParams - should params be ititialized. 1 - yes.
+ void FitBranchFast(const ca::Branch& t, TrackParamV& T, const bool dir, const fvec qp0,
+ const bool initParams = true);
+
+ private:
+ ///-------------------------------
+ /// Data members
+
+ const ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class
+ };
+
+} // namespace cbm::algo::ca
diff --git a/reco/L1/L1Algo/L1CaTrackFinder.cxx b/algo/ca/core/tracking/CaTrackFinder.cxx
similarity index 53%
rename from reco/L1/L1Algo/L1CaTrackFinder.cxx
rename to algo/ca/core/tracking/CaTrackFinder.cxx
index fc80be5d5d..5e9732bd68 100644
--- a/reco/L1/L1Algo/L1CaTrackFinder.cxx
+++ b/algo/ca/core/tracking/CaTrackFinder.cxx
@@ -20,14 +20,27 @@
#include <chrono>
+#include "CaFramework.h"
#include "CaTrack.h"
-#include "L1Algo.h"
+
+using namespace cbm::algo::ca;
using cbm::algo::ca::ECounter;
using cbm::algo::ca::ETimer;
using cbm::algo::ca::Track;
-void L1Algo::CaTrackFinder()
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+TrackFinder::TrackFinder(ca::Framework& algo) : frAlgo(algo) {}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+TrackFinder::~TrackFinder() {}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TrackFinder::FindTracks()
{
//
// The main CaTrackFinder routine
@@ -36,60 +49,59 @@ void L1Algo::CaTrackFinder()
//
// Reset monitor
- fMonitor.Reset();
- fMonitor.IncrementCounter(ECounter::RecoHit, fInputData.GetNhits());
- fMonitor.StartTimer(ETimer::Tracking);
+ frAlgo.fMonitor.Reset();
+ frAlgo.fMonitor.IncrementCounter(ECounter::RecoHit, frAlgo.fInputData.GetNhits());
+ frAlgo.fMonitor.StartTimer(ETimer::Tracking);
auto timerStart = std::chrono::high_resolution_clock::now();
-
// ----- Reset data arrays -------------------------------------------------------------------------------------------
- fvHitKeyFlags.reset(fInputData.GetNhitKeys(), 0);
+ frAlgo.fvHitKeyFlags.reset(frAlgo.fInputData.GetNhitKeys(), 0);
- fRecoTracks.clear();
- fRecoHits.clear();
- fRecoHits.reserve(2 * fInputData.GetNhits());
- fRecoTracks.reserve(2 * fInputData.GetNhits() / fParameters.GetNstationsActive());
+ frAlgo.fRecoTracks.clear();
+ frAlgo.fRecoHits.clear();
+ frAlgo.fRecoHits.reserve(2 * frAlgo.fInputData.GetNhits());
+ frAlgo.fRecoTracks.reserve(2 * frAlgo.fInputData.GetNhits() / frAlgo.GetParameters().GetNstationsActive());
- for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
- fSliceHitIds[iS].clear();
- fSliceHitIds[iS].reserve(fInputData.GetNhits());
+ for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); ++iS) {
+ frAlgo.fSliceHitIds[iS].clear();
+ frAlgo.fSliceHitIds[iS].reserve(frAlgo.fInputData.GetNhits());
}
- fvHitKeyFlags.reset(fInputData.GetNhitKeys(), 0);
- fHitTimeInfo.reset(fInputData.GetNhits());
+ frAlgo.fvHitKeyFlags.reset(frAlgo.fInputData.GetNhitKeys(), 0);
+ frAlgo.fHitTimeInfo.reset(frAlgo.fInputData.GetNhits());
- // TODO: move these values to CbmL1Parameters namespace (S.Zharko)
+ // TODO: move these values to Parameters namespace (S.Zharko)
const fscal tsLength = 10000; // length of sub-TS
const fscal tsOverlap = 15; // min length of overlap region
fscal minProtonMomentum = 0.1;
- fscal targX = fParameters.GetTargetPositionX()[0];
- fscal targY = fParameters.GetTargetPositionY()[0];
- fscal targZ = fParameters.GetTargetPositionZ()[0];
+ fscal targX = frAlgo.GetParameters().GetTargetPositionX()[0];
+ fscal targY = frAlgo.GetParameters().GetTargetPositionY()[0];
+ fscal targZ = frAlgo.GetParameters().GetTargetPositionZ()[0];
fscal tsStart = std::numeric_limits<fscal>::max(); // starting time of sub-TS
- // calculate possible event time for the hits (fHitTimeInfo array)
+ // calculate possible event time for the hits (frAlgo.fHitTimeInfo array)
- const int nDataStreams = fInputData.GetNdataStreams();
+ const int nDataStreams = frAlgo.fInputData.GetNdataStreams();
for (int iStream = 0; iStream < nDataStreams; ++iStream) {
fscal maxTimeBeforeHit = std::numeric_limits<fscal>::lowest();
- int nStreamHits = fInputData.GetStreamNhits(iStream);
+ int nStreamHits = frAlgo.fInputData.GetStreamNhits(iStream);
for (int ih = 0; ih < nStreamHits; ++ih) {
- ca::HitIndex_t caHitId = fInputData.GetStreamStartIndex(iStream) + ih;
- const ca::Hit& h = fInputData.GetHit(caHitId);
- const ca::Station& st = fParameters.GetStation(h.Station());
+ ca::HitIndex_t caHitId = frAlgo.fInputData.GetStreamStartIndex(iStream) + ih;
+ const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId);
+ const ca::Station& st = frAlgo.GetParameters().GetStation(h.Station());
- fscal dx = h.X() - targX;
- fscal dy = h.Y() - targY;
- fscal dz = h.Z() - targZ;
- fscal l = sqrt(dx * dx + dy * dy + dz * dz);
+ fscal dx = h.X() - targX;
+ fscal dy = h.Y() - targY;
+ fscal dz = h.Z() - targZ;
+ fscal l = sqrt(dx * dx + dy * dy + dz * dz);
fscal timeOfFlightMin = l * constants::phys::SpeedOfLightInv;
fscal timeOfFlightMax =
@@ -98,7 +110,7 @@ void L1Algo::CaTrackFinder()
* constants::phys::SpeedOfLightInvD;
fscal dt = h.RangeT();
- L1HitTimeInfo& info = fHitTimeInfo[caHitId];
+ CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId];
info.fEventTimeMin = h.T() - dt - timeOfFlightMax;
info.fEventTimeMax = h.T() + dt - timeOfFlightMin;
@@ -116,11 +128,11 @@ void L1Algo::CaTrackFinder()
}
fscal minTimeAfterHit = std::numeric_limits<fscal>::max();
- // loop in the reverse order to fill L1HitTimeInfo::fMinTimeAfterHit fields
+ // loop in the reverse order to fill CaHitTimeInfo::fMinTimeAfterHit fields
for (int ih = nStreamHits - 1; ih >= 0; --ih) {
- ca::HitIndex_t caHitId = fInputData.GetStreamStartIndex(iStream) + ih;
- L1HitTimeInfo& info = fHitTimeInfo[caHitId];
+ ca::HitIndex_t caHitId = frAlgo.fInputData.GetStreamStartIndex(iStream) + ih;
+ CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId];
if (minTimeAfterHit > info.fEventTimeMin) { minTimeAfterHit = info.fEventTimeMin; }
info.fMinTimeAfterHit = minTimeAfterHit;
}
@@ -133,15 +145,16 @@ void L1Algo::CaTrackFinder()
ca::HitIndex_t sliceFirstHit[nDataStreams];
for (int iStream = 0; iStream < nDataStreams; ++iStream) {
- sliceFirstHit[iStream] = fInputData.GetStreamStartIndex(iStream);
+ sliceFirstHit[iStream] = frAlgo.fInputData.GetStreamStartIndex(iStream);
}
while (areDataLeft) {
- fMonitor.IncrementCounter(ECounter::SubTS);
+
+ frAlgo.fMonitor.IncrementCounter(ECounter::SubTS);
// select the sub-slice hits
- for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
- fSliceHitIds[iS].clear();
+ for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); ++iS) {
+ frAlgo.fSliceHitIds[iS].clear();
}
areDataLeft = false;
@@ -150,11 +163,12 @@ void L1Algo::CaTrackFinder()
for (int iStream = 0; iStream < nDataStreams; ++iStream) {
- for (ca::HitIndex_t caHitId = sliceFirstHit[iStream]; caHitId < fInputData.GetStreamStopIndex(iStream);
+ for (ca::HitIndex_t caHitId = sliceFirstHit[iStream]; caHitId < frAlgo.fInputData.GetStreamStopIndex(iStream);
++caHitId) {
- L1HitTimeInfo& info = fHitTimeInfo[caHitId];
- const ca::Hit& h = fInputData.GetHit(caHitId);
- if (fvHitKeyFlags[h.FrontKey()] || fvHitKeyFlags[h.BackKey()]) { // the hit is already reconstructed
+ CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId];
+ const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId);
+ if (frAlgo.fvHitKeyFlags[h.FrontKey()]
+ || frAlgo.fvHitKeyFlags[h.BackKey()]) { // the hit is already reconstructed
continue;
}
if (info.fEventTimeMin > tsStart + tsLength) { // the hit is too late for the sub slice
@@ -166,7 +180,7 @@ void L1Algo::CaTrackFinder()
}
else {
if (tsStart <= info.fEventTimeMax) { // the hit belongs to the sub-slice
- fSliceHitIds[h.Station()].push_back(caHitId);
+ frAlgo.fSliceHitIds[h.Station()].push_back(caHitId);
if (info.fMaxTimeBeforeHit < tsStart + tsLength - tsOverlap) {
// this hit and all hits before are before the overlap
sliceFirstHit[iStream] = caHitId + 1;
@@ -176,14 +190,13 @@ void L1Algo::CaTrackFinder()
}
}
- fMonitor.StartTimer(ETimer::TrackFinder);
- CaTrackFinderSlice();
- fMonitor.StopTimer(ETimer::TrackFinder);
-
+ frAlgo.fMonitor.StartTimer(ETimer::TrackFinder);
+ frAlgo.fTrackFinderWindow.CaTrackFinderSlice();
+ frAlgo.fMonitor.StopTimer(ETimer::TrackFinder);
- fMonitor.StartTimer(ETimer::TrackFitter);
- L1KFTrackFitter();
- fMonitor.StopTimer(ETimer::TrackFitter);
+ frAlgo.fMonitor.StartTimer(ETimer::TrackFitter);
+ frAlgo.fTrackFitter.FitCaTracks();
+ frAlgo.fMonitor.StopTimer(ETimer::TrackFitter);
// save reconstructed tracks with no hits in the overlap region
@@ -194,11 +207,12 @@ void L1Algo::CaTrackFinder()
// TODO: only add those hits from the region before tsStartNew that belong to the not stored tracks
int trackFirstHit = 0;
- for (auto it = fSliceRecoTracks.begin(); it != fSliceRecoTracks.end(); trackFirstHit += it->fNofHits, it++) {
+ for (auto it = frAlgo.fSliceRecoTracks.begin(); it != frAlgo.fSliceRecoTracks.end();
+ trackFirstHit += it->fNofHits, it++) {
bool isTrackCompletelyInOverlap = true;
for (int ih = 0; ih < it->fNofHits; ih++) {
- int caHitId = fSliceRecoHits[trackFirstHit + ih];
- L1HitTimeInfo& info = fHitTimeInfo[caHitId];
+ int caHitId = frAlgo.fSliceRecoHits[trackFirstHit + ih];
+ CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId];
if (info.fEventTimeMax < tsStart) { // this hit is before the overlap
isTrackCompletelyInOverlap = false;
break;
@@ -216,21 +230,21 @@ void L1Algo::CaTrackFinder()
// release the track hits
for (int i = 0; i < it->fNofHits; i++) {
- int caHitId = fSliceRecoHits[trackFirstHit + i];
- const auto& h = fInputData.GetHit(caHitId);
- fvHitKeyFlags[h.FrontKey()] = 0;
- fvHitKeyFlags[h.BackKey()] = 0;
+ int caHitId = frAlgo.fSliceRecoHits[trackFirstHit + i];
+ const auto& h = frAlgo.fInputData.GetHit(caHitId);
+ frAlgo.fvHitKeyFlags[h.FrontKey()] = 0;
+ frAlgo.fvHitKeyFlags[h.BackKey()] = 0;
}
}
else { // save the track
- fRecoTracks.push_back(*it);
+ frAlgo.fRecoTracks.push_back(*it);
// mark the track hits as used
for (int i = 0; i < it->fNofHits; i++) {
- int caHitId = fSliceRecoHits[trackFirstHit + i];
- const auto& h = fInputData.GetHit(caHitId);
- fvHitKeyFlags[h.FrontKey()] = 1;
- fvHitKeyFlags[h.BackKey()] = 1;
- fRecoHits.push_back(caHitId);
+ int caHitId = frAlgo.fSliceRecoHits[trackFirstHit + i];
+ const auto& h = frAlgo.fInputData.GetHit(caHitId);
+ frAlgo.fvHitKeyFlags[h.FrontKey()] = 1;
+ frAlgo.fvHitKeyFlags[h.BackKey()] = 1;
+ frAlgo.fRecoHits.push_back(caHitId);
}
}
} // sub-timeslice tracks
@@ -238,10 +252,10 @@ void L1Algo::CaTrackFinder()
tsStart -= 5; // do 5 ns margin
}
- auto timerEnd = std::chrono::high_resolution_clock::now();
- fMonitor.StopTimer(ETimer::Tracking);
- fMonitor.IncrementCounter(ECounter::RecoTrack, fRecoTracks.size());
- fMonitor.IncrementCounter(ECounter::RecoHitUsed, fRecoHits.size());
+ frAlgo.fMonitor.StopTimer(ETimer::Tracking);
+ frAlgo.fMonitor.IncrementCounter(ECounter::RecoTrack, frAlgo.fRecoTracks.size());
+ frAlgo.fMonitor.IncrementCounter(ECounter::RecoHitUsed, frAlgo.fRecoHits.size());
- fCaRecoTime = (double) (std::chrono::duration<double>(timerEnd - timerStart).count());
+ auto timerEnd = std::chrono::high_resolution_clock::now();
+ frAlgo.fCaRecoTime = (double) (std::chrono::duration<double>(timerEnd - timerStart).count());
}
diff --git a/algo/ca/core/tracking/CaTrackFinder.h b/algo/ca/core/tracking/CaTrackFinder.h
new file mode 100644
index 0000000000..4e2db5dccc
--- /dev/null
+++ b/algo/ca/core/tracking/CaTrackFinder.h
@@ -0,0 +1,66 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer], Maksym Zyzak */
+
+/// \file CaTrackFinder.h
+
+#pragma once // include this header only once per compilation unit
+
+#include "CaBranch.h"
+#include "CaHit.h"
+#include "CaSimd.h"
+#include "CaTrackParam.h"
+#include "CaVector.h"
+
+
+namespace cbm::algo::ca
+{
+
+ class Track;
+ class Triplet;
+ class Framework;
+
+ namespace
+ {
+ using namespace cbm::algo; // to keep ca:: prefices in the code
+ }
+
+ /// Class implements a clones merger algorithm for the CA track finder
+ ///
+ class TrackFinder {
+ public:
+ /// Default constructor
+ TrackFinder(ca::Framework& algo);
+
+ /// Destructor
+ ~TrackFinder();
+
+ /// Copy constructor
+ TrackFinder(const TrackFinder&) = delete;
+
+ /// Move constructor
+ TrackFinder(TrackFinder&&) = delete;
+
+ /// Copy assignment operator
+ TrackFinder& operator=(const TrackFinder&) = delete;
+
+ /// Move assignment operator
+ TrackFinder& operator=(TrackFinder&&) = delete;
+
+ void FindTracks();
+
+ private:
+ ///-------------------------------
+ /// Private methods
+
+ bool checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2) const;
+
+
+ private:
+ ///-------------------------------
+ /// Data members
+
+ ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class
+ };
+
+} // namespace cbm::algo::ca
diff --git a/algo/ca/core/tracking/CaTrackFinderWindow.cxx b/algo/ca/core/tracking/CaTrackFinderWindow.cxx
new file mode 100644
index 0000000000..25b4dcdec6
--- /dev/null
+++ b/algo/ca/core/tracking/CaTrackFinderWindow.cxx
@@ -0,0 +1,771 @@
+/* Copyright (C) 2009-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Valentina Akishina, Ivan Kisel, Sergey Gorbunov [committer], Igor Kulakov, Maksym Zyzak */
+
+/*
+ *=====================================================
+ *
+ * CBM Level 1 4D Reconstruction
+ *
+ * Authors: V.Akishina, I.Kisel, S.Gorbunov, I.Kulakov, M.Zyzak
+ * Documentation: V.Akishina
+ *
+ * e-mail : v.akishina@gsi.de
+ *
+ *=====================================================
+ *
+ * Finds tracks using the Cellular Automaton algorithm
+ *
+ */
+
+
+#include "CaTrackFinderWindow.h"
+
+#include <algorithm>
+#include <fstream>
+#include <iostream>
+#include <list>
+#include <map>
+
+#include <stdio.h>
+
+#include "CaBranch.h"
+#include "CaFramework.h"
+#include "CaGrid.h"
+#include "CaGridEntry.h"
+#include "CaTrack.h"
+#include "CaTrackParam.h"
+#include "CaTripletConstructor.h"
+
+// #include "CaToolsDebugger.h"
+
+
+// using namespace std;
+using std::cout;
+using std::endl;
+using Track = cbm::algo::ca::Track;
+
+using namespace cbm::algo::ca;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+TrackFinderWindow::TrackFinderWindow(ca::Framework& algo) : frAlgo(algo) {}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+TrackFinderWindow::~TrackFinderWindow() {}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+
+
+bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2) const
+{
+ dchi2 = 1.e20;
+
+ if (r.GetMHit() != l.GetRHit()) return false;
+ if (r.GetLHit() != l.GetMHit()) return false;
+
+ if (r.GetMSta() != l.GetRSta()) return false;
+ if (r.GetLSta() != l.GetMSta()) return false;
+
+
+ if (r.IsMomentumFitted()) {
+ assert(l.IsMomentumFitted());
+
+ fscal dqp = l.GetQp() - r.GetQp();
+ fscal Cqp = l.GetCqp() + r.GetCqp();
+
+ if (!std::isfinite(dqp)) return false;
+ if (!std::isfinite(Cqp)) return false;
+
+ if (dqp * dqp > frAlgo.fTripletLinkChi2 * Cqp) {
+ return false; // bad neighbour // CHECKME why do we need recheck it?? (it really change result)
+ }
+ dchi2 = dqp * dqp / Cqp;
+ }
+ else {
+ fscal dtx = l.GetTx() - r.GetTx();
+ fscal Ctx = l.GetCtx() + r.GetCtx();
+
+ fscal dty = l.GetTy() - r.GetTy();
+ fscal Cty = l.GetCty() + r.GetCty();
+
+ // it shouldn't happen, but happens sometimes
+
+ if (!std::isfinite(dtx)) return false;
+ if (!std::isfinite(dty)) return false;
+ if (!std::isfinite(Ctx)) return false;
+ if (!std::isfinite(Cty)) return false;
+
+ if (dty * dty > frAlgo.fTripletLinkChi2 * Cty) return false;
+ if (dtx * dtx > frAlgo.fTripletLinkChi2 * Ctx) return false;
+
+ //dchi2 = 0.5f * (dtx * dtx / Ctx + dty * dty / Cty);
+ dchi2 = 0.;
+ }
+
+ if (!std::isfinite(dchi2)) return false;
+
+ return true;
+}
+
+
+// **************************************************************************************************
+// * *
+// * ------ CATrackFinder procedure ------ *
+// * *
+// **************************************************************************************************
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TrackFinderWindow::ReadWindowData()
+{
+ int nHits = 0;
+ for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); iS++) {
+ frAlgo.fStationHitsStartIndex[iS] = nHits;
+ frAlgo.fStationNhits[iS] = frAlgo.fSliceHitIds[iS].size();
+ nHits += frAlgo.fStationNhits[iS];
+ }
+ frAlgo.fStationHitsStartIndex[frAlgo.GetParameters().GetNstationsActive()] = nHits;
+
+ frAlgo.fWindowHits.reset(nHits);
+ frAlgo.fIsWindowHitSuppressed.reset(nHits, 0);
+
+ for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); iS++) {
+ for (ca::HitIndex_t ih = 0; ih < frAlgo.fSliceHitIds[iS].size(); ++ih) {
+ ca::Hit h = frAlgo.fInputData.GetHit(frAlgo.fSliceHitIds[iS][ih]);
+ h.SetId(frAlgo.fSliceHitIds[iS][ih]);
+
+ frAlgo.fWindowHits[frAlgo.fStationHitsStartIndex[iS] + ih] = h;
+ }
+ }
+
+ frAlgo.fHitFirstTriplet.reset(nHits, 0);
+ frAlgo.fHitNtriplets.reset(nHits, 0);
+
+ // TODO: SG: reduce the array size to the number of keys in subslice
+
+ frAlgo.fStripToTrack.reset(frAlgo.fInputData.GetNhitKeys(), -1);
+
+ frAlgo.fSliceRecoTracks.clear();
+ frAlgo.fSliceRecoTracks.reserve(2 * nHits / frAlgo.GetParameters().GetNstationsActive());
+
+ frAlgo.fSliceRecoHits.clear();
+ frAlgo.fSliceRecoHits.reserve(2 * nHits);
+
+ frAlgo.fTrackCandidates.clear();
+ frAlgo.fTrackCandidates.reserve(nHits / 10);
+ for (unsigned int iS = 0; iS < constants::size::MaxNstations; iS++) {
+ int nHitsStation = frAlgo.fSliceHitIds[iS].size();
+ frAlgo.fTriplets[iS].clear();
+ frAlgo.fTriplets[iS].reserve(2 * nHitsStation);
+ }
+}
+
+void TrackFinderWindow::CaTrackFinderSlice()
+{
+ frAlgo.fNFindIterations = frAlgo.GetParameters().GetNcaIterations();
+
+
+ /********************************/ /**
+ * CATrackFinder routine setting
+ ***********************************/
+
+ ReadWindowData();
+
+ for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); ++iS) {
+
+ bool timeUninitialised = 1;
+ fscal lasttime = 0;
+ fscal starttime = 0;
+
+
+ fscal gridMinX = -0.1;
+ fscal gridMaxX = 0.1;
+ fscal gridMinY = -0.1;
+ fscal gridMaxY = 0.1;
+
+ for (ca::HitIndex_t ih = 0; ih < frAlgo.fSliceHitIds[iS].size(); ++ih) {
+ const ca::Hit& h = frAlgo.fInputData.GetHit(frAlgo.fSliceHitIds[iS][ih]);
+
+ if (h.X() < gridMinX) { gridMinX = h.X(); }
+ if (h.X() > gridMaxX) { gridMaxX = h.X(); }
+ if (h.Y() < gridMinY) { gridMinY = h.Y(); }
+ if (h.Y() > gridMaxY) { gridMaxY = h.Y(); }
+
+ const fscal time = h.T();
+ assert(std::isfinite(time));
+ if (timeUninitialised || lasttime < time) { lasttime = time; }
+ if (timeUninitialised || starttime > time) { starttime = time; }
+ timeUninitialised = false;
+ }
+
+ // TODO: changing the grid also changes the result. Investigate why it happens.
+ // TODO: find the optimal grid size
+
+ int nSliceHits = frAlgo.fSliceHitIds[iS].size();
+ fscal sizeY = gridMaxY - gridMinY;
+ fscal sizeX = gridMaxX - gridMinX;
+ int nBins2D = 1 + nSliceHits;
+ fscal yStep = sizeY / sqrt(nBins2D);
+ fscal xStep = sizeX / sqrt(nBins2D);
+
+ fscal scale = frAlgo.GetParameters().GetStation(iS).GetZScal() - frAlgo.GetParameters().GetTargetPositionZ()[0];
+
+ if (yStep > 0.3 * scale) yStep = 0.3 * scale;
+ if (yStep < 0.01 * scale) yStep = 0.01 * scale;
+ if (xStep > 0.3 * scale) xStep = 0.3 * scale;
+ if (xStep < 0.01 * scale) xStep = 0.01 * scale;
+
+ frAlgo.vGrid[iS].BuildBins(gridMinX, gridMaxX, gridMinY, gridMaxY, xStep, yStep);
+ frAlgo.vGrid[iS].StoreHits(frAlgo.fWindowHits, frAlgo.fStationHitsStartIndex[iS], frAlgo.fStationNhits[iS],
+ frAlgo.fvHitKeyFlags);
+ }
+
+
+ // ---- Loop over Track Finder iterations ----------------------------------------------------------------//
+
+
+ assert(frAlgo.fNFindIterations == (int) frAlgo.GetParameters().GetCAIterations().size());
+
+ for (frAlgo.fCurrentIterationIndex = 0;
+ frAlgo.fCurrentIterationIndex < frAlgo.GetParameters().GetCAIterations().size();
+ ++frAlgo.fCurrentIterationIndex) {
+
+ // current iteration of the CA tracker
+ const auto& caIteration = frAlgo.GetParameters().GetCAIterations()[frAlgo.fCurrentIterationIndex];
+ frAlgo.fpCurrentIteration = &caIteration;
+
+ if (frAlgo.fCurrentIterationIndex > 0) {
+ for (int ista = 0; ista < frAlgo.GetParameters().GetNstationsActive(); ++ista) {
+ frAlgo.vGrid[ista].RemoveUsedHits(frAlgo.fWindowHits, frAlgo.fvHitKeyFlags);
+ }
+ }
+
+ frAlgo.fIsWindowHitSuppressed.reset(frAlgo.fWindowHits.size(), 0);
+
+ for (int j = 0; j < frAlgo.GetParameters().GetNstationsActive(); j++) {
+ frAlgo.fTriplets[j].clear();
+ }
+
+ frAlgo.fHitFirstTriplet.reset(frAlgo.fWindowHits.size(), 0);
+ frAlgo.fHitNtriplets.reset(frAlgo.fWindowHits.size(), 0);
+
+ {
+ // #pragma omp task
+ {
+ // --- SET PARAMETERS FOR THE ITERATION ---
+
+ frAlgo.fFirstCAstation = caIteration.GetFirstStationIndex();
+ frAlgo.fTrackChi2Cut = caIteration.GetTrackChi2Cut();
+ frAlgo.fDoubletChi2Cut = caIteration.GetDoubletChi2Cut(); //11.3449 * 2.f / 3.f; // prob = 0.1
+ frAlgo.fTripletChi2Cut = caIteration.GetTripletChi2Cut(); //21.1075; // prob = 0.01%
+ frAlgo.fTripletFinalChi2Cut = caIteration.GetTripletFinalChi2Cut();
+ frAlgo.fPickGather = caIteration.GetPickGather(); //3.0;
+ frAlgo.fTripletLinkChi2 = caIteration.GetTripletLinkChi2(); //5.0;
+ frAlgo.fMaxInvMom = caIteration.GetMaxQp(); //1.0 / 0.5; // max considered q/p
+ frAlgo.fMaxSlopePV = caIteration.GetMaxSlopePV(); //1.1;
+ frAlgo.fMaxSlope = caIteration.GetMaxSlope(); //2.748; // corresponds to 70 grad
+
+ // define the target
+ frAlgo.fTargX = frAlgo.GetParameters().GetTargetPositionX();
+ frAlgo.fTargY = frAlgo.GetParameters().GetTargetPositionY();
+ frAlgo.fTargZ = frAlgo.GetParameters().GetTargetPositionZ();
+
+ fscal SigmaTargetX = caIteration.GetTargetPosSigmaX();
+ fscal SigmaTargetY = caIteration.GetTargetPosSigmaY(); // target constraint [cm]
+
+ // Select magnetic field. For primary tracks - fVtxFieldValue, for secondary tracks - st.fieldSlice
+ if (caIteration.GetPrimaryFlag()) { frAlgo.fTargB = frAlgo.GetParameters().GetVertexFieldValue(); }
+ else {
+ frAlgo.fTargB = frAlgo.GetParameters().GetStation(0).fieldSlice.GetFieldValue(0, 0);
+ } // NOTE: calculates field frAlgo.fTargB in the center of 0th station
+
+ frAlgo.fIsTargetField = (fabs(frAlgo.fTargB.y[0]) > 0.001);
+
+ frAlgo.TargetMeasurement.SetX(frAlgo.fTargX);
+ frAlgo.TargetMeasurement.SetY(frAlgo.fTargY);
+ frAlgo.TargetMeasurement.SetDx2(SigmaTargetX * SigmaTargetX);
+ frAlgo.TargetMeasurement.SetDxy(0);
+ frAlgo.TargetMeasurement.SetDy2(SigmaTargetY * SigmaTargetY);
+ frAlgo.TargetMeasurement.SetNdfX(1);
+ frAlgo.TargetMeasurement.SetNdfY(1);
+
+ /// Set correction in order to take into account overlaping and iff z.
+ /// The reason is that low momentum tracks are too curved and goes not from target direction. That's why sort by hit_y/hit_z is not work idealy
+ /// If sort by y then it is max diff between same station's modules (~0.4cm)
+ frAlgo.fMaxDZ = caIteration.GetMaxDZ(); //0;
+ }
+ }
+
+
+ /// stage for triplets creation
+
+ frAlgo.fMonitor.StartTimer(ETimer::TripletConstruction);
+
+ ca::TripletConstructor constructor1(&frAlgo);
+ ca::TripletConstructor constructor2(&frAlgo);
+ ca::TripletConstructor constructor3(&frAlgo);
+
+ for (int istal = frAlgo.GetParameters().GetNstationsActive() - 2; istal >= frAlgo.fFirstCAstation;
+ istal--) { // start with downstream chambers
+ Tindex nGridEntriesL = frAlgo.vGrid[istal].GetEntries().size();
+
+ for (Tindex iel = 0; iel < nGridEntriesL; ++iel) {
+ ca::HitIndex_t ihitl = frAlgo.vGrid[istal].GetEntries()[iel].GetObjectId();
+ constructor1.CreateTripletsForHit(istal, istal + 1, istal + 2, ihitl);
+
+ if (frAlgo.fpCurrentIteration->GetJumpedFlag()) {
+ constructor2.CreateTripletsForHit(istal, istal + 1, istal + 3, ihitl);
+ constructor3.CreateTripletsForHit(istal, istal + 2, istal + 3, ihitl);
+ }
+
+ auto oldSize = frAlgo.fTriplets[istal].size();
+
+ frAlgo.fTriplets[istal].insert(frAlgo.fTriplets[istal].end(), constructor1.GetTriplets().begin(),
+ constructor1.GetTriplets().end());
+ frAlgo.fTriplets[istal].insert(frAlgo.fTriplets[istal].end(), constructor2.GetTriplets().begin(),
+ constructor2.GetTriplets().end());
+ frAlgo.fTriplets[istal].insert(frAlgo.fTriplets[istal].end(), constructor3.GetTriplets().begin(),
+ constructor3.GetTriplets().end());
+
+ frAlgo.fHitFirstTriplet[ihitl] = frAlgo.PackTripletId(istal, oldSize);
+ frAlgo.fHitNtriplets[ihitl] = frAlgo.fTriplets[istal].size() - oldSize;
+ }
+ } // istal
+
+ frAlgo.fMonitor.StopTimer(ETimer::TripletConstruction);
+
+ // search for neighbouring triplets
+ frAlgo.fMonitor.StartTimer(ETimer::NeighboringTripletSearch);
+ for (int istal = frAlgo.GetParameters().GetNstationsActive() - 2; istal >= frAlgo.fFirstCAstation;
+ istal--) { // start with downstream chambers
+
+ for (unsigned int it = 0; it < frAlgo.fTriplets[istal].size(); ++it) {
+
+ ca::Triplet& tr = frAlgo.fTriplets[istal][it];
+ tr.SetLevel(0);
+ tr.SetFNeighbour(0);
+ tr.SetNNeighbours(0);
+
+ if (istal > (frAlgo.GetParameters().GetNstationsActive() - 4)) break;
+
+ unsigned int nNeighbours = frAlgo.fHitNtriplets[tr.GetMHit()];
+
+ unsigned int neighLocation = frAlgo.fHitFirstTriplet[tr.GetMHit()];
+ unsigned int neighStation = frAlgo.TripletId2Station(neighLocation);
+ unsigned int neighTriplet = frAlgo.TripletId2Triplet(neighLocation);
+
+ if (nNeighbours > 0) { assert((int) neighStation == istal + 1 || (int) neighStation == istal + 2); }
+
+ unsigned char level = 0;
+
+ for (unsigned int iN = 0; iN < nNeighbours; ++iN, ++neighTriplet, ++neighLocation) {
+
+ ca::Triplet& neighbour = frAlgo.fTriplets[neighStation][neighTriplet];
+
+ fscal dchi2 = 0.;
+ if (!checkTripletMatch(tr, neighbour, dchi2)) continue;
+
+ if (tr.GetFNeighbour() == 0) { tr.SetFNeighbour(neighLocation); }
+
+ tr.SetNNeighbours(neighLocation - tr.GetFNeighbour() + 1);
+
+ const unsigned char neighbourLevel = neighbour.GetLevel();
+
+ if (level <= neighbourLevel) { level = neighbourLevel + 1; }
+ }
+ tr.SetLevel(level);
+ } // neighbour search
+
+ frAlgo.fMonitor.IncrementCounter(ECounter::Triplet, frAlgo.fTriplets[istal].size());
+
+ } // istal
+ frAlgo.fMonitor.StopTimer(ETimer::NeighboringTripletSearch);
+
+
+ ///====================================================================
+ ///= =
+ ///= Collect track candidates. CREATE TRACKS =
+ ///= =
+ ///====================================================================
+
+ // min level to start triplet. So min track length = min_level+3.
+ int min_level = std::min(caIteration.GetMinNhits(), caIteration.GetMinNhitsStation0()) - 3;
+
+ if (frAlgo.fpCurrentIteration->GetTrackFromTripletsFlag()) { min_level = 0; }
+
+ ca::Branch curr_tr;
+ ca::Branch new_tr[constants::size::MaxNstations];
+ ca::Branch best_tr;
+
+ int ndf = 1;
+
+
+ // collect consequtive: the longest tracks, shorter, more shorter and so on
+ for (int firstTripletLevel = frAlgo.GetParameters().GetNstationsActive() - 3; firstTripletLevel >= min_level;
+ firstTripletLevel--) {
+ // choose length in triplets number - firstTripletLevel - the maximum possible triplet level among all triplets in the searched candidate
+
+ // how many levels to check
+ int nlevel = (frAlgo.GetParameters().GetNstationsActive() - 2) - firstTripletLevel + 1;
+
+ const unsigned char min_best_l =
+ (firstTripletLevel > min_level) ? firstTripletLevel + 2 : min_level + 3; // loose maximum
+
+
+ frAlgo.fTrackCandidates.clear();
+
+ frAlgo.fStripToTrack.reset(frAlgo.fInputData.GetNhitKeys(), -1);
+
+ //== Loop over triplets with the required level, find and store track candidates
+
+ for (int istaF = frAlgo.fFirstCAstation;
+ istaF <= frAlgo.GetParameters().GetNstationsActive() - 3 - firstTripletLevel; ++istaF) {
+
+ if (--nlevel == 0) break; //TODO: SG: this is not needed
+
+ for (Tindex itrip = 0; itrip < (Tindex) frAlgo.fTriplets[istaF].size(); ++itrip) {
+
+ ca::Triplet& first_trip = (frAlgo.fTriplets[istaF][itrip]);
+ if (frAlgo.fvHitKeyFlags[frAlgo.fWindowHits[first_trip.GetLHit()].FrontKey()]
+ || frAlgo.fvHitKeyFlags[frAlgo.fWindowHits[first_trip.GetLHit()].BackKey()]) {
+ continue;
+ }
+ // ghost suppression !!!
+ //
+
+ if (!frAlgo.fpCurrentIteration->GetTrackFromTripletsFlag()) { // ghost suppression !!!
+ int nHits = 3 + first_trip.GetLevel();
+ if (frAlgo.fWindowHits[first_trip.GetLHit()].Station() == 0) {
+ if (nHits < frAlgo.fpCurrentIteration->GetMinNhitsStation0()) { continue; }
+ }
+ else {
+ if (nHits < frAlgo.fpCurrentIteration->GetMinNhits()) { continue; }
+ }
+ }
+
+ // Collect triplets, which can start a track with length equal to firstTipletLevel + 3. This cut suppresses
+ // ghost tracks, but does not affect the efficiency
+ if (first_trip.GetLevel() < firstTripletLevel) { continue; }
+
+ curr_tr.SetChi2(0.);
+ curr_tr.SetStation(0);
+ curr_tr.ResetHits();
+ curr_tr.AddHit(frAlgo.fWindowHits[first_trip.GetLHit()].Id());
+ curr_tr.SetChi2(first_trip.GetChi2());
+
+ best_tr = curr_tr;
+
+ CAFindTrack(istaF, best_tr, &first_trip, curr_tr, min_best_l,
+ new_tr); /// reqursive func to build a tree of possible track-candidates and choose the best
+
+ if (best_tr.NofHits() < firstTripletLevel + 2) continue; // loose maximum one hit
+
+ if (best_tr.NofHits() < min_level + 3) continue; // should find all hits for min_level
+
+ ndf = best_tr.NofHits() * 2 - 5;
+
+ // TODO: automatize the NDF calculation
+
+ if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode
+ || ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) {
+ ndf = best_tr.NofHits() * 2 - 4;
+ }
+
+ best_tr.SetChi2(best_tr.Chi2() / ndf);
+ if (frAlgo.fGhostSuppression) {
+ if (3 == best_tr.NofHits()) {
+ if (!frAlgo.fpCurrentIteration->GetPrimaryFlag() && (istaF != 0))
+ continue; // too /*short*/ non-MAPS track
+ if (frAlgo.fpCurrentIteration->GetPrimaryFlag() && (best_tr.Chi2() > 5.0)) continue;
+ }
+ }
+ frAlgo.fTrackCandidates.push_back(best_tr);
+ ca::Branch& tr = frAlgo.fTrackCandidates.back();
+ tr.SetStation(istaF);
+ tr.SetId(frAlgo.fTrackCandidates.size() - 1);
+ tr.SetAlive(true);
+ if (0) { // debug
+ cout << "track " << frAlgo.fTrackCandidates.size() - 1 << " found, L = " << best_tr.NofHits()
+ << " chi2= " << best_tr.Chi2() << endl;
+ cout << " hits: ";
+ for (auto hitId : tr.Hits()) {
+ cout << frAlgo.GetMcTrackIdForCaHit(hitId) << " ";
+ }
+ cout << endl;
+ }
+ } // itrip
+ } // istaF
+
+ for (ca::Branch& tr : frAlgo.fTrackCandidates) {
+ tr.SetAlive(false);
+ }
+
+ bool repeatCompetition = true;
+ int iComp = 0;
+ for (iComp = 0; (iComp < 100) && repeatCompetition; ++iComp) {
+
+ // == Loop over track candidates and mark their strips
+
+ repeatCompetition = false;
+
+ for (ca::Branch& tr : frAlgo.fTrackCandidates) {
+ if (tr.IsAlive()) { continue; }
+ for (auto& hitId : tr.Hits()) {
+ const ca::Hit& h = frAlgo.fInputData.GetHit(hitId);
+ bool isAlive = true;
+ { // front strip
+ auto& stripF = (frAlgo.fStripToTrack)[h.FrontKey()];
+ if ((stripF >= 0) && (stripF != tr.Id())) { // strip is used by other candidate
+ const auto& other = frAlgo.fTrackCandidates[stripF];
+ if (!other.IsAlive() && tr.IsBetterThan(other)) { stripF = tr.Id(); }
+ else {
+ isAlive = false;
+ }
+ }
+ else {
+ stripF = tr.Id();
+ }
+ if (!isAlive) { break; }
+ }
+
+ { // back strip
+ auto& stripB = (frAlgo.fStripToTrack)[h.BackKey()];
+ if ((stripB >= 0) && (stripB != tr.Id())) { // strip is used by other candidate
+ const auto& other = frAlgo.fTrackCandidates[stripB];
+ if (!other.IsAlive() && tr.IsBetterThan(other)) { stripB = tr.Id(); }
+ else {
+ isAlive = false;
+ }
+ }
+ else {
+ stripB = tr.Id();
+ }
+ if (!isAlive) { break; }
+ }
+ } // loop over hits
+ } // itrack
+
+ // == Check if some suppressed candidates are still alive
+
+ for (ca::Branch& tr : frAlgo.fTrackCandidates) {
+ if (tr.IsAlive()) { continue; }
+
+ tr.SetAlive(true);
+ for (int iHit = 0; tr.IsAlive() && (iHit < (int) tr.Hits().size()); ++iHit) {
+ const ca::Hit& h = frAlgo.fInputData.GetHit(tr.Hits()[iHit]);
+ tr.SetAlive(tr.IsAlive() && ((frAlgo.fStripToTrack)[h.FrontKey()] == tr.Id())
+ && ((frAlgo.fStripToTrack)[h.BackKey()] == tr.Id()));
+ }
+
+ if (!tr.IsAlive()) { // release strips
+ for (auto hitId : tr.Hits()) {
+ const ca::Hit& h = frAlgo.fInputData.GetHit(hitId);
+ if (frAlgo.fStripToTrack[h.FrontKey()] == tr.Id()) { frAlgo.fStripToTrack[h.FrontKey()] = -1; }
+ if (frAlgo.fStripToTrack[h.BackKey()] == tr.Id()) { frAlgo.fStripToTrack[h.BackKey()] = -1; }
+ }
+ }
+ else {
+ repeatCompetition = true;
+ }
+ } // itrack
+
+ } // competitions
+
+ // cout << " N Competitions: " << iComp << endl;
+
+ // ==
+
+ for (Tindex iCandidate = 0; iCandidate < (Tindex) frAlgo.fTrackCandidates.size(); ++iCandidate) {
+ ca::Branch& tr = frAlgo.fTrackCandidates[iCandidate];
+ if (!tr.IsAlive()) continue;
+ if (ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) {
+ if (tr.NofHits() <= 3) { continue; }
+ }
+ else if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode) {
+ if (tr.NofHits() < 3) { continue; }
+ }
+ else {
+ if (tr.NofHits() < 3) { continue; }
+ }
+ if (frAlgo.fpCurrentIteration->GetExtendTracksFlag()) {
+ if (tr.NofHits() < frAlgo.GetParameters().GetNstationsActive()) { frAlgo.fTrackExtender.ExtendBranch(tr); }
+ }
+
+ for (auto iHit : tr.Hits()) {
+ const ca::Hit& hit = frAlgo.fInputData.GetHit(iHit);
+
+ /// used strips are marked
+
+ frAlgo.fvHitKeyFlags[hit.FrontKey()] = 1;
+ frAlgo.fvHitKeyFlags[hit.BackKey()] = 1;
+
+ frAlgo.fSliceRecoHits.push_back(iHit);
+ }
+ Track t;
+ t.fNofHits = tr.NofHits();
+ frAlgo.fSliceRecoTracks.push_back(t);
+ if (0) { // SG debug
+ cout << "store track " << iCandidate << " chi2= " << tr.Chi2() << endl;
+ cout << " hits: ";
+ for (unsigned int ih = 0; ih < tr.Hits().size(); ih++) {
+ cout << frAlgo.GetMcTrackIdForCaHit(tr.Hits()[ih]) << " ";
+ }
+ cout << '\n';
+ }
+
+ } // tracks
+
+ } // firstTripletLevel
+
+
+ // suppress strips of suppressed hits
+ for (unsigned int ih = 0; ih < frAlgo.fWindowHits.size(); ih++) {
+ if (frAlgo.fIsWindowHitSuppressed[ih]) {
+ const ca::Hit& hit = frAlgo.fWindowHits[ih];
+ frAlgo.fvHitKeyFlags[hit.FrontKey()] = 1;
+ frAlgo.fvHitKeyFlags[hit.BackKey()] = 1;
+ }
+ }
+ } // ---- Loop over Track Finder iterations: END -----------------------------------------------------------//
+
+ // Fit tracks
+ frAlgo.fTrackFitter.FitCaTracks();
+
+ // Merge clones
+ frAlgo.fCloneMerger.Exec(frAlgo.fSliceRecoTracks, frAlgo.fSliceRecoHits);
+}
+
+
+/** *************************************************************
+ * *
+ * The routine performs recursive search for tracks *
+ * *
+ * I. Kisel 06.03.05 *
+ * I.Kulakov 2012 *
+ * *
+ ****************************************************************/
+
+void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Triplet* curr_trip, ca::Branch& curr_tr,
+ unsigned char min_best_l, ca::Branch* new_tr)
+/// recursive search for tracks
+/// input: @ista - station index, @&best_tr - best track for the privious call
+/// output: @&NCalls - number of function calls
+{
+
+ if (curr_trip->GetLevel() == 0) // the end of the track -> check and store
+ {
+ // -- finish with current track
+ // add rest of hits
+ const ca::HitIndex_t& ihitm = curr_trip->GetMHit();
+ const ca::HitIndex_t& ihitr = curr_trip->GetRHit();
+
+ if (!(frAlgo.fvHitKeyFlags[frAlgo.fWindowHits[ihitm].FrontKey()]
+ || frAlgo.fvHitKeyFlags[frAlgo.fWindowHits[ihitm].BackKey()])) {
+
+ // curr_tr.Hits.push_back(fGridHitIds[ihitm]);
+
+ curr_tr.AddHit(frAlgo.fWindowHits[ihitm].Id());
+ }
+
+ if (!(frAlgo.fvHitKeyFlags[frAlgo.fWindowHits[ihitr].FrontKey()]
+ || frAlgo.fvHitKeyFlags[frAlgo.fWindowHits[ihitr].BackKey()])) {
+
+ //curr_tr.Hits.push_back(fGridHitIds[ihitr]);
+ curr_tr.AddHit(frAlgo.fWindowHits[ihitr].Id());
+ }
+
+ //if( curr_tr.NofHits() < min_best_l - 1 ) return; // suppouse that only one hit can be added by extender
+ // TODO: SZh 21.08.2023: Replace hard-coded value with a parameter
+ int ndf = curr_tr.NofHits() * 2 - 5;
+
+ if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode
+ || ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) {
+ ndf = curr_tr.NofHits() * 2 - 4;
+ }
+
+ if (curr_tr.Chi2() > frAlgo.fTrackChi2Cut * ndf) { return; }
+
+
+ // -- select the best
+ if ((curr_tr.NofHits() > best_tr.NofHits())
+ || ((curr_tr.NofHits() == best_tr.NofHits()) && (curr_tr.Chi2() < best_tr.Chi2()))) {
+ best_tr = curr_tr;
+ }
+
+ return;
+ }
+ else //MEANS level ! = 0
+ {
+ int N_neighbour = (curr_trip->GetNNeighbours());
+
+ for (Tindex in = 0; in < N_neighbour; in++) {
+
+ unsigned int ID = curr_trip->GetFNeighbour() + in;
+
+ unsigned int Station = frAlgo.TripletId2Station(ID);
+ unsigned int Triplet = frAlgo.TripletId2Triplet(ID);
+
+ const ca::Triplet& new_trip = frAlgo.fTriplets[Station][Triplet];
+
+ fscal dchi2 = 0.;
+ if (!checkTripletMatch(*curr_trip, new_trip, dchi2)) continue;
+
+ if (frAlgo.fvHitKeyFlags[frAlgo.fWindowHits[new_trip.GetLHit()].FrontKey()]
+ || frAlgo.fvHitKeyFlags[frAlgo.fWindowHits[new_trip.GetLHit()].BackKey()]) { //hits are used
+ // no used hits allowed -> compare and store track
+ if ((curr_tr.NofHits() > best_tr.NofHits())
+ || ((curr_tr.NofHits() == best_tr.NofHits()) && (curr_tr.Chi2() < best_tr.Chi2()))) {
+ best_tr = curr_tr;
+ }
+ }
+ else { // hit is not used: add the left hit from the new triplet to the current track
+
+ unsigned char new_L = curr_tr.NofHits() + 1;
+ fscal new_chi2 = curr_tr.Chi2() + dchi2;
+
+
+ if (0) { //SGtrd2d debug!!
+ int mc01 = frAlgo.GetMcTrackIdForWindowHit(curr_trip->GetLHit());
+ int mc02 = frAlgo.GetMcTrackIdForWindowHit(curr_trip->GetMHit());
+ int mc03 = frAlgo.GetMcTrackIdForWindowHit(curr_trip->GetRHit());
+ int mc11 = frAlgo.GetMcTrackIdForWindowHit(new_trip.GetLHit());
+ int mc12 = frAlgo.GetMcTrackIdForWindowHit(new_trip.GetMHit());
+ int mc13 = frAlgo.GetMcTrackIdForWindowHit(new_trip.GetRHit());
+ if ((mc01 == mc02) && (mc02 == mc03)) {
+ cout << " sta " << ista << " mc0 " << mc01 << " " << mc02 << " " << mc03 << " mc1 " << mc11 << " " << mc12
+ << " " << mc13 << " chi2 " << curr_tr.Chi2() / (2 * (curr_tr.NofHits() + 2) - 4) << " new "
+ << new_chi2 / (2 * (new_L + 2) - 4) << endl;
+ cout << " hits " << curr_trip->GetLHit() << " " << curr_trip->GetMHit() << " " << curr_trip->GetRHit()
+ << " " << new_trip.GetLHit() << endl;
+ }
+ }
+
+ int ndf = 2 * (new_L + 2) - 5;
+
+ if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode) { //SGtrd2d!!!
+ ndf = 2 * (new_L + 2) - 4;
+ }
+ else if (ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) {
+ ndf = 2 * (new_L + 2) - 4;
+ }
+ else {
+ ndf = new_L; // TODO: SG: 2 * (new_L + 2) - 5
+ }
+
+ if (new_chi2 > frAlgo.fTrackChi2Cut * ndf) continue;
+
+ // add new hit
+ new_tr[ista] = curr_tr;
+ new_tr[ista].AddHit(frAlgo.fWindowHits[new_trip.GetLHit()].Id());
+
+ const int new_ista = ista + new_trip.GetMSta() - new_trip.GetLSta();
+ new_tr[ista].SetChi2(new_chi2);
+
+ CAFindTrack(new_ista, best_tr, &new_trip, new_tr[ista], min_best_l, new_tr);
+ } // add triplet to track
+ } // for neighbours
+ } // level = 0
+}
diff --git a/algo/ca/core/tracking/CaTrackFinderWindow.h b/algo/ca/core/tracking/CaTrackFinderWindow.h
new file mode 100644
index 0000000000..edf3846447
--- /dev/null
+++ b/algo/ca/core/tracking/CaTrackFinderWindow.h
@@ -0,0 +1,73 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer], Maksym Zyzak */
+
+/// \file CaTrackFinderWindow.h
+/// \authors S.Zharko <s.zharko@gsi.de> (interface), M.Zyzak (original algorithm)
+/// \brief A class wrapper over clones merger algorithm for the CA track finder (declaration)
+/// \since 22.07.2022
+
+#pragma once // include this header only once per compilation unit
+
+#include "CaBranch.h"
+#include "CaHit.h"
+#include "CaSimd.h"
+#include "CaTrackParam.h"
+#include "CaVector.h"
+
+
+namespace cbm::algo::ca
+{
+
+ class Track;
+ class Triplet;
+ class Framework;
+
+ namespace
+ {
+ using namespace cbm::algo; // to keep ca:: prefices in the code
+ }
+
+ /// Class implements a clones merger algorithm for the CA track finder
+ ///
+ class TrackFinderWindow {
+ public:
+ /// Default constructor
+ TrackFinderWindow(ca::Framework& algo);
+
+ /// Destructor
+ ~TrackFinderWindow();
+
+ /// Copy constructor
+ TrackFinderWindow(const TrackFinderWindow&) = delete;
+
+ /// Move constructor
+ TrackFinderWindow(TrackFinderWindow&&) = delete;
+
+ /// Copy assignment operator
+ TrackFinderWindow& operator=(const TrackFinderWindow&) = delete;
+
+ /// Move assignment operator
+ TrackFinderWindow& operator=(TrackFinderWindow&&) = delete;
+
+ void ReadWindowData();
+ void CaTrackFinderSlice();
+
+ void CAFindTrack(int ista, ca::Branch& best_tr, const ca::Triplet* curr_trip, ca::Branch& curr_tr,
+ unsigned char min_best_l, ca::Branch* new_tr);
+
+ private:
+ ///-------------------------------
+ /// Private methods
+
+ bool checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2) const;
+
+
+ private:
+ ///-------------------------------
+ /// Data members
+
+ ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class
+ };
+
+} // namespace cbm::algo::ca
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/algo/ca/core/tracking/CaTrackFitter.cxx
similarity index 81%
rename from reco/L1/L1Algo/L1TrackFitter.cxx
rename to algo/ca/core/tracking/CaTrackFitter.cxx
index 4d4e762724..b83eba8191 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/algo/ca/core/tracking/CaTrackFitter.cxx
@@ -2,29 +2,38 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Ivan Kisel, Sergey Gorbunov, Igor Kulakov [committer], Valentina Akishina, Maksym Zyzak */
+#include "CaTrackFitter.h"
+
#include <iostream>
#include <vector>
+#include "CaFramework.h"
#include "CaTrackFit.h"
#include "CaTrackParam.h"
-#include "L1Algo.h"
using std::cout;
using std::endl;
using std::vector;
using Track = cbm::algo::ca::Track;
-namespace NS_L1TrackFitter
-{
- const fvec c_light(0.000299792458), c_light_i(fvec(1.) / c_light);
-} // namespace NS_L1TrackFitter
-using namespace NS_L1TrackFitter;
+using namespace cbm::algo::ca;
+using namespace cbm::algo;
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+TrackFitter::TrackFitter(ca::Framework& algo) : frAlgo(algo) {}
-void L1Algo::L1KFTrackFitter()
+// ---------------------------------------------------------------------------------------------------------------------
+//
+TrackFitter::~TrackFitter() {}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void TrackFitter::FitCaTracks()
{
- // cout << " Start L1 Track Fitter " << endl;
- int start_hit = 0; // for interation in fSliceRecoHits[]
+ // cout << " Start CA Track Fitter " << endl;
+ int start_hit = 0; // for interation in frAlgo.fSliceRecoHits[]
// static ca::FieldValue fldB0, fldB1, fldB2 _fvecalignment;
// static ca::FieldRegion fld _fvecalignment;
@@ -35,17 +44,17 @@ void L1Algo::L1KFTrackFitter()
ca::FieldValue fldB01, fldB11, fldB21 _fvecalignment;
ca::FieldRegion fld1 _fvecalignment;
- const int nStations = fParameters.GetNstationsActive();
+ const int nStations = frAlgo.GetParameters().GetNstationsActive();
int nTracks_SIMD = fvec::size();
ca::TrackFit fit; // fit parameters coresponding to the current track
TrackParamV& tr = fit.Tr();
- fit.SetParticleMass(GetDefaultParticleMass());
+ fit.SetParticleMass(frAlgo.GetDefaultParticleMass());
fit.SetDoFitVelocity(true);
Track* t[fvec::size()] {nullptr};
- const ca::Station* sta = fParameters.GetStations().begin();
+ const ca::Station* sta = frAlgo.GetParameters().GetStations().begin();
// Spatial-time position of a hit vs. station and track in the portion
@@ -95,19 +104,19 @@ void L1Algo::L1KFTrackFitter()
mxy[ista].SetCov(1., 0., 1.);
}
- unsigned short N_vTracks = fSliceRecoTracks.size(); // number of tracks processed per one SSE register
+ unsigned short N_vTracks = frAlgo.fSliceRecoTracks.size(); // number of tracks processed per one SSE register
for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvec::size()) {
if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) nTracks_SIMD = N_vTracks - itrack;
for (int i = 0; i < nTracks_SIMD; i++) {
- t[i] = &fSliceRecoTracks[itrack + i];
+ t[i] = &frAlgo.fSliceRecoTracks[itrack + i];
}
// fill the rest of the SIMD vectors with something reasonable
for (uint i = nTracks_SIMD; i < fvec::size(); i++) {
- t[i] = &fSliceRecoTracks[itrack];
+ t[i] = &frAlgo.fSliceRecoTracks[itrack];
}
// get hits of current track
@@ -126,7 +135,7 @@ void L1Algo::L1KFTrackFitter()
for (int ih = 0; ih < nHitsTrack; ih++) {
- const ca::Hit& hit = fInputData.GetHit(fSliceRecoHits[start_hit++]);
+ const ca::Hit& hit = frAlgo.GetInputData().GetHit(frAlgo.fSliceRecoHits[start_hit++]);
const int ista = hit.Station();
//if (sta[ista].fieldStatus) { isFieldPresent[iVec] = true; }
@@ -187,14 +196,17 @@ void L1Algo::L1KFTrackFitter()
for (int ih = nHitsTrack - 1; ih >= 0; ih--) {
const int ista = iSta[ih];
- const ca::Station& st = fParameters.GetStation(ista);
+ const ca::Station& st = frAlgo.GetParameters().GetStation(ista);
By[ista][iVec] = st.fieldSlice.cy[0][0];
}
}
fit.GuessTrack(z_end, x, y, z, time, By, w, w_time, nStations);
- if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) { tr.Qp() = fvec(1. / 1.1); }
+ if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode
+ || ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) {
+ tr.Qp() = fvec(1. / 1.1);
+ }
// tr.Qp() = iif(isFieldPresent, tr.Qp(), fvec(1. / 0.25));
@@ -245,8 +257,8 @@ void L1Algo::L1KFTrackFitter()
fit.SetMask(initialised);
fit.Extrapolate(z[ista], fld1);
- fit.MultipleScattering(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()));
- fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()), fvec(1.f));
+ fit.MultipleScattering(frAlgo.GetParameters().GetMaterialThickness(ista, tr.X(), tr.Y()));
+ fit.EnergyLossCorrection(frAlgo.GetParameters().GetMaterialThickness(ista, tr.X(), tr.Y()), fvec(1.f));
fit.SetMask(initialised && w[ista]);
fit.FilterXY(mxy[ista]);
@@ -265,23 +277,23 @@ void L1Algo::L1KFTrackFitter()
{
fitpv.SetMask(fmask::One());
- ca::MeasurementXy<fvec> vtxInfo = TargetMeasurement;
+ ca::MeasurementXy<fvec> vtxInfo = frAlgo.TargetMeasurement;
vtxInfo.SetDx2(1.e-8);
vtxInfo.SetDxy(0.);
vtxInfo.SetDy2(1.e-8);
- if (kGlobal == fTrackingMode) {
+ if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode) {
for (int vtxIter = 0; vtxIter < 2; vtxIter++) {
fitpv.SetQp0(fitpv.Tr().Qp());
fitpv.Tr() = fit.Tr();
fitpv.Tr().Qp() = fitpv.Qp0();
- fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fld);
+ fitpv.Extrapolate(frAlgo.GetParameters().GetTargetPositionZ(), fld);
fitpv.FilterXY(vtxInfo);
}
}
else {
fitpv.SetQp0(fitpv.Tr().Qp());
- fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fld);
+ fitpv.Extrapolate(frAlgo.GetParameters().GetTargetPositionZ(), fld);
}
}
@@ -291,7 +303,7 @@ void L1Algo::L1KFTrackFitter()
//fit.FilterVi(TrackParamV::kClightNsInv);
TrackParamV Tf = fit.Tr();
- if (kGlobal == fTrackingMode) { Tf.Qp() = fitpv.Tr().Qp(); }
+ if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode) { Tf.Qp() = fitpv.Tr().Qp(); }
for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
t[iVec]->fParFirst.Set(Tf, iVec);
@@ -342,8 +354,8 @@ void L1Algo::L1KFTrackFitter()
fit.SetMask(initialised);
fit.Extrapolate(z[ista], fld);
- fit.MultipleScattering(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()));
- fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()), fvec(-1.f));
+ fit.MultipleScattering(frAlgo.GetParameters().GetMaterialThickness(ista, tr.X(), tr.Y()));
+ fit.EnergyLossCorrection(frAlgo.GetParameters().GetMaterialThickness(ista, tr.X(), tr.Y()), fvec(-1.f));
fit.SetMask(initialised && w[ista]);
fit.FilterXY(mxy[ista]);
fit.FilterTime(time[ista], dt2[ista], sta[ista].timeInfo);
@@ -360,7 +372,7 @@ void L1Algo::L1KFTrackFitter()
//fit.FilterVi(TrackParamV::kClightNsInv);
TrackParamV Tl = fit.Tr();
- if (kGlobal == fTrackingMode) { Tl.Qp() = fitpv.Tr().Qp(); }
+ if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode) { Tl.Qp() = fitpv.Tr().Qp(); }
for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
t[iVec]->fParLast.Set(Tl, iVec);
diff --git a/algo/ca/core/tracking/CaTrackFitter.h b/algo/ca/core/tracking/CaTrackFitter.h
new file mode 100644
index 0000000000..130a42903f
--- /dev/null
+++ b/algo/ca/core/tracking/CaTrackFitter.h
@@ -0,0 +1,66 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer], Maksym Zyzak */
+
+/// \file TrackFitter.h
+/// \brief A class wrapper over clones merger algorithm for the Ca track finder (declaration)
+/// \since 19.10.2023
+
+#pragma once // include this header only once per compilation unit
+
+#include "CaBranch.h"
+#include "CaHit.h"
+#include "CaSimd.h"
+#include "CaTrackParam.h"
+#include "CaVector.h"
+
+
+namespace cbm::algo::ca
+{
+ class Track;
+ class Framework;
+
+ namespace
+ {
+ using namespace cbm::algo; // to keep ca:: prefices in the code
+ }
+
+
+ /// Class implements a track fit the CA track finder
+ ///
+ class TrackFitter {
+ public:
+ /// Default constructor
+ TrackFitter(ca::Framework& algo);
+
+ /// Destructor
+ ~TrackFitter();
+
+ /// Copy constructor
+ TrackFitter(const TrackFitter&) = delete;
+
+ /// Move constructor
+ TrackFitter(TrackFitter&&) = delete;
+
+ /// Copy assignment operator
+ TrackFitter& operator=(const TrackFitter&) = delete;
+
+ /// Move assignment operator
+ TrackFitter& operator=(TrackFitter&&) = delete;
+
+
+ /// Fit tracks, found by the CA tracker
+ void FitCaTracks();
+
+ private:
+ ///-------------------------------
+ /// Private methods
+
+ private:
+ ///-------------------------------
+ /// Data members
+
+ ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class
+ };
+
+} // namespace cbm::algo::ca
diff --git a/reco/L1/L1Algo/L1TripletConstructor.cxx b/algo/ca/core/tracking/CaTripletConstructor.cxx
similarity index 93%
rename from reco/L1/L1Algo/L1TripletConstructor.cxx
rename to algo/ca/core/tracking/CaTripletConstructor.cxx
index 4ae4c1debb..2b0e2433f7 100644
--- a/reco/L1/L1Algo/L1TripletConstructor.cxx
+++ b/algo/ca/core/tracking/CaTripletConstructor.cxx
@@ -2,26 +2,25 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Igor Kulakov [committer], Maksym Zyzak, Valentina Akishina */
-#include "L1TripletConstructor.h"
-
-#include "CbmL1.h"
-#include "CbmL1MCTrack.h"
+#include "CaTripletConstructor.h"
#include <algorithm>
#include <iostream>
#include "CaDefines.h"
+#include "CaFramework.h"
#include "CaGridArea.h"
-#include "CaToolsDebugger.h"
+// #include "CaToolsDebugger.h"
#include "CaTrackFit.h"
#include "CaTrackParam.h"
-#include "L1Algo.h"
-using cbm::ca::tools::Debugger;
+// using cbm::ca::tools::Debugger;
+
+using namespace cbm::algo::ca;
-L1TripletConstructor::L1TripletConstructor(L1Algo* algo) { fAlgo = algo; }
+TripletConstructor::TripletConstructor(ca::Framework* algo) { fAlgo = algo; }
-void L1TripletConstructor::InitStations(int istal, int istam, int istar)
+void TripletConstructor::InitStations(int istal, int istam, int istar)
{
fIsInitialized = false;
@@ -74,7 +73,7 @@ void L1TripletConstructor::InitStations(int istal, int istam, int istar)
}
-void L1TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, ca::HitIndex_t iHitL)
+void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, ca::HitIndex_t iHitL)
{
InitStations(istal, istam, istar);
@@ -169,7 +168,7 @@ void L1TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar,
}
-void L1TripletConstructor::FindDoublets()
+void TripletConstructor::FindDoublets()
{
/// Find the doublets. Reformat data into portions of doublets.
@@ -183,12 +182,12 @@ void L1TripletConstructor::FindDoublets()
}
-void L1TripletConstructor::FitDoublets()
+void TripletConstructor::FitDoublets()
{
// ---- Add the middle hits to parameters estimation ----
- Vector<ca::HitIndex_t> hitsMtmp("L1TripletConstructor::hitsMtmp", fHitsM_2);
+ Vector<ca::HitIndex_t> hitsMtmp("TripletConstructor::hitsMtmp", fHitsM_2);
fHitsM_2.clear();
fTracks_2.clear();
@@ -227,7 +226,7 @@ void L1TripletConstructor::FitDoublets()
fFit.SetQp0(fFit.Tr().Qp());
// check if there are other hits close to the doublet on the same station
- if (L1Algo::kMcbm != fAlgo->fTrackingMode) {
+ if (ca::Framework::kMcbm != fAlgo->fTrackingMode) {
// TODO: SG: adjust cuts, put them to parameters
fscal tx = T2.Tx()[0];
@@ -268,7 +267,7 @@ void L1TripletConstructor::FitDoublets()
}
-void L1TripletConstructor::FindTriplets()
+void TripletConstructor::FindTriplets()
{
if (fIstaR >= fAlgo->fParameters.GetNstationsActive()) { return; }
FindRightHit();
@@ -277,7 +276,7 @@ void L1TripletConstructor::FindTriplets()
}
-void L1TripletConstructor::FindRightHit()
+void TripletConstructor::FindRightHit()
{
/// Add the middle hits to parameters estimation. Propagate to right station.
@@ -324,7 +323,7 @@ void L1TripletConstructor::FindRightHit()
CollectHits(T2, fIstaR, fAlgo->fTripletChi2Cut, iMC, collectedHits, fAlgo->fParameters.GetMaxTripletPerDoublets());
if (collectedHits.size() >= fAlgo->fParameters.GetMaxTripletPerDoublets()) {
- LOG(debug) << "L1: GetMaxTripletPerDoublets==" << fAlgo->fParameters.GetMaxTripletPerDoublets()
+ LOG(debug) << "Ca: GetMaxTripletPerDoublets==" << fAlgo->fParameters.GetMaxTripletPerDoublets()
<< " reached in findTripletsStep0()";
// reject already created triplets for this doublet
collectedHits.clear();
@@ -342,7 +341,7 @@ void L1TripletConstructor::FindRightHit()
}
-void L1TripletConstructor::FitTriplets()
+void TripletConstructor::FitTriplets()
{
constexpr bool dumpTriplets = 0;
@@ -356,8 +355,8 @@ void L1TripletConstructor::FitTriplets()
/// Refit Triplets
if (dumpTriplets) {
- cbm::ca::tools::Debugger::Instance().AddNtuple(
- "triplets", "ev:iter:i0:x0:y0:z0:i1:x1:y1:z1:i2:x2:y2:z2:mc:sta:p:vx:vy:vz:chi2:ndf:chi2time:ndfTime");
+ //cbm::ca::tools::Debugger::Instance().AddNtuple(
+ // "triplets", "ev:iter:i0:x0:y0:z0:i1:x1:y1:z1:i2:x2:y2:z2:mc:sta:p:vx:vy:vz:chi2:ndf:chi2time:ndfTime");
}
ca::TrackFit fit;
@@ -520,17 +519,17 @@ void L1TripletConstructor::FitTriplets()
int mc1 = fAlgo->GetMcTrackIdForCaHit(ih0);
int mc2 = fAlgo->GetMcTrackIdForCaHit(ih1);
int mc3 = fAlgo->GetMcTrackIdForCaHit(ih2);
-
- const ca::Hit& h0 = fAlgo->fInputData.GetHit(ih0);
- const ca::Hit& h1 = fAlgo->fInputData.GetHit(ih1);
- const ca::Hit& h2 = fAlgo->fInputData.GetHit(ih2);
-
if ((mc1 >= 0) && (mc1 == mc2) && (mc1 == mc3)) {
+ /*
+ const ca::Hit& h0 = fAlgo->fInputData.GetHit(ih0);
+ const ca::Hit& h1 = fAlgo->fInputData.GetHit(ih1);
+ const ca::Hit& h2 = fAlgo->fInputData.GetHit(ih2);
const CbmL1MCTrack& mctr = CbmL1::Instance()->GetMcTracks()[mc1];
cbm::ca::tools::Debugger::Instance().FillNtuple(
"triplets", mctr.iEvent, fAlgo->fCurrentIterationIndex, ih0, h0.X(), h0.Y(), h0.Z(), ih1, h1.X(), h1.Y(),
h1.Z(), ih2, h2.X(), h2.Y(), h2.Z(), mc1, fIstaL, mctr.p, mctr.x, mctr.y, mctr.z, (float) T.GetChiSq()[0],
(float) T.Ndf()[0], (float) T.ChiSqTime()[0], (float) T.NdfTime()[0]);
+ */
}
}
} //i3
@@ -538,7 +537,7 @@ void L1TripletConstructor::FitTriplets()
} // FitTriplets
-void L1TripletConstructor::StoreTriplets()
+void TripletConstructor::StoreTriplets()
{
/// Selects good triplets and saves them into fTriplets.
/// Finds neighbouring triplets at the next station.
@@ -587,8 +586,8 @@ void L1TripletConstructor::StoreTriplets()
}
-void L1TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, const double chi2Cut, const int iMC,
- Vector<ca::HitIndex_t>& collectedHits, int maxNhits)
+void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, const double chi2Cut, const int iMC,
+ Vector<ca::HitIndex_t>& collectedHits, int maxNhits)
{
/// Collect hits on a station
diff --git a/algo/ca/core/tracking/CaTripletConstructor.h b/algo/ca/core/tracking/CaTripletConstructor.h
new file mode 100644
index 0000000000..477bee185a
--- /dev/null
+++ b/algo/ca/core/tracking/CaTripletConstructor.h
@@ -0,0 +1,122 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov [committer] */
+
+/// \file CaTripletConstructor.h
+/// \brief Triplet constructor for the CA tracker
+/// \author Sergey Gorbunov
+
+#pragma once // include this header only once per compilation unit
+
+#include "CaField.h"
+#include "CaFramework.h"
+#include "CaGridEntry.h"
+#include "CaStation.h"
+#include "CaTrackFit.h"
+#include "CaTrackParam.h"
+#include "CaTriplet.h"
+#include "CaVector.h"
+
+namespace cbm::algo::ca
+{
+
+ namespace
+ {
+ using namespace cbm::algo; // to keep ca:: prefices in the code
+ }
+
+ /// Construction of triplets for the CA tracker
+ ///
+ class TripletConstructor {
+ public:
+ /// ------ Constructors and destructor ------
+
+ /// Constructor
+ /// \param nThreads Number of threads for multi-threaded mode
+ TripletConstructor(ca::Framework* algo);
+
+ /// Copy constructor
+ TripletConstructor(const TripletConstructor&) = delete;
+
+ /// Move constructor
+ TripletConstructor(TripletConstructor&&) = delete;
+
+ /// Copy assignment operator
+ TripletConstructor& operator=(const TripletConstructor&) = delete;
+
+ /// Move assignment operator
+ TripletConstructor& operator=(TripletConstructor&&) = delete;
+
+ /// Destructor
+ ~TripletConstructor() = default;
+
+
+ /// ------ FUNCTIONAL PART ------
+
+ void InitStations(int istal, int istam, int istar);
+
+ void CreateTripletsForHit(int istal, int istam, int istar, ca::HitIndex_t ihl);
+
+ const Vector<ca::Triplet>& GetTriplets() const { return fTriplets; }
+
+ /// Find the doublets. Reformat data in the portion of doublets.
+ void FindDoublets();
+ void FitDoublets();
+
+ /// Find triplets on station
+ void FindTriplets();
+
+ /// Add the middle hits to parameters estimation. Propagate to right station.
+ /// Find the triplets (right hit). Reformat data in the portion of triplets.
+ void FindRightHit();
+
+ /// Fit Triplets
+ void FitTriplets();
+
+ /// Select triplets. Save them into vTriplets.
+ void StoreTriplets();
+
+ void CollectHits(const TrackParamV& Tr, const int iSta, const double chi2Cut, const int iMC,
+ Vector<ca::HitIndex_t>& collectedHits, int maxNhits);
+
+ private:
+ /// left station
+ const ca::Station& staL() { return *fStaL; }
+ /// middle station
+ const ca::Station& staM() { return *fStaM; }
+ /// right station
+ const ca::Station& staR() { return *fStaR; }
+
+ private:
+ bool fIsInitialized {false}; ///< is initialized;
+ ca::Framework* fAlgo {nullptr}; ///< pointer to ca::Framework object
+
+ int fIstaL {-1}; ///< left station index
+ int fIstaM {-1}; ///< middle station index
+ int fIstaR {-1}; ///< right station index
+
+ const ca::Station* fStaL {nullptr}; ///< left station
+ const ca::Station* fStaM {nullptr}; ///< mid station
+ const ca::Station* fStaR {nullptr}; ///< right station
+
+ const ca::Station* fFld0Sta[2]; // two stations for approximating the field between the target and the left hit
+ const ca::Station* fFld1Sta[3]; // three stations for approximating the field between the left and the right hit
+
+ ca::HitIndex_t fIhitL; ///< index of the left hit in fAlgo->fWindowHits
+ TrackParamV fTrL;
+ ca::FieldRegion fFldL;
+
+ Vector<ca::HitIndex_t> fHitsM_2 {"TripletConstructor::fHitsM_2"};
+ Vector<TrackParamV> fTracks_2 {"TripletConstructor::fTracks_2"};
+
+ Vector<TrackParamV> fTracks_3 {"TripletConstructor::fTracks_3"};
+ Vector<ca::HitIndex_t> fHitsM_3 {"TripletConstructor::fHitsM_3"};
+ Vector<ca::HitIndex_t> fHitsR_3 {"TripletConstructor::fHitsR_3"};
+
+ Vector<ca::Triplet> fTriplets {"TripletConstructor::fTriplets"};
+
+ ca::TrackFit fFit;
+
+ } _fvecalignment;
+
+} // namespace cbm::algo::ca
diff --git a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
index b4a6162c7a..1b8b6c60b9 100644
--- a/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
+++ b/reco/KF/ParticleFitter/CbmL1PFFitter.cxx
@@ -27,7 +27,7 @@
#include "TClonesArray.h"
-//L1Algo tools
+//ca::Framework tools
#include "CbmKFVertex.h"
#include "FairRootManager.h"
@@ -35,12 +35,15 @@
#include "TDatabasePDG.h"
#include "CaField.h"
+#include "CaFramework.h"
+#include "CaSimd.h"
#include "CaStation.h"
#include "CaToolsField.h"
#include "CaTrackFit.h"
#include "CaTrackParam.h"
#include "KFParticleDatabase.h"
-#include "L1Algo.h" // Also defines L1Parameters
+
+using namespace cbm::algo::ca;
//typedef ca::TrackFit1 ca::TrackFit;
@@ -111,13 +114,11 @@ inline void CbmL1PFFitter::Initialize()
if (!manager) { LOG(fatal) << "CbmL1PFFitter: no FairRootManager"; }
- if (!CbmL1::Instance() || !CbmL1::Instance()->fpAlgo || !CbmL1::Instance()->fpAlgo->GetParameters()) {
- LOG(fatal) << "CbmL1PFFitter: no CbmL1 task initialised ";
- }
+ if (!CbmL1::Instance() || !CbmL1::Instance()->fpAlgo) { LOG(fatal) << "CbmL1PFFitter: no CbmL1 task initialised "; }
using cbm::algo::ca::EDetectorID;
- fNmvdStationsActive = CbmL1::Instance()->fpAlgo->GetParameters()->GetNstationsActive(EDetectorID::kMvd);
- fNstsStationsActive = CbmL1::Instance()->fpAlgo->GetParameters()->GetNstationsActive(EDetectorID::kSts);
+ fNmvdStationsActive = CbmL1::Instance()->fpAlgo->GetParameters().GetNstationsActive(EDetectorID::kMvd);
+ fNstsStationsActive = CbmL1::Instance()->fpAlgo->GetParameters().GetNstationsActive(EDetectorID::kSts);
if (fNmvdStationsActive > 0) { fMvdHitArray = static_cast<TClonesArray*>(manager->GetObject("MvdHit")); }
if (fNstsStationsActive > 0) { fStsHitArray = static_cast<TClonesArray*>(manager->GetObject("StsHit")); }
@@ -130,13 +131,13 @@ inline void CbmL1PFFitter::Initialize()
inline int CbmL1PFFitter::GetMvdStationIndex(const CbmMvdHit* hit)
{
using cbm::algo::ca::EDetectorID;
- return CbmL1::Instance()->fpAlgo->GetParameters()->GetStationIndexActive(hit->GetStationNr(), EDetectorID::kMvd);
+ return CbmL1::Instance()->fpAlgo->GetParameters().GetStationIndexActive(hit->GetStationNr(), EDetectorID::kMvd);
}
inline int CbmL1PFFitter::GetStsStationIndex(const CbmStsHit* hit)
{
using cbm::algo::ca::EDetectorID;
- return CbmL1::Instance()->fpAlgo->GetParameters()->GetStationIndexActive(
+ return CbmL1::Instance()->fpAlgo->GetParameters().GetStationIndexActive(
CbmStsSetup::Instance()->GetStationNumber(hit->GetAddress()), EDetectorID::kSts);
}
@@ -164,7 +165,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
ca::FieldRegion fld _fvecalignment;
// fld.SetUseOriginalField();
- static int nHits = CbmL1::Instance()->fpAlgo->GetParameters()->GetNstationsActive();
+ static int nHits = CbmL1::Instance()->fpAlgo->GetParameters().GetNstationsActive();
int nTracks_SIMD = fvec::size();
ca::TrackFit fit;
@@ -175,7 +176,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
CbmStsTrack* tr[fvec::size()] {nullptr};
int ista;
- const ca::Station* sta = CbmL1::Instance()->fpAlgo->GetParameters()->GetStations().begin();
+ const ca::Station* sta = CbmL1::Instance()->fpAlgo->GetParameters().GetStations().begin();
fvec x[nHits];
fvec y[nHits];
@@ -355,7 +356,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
fit.SetMask(initialised);
fit.Extrapolate(z[i], fld);
- auto radThick = CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, fit.Tr().X(), fit.Tr().Y());
+ auto radThick = CbmL1::Instance()->fpAlgo->GetParameters().GetMaterialThickness(i, fit.Tr().X(), fit.Tr().Y());
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, -fvec::One());
@@ -416,7 +417,7 @@ void CbmL1PFFitter::Fit(std::vector<CbmStsTrack>& Tracks, const std::vector<CbmM
fit.SetMask(initialised);
fit.Extrapolate(z[i], fld);
- auto radThick = CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, fit.Tr().X(), fit.Tr().Y());
+ auto radThick = CbmL1::Instance()->fpAlgo->GetParameters().GetMaterialThickness(i, fit.Tr().X(), fit.Tr().Y());
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, fvec::One());
@@ -493,7 +494,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
int nStations = fNmvdStationsActive + fNstsStationsActive;
- const ca::Station* sta = CbmL1::Instance()->fpAlgo->GetParameters()->GetStations().begin();
+ const ca::Station* sta = CbmL1::Instance()->fpAlgo->GetParameters().GetStations().begin();
fvec* zSta = new fvec[nStations];
for (int iSta = 0; iSta < nStations; iSta++) {
zSta[iSta] = sta[iSta].fZ;
@@ -565,8 +566,8 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
}
}
- fB[0] = CbmL1::Instance()->fpAlgo->GetParameters()->GetVertexFieldValue();
- zField[0] = CbmL1::Instance()->fpAlgo->GetParameters()->GetTargetPositionZ();
+ fB[0] = CbmL1::Instance()->fpAlgo->GetParameters().GetVertexFieldValue();
+ zField[0] = CbmL1::Instance()->fpAlgo->GetParameters().GetTargetPositionZ();
fld.Set(fB[2], zField[2], fB[1], zField[1], fB[0], zField[0]);
for (int i = 0; i < nTracks_SIMD; i++) {
field.emplace_back(fld, i);
@@ -577,7 +578,7 @@ void CbmL1PFFitter::GetChiToVertex(vector<CbmStsTrack>& Tracks, vector<PFFieldRe
for (int iSt = nStations - 4; iSt >= 0; iSt--) {
fit.SetMask(T.Z() > zSta[iSt] + fvec(2.5));
fit.Extrapolate(zSta[iSt], fld);
- auto radThick = CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(iSt, fit.Tr().X(), fit.Tr().Y());
+ auto radThick = CbmL1::Instance()->fpAlgo->GetParameters().GetMaterialThickness(iSt, fit.Tr().X(), fit.Tr().Y());
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, fvec::One());
}
@@ -645,7 +646,7 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<PFF
CbmStsTrack* tr[fvec::size()];
int ista;
- const ca::Station* sta = CbmL1::Instance()->fpAlgo->GetParameters()->GetStations().begin();
+ const ca::Station* sta = CbmL1::Instance()->fpAlgo->GetParameters().GetStations().begin();
ca::FieldValue fB[3], fB_temp _fvecalignment;
fvec zField[3];
@@ -693,8 +694,8 @@ void CbmL1PFFitter::CalculateFieldRegion(vector<CbmStsTrack>& Tracks, vector<PFF
}
}
- fB[0] = CbmL1::Instance()->fpAlgo->GetParameters()->GetVertexFieldValue();
- zField[0] = CbmL1::Instance()->fpAlgo->GetParameters()->GetTargetPositionZ();
+ fB[0] = CbmL1::Instance()->fpAlgo->GetParameters().GetVertexFieldValue();
+ zField[0] = CbmL1::Instance()->fpAlgo->GetParameters().GetTargetPositionZ();
fld.Set(fB[2], zField[2], fB[1], zField[1], fB[0], zField[0]);
for (int i = 0; i < nTracks_SIMD; i++) {
field.emplace_back(fld, i);
@@ -716,7 +717,7 @@ void CbmL1PFFitter::CalculateFieldRegionAtLastPoint(vector<CbmStsTrack>& Tracks,
CbmStsTrack* tr[fvec::size()];
int ista;
- const ca::Station* sta = CbmL1::Instance()->fpAlgo->GetParameters()->GetStations().begin();
+ const ca::Station* sta = CbmL1::Instance()->fpAlgo->GetParameters().GetStations().begin();
ca::FieldValue fB[3], fB_temp _fvecalignment;
fvec zField[3];
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 3ab91e31e3..773cd40127 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -29,13 +29,7 @@ set(SRCS
OffLineInterface/CbmL1StsTrackFinder.cxx
OffLineInterface/CbmL1GlobalTrackFinder.cxx
OffLineInterface/CbmL1GlobalFindTracksEvents.cxx
- L1Algo/utils/CaAlgoRandom.cxx
- L1Algo/L1Algo.cxx
- L1Algo/L1TripletConstructor.cxx
- L1Algo/L1CaTrackFinder.cxx
- L1Algo/L1CaTrackFinderSlice.cxx
- L1Algo/L1BranchExtender.cxx
- L1Algo/L1TrackFitter.cxx
+
CbmL1Util.cxx
CbmL1Performance.cxx
CbmL1ReadEvent.cxx
@@ -43,12 +37,13 @@ set(SRCS
CbmCaTimeSliceReader.cxx
CbmL1MCTrack.cxx
CbmL1Track.cxx
- L1Algo/L1CloneMerger.cxx
+
L1Algo/utils/L1AlgoDraw.cxx
L1Algo/utils/L1AlgoEfficiencyPerformance.cxx
L1Algo/utils/L1AlgoPulls.cxx
L1Algo/utils/CaUvConverter.cxx
-
+ L1Algo/utils/CaAlgoRandom.cxx
+
catools/CaToolsHitRecord.cxx
catools/CaToolsMCData.cxx
catools/CaToolsMCPoint.cxx
@@ -187,6 +182,3 @@ install(FILES CbmL1Counters.h
DESTINATION include
)
-install(FILES L1Algo/L1Algo.h
- DESTINATION include/L1Algo
-)
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 4cf9dec8e7..0a5df93def 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -64,9 +64,10 @@
#include <iostream>
#include <list>
+#include "CaFramework.h"
#include "CaHit.h"
+#include "CaToolsDebugger.h"
#include "CaToolsField.h"
-#include "L1Algo/L1Algo.h"
using namespace cbm::algo; // TODO: remove this line
@@ -79,7 +80,7 @@ using cbm::ca::tools::MaterialHelper;
ClassImp(CbmL1);
-static L1Algo gAlgo _fvecalignment; // TODO: Change coupling logic between L1Algo and CbmL1
+static ca::Framework gAlgo _fvecalignment; // TODO: Change coupling logic between ca::Framework and CbmL1
CbmL1* CbmL1::fpInstance = 0;
@@ -196,7 +197,7 @@ InitStatus CbmL1::Init()
FairRootManager* fairManager = FairRootManager::Instance();
- if (L1Algo::TrackingMode::kSts == fTrackingMode) {
+ if (ca::Framework::TrackingMode::kSts == fTrackingMode) {
fUseMVD = 1;
fUseSTS = 1;
fUseMUCH = 0;
@@ -207,7 +208,7 @@ InitStatus CbmL1::Init()
if (findTask) fUseMVD = findTask->MvdUsage();
}
- if (L1Algo::TrackingMode::kMcbm == fTrackingMode) {
+ if (ca::Framework::TrackingMode::kMcbm == fTrackingMode) {
fUseMVD = 1;
fUseSTS = 1;
fUseMUCH = 1;
@@ -216,7 +217,7 @@ InitStatus CbmL1::Init()
fInitManager.DevSetIgnoreHitSearchAreas(true); // uncomment for debug
}
- if (L1Algo::TrackingMode::kGlobal == fTrackingMode) {
+ if (ca::Framework::TrackingMode::kGlobal == fTrackingMode) {
//at the moment trd2d tracking only
fUseMVD = 1;
fUseSTS = 1;
@@ -413,9 +414,9 @@ InitStatus CbmL1::Init()
std::string mainConfig = std::string(gSystem->Getenv("VMCWORKDIR")) + "/macro/L1/configs/";
switch (fTrackingMode) {
- case L1Algo::TrackingMode::kSts: mainConfig += "ca_params_sts.yaml"; break;
- case L1Algo::TrackingMode::kMcbm: mainConfig += "ca_params_mcbm.yaml"; break;
- case L1Algo::TrackingMode::kGlobal: mainConfig += "ca_params_global.yaml"; break;
+ case ca::Framework::TrackingMode::kSts: mainConfig += "ca_params_sts.yaml"; break;
+ case ca::Framework::TrackingMode::kMcbm: mainConfig += "ca_params_mcbm.yaml"; break;
+ case ca::Framework::TrackingMode::kGlobal: mainConfig += "ca_params_global.yaml"; break;
}
fInitManager.SetConfigMain(mainConfig);
}
@@ -502,7 +503,7 @@ InitStatus CbmL1::Init()
// TODO: SZh 15.08.2022: Replace station type with ca::EDetectorID
stationInfo.SetStationType(1); // MVD
stationInfo.SetTimeInfo(mvdInterface->IsTimeInfoProvided(iSt));
- stationInfo.SetFieldStatus(fTrackingMode == L1Algo::TrackingMode::kMcbm ? 0 : 1);
+ stationInfo.SetFieldStatus(fTrackingMode == ca::Framework::TrackingMode::kMcbm ? 0 : 1);
stationInfo.SetZref(mvdInterface->GetZref(iSt));
stationInfo.SetZmin(mvdInterface->GetZmin(iSt));
stationInfo.SetZmax(mvdInterface->GetZmax(iSt));
@@ -525,9 +526,9 @@ InitStatus CbmL1::Init()
// TODO: SZh 15.08.2022: Replace station type with ca::EDetectorID
stationInfo.SetStationType(0); // STS
stationInfo.SetTimeInfo(stsInterface->IsTimeInfoProvided(iSt));
- stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt)
- : false);
- stationInfo.SetFieldStatus(L1Algo::TrackingMode::kMcbm == fTrackingMode ? 0 : 1);
+ stationInfo.SetTimeInfo(
+ ca::Framework::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt) : false);
+ stationInfo.SetFieldStatus(ca::Framework::TrackingMode::kMcbm == fTrackingMode ? 0 : 1);
stationInfo.SetZref(stsInterface->GetZref(iSt));
stationInfo.SetZmin(stsInterface->GetZmin(iSt));
stationInfo.SetZmax(stsInterface->GetZmax(iSt));
@@ -551,8 +552,8 @@ InitStatus CbmL1::Init()
// TODO: SZh 15.08.2022: Replace station type with ca::EDetectorID
stationInfo.SetStationType(2); // MuCh
stationInfo.SetTimeInfo(muchInterface->IsTimeInfoProvided(iSt));
- stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt)
- : false);
+ stationInfo.SetTimeInfo(
+ ca::Framework::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt) : false);
stationInfo.SetFieldStatus(0);
stationInfo.SetZref(muchInterface->GetZref(iSt));
stationInfo.SetZmin(muchInterface->GetZmin(iSt));
@@ -577,8 +578,8 @@ InitStatus CbmL1::Init()
// TODO: SZh 15.08.2022: Replace station type with ca::EDetectorID
stationInfo.SetStationType((iSt == 1 || iSt == 3) ? 6 : 3); // MuCh
stationInfo.SetTimeInfo(trdInterface->IsTimeInfoProvided(iSt));
- stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt)
- : false);
+ stationInfo.SetTimeInfo(
+ ca::Framework::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt) : false);
stationInfo.SetFieldStatus(0);
stationInfo.SetZref(trdInterface->GetZref(iSt));
stationInfo.SetZmin(trdInterface->GetZmin(iSt));
@@ -586,7 +587,7 @@ InitStatus CbmL1::Init()
stationInfo.SetXmax(trdInterface->GetXmax(iSt));
stationInfo.SetYmax(trdInterface->GetYmax(iSt));
- if (L1Algo::TrackingMode::kGlobal == fTrackingMode) { stationInfo.SetTimeInfo(false); }
+ if (ca::Framework::TrackingMode::kGlobal == fTrackingMode) { stationInfo.SetTimeInfo(false); }
stationInfo.SetTrackingStatus(true);
if (fvmDisabledStationIDs[ca::EDetectorID::kTrd].find(iSt)
!= fvmDisabledStationIDs[ca::EDetectorID::kTrd].cend()) {
@@ -604,8 +605,8 @@ InitStatus CbmL1::Init()
// TODO: SZh 15.08.2022: Replace station type with ca::EDetectorID
stationInfo.SetStationType(4);
stationInfo.SetTimeInfo(tofInterface->IsTimeInfoProvided(iSt));
- stationInfo.SetTimeInfo(L1Algo::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt)
- : false);
+ stationInfo.SetTimeInfo(
+ ca::Framework::TrackingMode::kMcbm != fTrackingMode ? stsInterface->IsTimeInfoProvided(iSt) : false);
stationInfo.SetFieldStatus(0);
stationInfo.SetZref(tofInterface->GetZref(iSt));
stationInfo.SetZmin(tofInterface->GetZmin(iSt));
@@ -654,27 +655,27 @@ InitStatus CbmL1::Init()
// Init L1 algo core
// FIXME: SZh 22.08.2023:
- // Re-organize the the relation between CbmL1 and L1Algo. I belive, we don't need a global pointer to the L1Algo
+ // Re-organize the the relation between CbmL1 and ca::Framework. I belive, we don't need a global pointer to the ca::Framework
// instance.
fpAlgo = &gAlgo;
fpAlgo->Init(fTrackingMode);
//
- // ** Send formed parameters object to L1Algo instance **
+ // ** Send formed parameters object to ca::Framework instance **
//
fpAlgo->ReceiveParameters(fInitManager.TakeParameters());
/*** Get numbers of active stations ***/
- fNMvdStations = fpAlgo->GetParameters()->GetNstationsActive(ca::EDetectorID::kMvd);
- fNStsStations = fpAlgo->GetParameters()->GetNstationsActive(ca::EDetectorID::kSts);
- fNTrdStations = fpAlgo->GetParameters()->GetNstationsActive(ca::EDetectorID::kTrd);
- fNMuchStations = fpAlgo->GetParameters()->GetNstationsActive(ca::EDetectorID::kMuch);
- fNTofStations = fpAlgo->GetParameters()->GetNstationsActive(ca::EDetectorID::kTof);
- fNStations = fpAlgo->GetParameters()->GetNstationsActive();
+ fNMvdStations = fpAlgo->GetParameters().GetNstationsActive(ca::EDetectorID::kMvd);
+ fNStsStations = fpAlgo->GetParameters().GetNstationsActive(ca::EDetectorID::kSts);
+ fNTrdStations = fpAlgo->GetParameters().GetNstationsActive(ca::EDetectorID::kTrd);
+ fNMuchStations = fpAlgo->GetParameters().GetNstationsActive(ca::EDetectorID::kMuch);
+ fNTofStations = fpAlgo->GetParameters().GetNstationsActive(ca::EDetectorID::kTof);
+ fNStations = fpAlgo->GetParameters().GetNstationsActive();
- LOG(info) << fpAlgo->GetParameters()->ToString(1);
+ LOG(info) << fpAlgo->GetParameters().ToString(1);
fpIODataManager->SetNofActiveStations(fNStations);
LOG(info) << "----- Numbers of stations active in tracking -----";
@@ -690,7 +691,7 @@ InitStatus CbmL1::Init()
LOG(info) << "\033[31;1m-------------------- L1 material -----------------------------\033[0m";
fMaterialMonitor.clear();
for (int i = 0; i < fNStations; i++) {
- ca::MaterialMonitor m(&(fpAlgo->GetParameters()->GetThicknessMaps()[i]), Form("station %d", i));
+ ca::MaterialMonitor m(&(fpAlgo->GetParameters().GetThicknessMaps()[i]), Form("station %d", i));
LOG(info) << m.ToString();
fMaterialMonitor.push_back(m);
}
@@ -757,7 +758,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
if (fVerbose > 1) { cout << "\nCbmL1::Exec event " << fEventNo << " ...\n\n"; }
- // ----- Read data from branches and send data from IODataManager to L1Algo ----------------------------------------
+ // ----- Read data from branches and send data from IODataManager to ca::Framework ----------------------------------------
ReadEvent(event);
//if (fPerformance) {
@@ -794,7 +795,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
// FieldIntegralCheck();
if (fVerbose > 1) { LOG(info) << "L1 Track finder..."; }
- fpAlgo->CaTrackFinder();
+ fpAlgo->fTrackFinder.FindTracks();
// IdealTrackFinder();
fTrackingTime = fpAlgo->fCaRecoTime; // TODO: remove (not used)
@@ -958,6 +959,7 @@ void CbmL1::Finish()
gFile = currentFile;
gDirectory = curr;
fpAlgo->Finish();
+ cbm::ca::tools::Debugger::Instance().Write();
}
@@ -1043,7 +1045,7 @@ void CbmL1::IdealTrackFinder()
CaTrack algoTr;
algoTr.fNofHits = 0;
- ca::Vector<int> hitIndices("CbmL1::hitIndices", fpAlgo->GetParameters()->GetNstationsActive(), -1);
+ ca::Vector<int> hitIndices("CbmL1::hitIndices", fpAlgo->GetParameters().GetNstationsActive(), -1);
for (unsigned int iH = 0; iH < MC.Hits.size(); iH++) {
const int hitI = MC.Hits[iH];
@@ -1053,7 +1055,7 @@ void CbmL1::IdealTrackFinder()
}
- for (int iH = 0; iH < fpAlgo->GetParameters()->GetNstationsActive(); iH++) {
+ for (int iH = 0; iH < fpAlgo->GetParameters().GetNstationsActive(); iH++) {
const int hitI = hitIndices[iH];
if (hitI < 0) continue;
@@ -1064,7 +1066,7 @@ void CbmL1::IdealTrackFinder()
if (algoTr.fNofHits < 3) continue;
- for (int iH = 0; iH < fpAlgo->GetParameters()->GetNstationsActive(); iH++) {
+ for (int iH = 0; iH < fpAlgo->GetParameters().GetNstationsActive(); iH++) {
const int hitI = hitIndices[iH];
if (hitI < 0) continue;
fpAlgo->fRecoHits.push_back(hitI);
@@ -1283,7 +1285,7 @@ void CbmL1::DumpMaterialToFile(TString fileName)
{
TFile* f = new TFile(fileName, "RECREATE");
f->cd();
- const auto& param = *fpAlgo->GetParameters();
+ const auto& param = fpAlgo->GetParameters();
const auto& map = param.GetThicknessMaps();
for (int ist = 0; ist < param.GetNstationsActive(); ist++) {
//const ca::Station& st = param.GetStation(ist);
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index ec4bba1492..30d34cbd70 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -54,13 +54,13 @@
#include "AlgoFairloggerCompat.h"
#include "CaDataManager.h"
+#include "CaFramework.h"
#include "CaInitManager.h"
#include "CaMaterialMonitor.h"
#include "CaMonitor.h"
+#include "CaTrackParam.h"
#include "CaVector.h"
-#include "L1Algo/L1Algo.h"
-class L1Algo;
class CbmL1MCTrack;
class CbmMCDataObject;
class CbmEvent;
@@ -106,10 +106,10 @@ namespace std
// TODO: insert documentation! (S.Zh.)
//
-/// L1Algo runtime constants modification can be performed in run_reco.C. Example:
+/// ca::Framework runtime constants modification can be performed in run_reco.C. Example:
///
/// l1->GetInitManager()->GetParameters()->SetMaxDoubletsPerSinglet(149);
-/// TODO: L1InitManager - main interface of communication between cbmroot/bmnroot and L1Algo (S.Zharko)
+/// TODO: L1InitManager - main interface of communication between cbmroot/bmnroot and ca::Framework (S.Zharko)
///
class CbmL1 : public FairTask {
public:
@@ -254,9 +254,9 @@ public:
}
void SetExtrapolateToTheEndOfSTS(bool b) { fExtrapolateToTheEndOfSTS = b; }
- void SetStsOnlyMode() { fTrackingMode = L1Algo::TrackingMode::kSts; }
- void SetMcbmMode() { fTrackingMode = L1Algo::TrackingMode::kMcbm; }
- void SetGlobalMode() { fTrackingMode = L1Algo::TrackingMode::kGlobal; }
+ void SetStsOnlyMode() { fTrackingMode = ca::Framework::TrackingMode::kSts; }
+ void SetMcbmMode() { fTrackingMode = ca::Framework::TrackingMode::kMcbm; }
+ void SetGlobalMode() { fTrackingMode = ca::Framework::TrackingMode::kGlobal; }
// Tracking monitor (prototype)
enum class EMonitorKey
@@ -408,7 +408,7 @@ private:
/// \note Should be called only after CbmL1::Performance()
void TrackFitPerformance();
- void FillFitHistos(TrackParamV& tr, const CbmL1MCPoint& mc, bool isTimeFitted, TH1F* h[]);
+ void FillFitHistos(cbm::algo::ca::TrackParamV& tr, const CbmL1MCPoint& mc, bool isTimeFitted, TH1F* h[]);
/// Fills performance histograms
void HistoPerformance();
@@ -488,9 +488,10 @@ private:
public:
// ** Basic data members **
- L1Algo* fpAlgo = nullptr; ///< Pointer to the L1 track finder algorithm
+ ca::Framework* fpAlgo = nullptr; ///< Pointer to the L1 track finder algorithm
- L1Algo::TrackingMode fTrackingMode = L1Algo::TrackingMode::kSts; ///< Tracking mode: kSts, kMcbm or kGlobal
+ ca::Framework::TrackingMode fTrackingMode =
+ ca::Framework::TrackingMode::kSts; ///< Tracking mode: kSts, kMcbm or kGlobal
DFSET fvSelectedMcEvents {}; ///< Set of selected MC events with fileID and eventId
diff --git a/reco/L1/CbmL1MCTrack.cxx b/reco/L1/CbmL1MCTrack.cxx
index 0b19922af3..331177e256 100644
--- a/reco/L1/CbmL1MCTrack.cxx
+++ b/reco/L1/CbmL1MCTrack.cxx
@@ -28,7 +28,7 @@
#include "CaConstants.h"
#include "CaHit.h"
-#include "L1Algo/L1Algo.h"
+
using cbm::algo::ca::constants::size::MaxNstations;
@@ -123,7 +123,7 @@ void CbmL1MCTrack::CountHitStations()
nHitContStations = 0;
int ncont = 0;
- for (int ista = 0; ista < L1->fpAlgo->GetParameters()->GetNstationsActive(); ista++) {
+ for (int ista = 0; ista < L1->fpAlgo->GetParameters().GetNstationsActive(); ista++) {
if (maxNStaHits < stationNhits[ista]) { maxNStaHits = stationNhits[ista]; }
if (stationNhits[ista] > 0) {
nStations++;
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 8cf2daa8f5..182f670d4e 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -55,9 +55,12 @@
#include <cmath>
+#include "CaFramework.h"
#include "CaToolsDebugger.h"
#include "CaTrackFit.h"
-#include "L1Algo.h"
+#include "CaTrackParam.h"
+
+using namespace cbm::algo::ca;
using cbm::ca::tools::Debugger;
using std::cout;
@@ -335,7 +338,7 @@ void CbmL1::EfficienciesPerformance()
cout << " n mc stations: " << t.NMCStations() << endl;
}
for (unsigned int i = 0; i < rtraIt->Hits.size(); i++) {
- const ca::Hit& h = fpAlgo->fInputData.GetHit(rtraIt->Hits[i]);
+ const ca::Hit& h = fpAlgo->GetInputData().GetHit(rtraIt->Hits[i]);
const CbmL1HitDebugInfo& s = fvHitDebugInfo[rtraIt->Hits[i]];
cout << " x y z t " << s.x << " " << s.y << " " << h.Z() << " dx " << s.dx << " dy " << s.dy << std::endl;
cbm::ca::tools::Debugger::Instance().FillNtuple("ghost", statNghost, i, fabs(1. / tr.GetQp()[0]), s.x, s.y,
@@ -346,10 +349,10 @@ void CbmL1::EfficienciesPerformance()
}
}
- int sta_nhits[fpAlgo->GetParameters()->GetNstationsActive()];
- int sta_nfakes[fpAlgo->GetParameters()->GetNstationsActive()];
+ int sta_nhits[fpAlgo->GetParameters().GetNstationsActive()];
+ int sta_nfakes[fpAlgo->GetParameters().GetNstationsActive()];
- for (int i = 0; i < fpAlgo->GetParameters()->GetNstationsActive(); i++) {
+ for (int i = 0; i < fpAlgo->GetParameters().GetNstationsActive(); i++) {
sta_nhits[i] = 0;
sta_nfakes[i] = 0;
}
@@ -486,7 +489,7 @@ void CbmL1::EfficienciesPerformance()
if (fVerbose > 1) {
ntra.PrintEff(true, true);
cout << "Number of true and fake hits in stations: " << endl;
- for (int i = 0; i < fpAlgo->GetParameters()->GetNstationsActive(); i++) {
+ for (int i = 0; i < fpAlgo->GetParameters().GetNstationsActive(); i++) {
cout << sta_nhits[i] - sta_nfakes[i] << "+" << sta_nfakes[i] << " ";
}
cout << endl;
@@ -839,16 +842,16 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
int iFstHit = prtra->GetFirstHitIndex();
auto& fstHit = fvHitDebugInfo[iFstHit];
- h2_fst_hit_yz->Fill(fpAlgo->GetParameters()->GetStation(fstHit.iStation).fZ[0], fstHit.GetY());
+ h2_fst_hit_yz->Fill(fpAlgo->GetParameters().GetStation(fstHit.iStation).fZ[0], fstHit.GetY());
int iLstHit = prtra->GetLastHitIndex();
auto& lstHit = fvHitDebugInfo[iLstHit];
- h2_lst_hit_yz->Fill(fpAlgo->GetParameters()->GetStation(lstHit.iStation).fZ[0], lstHit.GetY());
+ h2_lst_hit_yz->Fill(fpAlgo->GetParameters().GetStation(lstHit.iStation).fZ[0], lstHit.GetY());
for (int iH : prtra->Hits) {
const auto& hit = fvHitDebugInfo[iH];
int y = hit.GetY();
- int z = fpAlgo->GetParameters()->GetStation(hit.iStation).fZ[0];
+ int z = fpAlgo->GetParameters().GetStation(hit.iStation).fZ[0];
h2_all_hit_yz->Fill(z, y);
}
}
@@ -888,8 +891,8 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
CbmL1HitDebugInfo& h2 = fvHitDebugInfo[prtra->Hits[1]];
h_ghost_fstation->Fill(h1.iStation);
h_ghost_r->Fill(sqrt(fabs(h1.x * h1.x + h1.y * h1.y)));
- double z1 = fpAlgo->GetParameters()->GetStation(h1.iStation).fZ[0];
- double z2 = fpAlgo->GetParameters()->GetStation(h2.iStation).fZ[0];
+ double z1 = fpAlgo->GetParameters().GetStation(h1.iStation).fZ[0];
+ double z2 = fpAlgo->GetParameters().GetStation(h2.iStation).fZ[0];
if (fabs(z2 - z1) > 1.e-4) {
h_ghost_tx->Fill((h2.x - h1.x) / (z2 - z1));
h_ghost_ty->Fill((h2.y - h1.y) / (z2 - z1));
@@ -1092,8 +1095,8 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
// CbmL1HitDebugInfo &ph21 = fvHitDebugInfo[mtra.Hits[0]];
// CbmL1HitDebugInfo &ph22 = fvHitDebugInfo[mtra.Hits[1]];
- // double z21 = fpAlgo->GetParameters()->GetStation(ph21.iStation).fZ[0];
- // double z22 = fpAlgo->GetParameters()->GetStation(ph22.iStation).fZ[0];
+ // double z21 = fpAlgo->GetParameters().GetStation(ph21.iStation).fZ[0];
+ // double z22 = fpAlgo->GetParameters().GetStation(ph22.iStation).fZ[0];
// if( fabs(z22-z21)>1.e-4 ){
// h_notfound_tx->Fill((ph22.x-ph21.x)/(z22-z21));
// h_notfound_ty->Fill((ph22.y-ph21.y)/(z22-z21));
@@ -1221,7 +1224,7 @@ void CbmL1::TrackFitPerformance()
{"pvi", "Pull Vi [residual/estimated_error]", 100, -10., 10.},
{"distrVi", "Vi distribution [1/c]", 100, 0., 4.}};
- if (L1Algo::kGlobal == fpAlgo->fTrackingMode) {
+ if (ca::Framework::kGlobal == fpAlgo->fTrackingMode) {
Table[4].l = -1.;
Table[4].r = 1.;
}
@@ -1254,7 +1257,7 @@ void CbmL1::TrackFitPerformance()
{"pvi", "Pull Vi [residual/estimated_error]", 100, -10., 10.},
{"distrVi", "Vi distribution [1/c]", 100, -10., 10.}};
- if (L1Algo::kGlobal == fpAlgo->fTrackingMode) {
+ if (ca::Framework::kGlobal == fpAlgo->fTrackingMode) {
TableVertex[4].l = -1.;
TableVertex[4].r = 1.;
}
@@ -1296,7 +1299,7 @@ void CbmL1::TrackFitPerformance()
int nTimeMeasurenments = 0;
for (unsigned int ih = 0; ih < it->Hits.size(); ih++) {
int ista = fvHitDebugInfo[it->Hits[ih]].iStation;
- if (fpAlgo->GetParameters()->GetStation(ista).timeInfo) { nTimeMeasurenments++; }
+ if (fpAlgo->GetParameters().GetStation(ista).timeInfo) { nTimeMeasurenments++; }
}
isTimeFitted = (nTimeMeasurenments > 1);
}
@@ -1382,14 +1385,14 @@ void CbmL1::TrackFitPerformance()
fit.Extrapolate(mc.z, fld);
// add material
const int fSta = fvHitDebugInfo[it->Hits[0]].iStation;
- const int dir = int((mc.z - fpAlgo->GetParameters()->GetStation(fSta).fZ[0])
- / fabs(mc.z - fpAlgo->GetParameters()->GetStation(fSta).fZ[0]));
- // if (abs(mc.z - fpAlgo->GetParameters()->GetStation(fSta).fZ[0]) > 10.) continue; // can't extrapolate on large distance
+ const int dir = int((mc.z - fpAlgo->GetParameters().GetStation(fSta).fZ[0])
+ / fabs(mc.z - fpAlgo->GetParameters().GetStation(fSta).fZ[0]));
+ // if (abs(mc.z - fpAlgo->GetParameters().GetStation(fSta).fZ[0]) > 10.) continue; // can't extrapolate on large distance
for (int iSta = fSta /*+dir*/; (iSta >= 0) && (iSta < fNStations)
- && (dir * (mc.z - fpAlgo->GetParameters()->GetStation(iSta).fZ[0]) > 0);
+ && (dir * (mc.z - fpAlgo->GetParameters().GetStation(iSta).fZ[0]) > 0);
iSta += dir) {
// cout << iSta << " " << dir << endl;
- auto radThick = fpAlgo->GetParameters()->GetMaterialThickness(iSta, fit.Tr().GetX(), fit.Tr().GetY());
+ auto radThick = fpAlgo->GetParameters().GetMaterialThickness(iSta, fit.Tr().GetX(), fit.Tr().GetY());
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, fvec::One());
}
@@ -1446,16 +1449,16 @@ void CbmL1::TrackFitPerformance()
// add material
const int fSta = fvHitDebugInfo[it->Hits[0]].iStation;
- const int dir = (mc.z - fpAlgo->GetParameters()->GetStation(fSta).fZ[0])
- / abs(mc.z - fpAlgo->GetParameters()->GetStation(fSta).fZ[0]);
- // if (abs(mc.z - fpAlgo->GetParameters()->GetStation(fSta].fZ[0]) > 10.) continue; // can't extrapolate on large distance
+ const int dir = (mc.z - fpAlgo->GetParameters().GetStation(fSta).fZ[0])
+ / abs(mc.z - fpAlgo->GetParameters().GetStation(fSta).fZ[0]);
+ // if (abs(mc.z - fpAlgo->GetParameters().GetStation(fSta].fZ[0]) > 10.) continue; // can't extrapolate on large distance
for (int iSta = fSta + dir; (iSta >= 0) && (iSta < fNStations)
- && (dir * (mc.z - fpAlgo->GetParameters()->GetStation(iSta).fZ[0]) > 0);
+ && (dir * (mc.z - fpAlgo->GetParameters().GetStation(iSta).fZ[0]) > 0);
iSta += dir) {
- fit.Extrapolate(fpAlgo->GetParameters()->GetStation(iSta).fZ, fld);
- auto radThick = fpAlgo->GetParameters()->GetMaterialThickness(iSta, fit.Tr().GetX(), fit.Tr().GetY());
+ fit.Extrapolate(fpAlgo->GetParameters().GetStation(iSta).fZ, fld);
+ auto radThick = fpAlgo->GetParameters().GetMaterialThickness(iSta, fit.Tr().GetX(), fit.Tr().GetY());
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, fvec::One());
}
@@ -1632,9 +1635,9 @@ void CbmL1::FieldApproxCheck()
FairField* MF = FairRunAna::Instance()->GetField();
assert(MF);
- for (int ist = 0; ist < fpAlgo->GetParameters()->GetNstationsActive(); ist++) {
+ for (int ist = 0; ist < fpAlgo->GetParameters().GetNstationsActive(); ist++) {
- const ca::Station& st = fpAlgo->GetParameters()->GetStation(ist);
+ const ca::Station& st = fpAlgo->GetParameters().GetStation(ist);
double z = st.fZ[0];
double Xmax = st.Xmax[0];
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index c01b291e0d..f37f73ec4f 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -48,8 +48,8 @@
#include <iostream>
+#include "CaFramework.h"
#include "CaToolsDebugger.h"
-#include "L1Algo/L1Algo.h"
using std::cout;
using std::endl;
@@ -383,7 +383,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
CbmMvdHit* h = dynamic_cast<CbmMvdHit*>(fpMvdHits->At(hitIndex));
{
th.ExtIndex = hitIndex;
- th.iStation = fpAlgo->GetParameters()->GetStationIndexActive(
+ th.iStation = fpAlgo->GetParameters().GetStationIndexActive(
CbmMvdTrackingInterface::Instance()->GetTrackingStationIndex(h), ca::EDetectorID::kMvd);
if (th.iStation < 0) continue;
@@ -449,7 +449,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
th.ExtIndex = hitIndex;
th.Det = 1;
- th.iStation = fpAlgo->GetParameters()->GetStationIndexActive(
+ th.iStation = fpAlgo->GetParameters().GetStationIndexActive(
CbmStsTrackingInterface::Instance()->GetTrackingStationIndex(h), ca::EDetectorID::kSts);
if (th.iStation < 0) continue;
@@ -511,7 +511,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
th.ExtIndex = hitIndex;
th.Det = 2;
- th.iStation = fpAlgo->GetParameters()->GetStationIndexActive(
+ th.iStation = fpAlgo->GetParameters().GetStationIndexActive(
CbmMuchTrackingInterface::Instance()->GetTrackingStationIndex(h), ca::EDetectorID::kMuch);
if (th.iStation < 0) continue;
@@ -567,7 +567,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
Int_t hitIndex = (event ? event->GetIndex(ECbmDataType::kTrdHit, iHit) : iHit);
CbmTrdHit* h = dynamic_cast<CbmTrdHit*>(fpTrdHits->At(hitIndex));
- if ((L1Algo::TrackingMode::kGlobal == fTrackingMode) && (int) h->GetClassType() != 1) {
+ if ((ca::Framework::TrackingMode::kGlobal == fTrackingMode) && (int) h->GetClassType() != 1) {
// SGtrd2d!! skip TRD-1D hit
continue;
}
@@ -575,7 +575,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
th.ExtIndex = hitIndex;
th.Det = 3;
- th.iStation = fpAlgo->GetParameters()->GetStationIndexActive(
+ th.iStation = fpAlgo->GetParameters().GetStationIndexActive(
CbmTrdTrackingInterface::Instance()->GetTrackingStationIndex(h), ca::EDetectorID::kTrd);
if (th.iStation < 0) continue;
@@ -645,7 +645,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
if (5 == CbmTofAddress::GetSmType(h->GetAddress())) continue; // Skip T0 hits from TOF
- th.iStation = fpAlgo->GetParameters()->GetStationIndexActive(
+ th.iStation = fpAlgo->GetParameters().GetStationIndexActive(
CbmTofTrackingInterface::Instance()->GetTrackingStationIndex(h), ca::EDetectorID::kTof);
if (th.iStation < 0) continue;
@@ -673,7 +673,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmTofTrackingInterface::Instance()->GetHitRanges(*h);
//TODO: is it still needed here? (S.Zharko)
- if (L1Algo::TrackingMode::kMcbm == fTrackingMode && th.z > 400) continue;
+ if (ca::Framework::TrackingMode::kMcbm == fTrackingMode && th.z > 400) continue;
if (fPerformance) { std::tie(th.iBestMc, th.vMc) = MatchHitWithMc<ca::EDetectorID::kTof>(hitIndex); }
@@ -781,7 +781,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
if (fVerbose >= 2) cout << "ReadEvent: mvd and sts are saved." << endl;
- // ----- Send data from IODataManager to L1Algo --------------------------------------------------------------------
+ // ----- Send data from IODataManager to ca::Framework --------------------------------------------------------------------
if (1 == fSTAPDataMode) { WriteSTAPAlgoInputData(nCalls); }
if (2 == fSTAPDataMode) { ReadSTAPAlgoInputData(nCalls); }
// TODO: SZh: If we read data from file, we don't need to collect them above. This should be addressed
@@ -932,7 +932,7 @@ void CbmL1::ReadMCPoint(ca::EDetectorID iDet, int file, int event, int iPoint)
CbmL1MCPoint MC;
- MC.iStation = fpAlgo->GetParameters()->GetStationIndexActive(iStLoc, iDet);
+ MC.iStation = fpAlgo->GetParameters().GetStationIndexActive(iStLoc, iDet);
if (MC.iStation < 0) { return; }
TVector3 xyzI, PI;
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
deleted file mode 100644
index 41bd91f5c3..0000000000
--- a/reco/L1/L1Algo/L1Algo.h
+++ /dev/null
@@ -1,435 +0,0 @@
-/* Copyright (C) 2007-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Maksym Zyzak, Valentina Akishina, Igor Kulakov [committer], Sergei Zharko */
-
-#ifndef L1Algo_h
-#define L1Algo_h
-
-/// Debug features
-// #define PULLS // triplets pulls
-// #define TRIP_PERFORMANCE // triplets efficiencies
-// #define DOUB_PERFORMANCE // doublets efficiencies
-// #define DRAW // event display
-#ifdef DRAW
-class L1AlgoDraw;
-#endif
-//#define XXX // time debug
-//#define COUNTERS // diff counters (hits, doublets, ... )
-
-//#define GLOBAL
-//#define mCBM
-
-#include <array>
-#include <iomanip>
-#include <iostream>
-#include <limits>
-#include <map>
-
-#include "CaBranch.h"
-#include "CaConstants.h"
-#include "CaField.h"
-#include "CaGrid.h"
-#include "CaGridEntry.h"
-#include "CaHit.h"
-#include "CaInputData.h"
-#include "CaMeasurementXy.h"
-#include "CaParameters.h"
-#include "CaStation.h"
-#include "CaTimer.h"
-#include "CaTrack.h"
-#include "CaTrackParam.h"
-#include "CaTrackingMonitor.h"
-#include "CaTriplet.h"
-#include "CaVector.h"
-#include "L1CloneMerger.h"
-
-using namespace cbm::algo::ca; //TODO: remove
-
-class CbmL1MCTrack;
-
-// *******************************
-// ** Types definition (global) **
-// *******************************
-
-using L1StationsArray_t = std::array<ca::Station, constants::size::MaxNstations>;
-using L1MaterialArray_t = std::array<ca::MaterialMap, constants::size::MaxNstations>;
-using Tindex = int; // TODO: Replace with ca::HitIndex_t, if suitable
-
-#ifdef PULLS
-#define TRIP_PERFORMANCE
-class L1AlgoPulls;
-#endif
-#ifdef TRIP_PERFORMANCE
-template<Tindex NHits>
-class L1AlgoEfficiencyPerformance;
-#endif
-#ifdef DOUB_PERFORMANCE
-template<Tindex NHits>
-class L1AlgoEfficiencyPerformance;
-#endif
-
-struct L1HitTimeInfo {
- fscal fEventTimeMin {-std::numeric_limits<fscal>::max() / 2.};
- fscal fEventTimeMax {std::numeric_limits<fscal>::max() / 2.};
- fscal fMaxTimeBeforeHit {0.}; //< max event time for hits [0 .. hit] in the station hit array
- fscal fMinTimeAfterHit {0.}; //< min event time for hits [hit ... ] in the station hit array
-};
-
-/// Main class of L1 CA track finder algorithm
-///
-class L1Algo {
-public:
- // *************************
- // ** Friend classes list **
- // *************************
-
- friend class CbmL1; /// TODO: Remove this dependency
- friend class L1TripletConstructor;
-#ifdef DRAW
- friend class L1AlgoDraw;
-#endif
- // ***************************
- // ** Internal enumerations **
- // ***************************
-
- /// TODO: Move to L1
- enum TrackingMode
- {
- kSts,
- kGlobal,
- kMcbm
- };
-
-
- // **********************************
- // ** Member functions declaration **
- // **********************************
-
- // ** Constructors and destructor
-
- /// Constructor
- L1Algo();
-
- /// Copy constructor
- L1Algo(const L1Algo&) = delete;
-
- /// Move constructor
- L1Algo(L1Algo&&) = delete;
-
- /// Copy assignment operator
- L1Algo& operator=(const L1Algo&) = delete;
-
- /// Move assignment operator
- L1Algo& operator=(L1Algo&&) = delete;
-
- /// Destructor
- ~L1Algo() = default;
-
-
- // ** Functions, which pack and unpack indexes of station and triplet **
-
- // TODO: move to ca::Triplet class (S.Zharko)
- /// Packs station and triplet indices to an unique triplet ID
- /// \param iStation Index of station in the active stations array
- /// \param iTriplet Index of triplet
- static unsigned int PackTripletId(unsigned int iStation, unsigned int iTriplet);
-
- /// Unpacks the triplet ID to its station index
- /// \param id Unique triplet ID
- static unsigned int TripletId2Station(unsigned int id);
-
- /// Unpacks the triplet ID to its triplet index
- /// \param id Unique triplet ID
- static unsigned int TripletId2Triplet(unsigned int id);
-
-
- /// Sets L1Algo parameters object
- /// \param other - reference to the Parameters object
- void SetParameters(const Parameters& other) { fParameters = other; }
- // TODO: remove it (S.Zharko)
-
- /// Gets a pointer to the L1Algo parameters object
- const Parameters* GetParameters() const { return &fParameters; }
-
- /// Receives input data
- void ReceiveInputData(InputData&& inputData);
-
- /// Receives tracking parameters
- void ReceiveParameters(Parameters&& parameters);
-
- /// Gets pointer to input data object for external access
- const InputData& GetInputData() const { return fInputData; }
-
- int PackIndex(const int& a, const int& b, const int& c);
-
- int UnPackIndex(const int& i, int& a, int& b, int& c);
-
- /// \brief Sets a default particle mass for the track fit
- /// \param mass Default particle mass
- ///
- /// The function is used during the reconstruction in order to estimate the multiple scattering and energy loss
- void SetDefaultParticleMass(float mass) { fDefaultMass = mass; }
-
- /// \brief Gets default particle mass
- /// \return particle mass
- float GetDefaultParticleMass() const { return fDefaultMass; }
-
- /// \brief Gets default particle mass squared
- /// \return particle mass squared
- float GetDefaultParticleMass2() const { return fDefaultMass * fDefaultMass; }
-
-
- /*********************************************************************************************/ /**
- * ------ FUNCTIONAL PART ------
- ************************************************************************************************/
-
- /// ----- Subroutines used by L1Algo::CaTrackFinder() ------
-
- void ReadWindowData();
-
- void CAFindTrack(int ista, ca::Branch& best_tr, const ca::Triplet* curr_trip, ca::Branch& curr_tr,
- unsigned char min_best_l, ca::Branch* new_tr);
-
-
- /// \brief Fits track
- /// \param t - track with hits
- /// \param T - track parameters
- /// \param dir - false - forward, true - backward
- /// \param qp0 - momentum for extrapolation
- /// \param initParams - should be params ititialized. 1 - yes.
- void BranchFitterFast(const ca::Branch& t, TrackParamV& T, const bool dir, const fvec qp0 = 0.,
- const bool initParams = true);
-
- /// \brief Fits track. more precise than FitterFast
- /// \param t - track with hits
- /// \param T - track parameters
- /// \param dir - false - forward, true - backward
- /// \param qp0 - momentum for extrapolation
- /// \param initParams - should be params ititialized. 1 - yes.
- void BranchFitter(const ca::Branch& t, TrackParamV& T, const bool dir, const fvec qp0 = 0.,
- const bool initParams = true);
-
- /// \brief Finds additional hits for already found track
- /// \param t - track with hits
- /// \param T - track params
- /// \param dir - 0 - forward, 1 - backward
- /// \param qp0 - momentum for extrapolation
- void FindMoreHits(ca::Branch& t, TrackParamV& T, const bool dir, const fvec qp0 = 0.0);
-
- /// \brief Finds additional hits for existing track
- /// \return chi2
- fscal BranchExtender(ca::Branch& t);
-
-#ifdef DRAW
- L1AlgoDraw* draw {nullptr};
-#endif
-
-
- void Init(const TrackingMode mode);
-
- void Finish();
-
- void PrintHits();
-
- /// The main procedure - find tracks
- void CaTrackFinder();
-
- /// find tracks in a sub-timeslice
- void CaTrackFinderSlice();
-
- /// Track fitting procedures
-
- void L1KFTrackFitter(); // version from SIMD-KF benchmark
-
-
- float GetMaxInvMom() const { return fMaxInvMom[0]; }
-
- /// \brief Gets reference to the monitor
- const TrackingMonitor& GetMonitor() const { return fMonitor; }
-
-public:
- /// Gets number of stations before the pipe (MVD stations in CBM)
- int GetNstationsBeforePipe() const { return fNstationsBeforePipe; }
-
- /// Gets number of stations situated in field region (MVD + STS in CBM)
- int GetNfieldStations() const { return fNfieldStations; }
-
- /// Get mc track ID for a hit (debug tool)
- int GetMcTrackIdForCaHit(int iHit) const;
-
- /// Get mc track ID for a hit (debug tool)
- int GetMcTrackIdForWindowHit(int iGridHit) const;
-
- const CbmL1MCTrack* GetMcTrackForWindowHit(int iHit) const;
-
-private:
- bool checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2) const;
-
- int fNstationsBeforePipe {0}; ///< number of stations before pipe (MVD stations in CBM)
- int fNfieldStations {0}; ///< number of stations in the field region
- //alignas(16) L1StationsArray_t fStations {}; ///< array of ca::Station objects
- //alignas(16) L1MaterialArray_t fRadThick {}; ///< material for each station
- float fDefaultMass {constants::phys::MuonMass}; ///< mass of the propagated particle [GeV/c2]
-
-
- // ***************************
- // ** Member variables list **
- // ***************************
-
- Parameters fParameters; ///< Object of L1Algo parameters class
- InputData fInputData; ///< Tracking input data
-
- Vector<unsigned char> fvHitKeyFlags {
- "L1Algo::fvHitKeyFlags"}; ///< List of key flags: has been this hit or cluster already used
-
- ca::TrackingMonitor fMonitor {"CA Algo"}; ///< Tracking monitor
-
-public:
- Vector<L1HitTimeInfo> fHitTimeInfo;
-
- ca::Grid vGrid[constants::size::MaxNstations]; ///<
-
- /// ----- Data used during track finding -----
- ///
-
- ///\brief hits of the current time window
- /// it is a portion of fInputData to process in the current time window
- /// hit.Id is replaced by the hit index in fInputData
- Vector<ca::Hit> fWindowHits {"L1Algo::fWindowHits"};
-
- Vector<unsigned char> fIsWindowHitSuppressed {"L1Algo::fIsWindowHitSuppressed"};
-
- ///\brief first index of the station hits in the time window
- ca::HitIndex_t fStationHitsStartIndex[constants::size::MaxNstations + 1] {0};
-
- ///\brief number of station hits in the time window
- ca::HitIndex_t fStationNhits[constants::size::MaxNstations + 1] {0};
-
- /// ----------
-
- double fCaRecoTime {0.}; // time of the track finder + fitter
-
- Vector<Track> fRecoTracks {"L1Algo::fRecoTracks"}; ///< reconstructed tracks
- Vector<ca::HitIndex_t> fRecoHits {"L1Algo::fRecoHits"}; ///< packed hits of reconstructed tracks
-
- Vector<Track> fSliceRecoTracks {"L1Algo::fSliceRecoTracks"}; ///< reconstructed tracks in sub-timeslice
- Vector<ca::HitIndex_t> fSliceRecoHits {"L1Algo::fSliceRecoHits"}; ///< packed hits of reconstructed tracks
-
- /// Created triplets vs station index
- // FIXME: Better std::array<Vector<Triplet>, constants::size::MaxNstations>
- Vector<ca::Triplet> fTriplets[constants::size::MaxNstations] {{"L1Algo::fTriplets"}};
-
- /// Track candidates created out of adjacent triplets before the final track selection.
- /// The candidates may share any amount of hits.
- // FIXME: Better std::array<Vector<Branch>, constants::size::MaxNstations>
- Vector<ca::Branch> fTrackCandidates {"L1Algo::fTrackCandidates"};
-
- ///< indices of the sub-slice hits
- // FIXME: Better std::array<Vector<HitIndex_t>, constants::size::MaxNstations>
- Vector<ca::HitIndex_t> fSliceHitIds[constants::size::MaxNstations] {"L1Algo::fSliceHitIds"};
-
- Vector<int> fStripToTrack {"L1Algo::fStripToTrack"}; // strip to track pointers
-
- TrackingMode fTrackingMode {kSts};
-
- fvec EventTime {0.f};
- fvec Err {0.f};
-
- /// --- data used during finding iterations
- unsigned int fCurrentIterationIndex {0}; ///< index of the corrent iteration, needed for debug output
- const Iteration* fpCurrentIteration = nullptr; ///< pointer to the current CA track finder iteration
-
- Vector<int> fHitFirstTriplet {"L1Algo::fHitFirstTriplet"}; /// link hit -> first triplet { hit, *, *}
- Vector<int> fHitNtriplets {"L1Algo::fHitNtriplets"}; /// link hit ->n triplets { hit, *, *}
-
-
-private:
- L1CloneMerger fCloneMerger {*this}; ///< Object of L1Algo clones merger algorithm
-
- /// ================================= DATA PART =================================
-
- /// ----- Different parameters of CATrackFinder -----
-
- Tindex fFirstCAstation {0}; //first station used in CA
-
- // fNFindIterations - set number of interation for trackfinding
- // itetation of finding:
-
- int fNFindIterations = -1; // TODO investigate kAllPrimJumpIter & kAllSecJumpIter
-
- float fTrackChi2Cut {10.f};
- float fTripletFinalChi2Cut {10.f};
- float fTripletChi2Cut {5.f}; // cut for selecting triplets before collecting tracks.per one DoF
- float fDoubletChi2Cut {5.f};
- float fTimeCut1 {0.f}; // TODO: please, specify "1" and "2" (S.Zharko)
- float fTimeCut2 {0.f};
-
- /// correction in order to take into account overlaping and iff z. if sort by y then it is max diff between same station's modules (~0.4cm)
- fvec fMaxDZ {constants::Undef<fvec>};
-
- /// parameters which are different for different iterations. Set in the begin of CAL1TrackFinder
-
- float fPickGather {constants::Undef<fvec>}; ///< same for attaching additional hits to track
- float fTripletLinkChi2 {constants::Undef<fvec>}; ///< (dp2/dp_error2 < fTripletLinkChi2) => triplets are neighbours
- fvec fMaxInvMom {constants::Undef<fvec>}; ///< max considered q/p for tracks
- fvec fMaxSlopePV {constants::Undef<fvec>}; ///< max slope (tx\ty) in prim vertex
- float fMaxSlope {constants::Undef<fvec>}; ///< max slope (tx\ty) in 3d hit position of a triplet
-
- fvec fTargX {constants::Undef<fvec>}; ///< target position x coordinate for the current iteration (modifiable)
- fvec fTargY {constants::Undef<fvec>}; ///< target position y coordinate for the current iteration (modifiable)
- fvec fTargZ {constants::Undef<fvec>}; ///< target position z coordinate for the current iteration (modifiable)
-
- bool fIsTargetField {false}; ///< is the magnetic field present at the target
-
- ca::FieldValue fTargB _fvecalignment {}; // field in the target point (modifiable, do not touch!!)
- ca::MeasurementXy<fvec> TargetMeasurement _fvecalignment {}; // target constraint [cm]
-
- int fGhostSuppression {1}; // NOTE: Should be equal to 0 in TRACKS_FROM_TRIPLETS mode!
-
- /// ----- Debug features -----
-#ifdef PULLS
- L1AlgoPulls* fL1Pulls;
-#endif
-#ifdef TRIP_PERFORMANCE
- L1AlgoEfficiencyPerformance<3>* fL1Eff_triplets;
- L1AlgoEfficiencyPerformance<3>* fL1Eff_triplets2;
-#endif
-#ifdef DOUB_PERFORMANCE
- L1AlgoEfficiencyPerformance<2>* fL1Eff_doublets;
-#endif
-} _fvecalignment;
-
-
-// ********************************************
-// ** Inline member functions implementation **
-// ********************************************
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-[[gnu::always_inline]] inline unsigned int L1Algo::PackTripletId(unsigned int iStation, unsigned int iTriplet)
-{
-#ifndef FAST_CODE
- assert(iStation < constants::size::MaxNstations);
- assert(iTriplet < constants::size::MaxNtriplets);
-#endif
- constexpr unsigned int kMoveStation = constants::size::TripletBits;
- return (iStation << kMoveStation) + iTriplet;
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-[[gnu::always_inline]] inline unsigned int L1Algo::TripletId2Station(unsigned int id)
-{
- constexpr unsigned int kMoveStation = constants::size::TripletBits;
- return id >> kMoveStation;
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-[[gnu::always_inline]] inline unsigned int L1Algo::TripletId2Triplet(unsigned int id)
-{
- constexpr unsigned int kTripletMask = (1u << constants::size::TripletBits) - 1u;
- return id & kTripletMask;
-}
-
-#endif
diff --git a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
deleted file mode 100644
index c4b39e28e2..0000000000
--- a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
+++ /dev/null
@@ -1,871 +0,0 @@
-/* Copyright (C) 2009-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Valentina Akishina, Ivan Kisel, Sergey Gorbunov [committer], Igor Kulakov, Maksym Zyzak */
-
-/*
- *=====================================================
- *
- * CBM Level 1 4D Reconstruction
- *
- * Authors: V.Akishina, I.Kisel, S.Gorbunov, I.Kulakov, M.Zyzak
- * Documentation: V.Akishina
- *
- * e-mail : v.akishina@gsi.de
- *
- *=====================================================
- *
- * Finds tracks using the Cellular Automaton algorithm
- *
- */
-
-
-#include "CbmL1.h"
-#include "CbmL1MCTrack.h"
-
-#include "CaBranch.h"
-#include "CaGrid.h"
-#include "CaGridEntry.h"
-#include "CaTrack.h"
-#include "CaTrackParam.h"
-#include "L1Algo.h"
-#include "L1TripletConstructor.h"
-#ifdef DRAW
-#include "utils/L1AlgoDraw.h"
-#endif // DRAW
-#ifdef PULLS
-#include "utils/L1AlgoPulls.h"
-#endif // PULLS
-#ifdef TRIP_PERFORMANCE
-#include "utils/L1AlgoEfficiencyPerformance.h"
-#endif // TRIP_PERFORMANCE
-#ifdef DOUB_PERFORMANCE
-#include "utils/L1AlgoEfficiencyPerformance.h"
-#endif // DOUB_PERFORMANCE
-
-#include <algorithm>
-#include <fstream>
-#include <iostream>
-#include <list>
-#include <map>
-
-#include <stdio.h>
-
-#include "CaToolsDebugger.h"
-
-
-// using namespace std;
-using std::cout;
-using std::endl;
-using Track = cbm::algo::ca::Track;
-
-bool L1Algo::checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2) const
-{
- dchi2 = 1.e20;
-
- if (r.GetMHit() != l.GetRHit()) return false;
- if (r.GetLHit() != l.GetMHit()) return false;
-
- if (r.GetMSta() != l.GetRSta()) return false;
- if (r.GetLSta() != l.GetMSta()) return false;
-
-
- if (r.IsMomentumFitted()) {
- assert(l.IsMomentumFitted());
-
- fscal dqp = l.GetQp() - r.GetQp();
- fscal Cqp = l.GetCqp() + r.GetCqp();
-
- if (!std::isfinite(dqp)) return false;
- if (!std::isfinite(Cqp)) return false;
-
- if (dqp * dqp > fTripletLinkChi2 * Cqp) {
- return false; // bad neighbour // CHECKME why do we need recheck it?? (it really change result)
- }
- dchi2 = dqp * dqp / Cqp;
- }
- else {
- fscal dtx = l.GetTx() - r.GetTx();
- fscal Ctx = l.GetCtx() + r.GetCtx();
-
- fscal dty = l.GetTy() - r.GetTy();
- fscal Cty = l.GetCty() + r.GetCty();
-
- // it shouldn't happen, but happens sometimes
-
- if (!std::isfinite(dtx)) return false;
- if (!std::isfinite(dty)) return false;
- if (!std::isfinite(Ctx)) return false;
- if (!std::isfinite(Cty)) return false;
-
- if (dty * dty > fTripletLinkChi2 * Cty) return false;
- if (dtx * dtx > fTripletLinkChi2 * Ctx) return false;
-
- //dchi2 = 0.5f * (dtx * dtx / Ctx + dty * dty / Cty);
- dchi2 = 0.;
- }
-
- if (!std::isfinite(dchi2)) return false;
-
- return true;
-}
-
-
-// **************************************************************************************************
-// * *
-// * ------ CATrackFinder procedure ------ *
-// * *
-// **************************************************************************************************
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void L1Algo::ReadWindowData()
-{
- int nHits = 0;
- for (int iS = 0; iS < fParameters.GetNstationsActive(); iS++) {
- fStationHitsStartIndex[iS] = nHits;
- fStationNhits[iS] = fSliceHitIds[iS].size();
- nHits += fStationNhits[iS];
- }
- fStationHitsStartIndex[fParameters.GetNstationsActive()] = nHits;
-
- fWindowHits.reset(nHits);
- fIsWindowHitSuppressed.reset(nHits, 0);
-
- for (int iS = 0; iS < fParameters.GetNstationsActive(); iS++) {
- for (ca::HitIndex_t ih = 0; ih < fSliceHitIds[iS].size(); ++ih) {
- ca::Hit h = fInputData.GetHit(fSliceHitIds[iS][ih]);
- h.SetId(fSliceHitIds[iS][ih]);
-
- fWindowHits[fStationHitsStartIndex[iS] + ih] = h;
- }
- }
-
- fHitFirstTriplet.reset(nHits, 0);
- fHitNtriplets.reset(nHits, 0);
-
- // TODO: SG: reduce the array size to the number of keys in subslice
-
- fStripToTrack.reset(fInputData.GetNhitKeys(), -1);
-
- fSliceRecoTracks.clear();
- fSliceRecoTracks.reserve(2 * nHits / fParameters.GetNstationsActive());
-
- fSliceRecoHits.clear();
- fSliceRecoHits.reserve(2 * nHits);
-
- fTrackCandidates.clear();
- fTrackCandidates.reserve(nHits / 10);
- for (unsigned int iS = 0; iS < constants::size::MaxNstations; iS++) {
- int nHitsStation = fSliceHitIds[iS].size();
- fTriplets[iS].clear();
- fTriplets[iS].reserve(2 * nHitsStation);
- }
-}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void L1Algo::CaTrackFinderSlice()
-{
- fNFindIterations = fParameters.GetNcaIterations();
-
- /*******************************/ /**
- * Performance monitors setting
- **********************************/
-
-#ifdef PULLS
- static L1AlgoPulls* l1Pulls_ = new L1AlgoPulls();
- fL1Pulls = l1Pulls_;
- fL1Pulls->Init();
-#endif
-#ifdef TRIP_PERFORMANCE
- static L1AlgoEfficiencyPerformance<3>* l1Eff_triplets_ = new L1AlgoEfficiencyPerformance<3>();
- fL1Eff_triplets = l1Eff_triplets_;
- fL1Eff_triplets->Init();
- static L1AlgoEfficiencyPerformance<3>* l1Eff_triplets2_ = new L1AlgoEfficiencyPerformance<3>();
- fL1Eff_triplets2 = l1Eff_triplets2_;
- fL1Eff_triplets2->Init();
-#endif
-#ifdef DOUB_PERFORMANCE
- static L1AlgoEfficiencyPerformance<2>* l1Eff_doublets_ = new L1AlgoEfficiencyPerformance<2>();
- fL1Eff_doublets = l1Eff_doublets_;
- fL1Eff_doublets->Init();
-#endif
-
-#ifdef DRAW
- if (!draw) draw = new L1AlgoDraw;
- draw->InitL1Draw(this);
-#endif
-
-
-#if defined(COUNTERS)
- static unsigned int stat_N = 0; // number of events
- stat_N++;
-
- static Tindex stat_nStartHits = 0;
- static Tindex stat_nHits[fNFindIterations] = {0};
-
- static Tindex stat_nSinglets[fNFindIterations] = {0};
- // static Tindex stat_nDoublets[fNFindIterations] = {0};
- static Tindex stat_nTriplets[fNFindIterations] = {0};
-
- static Tindex stat_nLevels[constants::size::MaxNstations - 2][fNFindIterations] = {{0}};
- static Tindex stat_nCalls[fNFindIterations] = {0}; // n calls of CAFindTrack
- static Tindex stat_nTrCandidates[fNFindIterations] = {0};
-
- for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
- stat_nStartHits += fSliceHitIds[iS].size();
- }
-#endif
-
- /********************************/ /**
- * CATrackFinder routine setting
- ***********************************/
-
- ReadWindowData();
-
- for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
-
- bool timeUninitialised = 1;
- fscal lasttime = 0;
- fscal starttime = 0;
-
-
- fscal gridMinX = -0.1;
- fscal gridMaxX = 0.1;
- fscal gridMinY = -0.1;
- fscal gridMaxY = 0.1;
-
- for (ca::HitIndex_t ih = 0; ih < fSliceHitIds[iS].size(); ++ih) {
- const ca::Hit& h = fInputData.GetHit(fSliceHitIds[iS][ih]);
-
- if (h.X() < gridMinX) { gridMinX = h.X(); }
- if (h.X() > gridMaxX) { gridMaxX = h.X(); }
- if (h.Y() < gridMinY) { gridMinY = h.Y(); }
- if (h.Y() > gridMaxY) { gridMaxY = h.Y(); }
-
- const fscal time = h.T();
- assert(std::isfinite(time));
- if (timeUninitialised || lasttime < time) { lasttime = time; }
- if (timeUninitialised || starttime > time) { starttime = time; }
- timeUninitialised = false;
- }
-
- // TODO: changing the grid also changes the result. Investigate why it happens.
- // TODO: find the optimal grid size
-
- int nSliceHits = fSliceHitIds[iS].size();
- fscal sizeY = gridMaxY - gridMinY;
- fscal sizeX = gridMaxX - gridMinX;
- int nBins2D = 1 + nSliceHits;
- fscal yStep = sizeY / sqrt(nBins2D);
- fscal xStep = sizeX / sqrt(nBins2D);
-
- fscal scale = fParameters.GetStation(iS).GetZScal() - fParameters.GetTargetPositionZ()[0];
-
- if (yStep > 0.3 * scale) yStep = 0.3 * scale;
- if (yStep < 0.01 * scale) yStep = 0.01 * scale;
- if (xStep > 0.3 * scale) xStep = 0.3 * scale;
- if (xStep < 0.01 * scale) xStep = 0.01 * scale;
-
- vGrid[iS].BuildBins(gridMinX, gridMaxX, gridMinY, gridMaxY, xStep, yStep);
- vGrid[iS].StoreHits(fWindowHits, fStationHitsStartIndex[iS], fStationNhits[iS], fvHitKeyFlags);
- }
-
-
- // ---- Loop over Track Finder iterations ----------------------------------------------------------------//
-
-
- assert(fNFindIterations == (int) fParameters.GetCAIterations().size());
-
- for (fCurrentIterationIndex = 0; fCurrentIterationIndex < fParameters.GetCAIterations().size();
- ++fCurrentIterationIndex) {
-
- // current iteration of the CA tracker
- const auto& caIteration = fParameters.GetCAIterations()[fCurrentIterationIndex];
- fpCurrentIteration = &caIteration;
-
- if (fCurrentIterationIndex > 0) {
- for (int ista = 0; ista < fParameters.GetNstationsActive(); ++ista) {
- vGrid[ista].RemoveUsedHits(fWindowHits, fvHitKeyFlags);
- }
- }
-
- fIsWindowHitSuppressed.reset(fWindowHits.size(), 0);
-
- for (int j = 0; j < fParameters.GetNstationsActive(); j++) {
- fTriplets[j].clear();
- }
-
-#ifdef COUNTERS
- Tindex nSinglets = 0;
-#endif
-
- fHitFirstTriplet.reset(fWindowHits.size(), 0);
- fHitNtriplets.reset(fWindowHits.size(), 0);
-
- {
- // #pragma omp task
- {
- // --- SET PARAMETERS FOR THE ITERATION ---
-
- fFirstCAstation = caIteration.GetFirstStationIndex();
- fTrackChi2Cut = caIteration.GetTrackChi2Cut();
- fDoubletChi2Cut = caIteration.GetDoubletChi2Cut(); //11.3449 * 2.f / 3.f; // prob = 0.1
- fTripletChi2Cut = caIteration.GetTripletChi2Cut(); //21.1075; // prob = 0.01%
- fTripletFinalChi2Cut = caIteration.GetTripletFinalChi2Cut();
- fPickGather = caIteration.GetPickGather(); //3.0;
- fTripletLinkChi2 = caIteration.GetTripletLinkChi2(); //5.0;
- fMaxInvMom = caIteration.GetMaxQp(); //1.0 / 0.5; // max considered q/p
- fMaxSlopePV = caIteration.GetMaxSlopePV(); //1.1;
- fMaxSlope = caIteration.GetMaxSlope(); //2.748; // corresponds to 70 grad
-
- // define the target
- fTargX = fParameters.GetTargetPositionX();
- fTargY = fParameters.GetTargetPositionY();
- fTargZ = fParameters.GetTargetPositionZ();
-
- fscal SigmaTargetX = caIteration.GetTargetPosSigmaX();
- fscal SigmaTargetY = caIteration.GetTargetPosSigmaY(); // target constraint [cm]
-
- // Select magnetic field. For primary tracks - fVtxFieldValue, for secondary tracks - st.fieldSlice
- if (caIteration.GetPrimaryFlag()) { fTargB = fParameters.GetVertexFieldValue(); }
- else {
- fTargB = fParameters.GetStation(0).fieldSlice.GetFieldValue(0, 0);
- } // NOTE: calculates field fTargB in the center of 0th station
-
- fIsTargetField = (fabs(fTargB.y[0]) > 0.001);
-
- TargetMeasurement.SetX(fTargX);
- TargetMeasurement.SetY(fTargY);
- TargetMeasurement.SetDx2(SigmaTargetX * SigmaTargetX);
- TargetMeasurement.SetDxy(0);
- TargetMeasurement.SetDy2(SigmaTargetY * SigmaTargetY);
- TargetMeasurement.SetNdfX(1);
- TargetMeasurement.SetNdfY(1);
-
- /// Set correction in order to take into account overlaping and iff z.
- /// The reason is that low momentum tracks are too curved and goes not from target direction. That's why sort by hit_y/hit_z is not work idealy
- /// If sort by y then it is max diff between same station's modules (~0.4cm)
- fMaxDZ = caIteration.GetMaxDZ(); //0;
- }
- }
-
-
- /// stage for triplets creation
- fMonitor.StartTimer(ETimer::TripletConstruction);
- L1TripletConstructor constructor1(this);
- L1TripletConstructor constructor2(this);
- L1TripletConstructor constructor3(this);
-
- for (int istal = fParameters.GetNstationsActive() - 2; istal >= fFirstCAstation;
- istal--) { // start with downstream chambers
- Tindex nGridEntriesL = vGrid[istal].GetEntries().size();
-
- for (Tindex iel = 0; iel < nGridEntriesL; ++iel) {
- ca::HitIndex_t ihitl = vGrid[istal].GetEntries()[iel].GetObjectId();
- constructor1.CreateTripletsForHit(istal, istal + 1, istal + 2, ihitl);
-
- if (fpCurrentIteration->GetJumpedFlag()) {
- constructor2.CreateTripletsForHit(istal, istal + 1, istal + 3, ihitl);
- constructor3.CreateTripletsForHit(istal, istal + 2, istal + 3, ihitl);
- }
-
- auto oldSize = fTriplets[istal].size();
-
- fTriplets[istal].insert(fTriplets[istal].end(), constructor1.GetTriplets().begin(),
- constructor1.GetTriplets().end());
- fTriplets[istal].insert(fTriplets[istal].end(), constructor2.GetTriplets().begin(),
- constructor2.GetTriplets().end());
- fTriplets[istal].insert(fTriplets[istal].end(), constructor3.GetTriplets().begin(),
- constructor3.GetTriplets().end());
-
- fHitFirstTriplet[ihitl] = PackTripletId(istal, oldSize);
- fHitNtriplets[ihitl] = fTriplets[istal].size() - oldSize;
- }
- fMonitor.IncrementCounter(ECounter::Triplet, fTriplets[istal].size());
- } // istal
- fMonitor.StopTimer(ETimer::TripletConstruction);
-
- // search for neighbouring triplets
- fMonitor.StartTimer(ETimer::NeighboringTripletSearch);
- for (int istal = fParameters.GetNstationsActive() - 2; istal >= fFirstCAstation;
- istal--) { // start with downstream chambers
-
- for (unsigned int it = 0; it < fTriplets[istal].size(); ++it) {
-
- ca::Triplet& tr = fTriplets[istal][it];
- tr.SetLevel(0);
- tr.SetFNeighbour(0);
- tr.SetNNeighbours(0);
-
- if (istal > (fParameters.GetNstationsActive() - 4)) break;
-
- unsigned int nNeighbours = fHitNtriplets[tr.GetMHit()];
-
- unsigned int neighLocation = fHitFirstTriplet[tr.GetMHit()];
- unsigned int neighStation = TripletId2Station(neighLocation);
- unsigned int neighTriplet = TripletId2Triplet(neighLocation);
-
- if (nNeighbours > 0) { assert((int) neighStation == istal + 1 || (int) neighStation == istal + 2); }
-
- unsigned char level = 0;
-
- for (unsigned int iN = 0; iN < nNeighbours; ++iN, ++neighTriplet, ++neighLocation) {
-
- ca::Triplet& neighbour = fTriplets[neighStation][neighTriplet];
-
- fscal dchi2 = 0.;
- if (!checkTripletMatch(tr, neighbour, dchi2)) continue;
-
- if (tr.GetFNeighbour() == 0) { tr.SetFNeighbour(neighLocation); }
-
- tr.SetNNeighbours(neighLocation - tr.GetFNeighbour() + 1);
-
- const unsigned char neighbourLevel = neighbour.GetLevel();
-
- if (level <= neighbourLevel) { level = neighbourLevel + 1; }
- }
- tr.SetLevel(level);
- } // neighbour search
- } // istal
- fMonitor.StopTimer(ETimer::NeighboringTripletSearch);
-
- ///====================================================================
- ///= =
- ///= Collect track candidates. CREATE TRACKS =
- ///= =
- ///====================================================================
-
- // min level to start triplet. So min track length = min_level+3.
- int min_level = std::min(caIteration.GetMinNhits(), caIteration.GetMinNhitsStation0()) - 3;
-
- // TODO: SZh 04.10.2022: Why this fatal error is called?
- // NOTE: This line was wrapped into TRACKS_FROM_TRIPLETS ifdef
- //LOG(FATAL) << "L1CATrackFinder: min_level undefined in " << __FILE__ << " : " << __LINE__;
- if (fpCurrentIteration->GetTrackFromTripletsFlag()) { min_level = 0; }
-
- ca::Branch curr_tr;
- ca::Branch new_tr[constants::size::MaxNstations];
- ca::Branch best_tr;
-
- int ndf = 1;
-
-
- // collect consequtive: the longest tracks, shorter, more shorter and so on
- for (int firstTripletLevel = fParameters.GetNstationsActive() - 3; firstTripletLevel >= min_level;
- firstTripletLevel--) {
- // choose length in triplets number - firstTripletLevel - the maximum possible triplet level among all triplets in the searched candidate
-
- // how many levels to check
- int nlevel = (fParameters.GetNstationsActive() - 2) - firstTripletLevel + 1;
-
- const unsigned char min_best_l =
- (firstTripletLevel > min_level) ? firstTripletLevel + 2 : min_level + 3; // loose maximum
-
-
- fTrackCandidates.clear();
-
- fStripToTrack.reset(fInputData.GetNhitKeys(), -1);
-
- //== Loop over triplets with the required level, find and store track candidates
-
- for (int istaF = fFirstCAstation; istaF <= fParameters.GetNstationsActive() - 3 - firstTripletLevel; ++istaF) {
-
- if (--nlevel == 0) break; //TODO: SG: this is not needed
-
- for (Tindex itrip = 0; itrip < (Tindex) fTriplets[istaF].size(); ++itrip) {
-
- ca::Triplet& first_trip = (fTriplets[istaF][itrip]);
- if (fvHitKeyFlags[fWindowHits[first_trip.GetLHit()].FrontKey()]
- || fvHitKeyFlags[fWindowHits[first_trip.GetLHit()].BackKey()]) {
- continue;
- }
- // ghost suppression !!!
- //
-
- if (!fpCurrentIteration->GetTrackFromTripletsFlag()) { // ghost suppression !!!
- int nHits = 3 + first_trip.GetLevel();
- if (fWindowHits[first_trip.GetLHit()].Station() == 0) {
- if (nHits < fpCurrentIteration->GetMinNhitsStation0()) { continue; }
- }
- else {
- if (nHits < fpCurrentIteration->GetMinNhits()) { continue; }
- }
- }
-
- // Collect triplets, which can start a track with length equal to firstTipletLevel + 3. This cut suppresses
- // ghost tracks, but does not affect the efficiency
- if (first_trip.GetLevel() < firstTripletLevel) { continue; }
-
- curr_tr.SetChi2(0.);
- curr_tr.SetStation(0);
- curr_tr.ResetHits();
- curr_tr.AddHit(fWindowHits[first_trip.GetLHit()].Id());
- curr_tr.SetChi2(first_trip.GetChi2());
-
- best_tr = curr_tr;
-
- CAFindTrack(istaF, best_tr, &first_trip, curr_tr, min_best_l,
- new_tr); /// reqursive func to build a tree of possible track-candidates and choose the best
-
- if (best_tr.NofHits() < firstTripletLevel + 2) continue; // loose maximum one hit
-
- if (best_tr.NofHits() < min_level + 3) continue; // should find all hits for min_level
-
- ndf = best_tr.NofHits() * 2 - 5;
-
- // TODO: automatize the NDF calculation
-
- if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) { ndf = best_tr.NofHits() * 2 - 4; }
-
- best_tr.SetChi2(best_tr.Chi2() / ndf);
- if (fGhostSuppression) {
- if (3 == best_tr.NofHits()) {
- if (!fpCurrentIteration->GetPrimaryFlag() && (istaF != 0)) continue; // too /*short*/ non-MAPS track
- if (fpCurrentIteration->GetPrimaryFlag() && (best_tr.Chi2() > 5.0)) continue;
- }
- }
- fTrackCandidates.push_back(best_tr);
- ca::Branch& tr = fTrackCandidates.back();
- tr.SetStation(istaF);
- tr.SetId(fTrackCandidates.size() - 1);
- tr.SetAlive(true);
- if (0) { // debug
- cout << "track " << fTrackCandidates.size() - 1 << " found, L = " << best_tr.NofHits()
- << " chi2= " << best_tr.Chi2() << endl;
- cout << " hits: ";
- for (auto hitId : tr.Hits()) {
- cout << GetMcTrackIdForCaHit(hitId) << " ";
- }
- cout << endl;
- }
- } // itrip
- } // istaF
-
- for (ca::Branch& tr : fTrackCandidates) {
- tr.SetAlive(false);
- }
-
- bool repeatCompetition = true;
- int iComp = 0;
- for (iComp = 0; (iComp < 100) && repeatCompetition; ++iComp) {
-
- // == Loop over track candidates and mark their strips
-
- repeatCompetition = false;
-
- for (ca::Branch& tr : fTrackCandidates) {
- if (tr.IsAlive()) { continue; }
- for (auto& hitId : tr.Hits()) {
- const ca::Hit& h = fInputData.GetHit(hitId);
- bool isAlive = true;
- { // front strip
- auto& stripF = (fStripToTrack)[h.FrontKey()];
- if ((stripF >= 0) && (stripF != tr.Id())) { // strip is used by other candidate
- const auto& other = fTrackCandidates[stripF];
- if (!other.IsAlive() && tr.IsBetterThan(other)) { stripF = tr.Id(); }
- else {
- isAlive = false;
- }
- }
- else {
- stripF = tr.Id();
- }
- if (!isAlive) { break; }
- }
-
- { // back strip
- auto& stripB = (fStripToTrack)[h.BackKey()];
- if ((stripB >= 0) && (stripB != tr.Id())) { // strip is used by other candidate
- const auto& other = fTrackCandidates[stripB];
- if (!other.IsAlive() && tr.IsBetterThan(other)) { stripB = tr.Id(); }
- else {
- isAlive = false;
- }
- }
- else {
- stripB = tr.Id();
- }
- if (!isAlive) { break; }
- }
- } // loop over hits
- } // itrack
-
- // == Check if some suppressed candidates are still alive
-
- for (ca::Branch& tr : fTrackCandidates) {
- if (tr.IsAlive()) { continue; }
-
- tr.SetAlive(true);
- for (int iHit = 0; tr.IsAlive() && (iHit < (int) tr.Hits().size()); ++iHit) {
- const ca::Hit& h = fInputData.GetHit(tr.Hits()[iHit]);
- tr.SetAlive(tr.IsAlive() && ((fStripToTrack)[h.FrontKey()] == tr.Id())
- && ((fStripToTrack)[h.BackKey()] == tr.Id()));
- }
-
- if (!tr.IsAlive()) { // release strips
- for (auto hitId : tr.Hits()) {
- const ca::Hit& h = fInputData.GetHit(hitId);
- if (fStripToTrack[h.FrontKey()] == tr.Id()) { fStripToTrack[h.FrontKey()] = -1; }
- if (fStripToTrack[h.BackKey()] == tr.Id()) { fStripToTrack[h.BackKey()] = -1; }
- }
- }
- else {
- repeatCompetition = true;
- }
- } // itrack
-
- } // competitions
-
- // cout << " N Competitions: " << iComp << endl;
-
- // ==
-
- for (Tindex iCandidate = 0; iCandidate < (Tindex) fTrackCandidates.size(); ++iCandidate) {
- ca::Branch& tr = fTrackCandidates[iCandidate];
- if (!tr.IsAlive()) continue;
- if (kMcbm == fTrackingMode) {
- if (tr.NofHits() <= 3) { continue; }
- }
- else if (kGlobal == fTrackingMode) {
- if (tr.NofHits() < 3) { continue; }
- }
- else {
- if (tr.NofHits() < 3) { continue; }
- }
- if (fpCurrentIteration->GetExtendTracksFlag()) {
- if (tr.NofHits() < fParameters.GetNstationsActive()) { BranchExtender(tr); }
- }
-
- for (auto iHit : tr.Hits()) {
- const ca::Hit& hit = fInputData.GetHit(iHit);
-
- /// used strips are marked
-
- fvHitKeyFlags[hit.FrontKey()] = 1;
- fvHitKeyFlags[hit.BackKey()] = 1;
-
- fSliceRecoHits.push_back(iHit);
- }
- Track t;
- t.fNofHits = tr.NofHits();
- fSliceRecoTracks.push_back(t);
- if (0) { // SG debug
- cout << "store track " << iCandidate << " chi2= " << tr.Chi2() << endl;
- cout << " hits: ";
- for (unsigned int ih = 0; ih < tr.Hits().size(); ih++) {
- cout << GetMcTrackIdForCaHit(tr.Hits()[ih]) << " ";
- }
- cout << '\n';
- }
-
- } // tracks
-
- } // firstTripletLevel
-
-
- // suppress strips of suppressed hits
- for (unsigned int ih = 0; ih < fWindowHits.size(); ih++) {
- if (fIsWindowHitSuppressed[ih]) {
- const ca::Hit& hit = fWindowHits[ih];
- fvHitKeyFlags[hit.FrontKey()] = 1;
- fvHitKeyFlags[hit.BackKey()] = 1;
- }
- }
-
- {
- // #ifdef DRAW
- // draw->ClearVeiw();
- // // draw->DrawInfo();
- // draw->DrawRestHits(fGridHitStartIndex, fGridHitStopIndex, RealIHit);
- // draw->DrawRecoTracks();
- // draw->SaveCanvas("Reco_"+iCaIteration+"_");
- // draw->DrawAsk();
- // #endif
-
- // #ifdef PULLS
- // fL1Pulls->Build(1);
- // #endif
-#ifdef COUNTERS
- // stat_nHits[iCaIteration] += (fGridHits*)->Size();
- cout << "iter = " << iCaIteration << endl;
- cout << " NSinglets = " << stat_nSinglets[iCaIteration] / stat_N << endl;
- // cout << " NDoublets = " << stat_nDoublets[iCaIteration]/stat_N << endl;
- cout << " NTriplets = " << stat_nTriplets[iCaIteration] / stat_N << endl;
- cout << " NfGridHit = " << stat_nHits[iCaIteration] / stat_N << endl;
-#endif // COUNTERS
- }
-
- } // ---- Loop over Track Finder iterations: END -----------------------------------------------------------//
-
- // Fit tracks
- this->L1KFTrackFitter();
-
- // Merge clones
- fCloneMerger.Exec(fSliceRecoTracks, fSliceRecoHits);
-
- //==================================
-
-#ifdef DRAW
- draw->ClearVeiw();
- // draw->DrawInputHits();
- // draw->DrawInfo();
- draw->DrawRecoTracks();
-
- draw->SaveCanvas("Reco_");
- draw->DrawAsk();
-#endif
-#ifdef PULLS
- static int iEvee = 0;
- iEvee++;
- if (iEvee % 1 == 0) fL1Pulls->Build(1);
-#endif
-#ifdef DOUB_PERFORMANCE
- fL1Eff_doublets->CalculateEff();
- fL1Eff_doublets->Print("Doublets performance.", 1);
-#endif
-#ifdef TRIP_PERFORMANCE
- fL1Eff_triplets->CalculateEff();
- fL1Eff_triplets->Print("Triplet performance", 1);
-// fL1Eff_triplets2->CalculateEff();
-// fL1Eff_triplets2->Print("Triplet performance. After chi2 cut");
-#endif
-}
-
-
-// ---------------------------------------------------------------------------------------------------------------------
-// ****************************************************************
-// * *
-// * The routine performs recursive search for tracks *
-// * *
-// * I. Kisel 06.03.05 *
-// * I.Kulakov 2012 *
-// * *
-// ****************************************************************
-void L1Algo::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Triplet* curr_trip, ca::Branch& curr_tr,
- unsigned char min_best_l, ca::Branch* new_tr)
-/// recursive search for tracks
-/// input: @ista - station index, @&best_tr - best track for the privious call
-/// output: @&NCalls - number of function calls
-{
-
- if (curr_trip->GetLevel() == 0) // the end of the track -> check and store
- {
- // -- finish with current track
- // add rest of hits
- const ca::HitIndex_t& ihitm = curr_trip->GetMHit();
- const ca::HitIndex_t& ihitr = curr_trip->GetRHit();
-
-
- //if (!GetFUsed((*fStripFlag)[fGridHits[ihitm].f] | (*fStripFlag)[fGridHits[ihitm].b])) {
- //L1_SHOW(fInputData.GetNhitKeys());
- //L1_SHOW(fvHitKeyFlags.size());
- //L1_SHOW(fGridHits[ihitm].f);
- //L1_SHOW(fGridHits[ihitm].b);
- if (!(fvHitKeyFlags[fWindowHits[ihitm].FrontKey()] || fvHitKeyFlags[fWindowHits[ihitm].BackKey()])) {
-
- // curr_tr.Hits.push_back(fGridHitIds[ihitm]);
-
- curr_tr.AddHit(fWindowHits[ihitm].Id());
- }
-
- //if (!GetFUsed((*fStripFlag)[fGridHits[ihitr].f] | (*fStripFlag)[fGridHits[ihitr].b])) {
- //L1_SHOW(fInputData.GetNhitKeys());
- //L1_SHOW(fvHitKeyFlags.size());
- //L1_SHOW(fGridHits[ihitr].f);
- //L1_SHOW(fGridHits[ihitr].b);
- if (!(fvHitKeyFlags[fWindowHits[ihitr].FrontKey()] || fvHitKeyFlags[fWindowHits[ihitr].BackKey()])) {
-
- //curr_tr.Hits.push_back(fGridHitIds[ihitr]);
- curr_tr.AddHit(fWindowHits[ihitr].Id());
- }
-
- //if( curr_tr.NofHits() < min_best_l - 1 ) return; // suppouse that only one hit can be added by extender
- // TODO: SZh 21.08.2023: Replace hard-coded value with a parameter
- int ndf = curr_tr.NofHits() * 2 - 5;
-
- if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) { ndf = curr_tr.NofHits() * 2 - 4; }
-
- if (curr_tr.Chi2() > fTrackChi2Cut * ndf) { return; }
-
-
- // -- select the best
- if ((curr_tr.NofHits() > best_tr.NofHits())
- || ((curr_tr.NofHits() == best_tr.NofHits()) && (curr_tr.Chi2() < best_tr.Chi2()))) {
- best_tr = curr_tr;
- }
-
- return;
- }
- else //MEANS level ! = 0
- {
- int N_neighbour = (curr_trip->GetNNeighbours());
-
- for (Tindex in = 0; in < N_neighbour; in++) {
-
- unsigned int ID = curr_trip->GetFNeighbour() + in;
-
- unsigned int Station = TripletId2Station(ID);
- unsigned int Triplet = TripletId2Triplet(ID);
-
- const ca::Triplet& new_trip = fTriplets[Station][Triplet];
-
- fscal dchi2 = 0.;
- if (!checkTripletMatch(*curr_trip, new_trip, dchi2)) continue;
-
- if (fvHitKeyFlags[fWindowHits[new_trip.GetLHit()].FrontKey()]
- || fvHitKeyFlags[fWindowHits[new_trip.GetLHit()].BackKey()]) { //hits are used
- // no used hits allowed -> compare and store track
- if ((curr_tr.NofHits() > best_tr.NofHits())
- || ((curr_tr.NofHits() == best_tr.NofHits()) && (curr_tr.Chi2() < best_tr.Chi2()))) {
- best_tr = curr_tr;
- }
- }
- else { // hit is not used: add the left hit from the new triplet to the current track
-
- unsigned char new_L = curr_tr.NofHits() + 1;
- fscal new_chi2 = curr_tr.Chi2() + dchi2;
-
-
- if (0) { //SGtrd2d debug!!
- int mc01 = GetMcTrackIdForWindowHit(curr_trip->GetLHit());
- int mc02 = GetMcTrackIdForWindowHit(curr_trip->GetMHit());
- int mc03 = GetMcTrackIdForWindowHit(curr_trip->GetRHit());
- int mc11 = GetMcTrackIdForWindowHit(new_trip.GetLHit());
- int mc12 = GetMcTrackIdForWindowHit(new_trip.GetMHit());
- int mc13 = GetMcTrackIdForWindowHit(new_trip.GetRHit());
- if ((mc01 == mc02) && (mc02 == mc03)) {
- cout << " sta " << ista << " mc0 " << mc01 << " " << mc02 << " " << mc03 << " mc1 " << mc11 << " " << mc12
- << " " << mc13 << " chi2 " << curr_tr.Chi2() / (2 * (curr_tr.NofHits() + 2) - 4) << " new "
- << new_chi2 / (2 * (new_L + 2) - 4) << endl;
- cout << " hits " << curr_trip->GetLHit() << " " << curr_trip->GetMHit() << " " << curr_trip->GetRHit()
- << " " << new_trip.GetLHit() << endl;
- }
- }
-
- int ndf = 2 * (new_L + 2) - 5;
-
- if (kGlobal == fTrackingMode) { //SGtrd2d!!!
- ndf = 2 * (new_L + 2) - 4;
- }
- else if (kMcbm == fTrackingMode) {
- ndf = 2 * (new_L + 2) - 4;
- }
- else {
- ndf = new_L; // TODO: SG: 2 * (new_L + 2) - 5
- }
-
- if (new_chi2 > fTrackChi2Cut * ndf) continue;
-
- // add new hit
- new_tr[ista] = curr_tr;
- new_tr[ista].AddHit(fWindowHits[new_trip.GetLHit()].Id());
-
- const int new_ista = ista + new_trip.GetMSta() - new_trip.GetLSta();
- new_tr[ista].SetChi2(new_chi2);
-
- CAFindTrack(new_ista, best_tr, &new_trip, new_tr[ista], min_best_l, new_tr);
- } // add triplet to track
- } // for neighbours
- } // level = 0
-}
diff --git a/reco/L1/L1Algo/L1CloneMerger.h b/reco/L1/L1Algo/L1CloneMerger.h
deleted file mode 100644
index d115dac8c5..0000000000
--- a/reco/L1/L1Algo/L1CloneMerger.h
+++ /dev/null
@@ -1,133 +0,0 @@
-/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergei Zharko [committer], Maksym Zyzak */
-
-/// \file L1CloneMerger.h
-/// \authors S.Zharko <s.zharko@gsi.de> (interface), M.Zyzak (original algorithm)
-/// \brief A class wrapper over clones merger algorithm for the L1 track finder (declaration)
-/// \since 22.07.2022
-
-#ifndef L1CloneMerger_h
-#define L1CloneMerger_h 1
-
-#include "CaHit.h" // For ca::HitIndex_t
-#include "CaSimd.h" // TEMPORARY FOR fvec, fscal
-#include "CaVector.h"
-
-
-using namespace cbm::algo::ca; //TODO: remove
-using namespace cbm::algo; //TODO: remove
-
-namespace cbm::algo::ca
-{
- class Track;
-}
-
-class L1Algo;
-
-namespace
-{ // TMP
- using cbm::algo::ca::Track;
- using cbm::algo::ca::Vector;
-} // namespace
-
-/// Class implements a clones merger algorithm for the CA (L1) track finder
-///
-class L1CloneMerger {
-public:
- /// Default constructor
- L1CloneMerger(const L1Algo& algo);
-
- /// Destructor
- ~L1CloneMerger();
-
- /// Copy constructor
- L1CloneMerger(const L1CloneMerger&) = delete;
-
- /// Move constructor
- L1CloneMerger(L1CloneMerger&&) = delete;
-
- /// Copy assignment operator
- L1CloneMerger& operator=(const L1CloneMerger&) = delete;
-
- /// Move assignment operator
- L1CloneMerger& operator=(L1CloneMerger&&) = delete;
-
- /// Registers
-
- /// Executes track clones merging algorithm and updates input containers
- /// \param extTracks Reference to the external container of reconstructed tracks
- /// \param extRecoHits Reference to the external container of reconstructed hit indexes
- void Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>&);
-
-private:
- // ***************
- // ** Functions **
- // ***************
-
- ///
- static void InvertCholesky(fvec a[15]);
-
- /// Multiplication of two symmetric matryces 5x5
- /// \param C Left matrix:
- /// C[0] C[1] C[3] C[6] C[10]
- /// C[1] C[2] C[4] C[7] C[11]
- /// C = C[3] C[4] C[5] C[8] C[12]
- /// C[6] C[7] C[8] C[9] C[13]
- /// C[10] C[11] C[12] C[13] C[14]
- /// \param V Right matrix:
- /// V[0] V[1] V[3] V[6] V[10]
- /// V[1] V[2] V[4] V[7] V[11]
- /// V = V[3] V[4] V[5] V[8] V[12]
- /// V[6] V[7] V[8] V[9] V[13]
- /// V[10] V[11] V[12] V[13] V[14]
- /// \param K Output: K = C * V
- static void MultiplySS(fvec const C[15], fvec const V[15], fvec K[5][5]);
-
- ///
- static void MultiplyMS(fvec const C[5][5], fvec const V[15], fvec K[15]);
-
- ///
- static void MultiplySR(fvec const C[15], fvec const r_in[5], fvec r_out[5]);
-
- ///
- static void FilterTracks(fvec const r[5], fvec const C[15], fvec const m[5], fvec const V[15], fvec R[5], fvec W[15],
- fvec* chi2);
-
-
- // ***************
- // ** Variables **
- // ***************
-
- /// First station of a track
- Vector<unsigned short> fTrackFirstStation {"L1CloneMerger::fTrackFirstStation"};
-
- /// Last station of a track
- Vector<unsigned short> fTrackLastStation {"L1CloneMerger::fTrackLastStation"};
-
- /// Index of the first hit of a track
- Vector<ca::HitIndex_t> fTrackFirstHit {"L1CloneMerger::fTrackFirstHit"};
-
- /// Index of the last hit of a track
- Vector<ca::HitIndex_t> fTrackLastHit {"L1CloneMerger::fTrackLastHit"};
-
- /// Index (TODO:??) of a track that can be merge with the given track
- Vector<unsigned short> fTrackNeighbour {"L1CloneMerger::fTrackNeighbour"};
-
- /// Chi2 value of the track merging procedure
- Vector<float> fTrackChi2 {"L1CloneMerger::fTrackChi2"};
-
- /// Flag: is the given track already stored to the output
- Vector<char> fTrackIsStored {"L1CloneMerger::fTrackIsStored"};
-
- /// Flag: is the track a downstream neighbour of another track
- Vector<char> fTrackIsDownstreamNeighbour {"L1CloneMerger::fTrackIsDownstreamNeighbour"};
-
- Vector<Track> fTracksNew {"L1CAMergerClones::fTracksNew"}; ///< vector of tracks after the merge
-
- Vector<ca::HitIndex_t> fRecoHitsNew {"L1CAMergerClones::fRecoHitsNew"}; ///< vector of track hits after the merge
-
- const L1Algo& frAlgo; ///< Reference to the main track finder algorithm class
-};
-
-#endif // L1CloneMerger_h
diff --git a/reco/L1/L1Algo/L1TripletConstructor.h b/reco/L1/L1Algo/L1TripletConstructor.h
deleted file mode 100644
index 89187ad476..0000000000
--- a/reco/L1/L1Algo/L1TripletConstructor.h
+++ /dev/null
@@ -1,117 +0,0 @@
-/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-#ifndef L1TripletConstructor_h
-#define L1TripletConstructor_h
-
-#include "CaField.h"
-#include "CaGridEntry.h"
-#include "CaStation.h"
-#include "CaTrackFit.h"
-#include "CaTrackParam.h"
-#include "CaTriplet.h"
-#include "CaVector.h"
-#include "L1Algo.h"
-
-namespace
-{ // TMP!!
- using cbm::algo::ca::Vector;
-}
-
-
-/// Construction of triplets for the CA tracker
-///
-class L1TripletConstructor {
-public:
- /// ------ Constructors and destructor ------
-
- /// Constructor
- /// \param nThreads Number of threads for multi-threaded mode
- L1TripletConstructor(L1Algo* algo);
-
- /// Copy constructor
- L1TripletConstructor(const L1TripletConstructor&) = delete;
-
- /// Move constructor
- L1TripletConstructor(L1TripletConstructor&&) = delete;
-
- /// Copy assignment operator
- L1TripletConstructor& operator=(const L1TripletConstructor&) = delete;
-
- /// Move assignment operator
- L1TripletConstructor& operator=(L1TripletConstructor&&) = delete;
-
- /// Destructor
- ~L1TripletConstructor() = default;
-
-
- /// ------ FUNCTIONAL PART ------
-
- void InitStations(int istal, int istam, int istar);
-
- void CreateTripletsForHit(int istal, int istam, int istar, ca::HitIndex_t ihl);
-
- const Vector<ca::Triplet>& GetTriplets() const { return fTriplets; }
-
- /// Find the doublets. Reformat data in the portion of doublets.
- void FindDoublets();
- void FitDoublets();
-
- /// Find triplets on station
- void FindTriplets();
-
- /// Add the middle hits to parameters estimation. Propagate to right station.
- /// Find the triplets (right hit). Reformat data in the portion of triplets.
- void FindRightHit();
-
- /// Fit Triplets
- void FitTriplets();
-
- /// Select triplets. Save them into vTriplets.
- void StoreTriplets();
-
- void CollectHits(const TrackParamV& Tr, const int iSta, const double chi2Cut, const int iMC,
- Vector<ca::HitIndex_t>& collectedHits, int maxNhits);
-
-private:
- /// left station
- const ca::Station& staL() { return *fStaL; }
- /// middle station
- const ca::Station& staM() { return *fStaM; }
- /// right station
- const ca::Station& staR() { return *fStaR; }
-
-private:
- bool fIsInitialized {false}; ///< is initialized;
- L1Algo* fAlgo {nullptr}; ///< pointer to L1Algo object
-
- int fIstaL {-1}; ///< left station index
- int fIstaM {-1}; ///< middle station index
- int fIstaR {-1}; ///< right station index
-
- const ca::Station* fStaL {nullptr}; ///< left station
- const ca::Station* fStaM {nullptr}; ///< mid station
- const ca::Station* fStaR {nullptr}; ///< right station
-
- const ca::Station* fFld0Sta[2]; // two stations for approximating the field between the target and the left hit
- const ca::Station* fFld1Sta[3]; // three stations for approximating the field between the left and the right hit
-
- ca::HitIndex_t fIhitL; ///< index of the left hit in fAlgo->fWindowHits
- TrackParamV fTrL;
- ca::FieldRegion fFldL;
-
- Vector<ca::HitIndex_t> fHitsM_2 {"L1TripletConstructor::fHitsM_2"};
- Vector<TrackParamV> fTracks_2 {"L1TripletConstructor::fTracks_2"};
-
- Vector<TrackParamV> fTracks_3 {"L1TripletConstructor::fTracks_3"};
- Vector<ca::HitIndex_t> fHitsM_3 {"L1TripletConstructor::fHitsM_3"};
- Vector<ca::HitIndex_t> fHitsR_3 {"L1TripletConstructor::fHitsR_3"};
-
- Vector<ca::Triplet> fTriplets {"L1TripletConstructor::fTriplets"};
-
- ca::TrackFit fFit;
-
-} _fvecalignment;
-
-#endif
diff --git a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
index b4e6b7850b..aab0a32588 100644
--- a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
+++ b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
@@ -22,7 +22,10 @@
#include "TText.h"
#include "TView3D.h"
-#include "L1Algo.h"
+#include "CaFramework.h"
+#include "CaSimd.h"
+
+using namespace cbm::algo::ca;
// NOTE: using std::cout, std::endl were removed from one of the headers, and it
// caused errors. Please, don't use "using" in headers
@@ -71,7 +74,7 @@ L1AlgoDraw::L1AlgoDraw()
ask = true;
}
-void L1AlgoDraw::InitL1Draw(L1Algo* algo_)
+void L1AlgoDraw::InitL1Draw(ca::Framework* algo_)
{
// algo = CbmL1::Instance()->algo;
algo = algo_;
@@ -81,11 +84,11 @@ void L1AlgoDraw::InitL1Draw(L1Algo* algo_)
for (unsigned int i = 0; i < algo->fWindowHits.size(); i++) {
vHits.push_back(algo->fWindowHits[i]);
}
- NStations = algo->GetParameters()->GetNstationsActive();
+ NStations = algo->GetParameters().GetNstationsActive();
for (int i = 0; i < NStations; i++) {
HitsStartIndex[i] = algo->fStationHitsStartIndex[i];
HitsStopIndex[i] = algo->fStationHitsStartIndex[i] + algo->fStationNhits[i];
- vStations[i] = algo->GetParameters()->GetStation(i);
+ vStations[i] = algo->GetParameters().GetStation(i);
}
}
diff --git a/reco/L1/L1Algo/utils/L1AlgoDraw.h b/reco/L1/L1Algo/utils/L1AlgoDraw.h
index dbb0d68b08..39f7f9110e 100644
--- a/reco/L1/L1Algo/utils/L1AlgoDraw.h
+++ b/reco/L1/L1Algo/utils/L1AlgoDraw.h
@@ -20,7 +20,11 @@ namespace
}
class TCanvas;
-class L1Algo;
+
+namespace cbm::algo::ca
+{
+ class Framework;
+}
class L1AlgoDraw {
struct Point {
@@ -32,7 +36,7 @@ class L1AlgoDraw {
public:
L1AlgoDraw();
- void InitL1Draw(L1Algo* algo_);
+ void InitL1Draw(ca::Framework* algo_);
void DrawMCTracks();
void DrawRecoTracks();
@@ -63,7 +67,7 @@ private:
void DrawTriplet(int il, int im, int ir);
void DrawDoublet(int il, int ir);
- L1Algo* algo {nullptr};
+ ca::Framework* algo {nullptr};
std::vector<ca::Hit> vHits {};
int HitsStartIndex[20] {0};
diff --git a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
index 07656bb9c1..4c6600bde4 100644
--- a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
+++ b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
@@ -14,12 +14,11 @@
#include "CbmL1.h"
#include "CbmL1Counters.h"
+#include "CbmL1MCTrack.h"
#include <iostream>
#include <vector>
-#include "L1Algo.h"
-
using std::cout;
using std::endl;
using std::vector;
diff --git a/reco/L1/L1Algo/utils/L1AlgoPulls.cxx b/reco/L1/L1Algo/utils/L1AlgoPulls.cxx
index 2f84229276..cf6133e965 100644
--- a/reco/L1/L1Algo/utils/L1AlgoPulls.cxx
+++ b/reco/L1/L1Algo/utils/L1AlgoPulls.cxx
@@ -4,6 +4,8 @@
#include "L1AlgoPulls.h"
+using namespace cbm::algo::ca;
+
const TString L1TrackParametersNames[TL1TrackParameters::NParameters] = {"x", "y", "tx", "ty", "qp"};
void L1AlgoPulls::Init()
diff --git a/reco/L1/L1Algo/utils/L1AlgoPulls.h b/reco/L1/L1Algo/utils/L1AlgoPulls.h
index e82bc5cf25..9ab379840e 100644
--- a/reco/L1/L1Algo/utils/L1AlgoPulls.h
+++ b/reco/L1/L1Algo/utils/L1AlgoPulls.h
@@ -24,7 +24,6 @@ const int NStations = 0;
#include <vector>
#include "CaTrackParam.h"
-#include "L1Algo.h"
using std::cout;
using std::endl;
@@ -37,7 +36,7 @@ struct TL1TrackParameters {
TL1TrackParameters() {};
- TL1TrackParameters(TrackParamV& T, int i)
+ TL1TrackParameters(cbm::algo::ca::TrackParamV& T, int i)
: x(T.GetX()[i])
, y(T.GetY()[i])
, tx(T.GetTx()[i])
@@ -90,7 +89,7 @@ public:
void Init();
// void AddVec(TrackParamV& T, ca::HitIndex_t ih);
- void AddOne(TrackParamV& T, int i, ca::HitIndex_t ih);
+ void AddOne(cbm::algo::ca::TrackParamV& T, int i, ca::HitIndex_t ih);
void Print(); // fast method to see pulls :)
void Build(bool draw = 1);
diff --git a/reco/L1/OffLineInterface/CbmL1GlobalTrackFinder.cxx b/reco/L1/OffLineInterface/CbmL1GlobalTrackFinder.cxx
index efea0e6fcf..de59be1823 100644
--- a/reco/L1/OffLineInterface/CbmL1GlobalTrackFinder.cxx
+++ b/reco/L1/OffLineInterface/CbmL1GlobalTrackFinder.cxx
@@ -39,7 +39,6 @@
#include <iostream>
#include <vector>
-#include "L1Algo/L1Algo.h"
using std::cout;
using std::endl;
diff --git a/reco/L1/OffLineInterface/CbmL1StsTrackFinder.cxx b/reco/L1/OffLineInterface/CbmL1StsTrackFinder.cxx
index a0c074d7ae..c24827ea0a 100644
--- a/reco/L1/OffLineInterface/CbmL1StsTrackFinder.cxx
+++ b/reco/L1/OffLineInterface/CbmL1StsTrackFinder.cxx
@@ -34,7 +34,6 @@
#include <iostream>
#include <vector>
-#include "L1Algo/L1Algo.h"
using std::cout;
using std::endl;
diff --git a/reco/L1/catools/CaToolsDebugger.cxx b/reco/L1/catools/CaToolsDebugger.cxx
index 9394c4993e..5783a6c29e 100644
--- a/reco/L1/catools/CaToolsDebugger.cxx
+++ b/reco/L1/catools/CaToolsDebugger.cxx
@@ -38,13 +38,13 @@ Debugger& Debugger::Instance()
Debugger::Debugger(const char* fileName)
{
/// Default constructor
- fFile = new TFile(fileName, "RECREATE");
+ fFileName = fileName;
}
void Debugger::Write()
{
/// Write data to file
- fIsWritten = true;
+ fIsWritten = true;
TDirectory* currDir = gDirectory;
if (fFile) {
fFile->cd();
@@ -58,12 +58,14 @@ void Debugger::Write()
Debugger::~Debugger()
{
/// Destructor
- if (!fIsWritten) { std::cerr << "CaDebugger: you forgot to Write()!!" << std::endl; }
+ if (!fIsWritten && fFile) { std::cerr << "CaDebugger: you forgot to Write()!!" << std::endl; }
}
void Debugger::AddNtuple(const char* name, const char* varlist)
{
/// add ntuple
+
+ if (!fFile) { fFile = new TFile(fFileName.c_str(), "RECREATE"); }
assert(fFile);
if (GetNtupleIndex(name) >= 0) return;
TDirectory* currDir = gDirectory;
diff --git a/reco/L1/catools/CaToolsDebugger.h b/reco/L1/catools/CaToolsDebugger.h
index fc9f2c8c31..8b7698110b 100644
--- a/reco/L1/catools/CaToolsDebugger.h
+++ b/reco/L1/catools/CaToolsDebugger.h
@@ -86,6 +86,7 @@ namespace cbm::ca::tools
}
private:
+ std::string fFileName {"CAdebug.root"};
bool fIsWritten {0};
TFile* fFile {nullptr};
std::vector<TNtuple*> fNtuples;
--
GitLab