From 0c41dedd3686733ae3cc4fa2b988b840864220cc Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Mon, 25 Jul 2022 20:11:08 +0200
Subject: [PATCH] L1: CAMergeClones routine was incapsulated into
L1ClonesMerger class
---
reco/L1/CMakeLists.txt | 3 +-
reco/L1/L1Algo/L1Algo.cxx | 4 +-
reco/L1/L1Algo/L1Algo.h | 54 +-
reco/L1/L1Algo/L1Assert.h | 8 +-
reco/L1/L1Algo/L1CATrackFinder.cxx | 12 +-
...L1CAMergeClones.cxx => L1ClonesMerger.cxx} | 591 +++++++++---------
reco/L1/L1Algo/L1ClonesMerger.h | 119 ++++
reco/L1/L1Algo/L1Track.h | 32 +-
reco/L1/L1Algo/L1Utils.h | 74 +--
9 files changed, 472 insertions(+), 425 deletions(-)
rename reco/L1/L1Algo/{L1CAMergeClones.cxx => L1ClonesMerger.cxx} (62%)
create mode 100644 reco/L1/L1Algo/L1ClonesMerger.h
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 25d1c057ec..53f1bba9bd 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -129,7 +129,6 @@ L1Algo/L1Algo.cxx
L1Algo/L1CATrackFinder.cxx
L1Algo/L1TrackExtender.cxx
L1Algo/L1TrackFitter.cxx
-L1Algo/L1CAMergeClones.cxx
L1Algo/L1HitsSortHelper.cxx
L1Algo/L1Grid.cxx
CbmL1Performance.cxx
@@ -149,6 +148,7 @@ L1Algo/L1CAIteration.cxx
L1Algo/L1BaseStationInfo.cxx
L1Algo/L1InitManager.cxx
L1Algo/L1Parameters.cxx
+L1Algo/L1ClonesMerger.cxx
L1Algo/L1ConfigRW.cxx
L1Algo/utils/L1AlgoDraw.cxx
L1Algo/utils/L1AlgoEfficiencyPerformance.cxx
@@ -286,6 +286,7 @@ Install(FILES CbmL1Counters.h
L1Algo/L1InitManager.h
L1Algo/L1CAIteration.h
L1Algo/L1Parameters.h
+ L1Algo/L1ClonesMerger.h
L1Algo/L1ConfigRW.h
L1Algo/L1Constants.h
L1Algo/L1Utils.h
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 5a27570b36..a934089cbc 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -93,7 +93,7 @@ void L1Algo::Init(const bool UseHitErrors, const TrackingMode mode, const bool M
LOG(info) << fParameters.ToString(3);
// Get number of station
- fNstations = fInitManager.GetNstationsActive();
+ fNstations = fParameters.GetNstationsActive();
fTrackingLevel = fInitManager.GetTrackingLevel();
fGhostSuppression = fInitManager.GetGhostSuppression();
@@ -116,7 +116,7 @@ void L1Algo::SetData(L1Vector<L1Hit>& Hits_, int nStrips_, L1Vector<unsigned cha
int nHits = vHits->size();
- NHitsIsecAll = nHits;
+ NHitsIsecAll = nHits; // TODO: Is it needed?
vNotUsedHits_A.reset(nHits);
vNotUsedHits_B.reset(nHits);
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index ebb8d3bc47..b1f51f0d06 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -43,6 +43,7 @@ class L1AlgoDraw;
#include <map>
#include "L1Branch.h"
+#include "L1ClonesMerger.h"
#include "L1Field.h"
#include "L1Grid.h"
#include "L1Hit.h"
@@ -166,39 +167,6 @@ public:
Tindex fDupletPortionStopIndex[L1Constants::size::kMaxNstations] {0}; // end of the duplet portions for the station
L1Vector<Tindex> fDupletPortionSize {"L1Algo::fDupletPortionSize"}; // N duplets in a portion
- /********************************************************************************************/ /**
- * Temporary vectors used by the clone merger
- * TODO: Probably, the subclass L1TrackMerger for clones merger would help to improve
- * readability (S.Zharko)
- ***********************************************************************************************/
-
- //
- // Vectors that are parallel to fTracks
- //
- /// First station of a track
- L1Vector<unsigned short> fMergerTrackFirstStation {"L1Algo::fMergerTrackFirstStation"};
- /// Last station of a track
- L1Vector<unsigned short> fMergerTrackLastStation {"L1Algo::fMergerTrackLastStation"};
- /// Index of the first hit of a track
- L1Vector<L1HitIndex_t> fMergerTrackFirstHit {"L1Algo::fMergerTrackFirstHit"};
- /// Index of the last hit of a track
- L1Vector<L1HitIndex_t> fMergerTrackLastHit {"L1Algo::fMergerTrackLastHit"};
- /// Index (TODO:??) of a track that can be merge with the given track
- L1Vector<unsigned short> fMergerTrackNeighbour {"L1Algo::fMergerTrackNeighbour"};
- /// Chi2 value of the track merging procedure
- L1Vector<float> fMergerTrackChi2 {"L1Algo::fMergerTrackChi2"};
- /// Flag: is the given track already stored to the output
- L1Vector<char> fMergerTrackIsStored {"L1Algo::fMergerTrackIsStored"};
- /// Flag: is the track a downstream neighbour of another track
- L1Vector<char> fMergerTrackIsDownstreamNeighbour {"L1Algo::fMergerTrackIsDownstreamNeighbour"};
- //
- // Utility vectors
- //
- /// Tracks after the merging procedure
- L1Vector<L1Track> fMergerTracksNew {"L1Algo::fMergerTracksNew"}; // vector of tracks after the merge
- L1Vector<L1HitIndex_t> fMergerRecoHitsNew {"L1Algo::fMergerRecoHitsNew"}; // vector of track hits after the merge
-
-
#ifdef DRAW
L1AlgoDraw* draw {nullptr};
void DrawRecoTracksTime(const L1Vector<CbmL1Track>& tracks);
@@ -261,7 +229,7 @@ public:
public:
int fNstrips {0}; ///< number of strips
L1Vector<L1Hit>* vHits {nullptr}; ///< hits as a combination of front and back strips and z-position
- L1Grid vGrid[L1Constants::size::kMaxNstations]; ///< hits as a combination of front and back strips and z-position
+ L1Grid vGrid[L1Constants::size::kMaxNstations]; ///<
L1Grid vGridTime[L1Constants::size::kMaxNstations]; ///<
L1Vector<unsigned char>* fStripFlag {nullptr}; // information of hits station & using hits in tracks;
@@ -386,15 +354,15 @@ public:
// TODO: remove it (S.Zharko)
/// Gets a pointer to the L1Algo parameters object
- const L1Parameters* GetParameters() { return &fParameters; }
+ const L1Parameters* GetParameters() const { return &fParameters; }
/// Gets a pointer to the L1Algo initialization object
L1InitManager* GetInitManager() { return &fInitManager; }
private:
- L1Parameters fParameters {}; ///< Object of L1Algo parameters class
- L1InitManager fInitManager {}; ///< Object of L1Algo initialization manager class
-
+ L1Parameters fParameters {}; ///< Object of L1Algo parameters class
+ L1InitManager fInitManager {}; ///< Object of L1Algo initialization manager class
+ L1ClonesMerger fClonesMerger {*this}; ///< Object of L1Algo clones merger algorithm
/*********************************************************************************************/ /**
* ------ FUNCTIONAL PART ------
@@ -436,16 +404,6 @@ private:
/// \return chi2
fscal BranchExtender(L1Branch& t);
- /// ----- Subroutines used by L1Algo::CAMergeClones() ------
- void InvertCholetsky(fvec a[15]);
- void MultiplySS(fvec const C[15], fvec const V[15], fvec K[5][5]);
- void MultiplyMS(fvec const C[5][5], fvec const V[15], fvec K[15]);
- void MultiplySR(fvec const C[15], fvec const r_in[5], fvec r_out[5]);
- 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);
- ///
- void CAMergeClones();
-
inline __attribute__((always_inline)) void PackLocation(unsigned int& location, unsigned int& triplet,
unsigned int iStation, unsigned int& thread)
{
diff --git a/reco/L1/L1Algo/L1Assert.h b/reco/L1/L1Algo/L1Assert.h
index 9b3dc67257..a5e03567c5 100644
--- a/reco/L1/L1Algo/L1Assert.h
+++ b/reco/L1/L1Algo/L1Assert.h
@@ -38,7 +38,7 @@ namespace L1Assert
/// Basic template function. Usage: place "level <= L1Assert::kAssertionLevel" as a template parameter
template<bool IsAsserted>
- int DoAssertion(int level, bool condition, const char* msg, const char* fileName, int lineNo)
+ void DoAssertion(int level, bool condition, const char* msg, const char* fileName, int lineNo)
{
if (!condition) {
LOG(fatal) << '\n'
@@ -48,15 +48,13 @@ namespace L1Assert
<< " ***** line: " << lineNo;
std::abort(); // keep it here, because sometimes LOG(fatal) does not work (for example, in your unit testes)
}
- return 1;
}
/// Specialization in case of IsAsserted = false, i.e. the assertion is not made
template<>
- inline __attribute__((always_inline)) int DoAssertion<false>(int /*level*/, bool /*condition*/, const char* /*msg*/,
- const char* /*fileName*/, int /*lineNo*/)
+ inline __attribute__((always_inline)) void DoAssertion<false>(int /*level*/, bool /*condition*/, const char* /*msg*/,
+ const char* /*fileName*/, int /*lineNo*/)
{
- return 0;
}
}; // namespace L1Assert
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 1f34006ce5..1efe5236e2 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -754,6 +754,7 @@ inline void L1Algo::findTripletsStep0( // input
}
else if (kGlobal == fTrackingMode) {
fit.L1AddMaterial(T2, fParameters.GetMaterialThickness(istam, T2.x, T2.y), fMaxInvMom, 1);
+ }
else {
fit.L1AddMaterial(T2, fParameters.GetMaterialThickness(istam, T2.x, T2.y), T2.qp, 1);
}
@@ -827,6 +828,8 @@ inline void L1Algo::findTripletsStep0( // input
int indR = HitsUnusedStartIndex[iStaR] + irh;
if (GetMcTrackIdForUnusedHit(indM) != GetMcTrackIdForUnusedHit(indR)) { continue; }
}
+
+#ifdef USE_EVENT_NUMBER // NOTE:
if ((T2.n[i2_4] != hitr.n)) continue;
#endif // USE_EVENT_NUMBER
const fscal zr = hitr.Z();
@@ -2634,7 +2637,14 @@ void L1Algo::CATrackFinder()
c_timerG.Start();
#endif
- if constexpr (L1Constants::control::kIfMergeClones) { CAMergeClones(); }
+ if constexpr (L1Constants::control::kIfMergeClones) {
+ //CAMergeClones();
+ // Fit tracks
+ this->L1KFTrackFitter();
+
+ // Merge clones
+ fClonesMerger.Exec(fTracks, fRecoHits);
+ }
#ifdef XXX
c_timerG.Stop();
diff --git a/reco/L1/L1Algo/L1CAMergeClones.cxx b/reco/L1/L1Algo/L1ClonesMerger.cxx
similarity index 62%
rename from reco/L1/L1Algo/L1CAMergeClones.cxx
rename to reco/L1/L1Algo/L1ClonesMerger.cxx
index a5476902ae..26d1b3f94f 100644
--- a/reco/L1/L1Algo/L1CAMergeClones.cxx
+++ b/reco/L1/L1Algo/L1ClonesMerger.cxx
@@ -1,228 +1,52 @@
-/* Copyright (C) 2010-2018 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+/* Copyright (C) 2010-2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: Maksym Zyzak [committer] */
-
-/*
- *=====================================================
- *
- * CBM Level 1 Reconstruction
- *
- * Authors: M.Zyzak
- *
- * e-mail :
- *
- *=====================================================
- *
- * Merge Clones
- *
- */
+ Authors: Sergei Zharko, Maksym Zyzak [committer]*/
+
+/// \file L1ClonesMerger.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)
+/// \since 22.07.2022 (second version)
+
+#include "L1ClonesMerger.h"
#include <iostream>
#include "L1Algo.h"
#include "L1Extrapolation.h"
+#include "L1Parameters.h"
+#include "L1Track.h"
-// using namespace std;
-using std::cout;
-using std::endl;
-
-void L1Algo::InvertCholetsky(fvec a[15])
-{
- fvec d[5], uud, u[5][5];
- for (int i = 0; i < 5; i++) {
- d[i] = 0.f;
- for (int j = 0; j < 5; j++)
- u[i][j] = 0.f;
- }
-
- for (int i = 0; i < 5; i++) {
- uud = 0.f;
- for (int j = 0; j < i; j++)
- uud += u[j][i] * u[j][i] * d[j];
- uud = a[i * (i + 3) / 2] - uud;
-
- fvec smallval = 1.e-12;
- fvec initialised = fvec(uud < smallval);
- uud = ((!initialised) & uud) + (smallval & initialised);
-
- d[i] = uud / fabs(uud);
- u[i][i] = sqrt(fabs(uud));
-
- for (int j = i + 1; j < 5; j++) {
- uud = 0.f;
- for (int k = 0; k < i; k++)
- uud += u[k][i] * u[k][j] * d[k];
- uud = a[j * (j + 1) / 2 + i] /*a[i][j]*/ - uud;
- u[i][j] = d[i] / u[i][i] * uud;
- }
- }
-
- fvec u1[5];
-
- for (int i = 0; i < 5; i++) {
- u1[i] = u[i][i];
- u[i][i] = 1.f / u[i][i];
- }
- for (int i = 0; i < 4; i++) {
- u[i][i + 1] = -u[i][i + 1] * u[i][i] * u[i + 1][i + 1];
- }
- for (int i = 0; i < 3; i++) {
- u[i][i + 2] = u[i][i + 1] * u1[i + 1] * u[i + 1][i + 2] - u[i][i + 2] * u[i][i] * u[i + 2][i + 2];
- }
- for (int i = 0; i < 2; i++) {
- u[i][i + 3] = u[i][i + 2] * u1[i + 2] * u[i + 2][i + 3] - u[i][i + 3] * u[i][i] * u[i + 3][i + 3];
- u[i][i + 3] -= u[i][i + 1] * u1[i + 1] * (u[i + 1][i + 2] * u1[i + 2] * u[i + 2][i + 3] - u[i + 1][i + 3]);
- }
- u[0][4] = u[0][2] * u1[2] * u[2][4] - u[0][4] * u[0][0] * u[4][4];
- u[0][4] += u[0][1] * u1[1] * (u[1][4] - u[1][3] * u1[3] * u[3][4] - u[1][2] * u1[2] * u[2][4]);
- u[0][4] += u[3][4] * u1[3] * (u[0][3] - u1[2] * u[2][3] * (u[0][2] - u[0][1] * u1[1] * u[1][2]));
-
- for (int i = 0; i < 5; i++)
- a[i + 10] = u[i][4] * d[4] * u[4][4];
- for (int i = 0; i < 4; i++)
- a[i + 6] = u[i][3] * u[3][3] * d[3] + u[i][4] * u[3][4] * d[4];
- for (int i = 0; i < 3; i++)
- a[i + 3] = u[i][2] * u[2][2] * d[2] + u[i][3] * u[2][3] * d[3] + u[i][4] * u[2][4] * d[4];
- for (int i = 0; i < 2; i++)
- a[i + 1] =
- u[i][1] * u[1][1] * d[1] + u[i][2] * u[1][2] * d[2] + u[i][3] * u[1][3] * d[3] + u[i][4] * u[1][4] * d[4];
- a[0] = u[0][0] * u[0][0] * d[0] + u[0][1] * u[0][1] * d[1] + u[0][2] * u[0][2] * d[2] + u[0][3] * u[0][3] * d[3]
- + u[0][4] * u[0][4] * d[4];
-}
-
-void L1Algo::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];
- K[0][2] = C[0] * V[3] + C[1] * V[4] + C[3] * V[5] + C[6] * V[8] + C[10] * V[12];
- K[0][3] = C[0] * V[6] + C[1] * V[7] + C[3] * V[8] + C[6] * V[9] + C[10] * V[13];
- K[0][4] = C[0] * V[10] + C[1] * V[11] + C[3] * V[12] + C[6] * V[13] + C[10] * V[14];
-
- K[1][0] = C[1] * V[0] + C[2] * V[1] + C[4] * V[3] + C[7] * V[6] + C[11] * V[10];
- K[1][1] = C[1] * V[1] + C[2] * V[2] + C[4] * V[4] + C[7] * V[7] + C[11] * V[11];
- K[1][2] = C[1] * V[3] + C[2] * V[4] + C[4] * V[5] + C[7] * V[8] + C[11] * V[12];
- K[1][3] = C[1] * V[6] + C[2] * V[7] + C[4] * V[8] + C[7] * V[9] + C[11] * V[13];
- K[1][4] = C[1] * V[10] + C[2] * V[11] + C[4] * V[12] + C[7] * V[13] + C[11] * V[14];
-
- K[2][0] = C[3] * V[0] + C[4] * V[1] + C[5] * V[3] + C[8] * V[6] + C[12] * V[10];
- K[2][1] = C[3] * V[1] + C[4] * V[2] + C[5] * V[4] + C[8] * V[7] + C[12] * V[11];
- K[2][2] = C[3] * V[3] + C[4] * V[4] + C[5] * V[5] + C[8] * V[8] + C[12] * V[12];
- K[2][3] = C[3] * V[6] + C[4] * V[7] + C[5] * V[8] + C[8] * V[9] + C[12] * V[13];
- K[2][4] = C[3] * V[10] + C[4] * V[11] + C[5] * V[12] + C[8] * V[13] + C[12] * V[14];
+// ---------------------------------------------------------------------------------------------------------------------
+//
+L1ClonesMerger::L1ClonesMerger(const L1Algo& algo) : frAlgo(algo) {}
- K[3][0] = C[6] * V[0] + C[7] * V[1] + C[8] * V[3] + C[9] * V[6] + C[13] * V[10];
- K[3][1] = C[6] * V[1] + C[7] * V[2] + C[8] * V[4] + C[9] * V[7] + C[13] * V[11];
- K[3][2] = C[6] * V[3] + C[7] * V[4] + C[8] * V[5] + C[9] * V[8] + C[13] * V[12];
- K[3][3] = C[6] * V[6] + C[7] * V[7] + C[8] * V[8] + C[9] * V[9] + C[13] * V[13];
- K[3][4] = C[6] * V[10] + C[7] * V[11] + C[8] * V[12] + C[9] * V[13] + C[13] * V[14];
-
- K[4][0] = C[10] * V[0] + C[11] * V[1] + C[12] * V[3] + C[13] * V[6] + C[14] * V[10];
- K[4][1] = C[10] * V[1] + C[11] * V[2] + C[12] * V[4] + C[13] * V[7] + C[14] * V[11];
- K[4][2] = C[10] * V[3] + C[11] * V[4] + C[12] * V[5] + C[13] * V[8] + C[14] * V[12];
- K[4][3] = C[10] * V[6] + C[11] * V[7] + C[12] * V[8] + C[13] * V[9] + C[14] * V[13];
- K[4][4] = C[10] * V[10] + C[11] * V[11] + C[12] * V[12] + C[13] * V[13] + C[14] * V[14];
-}
-
-void L1Algo::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];
-
- K[1] = C[1][0] * V[0] + C[1][1] * V[1] + C[1][2] * V[3] + C[1][3] * V[6] + C[1][4] * V[10];
- K[2] = C[1][0] * V[1] + C[1][1] * V[2] + C[1][2] * V[4] + C[1][3] * V[7] + C[1][4] * V[11];
-
- K[3] = C[2][0] * V[0] + C[2][1] * V[1] + C[2][2] * V[3] + C[2][3] * V[6] + C[2][4] * V[10];
- K[4] = C[2][0] * V[1] + C[2][1] * V[2] + C[2][2] * V[4] + C[2][3] * V[7] + C[2][4] * V[11];
- K[5] = C[2][0] * V[3] + C[2][1] * V[4] + C[2][2] * V[5] + C[2][3] * V[8] + C[2][4] * V[12];
-
- K[6] = C[3][0] * V[0] + C[3][1] * V[1] + C[3][2] * V[3] + C[3][3] * V[6] + C[3][4] * V[10];
- K[7] = C[3][0] * V[1] + C[3][1] * V[2] + C[3][2] * V[4] + C[3][3] * V[7] + C[3][4] * V[11];
- K[8] = C[3][0] * V[3] + C[3][1] * V[4] + C[3][2] * V[5] + C[3][3] * V[8] + C[3][4] * V[12];
- K[9] = C[3][0] * V[6] + C[3][1] * V[7] + C[3][2] * V[8] + C[3][3] * V[9] + C[3][4] * V[13];
-
- K[10] = C[4][0] * V[0] + C[4][1] * V[1] + C[4][2] * V[3] + C[4][3] * V[6] + C[4][4] * V[10];
- K[11] = C[4][0] * V[1] + C[4][1] * V[2] + C[4][2] * V[4] + C[4][3] * V[7] + C[4][4] * V[11];
- K[12] = C[4][0] * V[3] + C[4][1] * V[4] + C[4][2] * V[5] + C[4][3] * V[8] + C[4][4] * V[12];
- K[13] = C[4][0] * V[6] + C[4][1] * V[7] + C[4][2] * V[8] + C[4][3] * V[9] + C[4][4] * V[13];
- K[14] = C[4][0] * V[10] + C[4][1] * V[11] + C[4][2] * V[12] + C[4][3] * V[13] + C[4][4] * V[14];
-}
-
-void L1Algo::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];
- r_out[2] = r_in[0] * C[3] + r_in[1] * C[4] + r_in[2] * C[5] + r_in[3] * C[8] + r_in[4] * C[12];
- r_out[3] = r_in[0] * C[6] + r_in[1] * C[7] + r_in[2] * C[8] + r_in[3] * C[9] + r_in[4] * C[13];
- r_out[4] = r_in[0] * C[10] + r_in[1] * C[11] + r_in[2] * C[12] + r_in[3] * C[13] + r_in[4] * C[14];
-}
-
-void L1Algo::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++) {
- if (W) W[i] = C[i];
- S[i] = C[i] + V[i];
- }
-
- InvertCholetsky(S);
-
- fvec dzeta[5];
-
- for (int i = 0; i < 5; i++)
- dzeta[i] = m[i] - r[i];
-
- if (W && R) {
- for (int i = 0; i < 5; i++)
- R[i] = r[i];
-
- fvec K[5][5];
- MultiplySS(C, S, K);
-
- fvec KC[15];
- MultiplyMS(K, C, KC);
- for (int i = 0; i < 15; i++)
- W[i] -= KC[i];
-
- fvec kd;
- for (int i = 0; i < 5; i++) {
- kd = 0.f;
- for (int j = 0; j < 5; j++)
- kd += K[i][j] * dzeta[j];
- R[i] += kd;
- }
- }
- if (chi2) {
- fvec S_dzeta[5];
- MultiplySR(S, dzeta, S_dzeta);
- *chi2 = dzeta[0] * S_dzeta[0] + dzeta[1] * S_dzeta[1] + dzeta[2] * S_dzeta[2] + dzeta[3] * S_dzeta[3]
- + dzeta[4] * S_dzeta[4];
- }
-}
+// ---------------------------------------------------------------------------------------------------------------------
+//
+L1ClonesMerger::~L1ClonesMerger() {}
-void L1Algo::CAMergeClones()
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void L1ClonesMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>& extRecoHits)
{
+ L1Vector<unsigned short>& firstStation = fTrackFirstStation;
+ L1Vector<unsigned short>& lastStation = fTrackLastStation;
+ L1Vector<L1HitIndex_t>& firstHit = fTrackFirstHit;
+ L1Vector<L1HitIndex_t>& lastHit = fTrackLastHit;
+ L1Vector<unsigned short>& neighbour = fTrackNeighbour;
+ L1Vector<float>& trackChi2 = fTrackChi2;
+ L1Vector<char>& isStored = fTrackIsStored;
+ L1Vector<char>& isDownstreamNeighbour = fTrackIsDownstreamNeighbour;
- L1Vector<unsigned short>& firstStation = fMergerTrackFirstStation;
- L1Vector<unsigned short>& lastStation = fMergerTrackLastStation;
- L1Vector<L1HitIndex_t>& firstHit = fMergerTrackFirstHit;
- L1Vector<L1HitIndex_t>& lastHit = fMergerTrackLastHit;
- L1Vector<unsigned short>& neighbour = fMergerTrackNeighbour;
- L1Vector<float>& trackChi2 = fMergerTrackChi2;
- L1Vector<char>& isStored = fMergerTrackIsStored;
- L1Vector<char>& isDownstreamNeighbour = fMergerTrackIsDownstreamNeighbour;
-
- int nTracks = fTracks.size();
+ int nTracks = extTracks.size();
assert(nTracks < std::numeric_limits<unsigned short>::max());
constexpr unsigned short kNoNeighbour = std::numeric_limits<unsigned short>::max();
- // vector< L1Track > vTracksNew;
- fMergerTracksNew.clear();
- fMergerTracksNew.reserve(nTracks);
- // vector< L1HitIndex_t > fRecoHitsNew;
- fMergerRecoHitsNew.clear();
- fMergerRecoHitsNew.reserve(fRecoHits.size());
+ fTracksNew.clear();
+ fTracksNew.reserve(nTracks);
+ fRecoHitsNew.clear();
+ fRecoHitsNew.reserve(extRecoHits.size());
firstStation.reset(nTracks);
lastStation.reset(nTracks);
@@ -240,10 +64,10 @@ void L1Algo::CAMergeClones()
#endif
for (int iTr = 0; iTr < nTracks; iTr++) {
firstHit[iTr] = start_hit;
- firstStation[iTr] = (*fStripFlag)[(*vHits)[fRecoHits[start_hit]].f] / 4;
- start_hit += fTracks[iTr].NHits - 1;
+ firstStation[iTr] = (*frAlgo.fStripFlag)[(*frAlgo.vHits)[extRecoHits[start_hit]].f] / 4;
+ start_hit += extTracks[iTr].NHits - 1;
lastHit[iTr] = start_hit;
- lastStation[iTr] = (*fStripFlag)[(*vHits)[fRecoHits[start_hit]].f] / 4;
+ lastStation[iTr] = (*frAlgo.fStripFlag)[(*frAlgo.vHits)[extRecoHits[start_hit]].f] / 4;
start_hit++;
isStored[iTr] = false;
@@ -252,103 +76,92 @@ void L1Algo::CAMergeClones()
isDownstreamNeighbour[iTr] = false;
}
- L1KFTrackFitter();
- //KFTrackFitter_simple();
-
L1TrackPar Tb;
L1TrackPar Tf;
L1FieldValue fBm, fBb, fBf _fvecalignment;
L1FieldRegion fld _fvecalignment;
- unsigned char maxLengthForMerge = static_cast<unsigned char>(fNstations - 3); // max length for merge
+ // Max length for merging
+ unsigned char maxLengthForMerge = static_cast<unsigned char>(frAlgo.GetParameters()->GetNstationsActive() - 3);
#ifdef OMP
#pragma omp parallel for
#endif
for (int iTr = 0; iTr < nTracks; iTr++) {
- if (fTracks[iTr].NHits > maxLengthForMerge) continue;
+ if (extTracks[iTr].NHits > maxLengthForMerge) continue;
for (int jTr = 0; jTr < nTracks; jTr++) {
- if (fTracks[jTr].NHits > maxLengthForMerge) continue;
+ if (extTracks[jTr].NHits > maxLengthForMerge) continue;
if (iTr == jTr) continue;
if (firstStation[iTr] <= lastStation[jTr]) continue;
- // if(vTracks[iTr].n != vTracks[jTr].n) continue;
- // if (fabs(vTracks[iTr].fTrackTime - vTracks[jTr].fTrackTime) > 3) continue;
-
- //if((vTracks[iTr].TFirst[4] - vTracks[jTr].TFirst[4])*(vTracks[iTr].TFirst[4] - vTracks[jTr].TFirst[4])
- // > 9*(vTracks[iTr].CFirst[14]+vTracks[jTr].CFirst[14]) ) continue;
-
unsigned short stab = firstStation[iTr];
- Tb.x = fTracks[iTr].TFirst[0];
- Tb.y = fTracks[iTr].TFirst[1];
- Tb.tx = fTracks[iTr].TFirst[2];
- Tb.ty = fTracks[iTr].TFirst[3];
- Tb.qp = fTracks[iTr].TFirst[4];
- Tb.z = fTracks[iTr].TFirst[5];
- Tb.t = fTracks[iTr].TFirst[6];
- Tb.C00 = fTracks[iTr].CFirst[0];
- Tb.C10 = fTracks[iTr].CFirst[1];
- Tb.C11 = fTracks[iTr].CFirst[2];
- Tb.C20 = fTracks[iTr].CFirst[3];
- Tb.C21 = fTracks[iTr].CFirst[4];
- Tb.C22 = fTracks[iTr].CFirst[5];
- Tb.C30 = fTracks[iTr].CFirst[6];
- Tb.C31 = fTracks[iTr].CFirst[7];
- Tb.C32 = fTracks[iTr].CFirst[8];
- Tb.C33 = fTracks[iTr].CFirst[9];
- Tb.C40 = fTracks[iTr].CFirst[10];
- Tb.C41 = fTracks[iTr].CFirst[11];
- Tb.C42 = fTracks[iTr].CFirst[12];
- Tb.C43 = fTracks[iTr].CFirst[13];
- Tb.C44 = fTracks[iTr].CFirst[14];
- Tb.C50 = fTracks[iTr].CFirst[15];
- Tb.C51 = fTracks[iTr].CFirst[16];
- Tb.C52 = fTracks[iTr].CFirst[17];
- Tb.C53 = fTracks[iTr].CFirst[18];
- Tb.C54 = fTracks[iTr].CFirst[19];
- Tb.C55 = fTracks[iTr].CFirst[20];
+ Tb.x = extTracks[iTr].TFirst[0];
+ Tb.y = extTracks[iTr].TFirst[1];
+ Tb.tx = extTracks[iTr].TFirst[2];
+ Tb.ty = extTracks[iTr].TFirst[3];
+ Tb.qp = extTracks[iTr].TFirst[4];
+ Tb.z = extTracks[iTr].TFirst[5];
+ Tb.t = extTracks[iTr].TFirst[6];
+ Tb.C00 = extTracks[iTr].CFirst[0];
+ Tb.C10 = extTracks[iTr].CFirst[1];
+ Tb.C11 = extTracks[iTr].CFirst[2];
+ Tb.C20 = extTracks[iTr].CFirst[3];
+ Tb.C21 = extTracks[iTr].CFirst[4];
+ Tb.C22 = extTracks[iTr].CFirst[5];
+ Tb.C30 = extTracks[iTr].CFirst[6];
+ Tb.C31 = extTracks[iTr].CFirst[7];
+ Tb.C32 = extTracks[iTr].CFirst[8];
+ Tb.C33 = extTracks[iTr].CFirst[9];
+ Tb.C40 = extTracks[iTr].CFirst[10];
+ Tb.C41 = extTracks[iTr].CFirst[11];
+ Tb.C42 = extTracks[iTr].CFirst[12];
+ Tb.C43 = extTracks[iTr].CFirst[13];
+ Tb.C44 = extTracks[iTr].CFirst[14];
+ Tb.C50 = extTracks[iTr].CFirst[15];
+ Tb.C51 = extTracks[iTr].CFirst[16];
+ Tb.C52 = extTracks[iTr].CFirst[17];
+ Tb.C53 = extTracks[iTr].CFirst[18];
+ Tb.C54 = extTracks[iTr].CFirst[19];
+ Tb.C55 = extTracks[iTr].CFirst[20];
unsigned short staf = lastStation[jTr];
- Tf.x = fTracks[jTr].TLast[0];
- Tf.y = fTracks[jTr].TLast[1];
- Tf.tx = fTracks[jTr].TLast[2];
- Tf.ty = fTracks[jTr].TLast[3];
- Tf.qp = fTracks[jTr].TLast[4];
- Tf.z = fTracks[jTr].TLast[5];
- Tf.t = fTracks[jTr].TLast[6];
- Tf.C00 = fTracks[jTr].CLast[0];
- Tf.C10 = fTracks[jTr].CLast[1];
- Tf.C11 = fTracks[jTr].CLast[2];
- Tf.C20 = fTracks[jTr].CLast[3];
- Tf.C21 = fTracks[jTr].CLast[4];
- Tf.C22 = fTracks[jTr].CLast[5];
- Tf.C30 = fTracks[jTr].CLast[6];
- Tf.C31 = fTracks[jTr].CLast[7];
- Tf.C32 = fTracks[jTr].CLast[8];
- Tf.C33 = fTracks[jTr].CLast[9];
- Tf.C40 = fTracks[jTr].CLast[10];
- Tf.C41 = fTracks[jTr].CLast[11];
- Tf.C42 = fTracks[jTr].CLast[12];
- Tf.C43 = fTracks[jTr].CLast[13];
- Tf.C44 = fTracks[jTr].CLast[14];
- Tf.C50 = fTracks[jTr].CLast[15];
- Tf.C51 = fTracks[jTr].CLast[16];
- Tf.C52 = fTracks[jTr].CLast[17];
- Tf.C53 = fTracks[jTr].CLast[18];
- Tf.C54 = fTracks[jTr].CLast[19];
- Tf.C55 = fTracks[jTr].CLast[20];
-
- //std::cout << "!!!!!!! Chi2 !!!!!! "<<fTracks[iTr].TFirst[0]<<" "<<fTracks[jTr].TLast[0]<<std::endl;
-
- //if(((Tf.qp - Tb.qp)*(Tf.qp - Tb.qp)/(Tb.C44+Tf.C44))[0] > 25*10*7) continue;
+ Tf.x = extTracks[jTr].TLast[0];
+ Tf.y = extTracks[jTr].TLast[1];
+ Tf.tx = extTracks[jTr].TLast[2];
+ Tf.ty = extTracks[jTr].TLast[3];
+ Tf.qp = extTracks[jTr].TLast[4];
+ Tf.z = extTracks[jTr].TLast[5];
+ Tf.t = extTracks[jTr].TLast[6];
+ Tf.C00 = extTracks[jTr].CLast[0];
+ Tf.C10 = extTracks[jTr].CLast[1];
+ Tf.C11 = extTracks[jTr].CLast[2];
+ Tf.C20 = extTracks[jTr].CLast[3];
+ Tf.C21 = extTracks[jTr].CLast[4];
+ Tf.C22 = extTracks[jTr].CLast[5];
+ Tf.C30 = extTracks[jTr].CLast[6];
+ Tf.C31 = extTracks[jTr].CLast[7];
+ Tf.C32 = extTracks[jTr].CLast[8];
+ Tf.C33 = extTracks[jTr].CLast[9];
+ Tf.C40 = extTracks[jTr].CLast[10];
+ Tf.C41 = extTracks[jTr].CLast[11];
+ Tf.C42 = extTracks[jTr].CLast[12];
+ Tf.C43 = extTracks[jTr].CLast[13];
+ Tf.C44 = extTracks[jTr].CLast[14];
+ Tf.C50 = extTracks[jTr].CLast[15];
+ Tf.C51 = extTracks[jTr].CLast[16];
+ Tf.C52 = extTracks[jTr].CLast[17];
+ Tf.C53 = extTracks[jTr].CLast[18];
+ Tf.C54 = extTracks[jTr].CLast[19];
+ Tf.C55 = extTracks[jTr].CLast[20];
+
+
if (fabs(Tf.t[0] - Tb.t[0]) > 3 * sqrt(Tf.C55[0] + Tb.C55[0])) continue;
- // if (fabs (Tf.time[0] - Tb.time[0]) > 500000) continue;
unsigned short stam;
- fParameters.GetStation(staf).fieldSlice.GetFieldValue(Tf.x, Tf.y, fBf);
- fParameters.GetStation(stab).fieldSlice.GetFieldValue(Tb.x, Tb.y, fBb);
+ frAlgo.GetParameters()->GetStation(staf).fieldSlice.GetFieldValue(Tf.x, Tf.y, fBf);
+ frAlgo.GetParameters()->GetStation(stab).fieldSlice.GetFieldValue(Tb.x, Tb.y, fBb);
unsigned short dist = firstStation[iTr] - lastStation[jTr];
@@ -356,10 +169,10 @@ void L1Algo::CAMergeClones()
else
stam = staf - 1;
- fvec zm = fParameters.GetStation(stam).z;
+ fvec zm = frAlgo.GetParameters()->GetStation(stam).z;
fvec xm = 0.5 * (Tf.x + Tf.tx * (zm - Tf.z) + Tb.x + Tb.tx * (zm - Tb.z));
fvec ym = 0.5 * (Tf.y + Tf.ty * (zm - Tf.z) + Tb.y + Tb.ty * (zm - Tb.z));
- fParameters.GetStation(stam).fieldSlice.GetFieldValue(xm, ym, fBm);
+ frAlgo.GetParameters()->GetStation(stam).fieldSlice.GetFieldValue(xm, ym, fBm);
fld.Set(fBb, Tb.z, fBm, zm, fBf, Tf.z);
fvec zMiddle = 0.5 * (Tb.z + Tf.z);
@@ -395,39 +208,211 @@ void L1Algo::CAMergeClones()
for (int iTr = 0; iTr < nTracks; iTr++) {
if (isStored[iTr]) continue;
- fMergerTracksNew.push_back(fTracks[iTr]);
+ fTracksNew.push_back(extTracks[iTr]);
if (!isDownstreamNeighbour[iTr]) {
for (L1HitIndex_t HI = firstHit[iTr]; HI <= lastHit[iTr]; HI++) {
- fMergerRecoHitsNew.push_back(fRecoHits[HI]);
+ fRecoHitsNew.push_back(extRecoHits[HI]);
}
}
if (neighbour[iTr] < kNoNeighbour) {
isStored[neighbour[iTr]] = true;
- fMergerTracksNew.back().NHits += fTracks[neighbour[iTr]].NHits;
+ fTracksNew.back().NHits += extTracks[neighbour[iTr]].NHits;
for (L1HitIndex_t HI = firstHit[neighbour[iTr]]; HI <= lastHit[neighbour[iTr]]; HI++)
- fMergerRecoHitsNew.push_back(fRecoHits[HI]);
+ fRecoHitsNew.push_back(extRecoHits[HI]);
}
if (isDownstreamNeighbour[iTr]) {
for (L1HitIndex_t HI = firstHit[iTr]; HI <= lastHit[iTr]; HI++) {
- fMergerRecoHitsNew.push_back(fRecoHits[HI]);
+ fRecoHitsNew.push_back(extRecoHits[HI]);
}
}
}
- fTracks.reset(fMergerTracksNew.size());
- for (unsigned int iTr = 0; iTr < fMergerTracksNew.size(); iTr++) {
- fTracks[iTr] = fMergerTracksNew[iTr];
+ // Save the merging results into external vectors
+ extTracks = std::move(fTracksNew);
+ extRecoHits = std::move(fRecoHitsNew);
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void L1ClonesMerger::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++) {
+ if (W) W[i] = C[i];
+ S[i] = C[i] + V[i];
+ }
+
+ InvertCholesky(S);
+
+ fvec dzeta[5];
+
+ for (int i = 0; i < 5; i++)
+ dzeta[i] = m[i] - r[i];
+
+ if (W && R) {
+ for (int i = 0; i < 5; i++)
+ R[i] = r[i];
+
+ fvec K[5][5];
+ MultiplySS(C, S, K);
+
+ fvec KC[15];
+ MultiplyMS(K, C, KC);
+ for (int i = 0; i < 15; i++)
+ W[i] -= KC[i];
+
+ fvec kd;
+ for (int i = 0; i < 5; i++) {
+ kd = 0.f;
+ for (int j = 0; j < 5; j++)
+ kd += K[i][j] * dzeta[j];
+ R[i] += kd;
+ }
+ }
+ if (chi2) {
+ fvec S_dzeta[5];
+ MultiplySR(S, dzeta, S_dzeta);
+ *chi2 = dzeta[0] * S_dzeta[0] + dzeta[1] * S_dzeta[1] + dzeta[2] * S_dzeta[2] + dzeta[3] * S_dzeta[3]
+ + dzeta[4] * S_dzeta[4];
+ }
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void L1ClonesMerger::InvertCholesky(fvec a[15])
+{
+ fvec d[5], uud, u[5][5];
+ for (int i = 0; i < 5; i++) {
+ d[i] = 0.f;
+ for (int j = 0; j < 5; j++)
+ u[i][j] = 0.f;
}
- assert(fRecoHits.size() == fMergerRecoHitsNew.size());
- fRecoHits.reset(fMergerRecoHitsNew.size());
- for (L1HitIndex_t iHi = 0; iHi < fMergerRecoHitsNew.size(); iHi++) {
- fRecoHits[iHi] = fMergerRecoHitsNew[iHi];
+ for (int i = 0; i < 5; i++) {
+ uud = 0.f;
+ for (int j = 0; j < i; j++)
+ uud += u[j][i] * u[j][i] * d[j];
+ uud = a[i * (i + 3) / 2] - uud;
+
+ fvec smallval = 1.e-12;
+ fvec initialised = fvec(uud < smallval);
+ uud = ((!initialised) & uud) + (smallval & initialised);
+
+ d[i] = uud / fabs(uud);
+ u[i][i] = sqrt(fabs(uud));
+
+ for (int j = i + 1; j < 5; j++) {
+ uud = 0.f;
+ for (int k = 0; k < i; k++)
+ uud += u[k][i] * u[k][j] * d[k];
+ uud = a[j * (j + 1) / 2 + i] /*a[i][j]*/ - uud;
+ u[i][j] = d[i] / u[i][i] * uud;
+ }
}
- NHitsIsecAll = fMergerRecoHitsNew.size();
+ fvec u1[5];
- //std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!! new "<<vTracksNew.size()<<" old "<< vTracks.size()<<std::endl;
+ for (int i = 0; i < 5; i++) {
+ u1[i] = u[i][i];
+ u[i][i] = 1.f / u[i][i];
+ }
+ for (int i = 0; i < 4; i++) {
+ u[i][i + 1] = -u[i][i + 1] * u[i][i] * u[i + 1][i + 1];
+ }
+ for (int i = 0; i < 3; i++) {
+ u[i][i + 2] = u[i][i + 1] * u1[i + 1] * u[i + 1][i + 2] - u[i][i + 2] * u[i][i] * u[i + 2][i + 2];
+ }
+ for (int i = 0; i < 2; i++) {
+ u[i][i + 3] = u[i][i + 2] * u1[i + 2] * u[i + 2][i + 3] - u[i][i + 3] * u[i][i] * u[i + 3][i + 3];
+ u[i][i + 3] -= u[i][i + 1] * u1[i + 1] * (u[i + 1][i + 2] * u1[i + 2] * u[i + 2][i + 3] - u[i + 1][i + 3]);
+ }
+ u[0][4] = u[0][2] * u1[2] * u[2][4] - u[0][4] * u[0][0] * u[4][4];
+ u[0][4] += u[0][1] * u1[1] * (u[1][4] - u[1][3] * u1[3] * u[3][4] - u[1][2] * u1[2] * u[2][4]);
+ u[0][4] += u[3][4] * u1[3] * (u[0][3] - u1[2] * u[2][3] * (u[0][2] - u[0][1] * u1[1] * u[1][2]));
+
+ for (int i = 0; i < 5; i++)
+ a[i + 10] = u[i][4] * d[4] * u[4][4];
+ for (int i = 0; i < 4; i++)
+ a[i + 6] = u[i][3] * u[3][3] * d[3] + u[i][4] * u[3][4] * d[4];
+ for (int i = 0; i < 3; i++)
+ a[i + 3] = u[i][2] * u[2][2] * d[2] + u[i][3] * u[2][3] * d[3] + u[i][4] * u[2][4] * d[4];
+ for (int i = 0; i < 2; i++)
+ a[i + 1] =
+ u[i][1] * u[1][1] * d[1] + u[i][2] * u[1][2] * d[2] + u[i][3] * u[1][3] * d[3] + u[i][4] * u[1][4] * d[4];
+ a[0] = u[0][0] * u[0][0] * d[0] + u[0][1] * u[0][1] * d[1] + u[0][2] * u[0][2] * d[2] + u[0][3] * u[0][3] * d[3]
+ + u[0][4] * u[0][4] * d[4];
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void L1ClonesMerger::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];
+
+ K[1] = C[1][0] * V[0] + C[1][1] * V[1] + C[1][2] * V[3] + C[1][3] * V[6] + C[1][4] * V[10];
+ K[2] = C[1][0] * V[1] + C[1][1] * V[2] + C[1][2] * V[4] + C[1][3] * V[7] + C[1][4] * V[11];
+
+ K[3] = C[2][0] * V[0] + C[2][1] * V[1] + C[2][2] * V[3] + C[2][3] * V[6] + C[2][4] * V[10];
+ K[4] = C[2][0] * V[1] + C[2][1] * V[2] + C[2][2] * V[4] + C[2][3] * V[7] + C[2][4] * V[11];
+ K[5] = C[2][0] * V[3] + C[2][1] * V[4] + C[2][2] * V[5] + C[2][3] * V[8] + C[2][4] * V[12];
+
+ K[6] = C[3][0] * V[0] + C[3][1] * V[1] + C[3][2] * V[3] + C[3][3] * V[6] + C[3][4] * V[10];
+ K[7] = C[3][0] * V[1] + C[3][1] * V[2] + C[3][2] * V[4] + C[3][3] * V[7] + C[3][4] * V[11];
+ K[8] = C[3][0] * V[3] + C[3][1] * V[4] + C[3][2] * V[5] + C[3][3] * V[8] + C[3][4] * V[12];
+ K[9] = C[3][0] * V[6] + C[3][1] * V[7] + C[3][2] * V[8] + C[3][3] * V[9] + C[3][4] * V[13];
+
+ K[10] = C[4][0] * V[0] + C[4][1] * V[1] + C[4][2] * V[3] + C[4][3] * V[6] + C[4][4] * V[10];
+ K[11] = C[4][0] * V[1] + C[4][1] * V[2] + C[4][2] * V[4] + C[4][3] * V[7] + C[4][4] * V[11];
+ K[12] = C[4][0] * V[3] + C[4][1] * V[4] + C[4][2] * V[5] + C[4][3] * V[8] + C[4][4] * V[12];
+ K[13] = C[4][0] * V[6] + C[4][1] * V[7] + C[4][2] * V[8] + C[4][3] * V[9] + C[4][4] * V[13];
+ K[14] = C[4][0] * V[10] + C[4][1] * V[11] + C[4][2] * V[12] + C[4][3] * V[13] + C[4][4] * V[14];
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void L1ClonesMerger::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];
+ r_out[2] = r_in[0] * C[3] + r_in[1] * C[4] + r_in[2] * C[5] + r_in[3] * C[8] + r_in[4] * C[12];
+ r_out[3] = r_in[0] * C[6] + r_in[1] * C[7] + r_in[2] * C[8] + r_in[3] * C[9] + r_in[4] * C[13];
+ r_out[4] = r_in[0] * C[10] + r_in[1] * C[11] + r_in[2] * C[12] + r_in[3] * C[13] + r_in[4] * C[14];
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+void L1ClonesMerger::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];
+ K[0][2] = C[0] * V[3] + C[1] * V[4] + C[3] * V[5] + C[6] * V[8] + C[10] * V[12];
+ K[0][3] = C[0] * V[6] + C[1] * V[7] + C[3] * V[8] + C[6] * V[9] + C[10] * V[13];
+ K[0][4] = C[0] * V[10] + C[1] * V[11] + C[3] * V[12] + C[6] * V[13] + C[10] * V[14];
+
+ K[1][0] = C[1] * V[0] + C[2] * V[1] + C[4] * V[3] + C[7] * V[6] + C[11] * V[10];
+ K[1][1] = C[1] * V[1] + C[2] * V[2] + C[4] * V[4] + C[7] * V[7] + C[11] * V[11];
+ K[1][2] = C[1] * V[3] + C[2] * V[4] + C[4] * V[5] + C[7] * V[8] + C[11] * V[12];
+ K[1][3] = C[1] * V[6] + C[2] * V[7] + C[4] * V[8] + C[7] * V[9] + C[11] * V[13];
+ K[1][4] = C[1] * V[10] + C[2] * V[11] + C[4] * V[12] + C[7] * V[13] + C[11] * V[14];
+
+ K[2][0] = C[3] * V[0] + C[4] * V[1] + C[5] * V[3] + C[8] * V[6] + C[12] * V[10];
+ K[2][1] = C[3] * V[1] + C[4] * V[2] + C[5] * V[4] + C[8] * V[7] + C[12] * V[11];
+ K[2][2] = C[3] * V[3] + C[4] * V[4] + C[5] * V[5] + C[8] * V[8] + C[12] * V[12];
+ K[2][3] = C[3] * V[6] + C[4] * V[7] + C[5] * V[8] + C[8] * V[9] + C[12] * V[13];
+ K[2][4] = C[3] * V[10] + C[4] * V[11] + C[5] * V[12] + C[8] * V[13] + C[12] * V[14];
+
+ K[3][0] = C[6] * V[0] + C[7] * V[1] + C[8] * V[3] + C[9] * V[6] + C[13] * V[10];
+ K[3][1] = C[6] * V[1] + C[7] * V[2] + C[8] * V[4] + C[9] * V[7] + C[13] * V[11];
+ K[3][2] = C[6] * V[3] + C[7] * V[4] + C[8] * V[5] + C[9] * V[8] + C[13] * V[12];
+ K[3][3] = C[6] * V[6] + C[7] * V[7] + C[8] * V[8] + C[9] * V[9] + C[13] * V[13];
+ K[3][4] = C[6] * V[10] + C[7] * V[11] + C[8] * V[12] + C[9] * V[13] + C[13] * V[14];
+
+ K[4][0] = C[10] * V[0] + C[11] * V[1] + C[12] * V[3] + C[13] * V[6] + C[14] * V[10];
+ K[4][1] = C[10] * V[1] + C[11] * V[2] + C[12] * V[4] + C[13] * V[7] + C[14] * V[11];
+ K[4][2] = C[10] * V[3] + C[11] * V[4] + C[12] * V[5] + C[13] * V[8] + C[14] * V[12];
+ K[4][3] = C[10] * V[6] + C[11] * V[7] + C[12] * V[8] + C[13] * V[9] + C[14] * V[13];
+ K[4][4] = C[10] * V[10] + C[11] * V[11] + C[12] * V[12] + C[13] * V[13] + C[14] * V[14];
}
diff --git a/reco/L1/L1Algo/L1ClonesMerger.h b/reco/L1/L1Algo/L1ClonesMerger.h
new file mode 100644
index 0000000000..9b76aefdfa
--- /dev/null
+++ b/reco/L1/L1Algo/L1ClonesMerger.h
@@ -0,0 +1,119 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer], Maksym Zyzak */
+
+/// \file L1ClonesMerger.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 L1ClonesMerger_h
+#define L1ClonesMerger_h 1
+
+#include "L1Hit.h" // For L1HitIndex_t
+#include "L1Vector.h"
+#include "vectors/P4_F32vec4.h"
+
+class L1Track;
+class L1Algo;
+
+/// Class implements a clones merger algorithm for the CA (L1) track finder
+///
+class L1ClonesMerger {
+public:
+ /// Default constructor
+ L1ClonesMerger(const L1Algo& algo);
+
+ /// Destructor
+ ~L1ClonesMerger();
+
+ /// Copy constructor
+ L1ClonesMerger(const L1ClonesMerger&) = delete;
+
+ /// Move constructor
+ L1ClonesMerger(L1ClonesMerger&&) = delete;
+
+ /// Copy assignment operator
+ L1ClonesMerger& operator=(const L1ClonesMerger&) = delete;
+
+ /// Move assignment operator
+ L1ClonesMerger& operator=(L1ClonesMerger&&) = 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(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_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
+ L1Vector<unsigned short> fTrackFirstStation {"L1ClonesMerger::fTrackFirstStation"};
+
+ /// Last station of a track
+ L1Vector<unsigned short> fTrackLastStation {"L1ClonesMerger::fTrackLastStation"};
+
+ /// Index of the first hit of a track
+ L1Vector<L1HitIndex_t> fTrackFirstHit {"L1ClonesMerger::fTrackFirstHit"};
+
+ /// Index of the last hit of a track
+ L1Vector<L1HitIndex_t> fTrackLastHit {"L1ClonesMerger::fTrackLastHit"};
+
+ /// Index (TODO:??) of a track that can be merge with the given track
+ L1Vector<unsigned short> fTrackNeighbour {"L1ClonesMerger::fTrackNeighbour"};
+
+ /// Chi2 value of the track merging procedure
+ L1Vector<float> fTrackChi2 {"L1ClonesMerger::fTrackChi2"};
+
+ /// Flag: is the given track already stored to the output
+ L1Vector<char> fTrackIsStored {"L1ClonesMerger::fTrackIsStored"};
+
+ /// Flag: is the track a downstream neighbour of another track
+ L1Vector<char> fTrackIsDownstreamNeighbour {"L1ClonesMerger::fTrackIsDownstreamNeighbour"};
+
+ L1Vector<L1Track> fTracksNew {"L1CAMergerClones::fTracksNew"}; ///< vector of tracks after the merge
+
+ L1Vector<L1HitIndex_t> fRecoHitsNew {"L1CAMergerClones::fRecoHitsNew"}; ///< vector of track hits after the merge
+
+ const L1Algo& frAlgo; ///< Reference to the main track finder algorithm class
+};
+
+#endif // L1ClonesMerger_h
diff --git a/reco/L1/L1Algo/L1Track.h b/reco/L1/L1Algo/L1Track.h
index 584d4b2338..21c4bc93af 100644
--- a/reco/L1/L1Algo/L1Track.h
+++ b/reco/L1/L1Algo/L1Track.h
@@ -31,24 +31,24 @@ class L1Track {
public:
static constexpr float kNaN {std::numeric_limits<float>::signaling_NaN()};
- unsigned char NHits; ///> Number of hits in track
- unsigned char n; ///> TODO: ??
- float Momentum {kNaN}; ///> TODO: ??
- float fTrackTime {kNaN}; ///> Track time
- fscal TFirst[7] {kNaN}; ///> Track parameters on the first station
- fscal CFirst[21] {kNaN}; ///> Track parameter covariation matrix elemenst on the first station
- fscal TLast[7] {kNaN}; ///> Track parameters on the last station
- fscal CLast[21] {kNaN}; ///> Track parameter covatiation matrix elements on the second station
- fscal Tpv[7] {kNaN}; ///> Track parameters in the primary vertex
- fscal Cpv[21] {kNaN}; ///> Track parameter covariation matrix elements in the primary vertex
- fscal chi2 {kNaN}; ///> Track fit chi-square value
- short int NDF; ///> Track fit NDF value
+ unsigned char NHits; ///< Number of hits in track
+ unsigned char n; ///< TODO: ??
+ float Momentum {kNaN}; ///< TODO: ??
+ float fTrackTime {kNaN}; ///< Track time
+ fscal TFirst[7] {kNaN}; ///< Track parameters on the first station
+ fscal CFirst[21] {kNaN}; ///< Track parameter covariation matrix elemenst on the first station
+ fscal TLast[7] {kNaN}; ///< Track parameters on the last station
+ fscal CLast[21] {kNaN}; ///< Track parameter covatiation matrix elements on the second station
+ fscal Tpv[7] {kNaN}; ///< Track parameters in the primary vertex
+ fscal Cpv[21] {kNaN}; ///< Track parameter covariation matrix elements in the primary vertex
+ fscal chi2 {kNaN}; ///< Track fit chi-square value
+ short int NDF; ///< Track fit NDF value
// TODO: shouldn't it be the L1HitIndex_t type instead of int? (S.Zharko)
- int FirstHitIndex; ///> Track first hit index in the hits vector
- int LastHitIndex; ///> Track last hit index in the hits vector
- int index; ///> Index of the track
- int ista; ///> TODO: ??
+ int FirstHitIndex; ///< Track first hit index in the hits vector
+ int LastHitIndex; ///< Track last hit index in the hits vector
+ int index; ///< Index of the track
+ int ista; ///< TODO: ??
/// Provides comparison of two L1Track objects
/// If two tracks have different number of hits, the smallest track has the largest number of hits.
diff --git a/reco/L1/L1Algo/L1Utils.h b/reco/L1/L1Algo/L1Utils.h
index 5f499b28ba..fdd3ba8f9e 100644
--- a/reco/L1/L1Algo/L1Utils.h
+++ b/reco/L1/L1Algo/L1Utils.h
@@ -10,13 +10,16 @@
#ifndef L1Utils_h
#define L1Utils_h 1
-#include <type_traits>
+#include <FairLogger.h>
+
+#include <chrono>
#include <iomanip>
#include <limits>
#include <map>
#include <sstream>
#include <string>
+#include <type_traits>
#include <unordered_map>
#include <cmath>
@@ -28,24 +31,23 @@
#endif
/// Class contains some utility functions for L1Algo
-struct L1Utils {
-
+namespace L1Utils
+{
/// NaN value for float
- static constexpr float kNaN {std::numeric_limits<float>::signaling_NaN()};
+ constexpr float kNaN {std::numeric_limits<float>::signaling_NaN()};
/// Comparison method for floats
/// \param lhs Left floating point to compare
/// \param rhs Right floating point to compare
/// \return Comparison result: true - equals within epsilon
template<typename T, typename std::enable_if<std::is_floating_point<T>::value, T>::type* = nullptr>
- static bool CmpFloats(T lhs, T rhs)
+ [[gnu::always_inline]] inline bool CmpFloats(T lhs, T rhs)
{
return fabs(lhs - rhs) < 2. * std::numeric_limits<T>::epsilon() * fabs(lhs + rhs)
|| fabs(lhs - rhs) < std::numeric_limits<T>::min();
}
-
- static __attribute__((always_inline)) void CheckSimdVectorEquality(fvec v, const char* name)
+ [[gnu::always_inline]] inline void CheckSimdVectorEquality(fvec v, const char* name)
{
if (!v.IsHorizontallyEqual()) {
std::stringstream msg;
@@ -63,50 +65,24 @@ struct L1Utils {
}
};
- /// Template function, which sets a value to an element of the map with a particular key
- /// \param key Key of the element to be modified
- /// \param value New value of the element under the selected key
- /// \param aMap A reference to the map, which element is to be modified
- template<class Key, class T, class Hash = std::hash<Key>>
- static void SetSingleValueToMap(Key key, T value, std::unordered_map<Key, T, Hash>& aMap)
- {
- aMap[key] = value;
- }
-
- /// Template function, which sets a value to ALL elements of the map
- /// \param value New value of the element under the selected key
- /// \param aMap A reference to the map, which element is to be modified
- template<class Key, class T, class Hash = std::hash<Key>>
- static void SetSameValueToMap(T value, std::unordered_map<Key, T, Hash>& aMap)
- {
- for (auto it = aMap.begin(); it != aMap.end(); ++it) {
- it->second = value;
- }
- }
+ /// A time profiler for measuring performace of scopes
+ class TimeProfiler {
+ public:
+ /// Constructor
+ TimeProfiler(const char* scopeName) : fScopeName(scopeName), fStart(std::chrono::high_resolution_clock::now()) {}
- /// Template function, which resets the elements of one map with the values defined in another map
- /// \param inMap A constant reference to the map containing new parameters
- /// \param aMap A reference to the map, which is going to be modified
- template<class Key, class T, class Hash = std::hash<Key>>
- static void SetMappedValuesToMap(const std::unordered_map<Key, T, Hash>& inMap,
- std::unordered_map<Key, T, Hash>& aMap)
- {
- for (auto it = aMap.begin(); it != aMap.end(); ++it) {
- if (inMap.find(it->first) != inMap.end()) { it->second = inMap.at(it->first); }
+ /// Destructor
+ ~TimeProfiler()
+ {
+ auto stop = std::chrono::high_resolution_clock::now();
+ auto time = std::chrono::duration_cast<std::chrono::nanoseconds>(stop - fStart).count();
+ LOG(info) << "Scope \033[1;32m" << fScopeName << "\033[0m was executed in " << time << " ns";
}
- }
- /// Template function to represent mapped contents into std::string
- /// NOTE: operator<< must be defined for value of the map
- template<class Key, class T, class Hash = std::hash<Key>>
- static std::string RepresentMapWithString(const std::unordered_map<Key, T, Hash>& aMap, int entryWidth = 15)
- {
- std::stringstream token;
- for (auto it = aMap.begin(); it != aMap.end(); ++it) {
- token << std::setw(entryWidth) << std::setfill(' ') << it->second << ' ';
- }
- return token.str();
- }
-};
+ private:
+ const char* fScopeName {};
+ const std::chrono::high_resolution_clock::time_point fStart {};
+ };
+} // namespace L1Utils
#endif // L1Utils_h
--
GitLab