From e820c3edb57e8753028d41063792db8c452a4a0b Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Thu, 31 Aug 2023 19:44:35 +0000
Subject: [PATCH] CA: de-SIMDIze and clean up the triplet constructor
---
reco/L1/CMakeLists.txt | 2 +-
reco/L1/L1Algo/L1Algo.cxx | 31 -
reco/L1/L1Algo/L1Algo.h | 114 +-
reco/L1/L1Algo/L1BaseStationInfo.cxx | 3 +-
reco/L1/L1Algo/L1BaseStationInfo.h | 1 +
reco/L1/L1Algo/L1CaTrackFinderSlice.cxx | 1526 +----------------------
reco/L1/L1Algo/L1Constants.h | 2 -
reco/L1/L1Algo/L1Field.cxx | 15 +-
reco/L1/L1Algo/L1Field.h | 11 +
reco/L1/L1Algo/L1Portion.h | 166 ---
reco/L1/L1Algo/L1Triplet.h | 3 +-
reco/L1/L1Algo/L1TripletConstructor.cxx | 690 ++++++++++
reco/L1/L1Algo/L1TripletConstructor.h | 116 ++
13 files changed, 894 insertions(+), 1786 deletions(-)
delete mode 100644 reco/L1/L1Algo/L1Portion.h
create mode 100644 reco/L1/L1Algo/L1TripletConstructor.cxx
create mode 100644 reco/L1/L1Algo/L1TripletConstructor.h
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index a43c86b078..b35b111739 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -34,6 +34,7 @@ set(SRCS
L1Algo/utils/CaAlgoRandom.cxx
L1Algo/L1Algo.cxx
L1Algo/L1TrackPar.cxx
+ L1Algo/L1TripletConstructor.cxx
L1Algo/L1CaTrackFinder.cxx
L1Algo/L1CaTrackFinderSlice.cxx
L1Algo/L1BranchExtender.cxx
@@ -205,7 +206,6 @@ install(FILES CbmL1Counters.h
L1Algo/L1Branch.h
L1Algo/L1HitPoint.h
L1Algo/L1Hit.h
- L1Algo/L1Portion.h
L1Algo/L1Triplet.h
L1Algo/L1Event.h
L1Algo/L1EventMatch.h
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index fba7db9add..b0a6dcff61 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -34,31 +34,6 @@ void L1Algo::SetNThreads(unsigned int n)
fTrackCandidates[i].SetName(std::stringstream() << "L1Algo::fTrackCandidates[" << i << "]");
-
- fTripletPar[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-
- fTripletHitsL[i].SetName(std::stringstream() << "L1Algo::fhitsl_3[" << i << "]");
- fTripletHitsM[i].SetName(std::stringstream() << "L1Algo::fhitsm_3[" << i << "]");
- fTripletHitsR[i].SetName(std::stringstream() << "L1Algo::fhitsr_3[" << i << "]");
-
- fTripletHitsL[i].reserve(L1Constants::size::kMaxPortionTriplets);
- fTripletHitsM[i].reserve(L1Constants::size::kMaxPortionTriplets);
- fTripletHitsR[i].reserve(L1Constants::size::kMaxPortionTriplets);
-
- fTripletHitR_Ufront[i].reserve(L1Constants::size::kMaxPortionTripletsP);
- fTripletHitR_Uback[i].reserve(L1Constants::size::kMaxPortionTripletsP);
- fTripletHitR_Z[i].reserve(L1Constants::size::kMaxPortionTripletsP);
- fTripletHitR_Time[i].reserve(L1Constants::size::kMaxPortionTripletsP);
- fTripletHitR_TimeErr[i].reserve(L1Constants::size::kMaxPortionTripletsP);
- fTripletHitR_dUfront[i].reserve(L1Constants::size::kMaxPortionTripletsP);
- fTripletHitR_dUback[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-
- fTripletHitR_X[i].reserve(L1Constants::size::kMaxPortionTripletsP);
- fTripletHitR_Y[i].reserve(L1Constants::size::kMaxPortionTripletsP);
- fTripletHitR_dX[i].reserve(L1Constants::size::kMaxPortionTripletsP);
- fTripletHitR_dXY[i].reserve(L1Constants::size::kMaxPortionTripletsP);
- fTripletHitR_dY[i].reserve(L1Constants::size::kMaxPortionTripletsP);
-
for (unsigned int j = 0; j < L1Constants::size::kMaxNstations; j++) {
fTriplets[j][i].SetName(std::stringstream() << "L1Algo::fTriplets[" << i << "][" << j << "]");
}
@@ -128,12 +103,6 @@ void L1Algo::ResetSliceData()
fSliceRecoHits.clear();
fSliceRecoHits.reserve(2 * nHits);
- for (int iS = 0; iS < L1Constants::size::kMaxNstations; iS++) {
- int nHitsStation = fSliceHitIds[iS].size();
- fSingletPortionSize[iS].clear();
- fSingletPortionSize[iS].reserve(2 * nHitsStation);
- }
-
for (int iT = 0; iT < fNThreads; iT++) {
fSliceRecoTracks_thread[iT].clear();
fSliceRecoTracks_thread[iT].reserve(nHits / 10);
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index a44de8a34c..a641957c1f 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -37,7 +37,6 @@ class L1AlgoDraw;
#include "L1HitPoint.h"
#include "L1InputData.h"
#include "L1Parameters.h"
-#include "L1Portion.h"
#include "L1Station.h"
#include "L1Track.h"
#include "L1TrackPar.h"
@@ -88,8 +87,7 @@ public:
// *************************
friend class CbmL1; /// TODO: Remove this dependency
- friend class ParalleledDup;
- friend class ParalleledTrip;
+ friend class L1TripletConstructor;
#ifdef DRAW
friend class L1AlgoDraw;
#endif
@@ -261,82 +259,6 @@ public:
flag = iStation * 4 + (flag % 4);
}
- /// Prepare the portion of left hits data
- void findSingletsStep0( // input
- Tindex start_lh, Tindex n1_l, L1HitPoint* Hits_l,
- // output
- L1HitIndex_t* hitsl, fvec* x_l, fvec* y_l, fvec* z_l, fvec* t_l, fvec* dx_l, fvec* dxy_l, fvec* dy_l, fvec* dt2_l);
-
- /// Get the field approximation. Add the target to parameters estimation. Propagate to middle station.
- void findSingletsStep1( // input
- int istal, int istam, Tindex n1_V, fvec* x_l, fvec* y_l, fvec* zPos_l, fvec* t_l, fvec* dx_l, fvec* dxy_l,
- fvec* dy_l, fvec* dt2_l,
- // output
- L1TrackPar* T_1, L1FieldRegion* fld_1);
-
- /// Find the doublets. Reformat data in the portion of doublets.
- void findDoubletsStep0( // input
- Tindex n1, const L1Station& stal, L1HitPoint* vHits_l, const L1Station& stam, L1HitPoint* vHits_m, L1TrackPar* T_1,
- L1HitIndex_t* hitsl_1,
- // output
- Tindex& n2, L1Vector<L1HitIndex_t>& i1_2,
-#ifdef DOUB_PERFORMANCE
- L1Vector<L1HitIndex_t>& hitsl_2,
-#endif // DOUB_PERFORMANCE
- L1Vector<L1HitIndex_t>& hitsm_2);
-
- /// Add the middle hits to parameters estimation. Propagate to right station.
- /// Find the triplets (right hit). Reformat data in the portion of triplets.
- void findTripletsStep0( // input
- L1HitPoint* vHits_r, int istal, int istam, int istar, L1HitPoint* vHits_m, L1TrackPar* T_1, L1FieldRegion* fld_1,
- L1HitIndex_t* hitsl_1,
-
- Tindex n2, L1Vector<L1HitIndex_t>& hitsm_2, L1Vector<L1HitIndex_t>& i1_2,
-
- // output
- Tindex& n3, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3, L1Vector<L1HitIndex_t>& hitsm_3,
- L1Vector<L1HitIndex_t>& hitsr_3,
-
- L1Vector<fvec>& x_3, L1Vector<fvec>& y_3, L1Vector<fvec>& z_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dx_3,
- L1Vector<fvec>& dxy_3, L1Vector<fvec>& dy_3, L1Vector<fvec>& dt2_3);
-
- /// Add the right hits to parameters estimation.
- void findTripletsStep1( // input
- Tindex n3_V, const L1Station& star, L1Vector<fvec>& x_3, L1Vector<fvec>& y_3, L1Vector<fvec>& z_3,
- L1Vector<fvec>& t_3, L1Vector<fvec>& dx_3, L1Vector<fvec>& dxy_3, L1Vector<fvec>& dy_3, L1Vector<fvec>& dt2_3,
- // output
- L1Vector<L1TrackPar>& T_3);
-
- /// Refit Triplets.
- void findTripletsStep2( // input
- Tindex n3, int istal, int istam, int istar, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3,
- L1Vector<L1HitIndex_t>& hitsm_3, L1Vector<L1HitIndex_t>& hitsr_3, int nIterations = 0);
-
- /// Select triplets. Save them into vTriplets.
- void findTripletsStep3( // input
- Tindex n3, int istal, int istam, int istar, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3,
- L1Vector<L1HitIndex_t>& hitsm_3, L1Vector<L1HitIndex_t>& hitsr_3);
-
-
- /// Find neighbours of triplets. Calculate level of triplets.
- void f5( // input
- // output
- int* nlevel);
-
-
- /// Find doublets on station
- void CreatePortionOfDoublets( // input
- int istal, int istam, Tindex iSingletPortion, Tindex singletPortionSize,
- // output
- L1TrackPar* T_1, L1FieldRegion* fld_1, L1HitIndex_t* hitsl_1,
- //
- Tindex& n_2, L1Vector<L1HitIndex_t>& i1_2, L1Vector<L1HitIndex_t>& hitsm_2);
-
- /// Find triplets on station
- void CreatePortionOfTriplets( // input
- int istal, int istam, int istar, L1TrackPar* T_1, L1FieldRegion* fld_1, L1HitIndex_t* hitsl_1, Tindex& n_2,
- L1Vector<L1HitIndex_t>& i1_2, L1Vector<L1HitIndex_t>& hitsm_2);
-
#ifdef DRAW
L1AlgoDraw* draw {nullptr};
void DrawRecoTracksTime(const L1Vector<CbmL1Track>& tracks);
@@ -432,10 +354,6 @@ public:
/// The candidates may share any amount of hits.
L1Vector<L1Branch> fTrackCandidates[L1Constants::size::kMaxNthreads] {"L1Algo::fTrackCandidates"};
- L1Vector<Tindex> fSingletPortionSize[L1Constants::size::kMaxNstations] {
- "L1Algo::fSingletPortionSize"}; ///< Number of doublets in a portion
-
-
///< indices of the sub-slice hits
L1Vector<L1HitIndex_t> fSliceHitIds[L1Constants::size::kMaxNstations] {"L1Algo::fSliceHitIds"};
@@ -476,36 +394,6 @@ public:
L1Vector<int> fHitNtriplets {"L1Algo::fHitNtriplets"}; /// link hit ->n triplets { hit, *, *}
- // fvec u_front[Portion/fvecLen], u_back[Portion/fvecLen];
- // fvec zPos[Portion/fvecLen];
- // fvec fHitTime[Portion/fvecLen];
-
- L1Vector<L1TrackPar> fTripletPar[L1Constants::size::kMaxNthreads] {"fTripletPar"};
-
- L1Vector<L1HitIndex_t> fTripletHitsL[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitsL"};
- L1Vector<L1HitIndex_t> fTripletHitsM[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitsM"};
- L1Vector<L1HitIndex_t> fTripletHitsR[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitsR"};
-
- L1Vector<fvec> fTripletHitR_Ufront[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_Ufront"};
- L1Vector<fvec> fTripletHitR_Uback[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_Uback"};
- L1Vector<fvec> fTripletHitR_Z[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_Z"};
- L1Vector<fvec> fTripletHitR_Time[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_Time"};
- L1Vector<fvec> fTripletHitR_TimeErr[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_TimeErr"};
- L1Vector<fvec> fTripletHitR_dUfront[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dUfront"};
- L1Vector<fvec> fTripletHitR_dUback[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dUback"};
-
- L1Vector<fvec> fTripletHitR_X[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_X"};
- L1Vector<fvec> fTripletHitR_Y[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_Y"};
- L1Vector<fvec> fTripletHitR_dX[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dX"};
- L1Vector<fvec> fTripletHitR_dXY[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dXY"};
- L1Vector<fvec> fTripletHitR_dY[L1Constants::size::kMaxNthreads] {"L1Algo::fTripletHitR_dY"};
-
-
- // Tindex NHits_l[L1Constants::size::kMaxNstations];
- // Tindex NHits_l_P[L1Constants::size::kMaxNstations];
- /// ----- Output data -----
-
-
private:
L1CloneMerger fCloneMerger {*this}; ///< Object of L1Algo clones merger algorithm
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.cxx b/reco/L1/L1Algo/L1BaseStationInfo.cxx
index 044ea49e2c..c5b214ecb8 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.cxx
+++ b/reco/L1/L1Algo/L1BaseStationInfo.cxx
@@ -205,6 +205,7 @@ void L1BaseStationInfo::SetFieldFunction(
fL1Station.fieldSlice.cy[j] = A[j][N + 1] / A[j][j];
fL1Station.fieldSlice.cz[j] = A[j][N + 2] / A[j][j];
}
+ fL1Station.fieldSlice.z = fZref;
fInitController.SetFlag(EInitKey::kFieldSlice);
}
@@ -285,7 +286,7 @@ void L1BaseStationInfo::SetXmax(double aSize)
//
void L1BaseStationInfo::SetYmax(double aSize)
{
- fYmax = aSize;
+ fYmax = aSize;
fL1Station.Ymax = aSize;
fInitController.SetFlag(EInitKey::kYmax);
}
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.h b/reco/L1/L1Algo/L1BaseStationInfo.h
index 60df5ababc..a23cbb5e7a 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.h
+++ b/reco/L1/L1Algo/L1BaseStationInfo.h
@@ -99,6 +99,7 @@ public:
/// Gets array of the coefficients for the field z-component approximation
const fvec* GetFieldSliceCz() const { return fL1Station.fieldSlice.cz; }
+
/// Gets field status: 0 - station is outside the field, 1 - station is inside the field
int GetFieldStatus() const { return fL1Station.fieldStatus; }
diff --git a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
index e521b76627..d2e9a0d04f 100644
--- a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
+++ b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
@@ -29,9 +29,9 @@
#include "L1Grid.h"
#include "L1HitArea.h"
#include "L1HitPoint.h"
-#include "L1Portion.h"
#include "L1Track.h"
#include "L1TrackPar.h"
+#include "L1TripletConstructor.h"
#ifdef _OPENMP
#include "omp.h"
@@ -122,1355 +122,6 @@ bool L1Algo::checkTripletMatch(const L1Triplet& l, const L1Triplet& r, fscal& dc
return true;
}
-inline void L1Algo::findSingletsStep0( // input
- Tindex start_lh, Tindex n1_l, L1HitPoint* Hits_l,
- // output
- L1HitIndex_t* hitsl, fvec* x_l, fvec* y_l, fvec* z_l, fvec* t_l, fvec* dx_l, fvec* dxy_l, fvec* dy_l, fvec* dt2_l)
-{
-
- /// Prepare the portion of data of left hits of a triplet:
- /// all hits except the last and the second last station will be procesesed in the algorithm,
- /// the data is orginesed in order to be used by SIMD
-
- const Tindex end_lh = start_lh + n1_l;
- const int lastV = (n1_l - 1) / fvec::size();
- if (lastV >= 0) {
- // set some positive errors to unfilled part of vectors in order to avoid nans
- L1HitPoint& hitl = Hits_l[0];
- x_l[lastV] = hitl.X();
- y_l[lastV] = hitl.Y();
- z_l[lastV] = hitl.Z();
- t_l[lastV] = hitl.T();
- dx_l[lastV] = hitl.dX();
- dxy_l[lastV] = hitl.dXY();
- dy_l[lastV] = hitl.dY();
- dt2_l[lastV] = hitl.dT2();
- }
-
- for (Tindex ilh = start_lh, i1 = 0; ilh < end_lh; ++ilh, ++i1) {
- Tindex i1_V = i1 / fvec::size();
- Tindex i1_4 = i1 % fvec::size();
- L1HitPoint& hitl = Hits_l[ilh];
- hitsl[i1] = ilh;
- x_l[i1_V][i1_4] = hitl.X();
- y_l[i1_V][i1_4] = hitl.Y();
- z_l[i1_V][i1_4] = hitl.Z();
- t_l[i1_V][i1_4] = hitl.T();
- dx_l[i1_V][i1_4] = hitl.dX();
- dxy_l[i1_V][i1_4] = hitl.dXY();
- dy_l[i1_V][i1_4] = hitl.dY();
- dt2_l[i1_V][i1_4] = hitl.dT2();
- }
-}
-
-
-inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
- int istal,
- int istam, /// indexes of left and middle stations of a triplet
- Tindex n1_V, ///
- fvec* x_l, fvec* y_l, fvec* z_l, fvec* t_l, fvec* dx_l, fvec* dxy_l, fvec* dy_l, fvec* dt2_l,
- // output
- // L1TrackPar *T_1,
- L1TrackPar* T_1, L1FieldRegion* fld_1)
-{
-
- /// Get the field approximation. Add the target to parameters estimation.
- /// Propagaete to the middle station of a triplet.
-
- const L1Station& stal = fParameters.GetStation(istal);
- const L1Station& stam = fParameters.GetStation(istam);
-
- // stations for approximating the field between the target and the left hit
- int fld0Ista1 = istal;
- if (fld0Ista1 >= fNfieldStations) { fld0Ista1 = fNfieldStations - 1; }
- if (fld0Ista1 < 1) { fld0Ista1 = 1; }
- int fld0Ista0 = fld0Ista1 / 2; // station between istal and the target
-
- assert(0 <= fld0Ista0 && fld0Ista0 < fld0Ista1 && fld0Ista1 < fParameters.GetNstationsActive());
-
- // stations for approximating the field between the left and the right hit
- int fld1Ista0 = istal;
- int fld1Ista1 = istam;
- int fld1Ista2 = istam + 1;
- if (fld1Ista2 >= fParameters.GetNstationsActive()) {
- fld1Ista2 = istam;
- fld1Ista1 = istam - 1;
- if (fld1Ista0 >= fld1Ista1) { fld1Ista0 = fld1Ista1 - 1; }
- }
-
- assert(0 <= fld1Ista0 && fld1Ista0 < fld1Ista1 && fld1Ista1 < fld1Ista2
- && fld1Ista2 < fParameters.GetNstationsActive());
-
- const L1Station& fld0Sta0 = fParameters.GetStation(fld0Ista0);
- const L1Station& fld0Sta1 = fParameters.GetStation(fld0Ista1);
- const L1Station& fld1Sta0 = fParameters.GetStation(fld1Ista0);
- const L1Station& fld1Sta1 = fParameters.GetStation(fld1Ista1);
- const L1Station& fld1Sta2 = fParameters.GetStation(fld1Ista2);
-
- L1Fit fit;
- fit.SetParticleMass(GetDefaultParticleMass());
- fit.SetMask(fmask::One());
- fit.SetQp0(fMaxInvMom);
-
- if (fpCurrentIteration->GetElectronFlag()) { fit.SetParticleMass(L1Constants::phys::kElectronMass); }
- else {
- fit.SetParticleMass(fDefaultMass);
- }
-
- for (int i1_V = 0; i1_V < n1_V; i1_V++) {
-
- // field made by the left hit, the target and the station istac in-between.
- // is used for extrapolation to the target and to the middle hit
- L1FieldRegion fld0;
-
- // field, made by the left hit, the middle station and the right station
- // Will be used for extrapolation to the right hit
- L1FieldRegion& fld1 = fld_1[i1_V];
-
- L1FieldValue lB, mB, cB, rB _fvecalignment;
- L1FieldValue l_B_global, centB_global _fvecalignment;
-
- // get the magnetic field along the track.
- // (suppose track is straight line with origin in the target)
-
- fvec& xl = x_l[i1_V];
- fvec& yl = y_l[i1_V];
- fvec& zl = z_l[i1_V];
- fvec& time = t_l[i1_V];
- fvec& timeEr2 = dt2_l[i1_V];
- const fvec dzli = 1.f / (zl - fTargZ);
-
- const fvec tx = (xl - fTargX) * dzli;
- const fvec ty = (yl - fTargY) * dzli;
-
- L1FieldValue b00, b01, b10, b11, b12;
- fld0Sta0.fieldSlice.GetFieldValue(fTargX + tx * (fld0Sta0.fZ - fTargZ), fTargY + ty * (fld0Sta0.fZ - fTargZ), b00);
- fld0Sta1.fieldSlice.GetFieldValue(fTargX + tx * (fld0Sta1.fZ - fTargZ), fTargY + ty * (fld0Sta1.fZ - fTargZ), b01);
- fld1Sta0.fieldSlice.GetFieldValue(fTargX + tx * (fld1Sta0.fZ - fTargZ), fTargY + ty * (fld1Sta0.fZ - fTargZ), b10);
- fld1Sta1.fieldSlice.GetFieldValue(fTargX + tx * (fld1Sta1.fZ - fTargZ), fTargY + ty * (fld1Sta1.fZ - fTargZ), b11);
- fld1Sta2.fieldSlice.GetFieldValue(fTargX + tx * (fld1Sta2.fZ - fTargZ), fTargY + ty * (fld1Sta2.fZ - fTargZ), b12);
-
- fld0.Set(fTargB, fTargZ, b00, fld0Sta0.fZ, b01, fld0Sta1.fZ);
- fld1.Set(b10, fld1Sta0.fZ, b11, fld1Sta1.fZ, b12, fld1Sta2.fZ);
-
- L1TrackPar& T = fit.Tr();
-
- T.x = xl;
- T.y = yl;
- T.z = zl;
- T.tx = tx;
- T.ty = ty;
- T.qp = fvec(0.);
- T.t = time;
- T.vi = fvec(0.);
-
- fvec txErr2 = fMaxSlopePV * fMaxSlopePV / fvec(9.);
- fvec qpErr2 = fMaxInvMom * fMaxInvMom / fvec(9.);
-
- T.ResetErrors(1., 1., txErr2, txErr2, qpErr2, (stal.timeInfo ? timeEr2 : 1.e6), 1.e2);
- T.C00 = dx_l[i1_V] * dx_l[i1_V];
- T.C10 = dxy_l[i1_V];
- T.C11 = dy_l[i1_V] * dy_l[i1_V];
-
- T.chi2 = 0.;
- T.NDF = (fpCurrentIteration->GetPrimaryFlag()) ? fvec(2.) : fvec(0.);
-
- // TODO: iteration parameter: "Starting NDF of track parameters"
- // NDF = number of track parameters (6: x, y, tx, ty, qp, time) - number of measured parameters (3: x, y, time) on station or (2: x, y) on target
-
- // add the target constraint
-
- fit.AddTargetToLine(fTargX, fTargY, fTargZ, TargetXYInfo, fld0);
-
- fit.AddMsInMaterial(fParameters.GetMaterialThickness(istal, T.x, T.y));
-
- // extrapolate to the middle hit
-
- fit.ExtrapolateLine(stam.fZ, fld1);
-
- T_1[i1_V] = T;
- } // i1_V
-}
-
-
-inline void L1Algo::findDoubletsStep0(
- /// input
- Tindex n1, const L1Station& stal, L1HitPoint* vHits_l, const L1Station& stam, L1HitPoint* vHits_m, L1TrackPar* T_1,
- L1HitIndex_t* hitsl_1,
- /// output
- Tindex& n2, L1Vector<L1HitIndex_t>& i1_2,
-#ifdef DOUB_PERFORMANCE
- L1Vector<L1HitIndex_t>& hitsl_2,
-#endif // DOUB_PERFORMANCE
- L1Vector<L1HitIndex_t>& hitsm_2)
-{
- /// Find the doublets. Reformat data into portions of doublets.
-
- assert(i1_2.size() == 0);
-
- int iStaL = &stal - fParameters.GetStations().begin();
- int iStaM = &stam - fParameters.GetStations().begin();
-
- n2 = 0; // number of doublets
-
- L1Fit fit;
- fit.SetParticleMass(GetDefaultParticleMass());
-
- for (Tindex i1 = 0; i1 < n1; ++i1) // for each singlet
- {
- unsigned int Ndoublets = 0;
- const Tindex i1_V = i1 / fvec::size();
- const Tindex i1_4 = i1 % fvec::size();
-
- fit.SetTrack(T_1[i1_V]);
- L1TrackPar& T1 = fit.Tr();
-
- // assert(T1.IsEntryConsistent(true, i1_4));
- // if (!T1.IsEntryConsistent(false, i1_4)) continue;
-
- const fvec Pick_m22 = (fDoubletChi2Cut - T1.chi2);
- // if make it bigger the found hits will be rejected later because of the chi2 cut.
- // Pick_m22 is not used, search for mean squared, 2nd version
-
- unsigned int indFirstDoublet = hitsm_2.size();
-
- // -- collect possible doublets --
- const fscal iz = 1.f / (T1.z[i1_4] - fParameters.GetTargetPositionZ()[0]);
- const fscal timeError2 = T1.C55[i1_4];
- const fscal time = T1.t[i1_4];
-
- L1HitAreaTime areaTime(vGridTime[iStaM], T1.x[i1_4] * iz, T1.y[i1_4] * iz,
- (sqrt(Pick_m22 * T1.C00) + fMaxDx[iStaM] + fMaxDZ * abs(T1.tx))[i1_4] * iz,
- (sqrt(Pick_m22 * T1.C11) + fMaxDy[iStaM] + fMaxDZ * abs(T1.ty))[i1_4] * iz, time,
- sqrt(timeError2) * 5);
-
- for (L1HitIndex_t imh = -1; true;) { // loop over the hits in the area
- if (fParameters.DevIsIgnoreHitSearchAreas()) {
- imh++;
- if ((L1HitIndex_t) imh >= (fGridHitStopIndex[iStaM] - fGridHitStartIndex[iStaM])) { break; }
- }
- else {
- if (!areaTime.GetNext(imh)) { break; }
- }
-
- L1HitPoint& hitm = vHits_m[imh];
- if (hitm.IsSuppressed()) continue;
-
- const L1HitPoint& hitl = vHits_l[hitsl_1[i1]];
-
- if (fParameters.DevIsMatchDoubletsViaMc()) { // trd2d
- int indL = fGridHitStartIndex[iStaL] + hitsl_1[i1];
- int indM = fGridHitStartIndex[iStaM] + imh;
- int iMC = GetMcTrackIdForGridHit(indL);
- if (iMC < 0 || iMC != GetMcTrackIdForGridHit(indM)) { continue; }
- }
-
- // check y-boundaries
- //TODO: move hardcoded cuts to parameters
- if ((stam.timeInfo) && (stal.timeInfo)) {
- if (fabs(time - hitm.T()) > sqrt(timeError2 + hitm.dT2()) * 5) continue;
- if (fabs(time - hitm.T()) > 30) continue;
- }
-
- // - check whether hit belong to the window ( track position +\- errors ) -
- const fscal zm = hitm.Z();
- // L1TrackPar T1_new = T1;
- // fvec dz = fvec(zm) - T1.z;
- // L1ExtrapolateTime(T1, dz);
-
-
- // if (fabs(T1_new.t[i1_4]-hitm.time)>sqrt(T1_new.C55[i1_4]+hitm.timeEr*hitm.timeEr)*4) continue;
- // if (fabs(T1_new.t[i1_4]-hitm.time)>sqrt(2.9*2.9)*5) continue;
- // const fscal dt_est2 = Pick_m22[i1_4] * fabs(T1_new.C55[i1_4] + hitm.timeEr*hitm.timeEr);
- // const fscal dt = hitm.time - T1_new.t[i1_4];
- // if ( dt*dt > dt_est2 && dt < 0 ) continue;
-
- fvec y, C11;
- fit.ExtrapolateYC11Line(zm, y, C11);
-
- fscal dy_est2 = Pick_m22[i1_4] * fabs(C11[i1_4] + hitm.dY2());
-
- fscal xm = hitm.X();
- fscal ym = hitm.Y();
-
- const fscal dY = ym - y[i1_4];
-
- if (dY * dY > dy_est2) continue;
-
- // check x-boundaries
- fvec x, C00;
- fit.ExtrapolateXC00Line(zm, x, C00);
-
- fscal dx_est2 = Pick_m22[i1_4] * fabs(C00[i1_4] + hitm.dX2());
-
- const fscal dX = xm - x[i1_4];
-
- if (dX * dX > dx_est2) continue;
-
- // check chi2
-
- fvec C10;
- fit.ExtrapolateC10Line(zm, C10);
-
- auto [chi2x, chi2u] =
- L1Fit::GetChi2XChi2U(hitm.X(), hitm.Y(), hitm.dX2(), hitm.dXY(), hitm.dY2(), x, y, C00, C10, C11);
-
- // TODO: adjust the cut, cut on chi2x & chi2u separately
- if (!fpCurrentIteration->GetTrackFromTripletsFlag()) {
- if (chi2x[i1_4] > fDoubletChi2Cut) continue;
- if (chi2x[i1_4] + chi2u[i1_4] > fDoubletChi2Cut) continue;
- }
-
- // FilterTime(T1, hitm.T(), hitm.dT2());
-
- // check if there is a second hit on the same station
- {
- bool isOtherHit = 0;
- fscal dz = hitm.Z() - hitl.Z();
- fscal tx = (hitm.X() - hitl.X()) / dz;
- fscal ty = (hitm.Y() - hitl.Y()) / dz;
- fscal tt = (hitm.T() - hitl.T()) / dz;
-
- for (unsigned int imh1 = indFirstDoublet; imh1 < hitsm_2.size(); imh1++) {
- const L1HitPoint& hitm1 = vHits_m[hitsm_2[imh1]];
-
- if ((stam.timeInfo) && (stal.timeInfo)) {
- fscal dt = hitm.T() + tt * (hitm1.Z() - hitm.Z()) - hitm1.T();
- if (dt * dt > 30. * (hitm.dT2() + hitm1.dT2())) { continue; }
- }
-
- fscal dx = hitm.X() + tx * (hitm1.Z() - hitm.Z()) - hitm1.X();
- if (dx * dx > 20. * (hitm.dX2() + hitm1.dX2())) { continue; }
-
- fscal dy = hitm.Y() + ty * (hitm1.Z() - hitm.Z()) - hitm1.Y();
- if (dy * dy > 30. * (hitm.dY2() + hitm1.dY2())) { continue; }
-
- if (fParameters.DevIsSuppressOverlapHitsViaMc()) {
- int indL = fGridHitStartIndex[iStaL] + hitsl_1[i1];
- int indM = fGridHitStartIndex[iStaM] + imh;
- int indM1 = fGridHitStartIndex[iStaM] + hitsm_2[imh1];
- int iMC = GetMcTrackIdForGridHit(indL);
- if ((iMC != GetMcTrackIdForGridHit(indM)) || (iMC != GetMcTrackIdForGridHit(indM1))) { continue; }
- }
-
- isOtherHit = true;
- break;
- }
-
- if (isOtherHit) {
- hitm.SetIsSuppresed(1);
- continue;
- }
- }
-
- i1_2.push_back(i1);
-#ifdef DOUB_PERFORMANCE
- hitsl_2.push_back(hitsl_1[i1]);
-#endif // DOUB_PERFORMANCE
- hitsm_2.push_back(imh);
-
- Ndoublets++;
- n2++;
-
-#ifndef FAST_CODE
-//TODO: optimise triplet portion limit and put a warning
-// assert(Ndoublets < fParameters.GetMaxDoubletsPerSinglet());
-#endif // FAST_CODE
-
- if (Ndoublets >= fParameters.GetMaxDoubletsPerSinglet()) {
- n2 = n2 - Ndoublets;
- i1_2.reduce(i1_2.size() - Ndoublets);
- hitsm_2.reduce(hitsm_2.size() - Ndoublets);
- break;
- }
- } // loop over the hits in the area
-
- T_1[i1_V] = fit.Tr();
-
- } // for i1
-}
-
-
-/// Add the middle hits to parameters estimation. Propagate to right station.
-/// Find the triplets(right hit). Reformat data in the portion of triplets.
-inline void L1Algo::findTripletsStep0( // input
- L1HitPoint* vHits_r, int iStaL, int iStaM, int iStaR, L1HitPoint* vHits_m, L1TrackPar* T_1, L1FieldRegion* fld_1,
- L1HitIndex_t* hitsl_1, Tindex n2, L1Vector<L1HitIndex_t>& hitsm_2, L1Vector<L1HitIndex_t>& i1_2,
- // output
- Tindex& n3, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3, L1Vector<L1HitIndex_t>& hitsm_3,
- L1Vector<L1HitIndex_t>& hitsr_3, L1Vector<fvec>& x_3, L1Vector<fvec>& y_3, L1Vector<fvec>& z_3, L1Vector<fvec>& t_3,
- L1Vector<fvec>& dx_3, L1Vector<fvec>& dxy_3, L1Vector<fvec>& dy_3, L1Vector<fvec>& dt2_3)
-{
- const L1Station& stal = fParameters.GetStation(iStaL);
- const L1Station& stam = fParameters.GetStation(iStaM);
- const L1Station& star = fParameters.GetStation(iStaR);
-
- std::array<L1HitIndex_t, fvec::size()> hitsl_2 = {};
- std::array<L1HitIndex_t, fvec::size()> hitsm_2_tmp = {};
-
- std::fill(hitsl_2.begin(), hitsl_2.end(), undef::kU32);
- std::fill(hitsm_2_tmp.begin(), hitsm_2_tmp.end(), undef::kU32);
-
- L1TrackPar L1TrackPar_0;
- // SG!! to avoid nans in unfilled part
- //TODO: SG: investigate, it changes the results !!
- /*
- L1TrackPar_0.C00 = 1.f;
- L1TrackPar_0.C11 = 1.f;
- L1TrackPar_0.C22 = 1.f;
- L1TrackPar_0.C33 = 1.f;
- L1TrackPar_0.C44 = 1.f;
- L1TrackPar_0.C55 = 1.f;
- */
-
- bool isMomentumFitted = (fIsTargetField || (stal.fieldStatus != 0) || (stam.fieldStatus != 0));
- //bool isTimeFitted = ((stal.timeInfo != 0) || (stam.timeInfo != 0));
-
- L1Fit fit;
- fit.SetParticleMass(fDefaultMass);
- fit.SetMask(fmask::One());
- fit.SetQp0(fvec(0.));
-
- n3 = 0;
- Tindex n3_V = 0, n3_4 = 0;
-
- T_3.reset(1, L1TrackPar_0);
-
- x_3.reset(1, fvec::Zero());
- y_3.reset(1, fvec::Zero());
- z_3.reset(1, fvec::Zero());
- t_3.reset(1, fvec::Zero());
-
- dx_3.reset(1, fvec::One());
- dxy_3.reset(1, fvec::Zero());
- dy_3.reset(1, fvec::One());
- dt2_3.reset(1, fvec::One());
-
- assert(iStaR < fParameters.GetNstationsActive()); //TODO SG!!! check if it is needed
-
- if (iStaR >= fParameters.GetNstationsActive()) return;
-
- // ---- Add the middle hits to parameters estimation. Propagate to right station. ----
-
- for (Tindex i2 = 0; i2 < n2;) {
- fit.SetTrack(L1TrackPar_0);
- fit.SetQp0(fvec(0.));
- L1TrackPar& T2 = fit.Tr();
-
- L1FieldRegion f2;
- // pack the data
- fvec x_2 = 0.f;
- fvec y_2 = 0.f;
- fvec z_2 = 0.f;
- fvec t_2 = 0.f;
- L1XYMeasurementInfo cov_2;
- cov_2.C00 = 1.f;
- cov_2.C10 = 0.f;
- cov_2.C11 = 1.f;
- fvec dt2_2 = 1.f;
-
- size_t n2_4 = 0;
- for (; n2_4 < fvec::size() && i2 < n2; i2++, n2_4++) {
-
- const Tindex& i1 = i1_2[i2];
- const Tindex i1_V = i1 / fvec::size();
- const Tindex i1_4 = i1 % fvec::size();
-
- const L1TrackPar& T1 = T_1[i1_V];
-
- //assert(T1.IsEntryConsistent(true, i1_4));
- //if (!T1.IsEntryConsistent(false, i1_4)) continue;
-
- const L1FieldRegion& f1 = fld_1[i1_V];
- T2.SetOneEntry(n2_4, T1, i1_4);
- f2.SetOneEntry(n2_4, f1, i1_4);
-
- const Tindex& imh = hitsm_2[i2];
- const L1HitPoint& hitm = vHits_m[imh];
-
- x_2[n2_4] = hitm.X();
- y_2[n2_4] = hitm.Y();
- z_2[n2_4] = hitm.Z();
- t_2[n2_4] = hitm.T();
- cov_2.C00[n2_4] = hitm.dX2();
- cov_2.C10[n2_4] = hitm.dXY();
- cov_2.C11[n2_4] = hitm.dY2();
- dt2_2[n2_4] = hitm.dT2();
-
- hitsl_2[n2_4] = hitsl_1[i1];
- hitsm_2_tmp[n2_4] = hitsm_2[i2];
- } // n2_4
-
- // add the middle hit
-
- fit.ExtrapolateLineNoField(z_2);
-
- fit.FilterXY(cov_2, x_2, y_2);
-
- fit.FilterTime(t_2, dt2_2, stam.timeInfo);
-
- fit.SetQp0(isMomentumFitted ? fit.Tr().qp : fMaxInvMom);
-
- fit.AddMsInMaterial(fParameters.GetMaterialThickness(iStaM, T2.x, T2.y));
-
- fit.SetQp0(fit.Tr().qp);
-
- // extrapolate to the right hit station
-
- fit.Extrapolate(star.fZ, f2);
-
- // assert(T2.IsConsistent(true, n2_4));
-
- // ---- Find the triplets(right hit). Reformat data in the portion of triplets. ----
- for (size_t i2_4 = 0; i2_4 < n2_4; ++i2_4) {
-
- //if (!T2.IsEntryConsistent(false, i2_4)) { continue; }
- // FIXME: SZh 22.08.2023: Remove the hack!! (Why is it used only for kSts?)
- if (kSts == fTrackingMode && (T2.C44[i2_4] < 0)) { continue; }
- if (T2.C00[i2_4] < 0 || T2.C11[i2_4] < 0 || T2.C22[i2_4] < 0 || T2.C33[i2_4] < 0 || T2.C55[i2_4] < 0) continue;
-
- if (fabs(T2.tx[i2_4]) > fMaxSlope) continue;
- if (fabs(T2.ty[i2_4]) > fMaxSlope) continue;
-
- const fvec Pick_r22 = fTripletChi2Cut - T2.chi2 + (!fpCurrentIteration->GetTrackFromTripletsFlag() ? 0 : 1);
- const fscal timeError2 = T2.C55[i2_4];
- const fscal time = T2.t[i2_4];
- // find first possible hit
-
-
- const fscal iz = 1.f / (T2.z[i2_4] - fParameters.GetTargetPositionZ()[0]);
- L1HitAreaTime area(vGridTime[&star - fParameters.GetStations().begin()], T2.x[i2_4] * iz, T2.y[i2_4] * iz,
- (sqrt(Pick_r22 * T2.C00) + fMaxDx[iStaM] + fMaxDZ * abs(T2.tx))[i2_4] * iz,
- (sqrt(Pick_r22 * T2.C11) + fMaxDy[iStaM] + fMaxDZ * abs(T2.ty))[i2_4] * iz, time,
- sqrt(timeError2) * 5);
-
- L1HitIndex_t irh = 0;
- L1HitIndex_t doubletNtriplets = 0;
- int irh1 = -1;
- while (true) {
- if (fParameters.DevIsIgnoreHitSearchAreas()) {
- irh1++;
- if ((L1HitIndex_t) irh1 >= (fGridHitStopIndex[iStaR] - fGridHitStartIndex[iStaR])) break;
- irh = irh1;
- }
- else {
- if (!area.GetNext(irh)) break;
- }
-
- // while (area.GetNext(irh)) {
- //for (int irh = 0; irh < ( fGridHitStopIndex[iStaR] - fGridHitStartIndex[iStaR] ); irh++){
- const L1HitPoint& hitr = vHits_r[irh];
- if (hitr.IsSuppressed()) continue;
-
- if (fParameters.DevIsMatchTripletsViaMc()) {
- int indL = fGridHitStartIndex[iStaL] + hitsl_2[i2_4];
- int indM = fGridHitStartIndex[iStaM] + hitsm_2_tmp[i2_4];
- int indR = fGridHitStartIndex[iStaR] + irh;
- int mcL = GetMcTrackIdForGridHit(indL);
- if (mcL < 0 || mcL != GetMcTrackIdForGridHit(indM) || mcL != GetMcTrackIdForGridHit(indR)) { continue; }
- }
-
- const fscal zr = hitr.Z();
-
- fscal xr = hitr.X();
- fscal yr = hitr.Y();
-
- //auto [xr, yr] = star.ConvUVtoXY<fscal>(hitr.U(), hitr.V());
-
- L1Fit fit3;
- fit3.SetParticleMass(fDefaultMass);
- fit3.SetMask(fmask::One());
- fit3.SetTrack(T2);
- L1TrackPar& T_cur = fit3.Tr();
-
- fit3.ExtrapolateLineNoField(zr);
-
- if ((star.timeInfo) && (stam.timeInfo)) {
- auto dt = T_cur.t[i2_4] - hitr.T();
- if (dt > (T_cur.C55[i2_4] + hitr.dT2()) * 5) continue;
- }
-
- // TODO: SG: hardcoded cut of 30 ns
- if ((star.timeInfo) && (stam.timeInfo)) {
- if (fabs(T_cur.t[i2_4] - hitr.T()) > 30) continue;
- }
-
- // - check whether hit belong to the window ( track position +\- errors ) -
- // check lower boundary
- fvec y, C11;
- fit3.ExtrapolateYC11Line(zr, y, C11);
-
- fscal dy_est2 = (Pick_r22[i2_4] * (fabs(C11[i2_4] + hitr.dY2())));
-
- const fscal dY = yr - y[i2_4];
- const fscal dY2 = dY * dY;
-
- if (dY2 > dy_est2) continue; // if (yr > y_plus_new [i2_4] ) continue;
-
- // check x-boundaries
- fvec x, C00;
- fit3.ExtrapolateXC00Line(zr, x, C00);
-
- fscal dx_est2 = (Pick_r22[i2_4] * (fabs(C00[i2_4] + hitr.dX2())));
-
- const fscal dX = xr - x[i2_4];
- if (dX * dX > dx_est2) continue;
- // check chi2 // not effective
- fvec C10;
- fit3.ExtrapolateC10Line(zr, C10);
-
- auto [chi2x, chi2u] =
- L1Fit::GetChi2XChi2U(hitr.X(), hitr.Y(), hitr.dX2(), hitr.dXY(), hitr.dY2(), x, y, C00, C10, C11);
-
- //fit3.FilterTime(hitr.T(), hitr.dT2(), star.timeInfo);
-
- if (!fpCurrentIteration->GetTrackFromTripletsFlag()) {
- if (chi2x[i2_4] + chi2u[i2_4] > fTripletChi2Cut || T_cur.C55[i2_4] < 0) {
- continue; // chi2_triplet < CHI2_CUT
- }
- }
-
-#ifndef FAST_CODE
-//TODO: optimise triplet portion limit and put a warning
-// assert(doubletNtriplets < fParameters.GetMaxTripletPerDoublets());
-#endif // FAST_CODE
-
- if (doubletNtriplets >= fParameters.GetMaxTripletPerDoublets()) {
- LOG(debug) << "L1: GetMaxTripletPerDoublets==" << fParameters.GetMaxTripletPerDoublets()
- << " reached in findTripletsStep0()";
- // reject already created triplets for this doublet
- n3 = n3 - doubletNtriplets;
- n3_V = n3 / fvec::size();
- n3_4 = n3 % fvec::size();
- hitsl_3.reduce(n3);
- hitsm_3.reduce(n3);
- hitsr_3.reduce(n3);
- break;
- }
-
- // pack triplet
-
- doubletNtriplets++;
-
- L1TrackPar& T3 = T_3[n3_V];
-
- hitsl_3.push_back(hitsl_2[i2_4]);
- hitsm_3.push_back(hitsm_2_tmp[i2_4]);
- hitsr_3.push_back(irh);
-
-
- T3.SetOneEntry(n3_4, T2, i2_4);
-
- x_3[n3_V][n3_4] = hitr.X();
- y_3[n3_V][n3_4] = hitr.Y();
- z_3[n3_V][n3_4] = zr;
- t_3[n3_V][n3_4] = hitr.T();
- dx_3[n3_V][n3_4] = hitr.dX();
- dxy_3[n3_V][n3_4] = hitr.dXY();
- dy_3[n3_V][n3_4] = hitr.dY();
- dt2_3[n3_V][n3_4] = hitr.dT2();
-
- n3++;
- n3_V = n3 / fvec::size();
- n3_4 = n3 % fvec::size();
-
- assert(n3 == (int) hitsl_3.size());
- assert(n3 == (int) hitsm_3.size());
- assert(n3 == (int) hitsr_3.size());
-
- if (0 == n3_4) {
- T_3.push_back(L1TrackPar_0);
- x_3.push_back(fvec::Zero());
- y_3.push_back(fvec::Zero());
- z_3.push_back(fvec::Zero());
- t_3.push_back(fvec::Zero());
- dx_3.push_back(fvec::One());
- dxy_3.push_back(fvec::Zero());
- dy_3.push_back(fvec::One());
- dt2_3.push_back(fvec::One());
- }
- } // search area
-
- } // i2_4
- } // i2_V
-}
-
-/// Add the right hits to parameters estimation.
-inline void L1Algo::findTripletsStep1( // input
- Tindex n3_V, const L1Station& star, L1Vector<fvec>& x_3, L1Vector<fvec>& y_3, L1Vector<fvec>& z_3,
- L1Vector<fvec>& t_3, L1Vector<fvec>& dx_3, L1Vector<fvec>& dxy_3, L1Vector<fvec>& dy_3, L1Vector<fvec>& dt2_3,
- // output
- L1Vector<L1TrackPar>& T_3)
-{
- L1Fit fit;
- fit.SetParticleMass(fDefaultMass);
- fit.SetMask(fmask::One());
- for (Tindex i3_V = 0; i3_V < n3_V; ++i3_V) {
- fit.SetTrack(T_3[i3_V]);
- fit.ExtrapolateLineNoField(z_3[i3_V]);
- L1XYMeasurementInfo info;
- info.C00 = dx_3[i3_V] * dx_3[i3_V];
- info.C10 = dxy_3[i3_V];
- info.C11 = dy_3[i3_V] * dy_3[i3_V];
- fit.FilterXY(info, x_3[i3_V], y_3[i3_V]);
- if (kMcbm != fTrackingMode) { fit.FilterTime(t_3[i3_V], dt2_3[i3_V], star.timeInfo); }
- T_3[i3_V] = fit.Tr();
- }
-}
-
-
-inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar, L1Vector<L1TrackPar>& T_3,
- L1Vector<L1HitIndex_t>& hitsl_3, L1Vector<L1HitIndex_t>& hitsm_3,
- L1Vector<L1HitIndex_t>& hitsr_3, int nIterations)
-{
-
- nIterations = 2;
-
- /// Refit Triplets
-
- ca::tools::Debugger::Instance().AddNtuple("triplets", "ev:mc:sta:p:vx:vy:vz:chi2:ndf");
-
- L1Fit fit;
- fit.SetMask(fmask::One());
-
- if (fpCurrentIteration->GetElectronFlag()) { fit.SetParticleMass(L1Constants::phys::kElectronMass); }
- else {
- fit.SetParticleMass(GetDefaultParticleMass());
- }
- const int NHits = 3; // triplets
-
- // prepare data
- int ista[NHits] = {istal, istam, istar};
-
- L1Station sta[3];
- for (int is = 0; is < NHits; ++is) {
- sta[is] = fParameters.GetStation(ista[is]);
- };
-
- const L1Station& stal = fParameters.GetStation(istal);
- const L1Station& stam = fParameters.GetStation(istam);
- const L1Station& star = fParameters.GetStation(istar);
-
- bool isMomentumFitted = ((stal.fieldStatus != 0) || (stam.fieldStatus != 0) || (star.fieldStatus != 0));
- bool isTimeFitted = ((stal.timeInfo != 0) || (stam.timeInfo != 0) || (star.timeInfo != 0));
- fvec ndf = -4; // straight line
- ndf += isMomentumFitted ? -1 : 0;
- ndf += isTimeFitted ? -1 : 0;
-
- for (int i3 = 0; i3 < n3; ++i3) {
- int i3_V = i3 / fvec::size();
- int i3_4 = i3 % fvec::size();
-
- L1TrackPar& T3 = T_3[i3_V];
-
- //if (!T3.IsEntryConsistent(false, i3_4)) continue;
-
- // prepare data
- L1HitIndex_t ihit[NHits] = {fGridHitIds[hitsl_3[i3] + fGridHitStartIndex[ista[0]]],
- fGridHitIds[hitsm_3[i3] + fGridHitStartIndex[ista[1]]],
- fGridHitIds[hitsr_3[i3] + fGridHitStartIndex[ista[2]]]};
-
- if (fParameters.DevIsMatchTripletsViaMc()) {
- int mc1 = GetMcTrackIdForCaHit(ihit[0]);
- int mc2 = GetMcTrackIdForCaHit(ihit[1]);
- int mc3 = GetMcTrackIdForCaHit(ihit[2]);
- if ((mc1 != mc2) || (mc1 != mc3)) { continue; }
- }
-
- fscal x[NHits], y[NHits], z[NHits], t[NHits];
- fscal dt2[NHits];
- L1XYMeasurementInfo cov[NHits];
-
- for (int ih = 0; ih < NHits; ++ih) {
- const L1Hit& hit = fInputData.GetHit(ihit[ih]);
- x[ih] = hit.x;
- y[ih] = hit.y;
- z[ih] = hit.z;
- t[ih] = hit.t;
- cov[ih].C00 = hit.dx * hit.dx;
- cov[ih].C10 = hit.dxy;
- cov[ih].C11 = hit.dy * hit.dy;
- dt2[ih] = hit.dt2;
- };
-
- // find the field along the track
-
- L1FieldValue B[3] _fvecalignment;
- L1FieldRegion fld _fvecalignment;
- L1FieldRegion fldTarget _fvecalignment;
-
- fvec tx[3] = {(x[1] - x[0]) / (z[1] - z[0]), (x[2] - x[0]) / (z[2] - z[0]), (x[2] - x[1]) / (z[2] - z[1])};
- fvec ty[3] = {(y[1] - y[0]) / (z[1] - z[0]), (y[2] - y[0]) / (z[2] - z[0]), (y[2] - y[1]) / (z[2] - z[1])};
- for (int ih = 0; ih < NHits; ++ih) {
- fvec dz = (sta[ih].fZ - z[ih]);
- sta[ih].fieldSlice.GetFieldValue(x[ih] + tx[ih] * dz, y[ih] + ty[ih] * dz, B[ih]);
- };
-
- fld.Set(B[0], sta[0].fZ, B[1], sta[1].fZ, B[2], sta[2].fZ);
- fldTarget.Set(fTargB, fTargZ, B[0], sta[0].fZ, B[1], sta[1].fZ);
-
- L1TrackPar& T = fit.Tr();
-
- // initial parameters
- {
- fvec dz01 = 1. / (z[1] - z[0]);
- T.tx = (x[1] - x[0]) * dz01;
- T.ty = (y[1] - y[0]) * dz01;
- T.qp = 0.; //(fscal) T3.qp[i3_4];
- T.vi = 0.;
- }
-
- // repeat several times in order to improve the precision
- for (int iiter = 0; iiter < nIterations; ++iiter) {
- // fit forward
- {
- fit.SetQp0(T.qp);
- fit.Qp0()(fit.Qp0() > GetMaxInvMom()) = GetMaxInvMom();
- fit.Qp0()(fit.Qp0() < -GetMaxInvMom()) = -GetMaxInvMom();
-
- T.qp = 0.;
- T.vi = 0.;
- int ih0 = 0;
- T.x = x[ih0];
- T.y = y[ih0];
- T.z = z[ih0];
- T.t = t[ih0];
-
- T.ResetErrors(1., 1., 1., 1., 100., (sta[ih0].timeInfo ? dt2[ih0] : 1.e6), 1.e2);
- T.C00 = cov[ih0].C00;
- T.C10 = cov[ih0].C10;
- T.C11 = cov[ih0].C11;
-
- T.NDF = ndf;
-
- // add the target constraint
- fit.AddTargetToLine(fTargX, fTargY, fTargZ, TargetXYInfo, fldTarget);
-
- for (int ih = 1; ih < NHits; ++ih) {
- fit.Extrapolate(z[ih], fld);
- auto radThick = fParameters.GetMaterialThickness(ista[ih], T.x, T.y);
- fit.AddMsInMaterial(radThick);
- fit.EnergyLossCorrection(radThick, fvec(-1.f));
- fit.FilterXY(cov[ih], x[ih], y[ih]);
- fit.FilterTime(t[ih], dt2[ih], sta[ih].timeInfo);
- }
- }
-
- if (iiter == nIterations - 1) break;
-
- // fit backward
- {
- fit.SetQp0(T.qp);
- fit.Qp0()(fit.Qp0() > GetMaxInvMom()) = GetMaxInvMom();
- fit.Qp0()(fit.Qp0() < -GetMaxInvMom()) = -GetMaxInvMom();
-
- T.NDF = ndf;
- T.qp = 0.;
- T.vi = 0.;
- int ih0 = NHits - 1;
- T.x = x[ih0];
- T.y = y[ih0];
- T.z = z[ih0];
- T.t = t[ih0];
-
- T.ResetErrors(1., 1., 1., 1., 100., (sta[ih0].timeInfo ? dt2[ih0] : 1.e6), 1.e2);
- T.C00 = cov[ih0].C00;
- T.C10 = cov[ih0].C10;
- T.C11 = cov[ih0].C11;
-
- T.NDF = ndf;
-
- for (int ih = NHits - 2; ih >= 0; --ih) {
- fit.Extrapolate(z[ih], fld);
- auto radThick = fParameters.GetMaterialThickness(ista[ih], T.x, T.y);
- fit.AddMsInMaterial(radThick);
- fit.EnergyLossCorrection(radThick, fvec(1.f));
- fit.FilterXY(cov[ih], x[ih], y[ih]);
- fit.FilterTime(t[ih], dt2[ih], sta[ih].timeInfo);
- }
- }
- } // for iiter
-
- //cout << " chi2 " << T3.chi2[i3_4] << " " << T.chi2[0] << endl;
- T.chi2 = (fscal) T3.chi2[i3_4]; //SG!!!
- T3.SetOneEntry(i3_4, T, i3_4);
-
- {
- int mc1 = GetMcTrackIdForCaHit(ihit[0]);
- int mc2 = GetMcTrackIdForCaHit(ihit[1]);
- int mc3 = GetMcTrackIdForCaHit(ihit[2]);
- if ((mc1 >= 0) && (mc1 == mc2) && (mc1 == mc3)) {
- const CbmL1MCTrack& mctr = CbmL1::Instance()->GetMcTracks()[mc1];
- float ev = 0;
- float chi2 = T.chi2[i3_4];
- float nd = T.NDF[i3_4];
- ca::tools::Debugger::Instance().FillNtuple("triplets", ev, mc1, istal, mctr.p, mctr.x, mctr.y, mctr.z, chi2,
- nd);
- }
- }
- } //i3
-
-} // findTripletsStep2
-
-
-inline void L1Algo::findTripletsStep3( // input
- Tindex n3, int istal, int istam, int istar, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3,
- L1Vector<L1HitIndex_t>& hitsm_3, L1Vector<L1HitIndex_t>& hitsr_3)
-{
- /// Selects good triplets and saves them into fTriplets.
- /// Finds neighbouring triplets at the next station.
- ///
- /// \param n3 Number of triplets to process
- /// \param istal index of the left station
- /// \param istam index of the middle station
- /// \param istar index of the right station
- /// \param T_3 fitted trajectories of triplets
- /// \param hitsl_3 index of the left triplet hit in Grid hits on the left station
- /// \param hitsm_3 index of the middle triplet hit in Grid hits on the middle station
- /// \param hitsr_3 index of the right triplet hit in Grid hits on the right station
-
-#ifdef _OPENMP
- unsigned int Thread = omp_get_thread_num();
-#else
- unsigned int Thread = 0;
-#endif
-
- L1HitIndex_t ihitl_prev = 0;
-
- const L1Station& stal = fParameters.GetStation(istal);
- const L1Station& stam = fParameters.GetStation(istam);
- const L1Station& star = fParameters.GetStation(istar);
-
- bool isMomentumFitted = ((stal.fieldStatus != 0) || (stam.fieldStatus != 0) || (star.fieldStatus != 0));
-
- for (Tindex i3 = 0; i3 < n3; ++i3) {
- const Tindex i3_V = i3 / fvec::size();
- const Tindex i3_4 = i3 % fvec::size();
-
- L1TrackPar& T3 = T_3[i3_V];
-
- //if (!T3.IsEntryConsistent(false, i3_4)) continue;
-
- // select
- fscal chi2 = T3.chi2[i3_4]; // / T3.NDF[i3_4];
-
- const L1HitIndex_t ihitl = hitsl_3[i3] + fGridHitStartIndex[istal];
- const L1HitIndex_t ihitm = hitsm_3[i3] + fGridHitStartIndex[istam];
- const L1HitIndex_t ihitr = hitsr_3[i3] + fGridHitStartIndex[istar];
- L1_ASSERT(ihitl < fGridHitStopIndex[istal], ihitl << " < " << fGridHitStopIndex[istal]);
- L1_ASSERT(ihitm < fGridHitStopIndex[istam], ihitm << " < " << fGridHitStopIndex[istam]);
- L1_ASSERT(ihitr < fGridHitStopIndex[istar], ihitr << " < " << fGridHitStopIndex[istar]);
-
- unsigned int Location = PackTripletId(istal, Thread, fTriplets[istal][Thread].size());
-
- if (ihitl_prev != 1 + hitsl_3[i3]) {
- // TODO: SG: A bug! The code doesn't work for iterations with missing hits.
- // TODO: Triplets with the same left hit are not placed continiously in fTriplets vector.
- // TODO: Urgent!!!
- // assert(fHitNtriplets[ihitl] == 0);
- fHitFirstTriplet[ihitl] = Location;
- fHitNtriplets[ihitl] = 0;
- }
-
- ihitl_prev = 1 + hitsl_3[i3];
-
- if (!fpCurrentIteration->GetTrackFromTripletsFlag()) {
- if (chi2 > fTripletFinalChi2Cut) { continue; }
- }
-
- // assert(std::isfinite(chi2) && chi2 > 0);
-
- if (!std::isfinite(chi2) || chi2 < 0) { continue; }
- //if (!T3.IsEntryConsistent(false, i3_4)) { continue; }
-
- fscal qp = T3.qp[i3_4];
-
- //TODO: why sqrt's? Wouldn't it be faster to skip sqrt() here and
- //TODO: compare the squared differences dqr*dqp later?
-
- fscal Cqp = T3.C44[i3_4];
-
- // TODO: SG: a magic correction that comes from the legacy code
- // removing it leads to a higher ghost ratio
- Cqp += 0.001;
-
- fTriplets[istal][Thread].emplace_back(ihitl, ihitm, ihitr, istal, istam, istar, 0, 0, 0, chi2, qp, Cqp, T3.tx[i3_4],
- T3.C22[i3_4], T3.ty[i3_4], T3.C33[i3_4]);
-
- L1Triplet& tr1 = fTriplets[istal][Thread].back();
- tr1.SetIsMomentumFitted(isMomentumFitted);
- tr1.SetLevel(0);
- tr1.SetFNeighbour(0);
- tr1.SetNNeighbours(0);
-
- fHitNtriplets[ihitl]++;
-
- if (istal > (fParameters.GetNstationsActive() - 4)) continue;
-
- unsigned int nNeighbours = fHitNtriplets[ihitm];
-
- unsigned int neighLocation = fHitFirstTriplet[ihitm];
- unsigned int neighStation = TripletId2Station(neighLocation);
- unsigned int neighThread = TripletId2Thread(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) {
-
- L1Triplet& curNeighbour = fTriplets[neighStation][neighThread][neighTriplet];
-
- fscal dchi2 = 0.;
- if (!checkTripletMatch(tr1, curNeighbour, dchi2)) continue;
-
- if (tr1.GetFNeighbour() == 0) tr1.SetFNeighbour(neighLocation);
-
- tr1.SetNNeighbours(neighLocation - tr1.GetFNeighbour() + 1);
-
- const unsigned char jlevel = curNeighbour.GetLevel();
-
- if (level <= jlevel) level = jlevel + 1;
- }
-
- tr1.SetLevel(level);
- }
-}
-
-/// Find neighbours of triplets. Calculate level of triplets.
-inline void L1Algo::f5( // input
- // output
-
- int* nlevel)
-{
- if (!fpCurrentIteration->GetTrackFromTripletsFlag()) {
- for (int istal = fParameters.GetNstationsActive() - 4; istal >= fFirstCAstation; istal--) {
- for (int tripType = 0; tripType < 3; tripType++) // tT = 0 - 123triplet, tT = 1 - 124triplet, tT = 2 - 134triplet
- {
- // All iters - jump
- if ((!fpCurrentIteration->GetJumpedFlag() && tripType != 0)
- || (fpCurrentIteration->GetJumpedFlag() && tripType != 0
- && (fParameters.GetNstationsActive() - 4 == istal))) {
- continue;
- }
-
- int istam = istal + 1;
- int istar = istal + 2;
- switch (tripType) {
- case 1: istar++; break;
- case 2:
- istam++;
- istar++;
- break;
- }
-
- for (Tindex iThread = 0; iThread < fNThreads; ++iThread) {
- for (Tindex itrip = 0; itrip < (Tindex) fTriplets[istal][iThread].size(); ++itrip) {
- L1Triplet& trip = fTriplets[istal][iThread][itrip];
- if (istam != trip.GetMSta()) continue;
- if (istar != trip.GetRSta()) continue;
-
- unsigned char level = 0;
- // fscal chi2 = trip->GetChi2();
-
- L1HitIndex_t ihitl = trip.GetLHit();
- L1HitIndex_t ihitm = trip.GetMHit();
- L1HitIndex_t ihitr = trip.GetRHit();
- L1_ASSERT(ihitl < fGridHitStopIndex[istal], ihitl << " < " << fGridHitStopIndex[istal]);
- L1_ASSERT(ihitm < fGridHitStopIndex[istam], ihitm << " < " << fGridHitStopIndex[istam]);
- L1_ASSERT(ihitr < fGridHitStopIndex[istar], ihitr << " < " << fGridHitStopIndex[istar]);
-
- L1Vector<unsigned int> neighCands("L1CATrackFinder::neighCands"); // save neighbour candidates
- neighCands.reserve(8); // ~average is 1-2 for central, up to 5
-
- unsigned int nNeighbours = fHitNtriplets[ihitm];
-
- unsigned int neighLocation = fHitFirstTriplet[ihitm];
- unsigned int neighStation = TripletId2Station(neighLocation);
- unsigned int neighThread = TripletId2Thread(neighLocation);
- unsigned int neighTriplet = TripletId2Triplet(neighLocation);
-
- for (unsigned int iN = 0; iN < nNeighbours; ++iN, ++neighTriplet, ++neighLocation) {
- // for (iN = first_triplet; iN <= last_triplet; ++iN){
- L1Triplet& neigh = fTriplets[neighStation][neighThread][neighTriplet];
- // if (neigh.GetMSta() != istar) continue; // neighbours should have 2 common hits
- // if (neigh.GetMHit() != ihitr) continue; //!!!
-
- fscal dchi2 = 0.;
- if (!checkTripletMatch(trip, neigh, dchi2)) continue;
-
- // calculate level
- unsigned char jlevel = neigh.GetLevel();
- if (level <= jlevel) level = jlevel + 1;
- if (level == jlevel + 1) neighCands.push_back(neighLocation);
- }
-
- // trip->neighbours.reset(0);
- // for (unsigned int in = 0; in < neighCands.size(); in++) {
- // int ID = neighCands[in];
- // unsigned int Station = TripletId2Station(ID);
- // unsigned int Thread = TripletId2Thread(ID);
- // unsigned int Triplet = TripletId2Triplet(ID);
- // const int nLevel = fTriplets[Station][Thread][Triplet].GetLevel();
- // if (level == nLevel + 1) trip->neighbours.push_back(Location);
- // }
- nlevel[level]++;
- } // vTriplets
- }
- } // tripType
- } // istal
- } // if !fpCurrentIteration->GetTrackFromTripletsFlag()
-}
-
-/// ------------------- doublets on station ----------------------
-
-inline void L1Algo::CreatePortionOfDoublets(
- /// input:
- const int istal, const int istam, const Tindex iSingletPortion, const Tindex singletPortionSize,
- /// output:
- L1TrackPar* T_1, L1FieldRegion* fld_1, L1HitIndex_t* hitsl_1, Tindex& n_2, L1Vector<L1HitIndex_t>& i1_2,
- L1Vector<L1HitIndex_t>& hitsm_2
- ///
-)
-{
- /// creates doublets
- /// input:
- /// @istal - start station number
- /// @istam - last station number
- /// @iSingletPortion - index of portion of left hits
- /// @singletPortionSize - number of left hits in the portion
- /// output:
- /// @*T_1 - singlet parameters
- /// @*fld_1 - field aproximation for singlets
- /// @*hitsl_1- left hits of future triplets
- /// @&n_2 - number of doublets
- /// @&i1_2 - index of 1st hit in portion indexed by doublet index
- /// @&hitsm_2 - index of middle hit in hits array indexed by doublet index
-
- if (istam < fParameters.GetNstationsActive()) {
- const L1Station& stal = fParameters.GetStation(istal);
- const L1Station& stam = fParameters.GetStation(istam);
-
- // prepare data
- L1HitPoint* vHits_l = &(fGridPoints[0]) + fGridHitStartIndex[istal];
- L1HitPoint* vHits_m = &(fGridPoints[0]) + fGridHitStartIndex[istam];
-
- fvec x[L1Constants::size::kSingletPortionSizeVec];
- fvec y[L1Constants::size::kSingletPortionSizeVec];
- fvec z[L1Constants::size::kSingletPortionSizeVec];
- fvec t[L1Constants::size::kSingletPortionSizeVec];
- fvec dx[L1Constants::size::kSingletPortionSizeVec];
- fvec dxy[L1Constants::size::kSingletPortionSizeVec];
- fvec dy[L1Constants::size::kSingletPortionSizeVec];
- fvec dt2[L1Constants::size::kSingletPortionSizeVec];
-
- /// prepare the portion of left hits data
-
- findSingletsStep0( // input
- iSingletPortion * L1Constants::size::kSingletPortionSize, singletPortionSize, vHits_l,
- // output
- hitsl_1, x, y, z, t, dx, dxy, dy, dt2);
-
- for (Tindex i = 0; i < singletPortionSize; ++i)
- L1_ASSERT(hitsl_1[i] < fGridHitStopIndex[istal] - fGridHitStartIndex[istal],
- hitsl_1[i] << " < " << fGridHitStopIndex[istal] - fGridHitStartIndex[istal]);
-
- Tindex portionSize_V = (singletPortionSize + fvec::size() - 1) / fvec::size();
-
- /// Get the field approximation. Add the target to parameters estimation. Propagaete to middle station.
-
- findSingletsStep1(istal, istam, portionSize_V, x, y, z, t, dx, dxy, dy, dt2,
- // output
- T_1, fld_1);
-
- /// Find the doublets. Reformat data in the portion of doublets.
-
- // TODO: repalce with constexpr if (C++17) (S.Zharko)
-#ifdef DOUB_PERFORMANCE
- L1Vector<L1HitIndex_t> hitsl_2("L1CATrackFinder::hitsl_2");
-#endif // DOUB_PERFORMANCE
-
- findDoubletsStep0( // input
- singletPortionSize, stal, vHits_l, stam, vHits_m, T_1, hitsl_1,
- // output
- n_2, i1_2,
-#ifdef DOUB_PERFORMANCE
- hitsl_2,
-#endif // DOUB_PERFORMANCE
- hitsm_2);
-
- for (Tindex i = 0; i < static_cast<Tindex>(hitsm_2.size()); ++i)
- L1_ASSERT(hitsm_2[i] < fGridHitStopIndex[istam] - fGridHitStartIndex[istam],
- hitsm_2[i] << " " << fGridHitStopIndex[istam] - fGridHitStartIndex[istam]);
-
-#ifdef DOUB_PERFORMANCE
- L1HitIndex_t* RealIHitL = &(fGridHitIds[fGridHitStartIndex[istal]]);
- L1HitIndex_t* RealIHitM = &(fGridHitIds[fGridHitStartIndex[istam]]);
- for (Tindex i = 0; i < n2; ++i) {
- // int i_4 = i%4;
- // int i_V = i/4;
- L1HitIndex_t iHits[2] = {RealIHitL[hitsl_2[i]], RealIHitM[hitsm_2[i]]};
- fL1Eff_doublets->AddOne(iHits);
- }
-#endif // DOUB_PERFORMANCE
- }
-}
-
-
-/// ------------------- Triplets on station ----------------------
-
-inline void L1Algo::CreatePortionOfTriplets(
- /// input
- int istal, int istam, int istar,
-
- /// input / output
- L1TrackPar* T_1, L1FieldRegion* fld_1, L1HitIndex_t* hitsl_1,
-
- Tindex& n_2, L1Vector<L1HitIndex_t>& i1_2, L1Vector<L1HitIndex_t>& hitsm_2)
-{
-
- /// creates a portion of triplets:
- /// input:
- /// @istal - left station number,
- /// @istam - middle station number,
- /// @istar - right station number,
-
- /// input / output:
-
- /// @*T_1 - track parameters for singlets,
- /// @*fld_1 - field approximation for singlets,
- /// @&n_2 - number of doublets in portion
- /// @&n_2 - number of doublets,
- /// @&i1_2 - index of 1st hit in portion indexed by doublet index,
- /// @&hitsm_2 - index of middle hit in hits array indexed by doublet index
- /// output:
- // @*vTriplets_part - array of triplets,
- // @*TripStartIndexH,
- // @*TripStopIndexH - start/stop index of a triplet in the array
-
- if (istar < fParameters.GetNstationsActive()) {
- // prepare data
- const L1Station& star = fParameters.GetStation(istar);
-
- L1HitPoint* vHits_m = &(fGridPoints[0]) + fGridHitStartIndex[istam];
-
- L1HitPoint* vHits_r = 0;
- vHits_r = &(fGridPoints[0]) + fGridHitStartIndex[istar];
-
- Tindex n3 = 0;
-
- /// Add the middle hits to parameters estimation. Propagate to right station.
-
-
-#ifdef _OPENMP
- int Thread = omp_get_thread_num();
-#else
- int Thread = 0;
-#endif
-
- L1Vector<L1TrackPar>& T_3 = fTripletPar[Thread];
- L1Vector<L1HitIndex_t>& hitsl_3 = fTripletHitsL[Thread];
- L1Vector<L1HitIndex_t>& hitsm_3 = fTripletHitsM[Thread];
- L1Vector<L1HitIndex_t>& hitsr_3 = fTripletHitsR[Thread];
-
- L1Vector<fvec>& x_3 = fTripletHitR_X[Thread];
- L1Vector<fvec>& y_3 = fTripletHitR_Y[Thread];
- L1Vector<fvec>& z_3 = fTripletHitR_Z[Thread];
- L1Vector<fvec>& t_3 = fTripletHitR_Time[Thread];
-
- L1Vector<fvec>& dx_3 = fTripletHitR_dX[Thread];
- L1Vector<fvec>& dxy_3 = fTripletHitR_dXY[Thread];
- L1Vector<fvec>& dy_3 = fTripletHitR_dY[Thread];
- L1Vector<fvec>& dt2_3 = fTripletHitR_TimeErr[Thread];
-
- T_3.clear();
- hitsl_3.clear();
- hitsm_3.clear();
- hitsr_3.clear();
-
- x_3.clear();
- y_3.clear();
- z_3.clear();
- t_3.clear();
-
- dx_3.clear();
- dxy_3.clear();
- dy_3.clear();
- dt2_3.clear();
-
- /// Find the triplets(right hit). Reformat data in the portion of triplets.
-
- findTripletsStep0( // input
- vHits_r,
-
- istal, istam, istar, vHits_m, T_1, fld_1, hitsl_1,
-
- n_2, hitsm_2, i1_2,
-
- // output
- n3, T_3, hitsl_3, hitsm_3, hitsr_3, x_3, y_3, z_3, t_3, dx_3, dxy_3, dy_3, dt2_3);
-
-
- for (Tindex i = 0; i < static_cast<Tindex>(hitsl_3.size()); ++i)
- L1_assert(hitsl_3[i] < fGridHitStopIndex[istal] - fGridHitStartIndex[istal]);
- for (Tindex i = 0; i < static_cast<Tindex>(hitsm_3.size()); ++i)
- L1_assert(hitsm_3[i] < fGridHitStopIndex[istam] - fGridHitStartIndex[istam]);
- for (Tindex i = 0; i < static_cast<Tindex>(hitsr_3.size()); ++i)
- L1_assert(hitsr_3[i] < fGridHitStopIndex[istar] - fGridHitStartIndex[istar]);
-
- /// Add the right hits to parameters estimation.
-
-
- Tindex n3_V = (n3 + fvec::size() - 1) / fvec::size();
- findTripletsStep1( // input
- n3_V, star, x_3, y_3, z_3, t_3, dx_3, dxy_3, dy_3, dt2_3,
- // output
- T_3);
-
-
- /// refit
- findTripletsStep2(n3, istal, istam, istar, T_3, hitsl_3, hitsm_3, hitsr_3, 1);
-
-#ifdef TRIP_PERFORMANCE
- L1HitIndex_t* RealIHitL = &(fGridHitIds[fGridHitStartIndex[istal]]);
- L1HitIndex_t* RealIHitM = &(fGridHitIds[fGridHitStartIndex[istam]]);
- L1HitIndex_t* RealIHitR = &(fGridHitIds[fGridHitStartIndex[istar]]);
- for (Tindex i = 0; i < n3; ++i) {
- Tindex i_4 = i % 4;
- Tindex i_V = i / 4;
- L1HitIndex_t iHits[3] = {RealIHitL[hitsl_3[i]], RealIHitM[hitsm_3[i]], RealIHitR[hitsr_3[i]]};
-#ifdef PULLS
- if (fL1Eff_triplets->AddOne(iHits)) fL1Pulls->AddOne(T_3[i_V], i_4, iHits[2]);
-#else
- fL1Eff_triplets->AddOne(iHits);
-#endif
- if (T_3[i_V].chi2[i_4] < fTripletFinalChi2Cut) { fL1Eff_triplets2->AddOne(iHits); }
- }
-#endif // TRIP_PERFORMANCE
-
-
- /// Fill Triplets.
- findTripletsStep3( // input
- n3, istal, istam, istar, T_3, hitsl_3, hitsm_3, hitsr_3);
- }
-}
-
-
-// hitCheck::hitCheck()
-// {
-// omp_init_lock(&Occupied);
-// trackCandidateIndex = -1;
-// UsedByTrack=0;
-// Chi2Track = 100000000;
-// Length = 0;
-// ista = 1000;
-//
-// }
-// hitCheck::~hitCheck()
-// {
-// omp_destroy_lock(&Occupied);
-// }
// **************************************************************************************************
@@ -1761,142 +412,81 @@ void L1Algo::CaTrackFinderSlice()
}
- {
- /// possible left hits of triplets are splited in portions of 16 (4 SIMDs) to use memory faster
-
- for (int istal = fParameters.GetNstationsActive() - 2; istal >= fFirstCAstation;
- istal--) { //start downstream chambers
- fSingletPortionSize[istal].clear();
- int NHits_l = fGridHitStopIndex[istal] - fGridHitStartIndex[istal];
- int nPortions = NHits_l / L1Constants::size::kSingletPortionSize;
- int rest = NHits_l - nPortions * L1Constants::size::kSingletPortionSize;
- for (int ipp = 0; ipp < nPortions; ipp++) {
- fSingletPortionSize[istal].push_back(L1Constants::size::kSingletPortionSize);
- } // start_lh - portion of left hits
- if (rest > 0) fSingletPortionSize[istal].push_back(rest);
- } // lstations
-
-
-#ifdef COUNTERS
- stat_nSinglets[isec] += nSinglets;
-#endif
- }
-
- /* {
- /// possible left hits of triplets are splited in portions of 16 (4 SIMDs) to use memory faster
- fDoubletPortionStopIndex[fParameters.GetNstationsActive() - 1] = 0;
- unsigned int ip = 0; //index of curent portion
-
- for (int istal = fParameters.GetNstationsActive() - 2; istal >= fFirstCAstation; istal--) //start downstream chambers
- {
- int nHits = fGridHitStopIndex[istal] - fGridHitStartIndex[istal];
-
- int NHits_P = nHits/L1Constants::size::kSingletPortionSize;
-
- for( int ipp = 0; ipp < NHits_P; ipp++ )
- {
- n_g1[ip] = L1Constants::size::kSingletPortionSize;
- ip++;
- } // start_lh - portion of left hits
+ /// stage for triplets creation
- n_g1[ip] = nHits - NHits_P * L1Constants::size::kSingletPortionSize;
+ L1TripletConstructor constructor1(this);
+ L1TripletConstructor constructor2(this);
+ L1TripletConstructor constructor3(this);
- ip++;
- fDoubletPortionStopIndex[istal] = ip;
- }// lstations
- // nPortions = ip;
- } */
+ for (int istal = fParameters.GetNstationsActive() - 2; istal >= fFirstCAstation;
+ istal--) { // start with downstream chambers
+ Tindex NHits_l = fGridHitStopIndex[istal] - fGridHitStartIndex[istal];
+ for (Tindex ih = 0; ih < NHits_l; ++ih) {
+ constructor1.CreateTripletsForHit(istal, istal + 1, istal + 2, ih);
+ if (fpCurrentIteration->GetJumpedFlag()) {
+ constructor2.CreateTripletsForHit(istal, istal + 1, istal + 3, ih);
+ constructor3.CreateTripletsForHit(istal, istal + 2, istal + 3, ih);
+ }
- /// stage for triplets creation
+ auto oldSize = fTriplets[istal][0].size();
-#ifdef XXX
- TStopwatch c_timer;
- c_timer.Start();
-#endif
+ fTriplets[istal][0].insert(fTriplets[istal][0].end(), constructor1.GetTriplets().begin(),
+ constructor1.GetTriplets().end());
+ fTriplets[istal][0].insert(fTriplets[istal][0].end(), constructor2.GetTriplets().begin(),
+ constructor2.GetTriplets().end());
+ fTriplets[istal][0].insert(fTriplets[istal][0].end(), constructor3.GetTriplets().begin(),
+ constructor3.GetTriplets().end());
+ const L1HitIndex_t ihitl = ih + fGridHitStartIndex[istal];
+ fHitFirstTriplet[ihitl] = PackTripletId(istal, 0, oldSize);
+ fHitNtriplets[ihitl] = fTriplets[istal][0].size() - oldSize;
+ }
+ } // istal
- L1TrackPar T_1[L1Constants::size::kSingletPortionSizeVec];
- L1FieldRegion fld_1[L1Constants::size::kSingletPortionSizeVec];
- L1HitIndex_t hitsl_1[L1Constants::size::kSingletPortionSize];
- L1TrackPar TG_1[L1Constants::size::kSingletPortionSizeVec];
- L1FieldRegion fldG_1[L1Constants::size::kSingletPortionSizeVec];
- L1HitIndex_t hitslG_1[L1Constants::size::kSingletPortionSize];
- /// middle hits indexed by number of doublets in portion(i2)
- L1Vector<L1HitIndex_t> hitsm_2("L1CATrackFinder::hitsm_2");
+ // search for neighbouring triplets
- /// index in portion of singlets(i1) indexed by index in portion of doublets(i2)
- L1Vector<L1HitIndex_t> i1_2("L1CATrackFinder::i1_2");
+ for (int istal = fParameters.GetNstationsActive() - 2; istal >= fFirstCAstation;
+ istal--) { // start with downstream chambers
- /// middle hits indexed by number of doublets in portion(i2)
- L1Vector<L1HitIndex_t> hitsmG_2("L1CATrackFinder::hitsmG_2");
+ for (unsigned int it = 0; it < fTriplets[istal][0].size(); ++it) {
- /// index in portion of singlets(i1) indexed by index in portion of doublets(i2)
- L1Vector<L1HitIndex_t> i1G_2("L1CATrackFinder::i1G_2");
+ L1Triplet& tr = fTriplets[istal][0][it];
+ tr.SetLevel(0);
+ tr.SetFNeighbour(0);
+ tr.SetNNeighbours(0);
- hitsm_2.reserve(9000); // TODO: make reasonable cuts on n combinations, put them to the header
- i1_2.reserve(9000); // TODO: why that large numbers are needed even for mbias??? something goes wrong sometimes..
- hitsmG_2.reserve(9000);
- i1G_2.reserve(9000);
+ if (istal > (fParameters.GetNstationsActive() - 4)) break;
- for (int istal = fParameters.GetNstationsActive() - 2; istal >= fFirstCAstation;
- istal--) // start downstream chambers
- {
+ unsigned int nNeighbours = fHitNtriplets[tr.GetMHit()];
-#ifdef _OPENMP
-#pragma omp parallel for firstprivate(T_1, fld_1, hitsl_1, hitsm_2, i1_2, TG_1, fldG_1, hitslG_1, hitsmG_2, \
- i1G_2) //schedule(dynamic, 2)
-#endif
- for (Tindex ip = 0; ip < (Tindex) fSingletPortionSize[istal].size(); ++ip) {
- Tindex n_2 = 0; /// number of doublets in portion
+ unsigned int neighLocation = fHitFirstTriplet[tr.GetMHit()];
+ unsigned int neighStation = TripletId2Station(neighLocation);
+ unsigned int neighThread = TripletId2Thread(neighLocation);
+ unsigned int neighTriplet = TripletId2Triplet(neighLocation);
- hitsm_2.clear();
- i1_2.clear();
+ if (nNeighbours > 0) { assert((int) neighStation == istal + 1 || (int) neighStation == istal + 2); }
- CreatePortionOfDoublets(istal, istal + 1, ip, fSingletPortionSize[istal][ip],
- // output
- T_1, fld_1, hitsl_1, n_2, i1_2, hitsm_2);
+ unsigned char level = 0;
+ for (unsigned int iN = 0; iN < nNeighbours; ++iN, ++neighTriplet, ++neighLocation) {
- CreatePortionOfTriplets( // input
- istal, istal + 1, istal + 2, T_1, fld_1, hitsl_1, n_2, i1_2, hitsm_2);
+ L1Triplet& neighbour = fTriplets[neighStation][neighThread][neighTriplet];
- if (fpCurrentIteration->GetJumpedFlag()) {
- // All iterations are "jump"!
- Tindex nG_2;
- hitsmG_2.clear();
- i1G_2.clear();
+ fscal dchi2 = 0.;
+ if (!checkTripletMatch(tr, neighbour, dchi2)) continue;
- CreatePortionOfDoublets( // input
- istal, istal + 2, ip, fSingletPortionSize[istal][ip],
- // output
- TG_1, fldG_1, hitslG_1, nG_2, i1G_2, hitsmG_2);
+ if (tr.GetFNeighbour() == 0) { tr.SetFNeighbour(neighLocation); }
- CreatePortionOfTriplets( // input
- istal, istal + 1, istal + 3, T_1, fld_1, hitsl_1,
+ tr.SetNNeighbours(neighLocation - tr.GetFNeighbour() + 1);
- n_2, i1_2, hitsm_2);
+ const unsigned char neighbourLevel = neighbour.GetLevel();
- CreatePortionOfTriplets( // input
- istal, istal + 2, istal + 3, TG_1, fldG_1, hitslG_1, nG_2, i1G_2, hitsmG_2);
+ if (level <= neighbourLevel) { level = neighbourLevel + 1; }
}
- } //
- }
-
- // int nlevels[L1Constants::size::kMaxNstations]; // number of triplets with some number of neighbours.
- // for (int il = 0; il < fParameters.GetNstationsActive(); ++il) nlevels[il] = 0;
- //
- // f5( // input
- // // output
- // 0,
- // nlevels
- // );
-
-#ifdef XXX
- c_timer.Stop();
- ti[isec]["tripl1"] = c_timer;
- c_timer.Start();
-#endif
+ tr.SetLevel(level);
+ } // neighbour search
+ } // istal
///====================================================================
///= =
@@ -2408,8 +998,8 @@ void L1Algo::CaTrackFinderSlice()
#ifdef TRIP_PERFORMANCE
fL1Eff_triplets->CalculateEff();
fL1Eff_triplets->Print("Triplet performance", 1);
- // fL1Eff_triplets2->CalculateEff();
- // fL1Eff_triplets2->Print("Triplet performance. After chi2 cut");
+// fL1Eff_triplets2->CalculateEff();
+// fL1Eff_triplets2->Print("Triplet performance. After chi2 cut");
#endif
}
@@ -2423,9 +1013,9 @@ void L1Algo::CaTrackFinderSlice()
* *
****************************************************************/
-inline void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best_L, fscal& best_chi2,
- const L1Triplet* curr_trip, L1Branch& curr_tr, unsigned char& curr_L, fscal& curr_chi2,
- unsigned char min_best_l, L1Branch* new_tr)
+void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best_L, fscal& best_chi2,
+ const L1Triplet* curr_trip, L1Branch& curr_tr, unsigned char& curr_L, fscal& curr_chi2,
+ unsigned char min_best_l, L1Branch* new_tr)
/// recursive search for tracks
/// input: @ista - station index, @&best_tr - best track for the privious call, @&best_L -
/// output: @&NCalls - number of function calls
diff --git a/reco/L1/L1Algo/L1Constants.h b/reco/L1/L1Algo/L1Constants.h
index eed01d2eec..e228cbea1e 100644
--- a/reco/L1/L1Algo/L1Constants.h
+++ b/reco/L1/L1Algo/L1Constants.h
@@ -49,8 +49,6 @@ namespace L1Constants
// TODO: Clarify the meaning of these coefficients
constexpr int kCoeff = 64 / 4; ///< TODO:
- constexpr int kSingletPortionSize = 1024 / kCoeff; ///< portion of left hits
- constexpr int kSingletPortionSizeVec = 1024 / kCoeff / fvec::size(); ///< portion of left hits per one vector word
constexpr int kMaxPortionDoublets = 10000 / 5 * 64 / 2 / kCoeff; ///< Max size of the doublets portion
constexpr int kMaxPortionTriplets = 10000 * 5 * 64 / 2 / kCoeff; ///< Max size of the triplets portion
constexpr int kMaxPortionTripletsP = kMaxPortionTriplets / fvec::size(); ///< Max size of the triplets portion
diff --git a/reco/L1/L1Algo/L1Field.cxx b/reco/L1/L1Algo/L1Field.cxx
index 57172997aa..c50562a226 100644
--- a/reco/L1/L1Algo/L1Field.cxx
+++ b/reco/L1/L1Algo/L1Field.cxx
@@ -11,6 +11,7 @@
#include <iostream>
#include <sstream>
+#include "L1TrackPar.h"
#include "L1Utils.h"
//
@@ -60,10 +61,11 @@ std::ostream& operator<<(std::ostream& out, const L1FieldValue& B)
L1FieldSlice::L1FieldSlice()
{
for (int i = 0; i < L1Constants::size::kMaxNFieldApproxCoefficients; ++i) {
- cx[i] = L1Utils::kNaN;
- cy[i] = L1Utils::kNaN;
- cz[i] = L1Utils::kNaN;
+ cx[i] = undef::kFvc;
+ cy[i] = undef::kFvc;
+ cz[i] = undef::kFvc;
}
+ z = undef::kFvc;
}
//----------------------------------------------------------------------------------------------------------------------
@@ -76,6 +78,7 @@ void L1FieldSlice::CheckConsistency() const
L1Utils::CheckSimdVectorEquality(cy[i], "L1FieldSlice: cy");
L1Utils::CheckSimdVectorEquality(cz[i], "L1FieldSlice: cz");
}
+ L1Utils::CheckSimdVectorEquality(z, "L1FieldSlice: z");
}
//----------------------------------------------------------------------------------------------------------------------
@@ -117,6 +120,12 @@ void L1FieldSlice::GetFieldValue(const fvec& x, const fvec& y, L1FieldValue& B)
+ cz[16] * x4y + cz[17] * x3y2 + cz[18] * x2y3 + cz[19] * xy4 + cz[20] * y5;
}
+void L1FieldSlice::GetFieldValueForLine(const L1TrackPar& t, L1FieldValue& B) const
+{
+ fvec dz = z - t.z;
+ GetFieldValue(t.x + t.tx * dz, t.y + t.ty * dz, B);
+}
+
//----------------------------------------------------------------------------------------------------------------------
//
std::string L1FieldSlice::ToString(int indentLevel) const
diff --git a/reco/L1/L1Algo/L1Field.h b/reco/L1/L1Algo/L1Field.h
index 086dca9c1c..e8069ce082 100644
--- a/reco/L1/L1Algo/L1Field.h
+++ b/reco/L1/L1Algo/L1Field.h
@@ -10,6 +10,9 @@
#include "L1Constants.h"
#include "L1Def.h"
#include "L1SimdSerializer.h"
+#include "L1Undef.h"
+
+class L1TrackPar;
class L1FieldValue {
public:
@@ -76,6 +79,11 @@ public:
/// \param B the L1FieldValue output
void GetFieldValue(const fvec& x, const fvec& y, L1FieldValue& B) const;
+ /// Gets field value for the intersection with a straight track
+ /// \param t straight track
+ /// \param B the L1FieldValue output
+ void GetFieldValueForLine(const L1TrackPar& t, L1FieldValue& B) const;
+
/// String representation of class contents
/// \param indentLevel number of indent characters in the output
std::string ToString(int indentLevel = 0) const;
@@ -91,6 +99,8 @@ public:
fvec cz
[L1Constants::size::kMaxNFieldApproxCoefficients]; ///< Polynomial coefficients for z-component of the field value
+ fvec z {undef::kFvc}; ///< z coordinate of the slice
+
/// Serialization function
friend class boost::serialization::access;
template<class Archive>
@@ -99,6 +109,7 @@ public:
ar& cx;
ar& cy;
ar& cz;
+ ar& z;
}
} _fvecalignment;
diff --git a/reco/L1/L1Algo/L1Portion.h b/reco/L1/L1Algo/L1Portion.h
deleted file mode 100644
index 56c04ce973..0000000000
--- a/reco/L1/L1Algo/L1Portion.h
+++ /dev/null
@@ -1,166 +0,0 @@
-/* Copyright (C) 2010-2021 Frankfurt Institute for Advanced Studies, Goethe-Universität Frankfurt, Frankfurt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Igor Kulakov [committer] */
-
-#ifndef L1Portion_H
-#define L1Portion_H
-
-#include <vector>
-
-#include "L1TrackPar.h"
-
-using std::vector;
-// class for portions organization
-
-template<typename T>
-class L1Portion;
-
-template<>
-class L1Portion<L1TrackPar> {
-public:
- typedef L1TrackPar T;
- typedef L1Vector<T> vType;
-
- L1Portion() {};
- L1Portion(int size) { reserve(size); };
- L1Portion(int size, int size2) : dataSize(size2)
- {
- reserve(size);
- // reserve2(size2);
- };
- vType& operator[](int i) { return a[i]; }
- void resize(int size) { a.resize(size); };
- void reserve(int size) { a.reserve(size); };
- void reserve2(int size)
- {
- for (unsigned int i = 0; i < a.size(); i++)
- a[i].reserve(size);
- };
- void push_back(vType& v) { a.push_back(v); };
- void add_void()
- {
- vType v("L1Portion<L1TrackPar>::add_void");
- // v.resize(dataSize);
- a.push_back(v);
- a[a.size() - 1].reserve(dataSize);
- };
- void add_void(int i)
- {
- T t;
- a[i].push_back(t);
- };
-
- int CalcSize()
- {
- int size = 0;
- for (unsigned int i = 0; i < a.size(); i++)
- size += a[i].size();
- return size * sizeof(T);
- };
-
-private:
- vector<vType> a {"L1Portion<L1TrackPar>::a"};
- // int mainSize; // size of a
- int dataSize {0}; // size of vType
-};
-
-template<>
-class L1Portion<L1FieldRegion> {
-public:
- typedef L1FieldRegion T;
- typedef L1Vector<T> vType;
-
- L1Portion() {};
- L1Portion(int size) { reserve(size); };
- L1Portion(int size, int size2) : dataSize(size2)
- {
- reserve(size);
- // reserve2(size2);
- };
-
- vType& operator[](int i) { return a[i]; }
- void resize(int size) { a.resize(size); };
- void reserve(int size) { a.reserve(size); };
- void reserve2(int size)
- {
- for (unsigned int i = 0; i < a.size(); i++)
- a[i].reserve(size);
- };
- void push_back(vType& v) { a.push_back(v); };
- void add_void()
- {
- vType v("L1Portion<L1FieldRegion>::add_void");
- // v.resize(dataSize);
- a.push_back(v);
- a[a.size() - 1].reserve(dataSize);
- };
- void add_void(int i)
- {
- T t;
- a[i].push_back(t);
- };
-
- int CalcSize()
- {
- int size = 0;
- for (unsigned int i = 0; i < a.size(); i++)
- size += a[i].size();
- return size * sizeof(T);
- };
-
-private:
- vector<vType> a {"L1Portion<L1FieldRegion>::a"};
- // int mainSize; // size of a
- int dataSize {0}; // size of vType
-};
-
-
-template<typename T>
-class L1Portion {
-public:
- typedef vector<T> vType;
-
- L1Portion() {};
- L1Portion(int size) { reserve(size); };
- L1Portion(int size, int size2) : dataSize(size2)
- {
- reserve(size);
- // reserve2(size2);
- };
- vType& operator[](int i) { return a[i]; }
- void resize(int size) { a.resize(size); };
- void reserve(int size) { a.reserve(size); };
- void reserve2(int size)
- {
- for (int i = 0; i < a.size(); i++)
- a[i].reserve(size);
- };
- void push_back(vType& v) { a.push_back(v); };
- void add_void()
- {
- vType v("L1Portion<T>::add_void");
- // v.resize(dataSize);
- a.push_back(v);
- a[a.size() - 1].reserve(dataSize);
- };
- void add_void(int i)
- {
- T t;
- a[i].push_back(t);
- };
-
- int CalcSize()
- {
- int size = 0;
- for (int i = 0; i < a.size(); i++)
- size += a[i].size();
- return size * sizeof(T);
- };
-
-private:
- vector<vType> a {"L1Portion<T>::a"};
- // int mainSize; // size of a
- int dataSize {0}; // size of vType
-};
-
-#endif // L1Portion_H
diff --git a/reco/L1/L1Algo/L1Triplet.h b/reco/L1/L1Algo/L1Triplet.h
index 0e0f4adb1f..4731592963 100644
--- a/reco/L1/L1Algo/L1Triplet.h
+++ b/reco/L1/L1Algo/L1Triplet.h
@@ -22,7 +22,7 @@ public:
/// constructor
L1Triplet(unsigned int iHitL, unsigned int iHitM, unsigned int iHitR, unsigned int iStaL, unsigned int iStaM,
unsigned int iStaR, unsigned char Level, unsigned int firstNeighbour, char nNeighbours, fscal Chi2,
- fscal Qp, fscal Cqp, fscal tx, fscal Ctx, fscal ty, fscal Cty)
+ fscal Qp, fscal Cqp, fscal tx, fscal Ctx, fscal ty, fscal Cty, bool isMomentumFitted)
: fChi2(Chi2)
, fQp(Qp)
, fCqp(Cqp)
@@ -37,6 +37,7 @@ public:
, fNneighbours(nNeighbours)
, fLevel(Level)
, fSta((iStaL << 4) + ((iStaM - iStaL - 1) << 2) + (iStaR - iStaL - 2))
+ , fIsMomentumFitted(isMomentumFitted)
{
}
diff --git a/reco/L1/L1Algo/L1TripletConstructor.cxx b/reco/L1/L1Algo/L1TripletConstructor.cxx
new file mode 100644
index 0000000000..96f4b28f62
--- /dev/null
+++ b/reco/L1/L1Algo/L1TripletConstructor.cxx
@@ -0,0 +1,690 @@
+/* Copyright (C) 2010-2021 Frankfurt Institute for Advanced Studies, Goethe-Universität Frankfurt, Frankfurt
+ 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 <algorithm>
+#include <iostream>
+
+#include "CaToolsDebugger.h"
+#include "L1Algo.h"
+#include "L1Assert.h"
+#include "L1Fit.h"
+#include "L1Grid.h"
+#include "L1HitArea.h"
+#include "L1HitPoint.h"
+#include "L1TrackPar.h"
+
+L1TripletConstructor::L1TripletConstructor(L1Algo* algo) { fAlgo = algo; }
+
+void L1TripletConstructor::InitStations(int istal, int istam, int istar)
+{
+ fIsInitialized = false;
+
+ fIstaL = istal;
+ fIstaM = istam;
+ fIstaR = istar;
+
+ if (fIstaM >= fAlgo->fParameters.GetNstationsActive()) { return; }
+
+ fStaL = &fAlgo->fParameters.GetStation(fIstaL);
+ fStaM = &fAlgo->fParameters.GetStation(fIstaM);
+ fStaR = &fAlgo->fParameters.GetStation(fIstaR);
+
+ fHitsL = &(fAlgo->fGridPoints[0]) + fAlgo->fGridHitStartIndex[fIstaL];
+ fHitsM = &(fAlgo->fGridPoints[0]) + fAlgo->fGridHitStartIndex[fIstaM];
+ fHitsR = &(fAlgo->fGridPoints[0]) + fAlgo->fGridHitStartIndex[fIstaR];
+
+ { // two stations for approximating the field between the target and the left hit
+ int sta1 = fIstaL;
+ if (sta1 >= fAlgo->fNfieldStations) { sta1 = fAlgo->fNfieldStations - 1; }
+ if (sta1 < 1) { sta1 = 1; }
+ int sta0 = sta1 / 2; // station between fIstaL and the target
+
+ assert(0 <= sta0 && sta0 < sta1 && sta1 < fAlgo->fParameters.GetNstationsActive());
+
+ fFld0Sta[0] = &fAlgo->fParameters.GetStation(sta0);
+ fFld0Sta[1] = &fAlgo->fParameters.GetStation(sta1);
+ }
+
+ { // three stations for approximating the field between the left and the right hit
+
+ int sta0 = fIstaL;
+ int sta1 = fIstaM;
+ int sta2 = fIstaM + 1;
+ if (sta2 >= fAlgo->fParameters.GetNstationsActive()) {
+ sta2 = fIstaM;
+ sta1 = fIstaM - 1;
+ if (sta0 >= sta1) { sta0 = sta1 - 1; }
+ }
+
+ assert(0 <= sta0 && sta0 < sta1 && sta1 < sta2 && sta2 < fAlgo->fParameters.GetNstationsActive());
+
+ fFld1Sta[0] = &fAlgo->fParameters.GetStation(sta0);
+ fFld1Sta[1] = &fAlgo->fParameters.GetStation(sta1);
+ fFld1Sta[2] = &fAlgo->fParameters.GetStation(sta2);
+ }
+
+ fFit.SetParticleMass(fAlgo->GetDefaultParticleMass());
+ if (fAlgo->fpCurrentIteration->GetElectronFlag()) { fFit.SetParticleMass(L1Constants::phys::kElectronMass); }
+ fFit.SetMask(fmask::One());
+
+ fIsInitialized = true;
+}
+
+
+void L1TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, L1HitIndex_t ihl)
+{
+
+ InitStations(istal, istam, istar);
+
+ if (!fIsInitialized) return;
+
+ fIhitL = ihl;
+
+ L1HitPoint& hitl = fHitsL[fIhitL];
+
+ // fit a straight line through the target and the left hit.
+ {
+
+ /// Get the field approximation. Add the target to parameters estimation.
+ /// Propagaete to the middle station of a triplet.
+
+ fFit.SetQp0(fAlgo->fMaxInvMom);
+
+ L1FieldValue lB, mB, cB, rB _fvecalignment;
+ L1FieldValue l_B_global, centB_global _fvecalignment;
+
+ // get the magnetic field along the track.
+ // (suppose track is straight line with origin in the target)
+
+ fvec xl = hitl.X();
+ fvec yl = hitl.Y();
+ fvec zl = hitl.Z();
+ fvec time = hitl.T();
+ fvec timeEr2 = hitl.dT2();
+ const fvec dzli = 1.f / (zl - fAlgo->fTargZ);
+
+ const fvec tx = (xl - fAlgo->fTargX) * dzli;
+ const fvec ty = (yl - fAlgo->fTargY) * dzli;
+
+ L1TrackPar& T = fFit.Tr();
+
+ T.x = xl;
+ T.y = yl;
+ T.z = zl;
+ T.tx = tx;
+ T.ty = ty;
+ T.qp = fvec(0.);
+ T.t = time;
+ T.vi = fvec(0.);
+
+ fvec txErr2 = fAlgo->fMaxSlopePV * fAlgo->fMaxSlopePV / fvec(9.);
+ fvec qpErr2 = fAlgo->fMaxInvMom * fAlgo->fMaxInvMom / fvec(9.);
+
+ T.ResetErrors(1., 1., txErr2, txErr2, qpErr2, (staL().timeInfo ? timeEr2 : 1.e6), 1.e2);
+ T.C00 = hitl.dX2();
+ T.C10 = hitl.dXY();
+ T.C11 = hitl.dY2();
+
+ T.chi2 = 0.;
+ T.NDF = (fAlgo->fpCurrentIteration->GetPrimaryFlag()) ? fvec(2.) : fvec(0.);
+ T.nTimeMeasurements = (staL().timeInfo ? 1 : 0);
+
+ // TODO: iteration parameter: "Starting NDF of track parameters"
+ // NDF = number of track parameters (6: x, y, tx, ty, qp, time) - number of measured parameters (3: x, y, time) on station or (2: x, y) on target
+
+
+ // field made by the left hit, the target and the station istac in-between.
+ // is used for extrapolation to the target and to the middle hit
+ L1FieldRegion fld0;
+ {
+ L1FieldValue B0, B1;
+ fFld0Sta[0]->fieldSlice.GetFieldValueForLine(T, B0);
+ fFld0Sta[1]->fieldSlice.GetFieldValueForLine(T, B1);
+ fld0.Set(fAlgo->fTargB, fAlgo->fTargZ, B0, fFld0Sta[0]->fZ, B1, fFld0Sta[1]->fZ);
+ }
+
+ { // field, made by the left hit, the middle station and the right station
+ // Will be used for extrapolation to the right hit
+ L1FieldValue B0, B1, B2;
+ fFld1Sta[0]->fieldSlice.GetFieldValueForLine(T, B0);
+ fFld1Sta[1]->fieldSlice.GetFieldValueForLine(T, B1);
+ fFld1Sta[2]->fieldSlice.GetFieldValueForLine(T, B2);
+ fFldL.Set(B0, fFld1Sta[0]->fZ, B1, fFld1Sta[1]->fZ, B2, fFld1Sta[2]->fZ);
+ }
+
+ // add the target constraint
+
+ fFit.AddTargetToLine(fAlgo->fTargX, fAlgo->fTargY, fAlgo->fTargZ, fAlgo->TargetXYInfo, fld0);
+
+ fFit.AddMsInMaterial(fAlgo->fParameters.GetMaterialThickness(fIstaL, T.x, T.y));
+
+ // extrapolate to the middle hit
+
+ fFit.ExtrapolateLine(staM().fZ, fFldL);
+
+ fTrL = T;
+ }
+
+ FindDoublets();
+ FitDoublets();
+ FindTriplets();
+}
+
+
+void L1TripletConstructor::FindDoublets()
+{
+ /// Find the doublets. Reformat data into portions of doublets.
+
+ int iMC = -1;
+ if (fAlgo->fParameters.DevIsMatchDoubletsViaMc()) {
+ int indL = fAlgo->fGridHitStartIndex[fIstaL] + fIhitL;
+ iMC = fAlgo->GetMcTrackIdForGridHit(indL);
+ if (iMC < 0) { return; }
+ }
+
+ CollectHits(fTrL, fIstaM, fAlgo->fDoubletChi2Cut, iMC, fHitsM_2, fAlgo->fParameters.GetMaxDoubletsPerSinglet());
+}
+
+
+void L1TripletConstructor::FitDoublets()
+{
+
+ // ---- Add the middle hits to parameters estimation ----
+
+ L1Vector<L1HitIndex_t> hitsMtmp("L1TripletConstructor::hitsMtmp", fHitsM_2);
+
+ fHitsM_2.clear();
+ fTracks_2.clear();
+ fTracks_2.reserve(hitsMtmp.size());
+
+ bool isMomentumFitted = (fAlgo->fIsTargetField || (staL().fieldStatus != 0) || (staM().fieldStatus != 0));
+
+ for (unsigned int i2 = 0; i2 < hitsMtmp.size(); i2++) {
+
+ const L1HitIndex_t imh = hitsMtmp[i2];
+ const L1HitPoint& hitm = fHitsM[imh];
+
+ if (hitm.IsSuppressed()) continue;
+
+ L1TrackPar& T2 = fFit.Tr();
+ T2 = fTrL;
+ fFit.SetQp0(fvec(0.f));
+
+ fvec x_2 = hitm.X();
+ fvec y_2 = hitm.Y();
+ fvec z_2 = hitm.Z();
+ fvec t_2 = hitm.T();
+ L1XYMeasurementInfo cov_2;
+ cov_2.C00 = hitm.dX2();
+ cov_2.C10 = hitm.dXY();
+ cov_2.C11 = hitm.dY2();
+ fvec dt2_2 = hitm.dT2();
+
+ // add the middle hit
+
+ fFit.ExtrapolateLineNoField(z_2);
+
+ fFit.FilterXY(cov_2, x_2, y_2);
+
+ fFit.FilterTime(t_2, dt2_2, staM().timeInfo);
+
+ fFit.SetQp0(isMomentumFitted ? fFit.Tr().qp : fAlgo->fMaxInvMom);
+
+ fFit.AddMsInMaterial(fAlgo->fParameters.GetMaterialThickness(fIstaM, T2.x, T2.y));
+
+ fFit.SetQp0(fFit.Tr().qp);
+
+ // check if there are other hits close to the doublet on the same station
+ {
+ fscal tx = T2.tx[0];
+ fscal ty = T2.ty[0];
+ fscal tt = T2.vi[0] * sqrt(1. + tx * tx + ty * ty); // dt/dl * dl/dz
+
+ for (unsigned int j2 = i2 + 1; j2 < hitsMtmp.size(); j2++) {
+
+ L1HitPoint& hitm1 = fHitsM[hitsMtmp[j2]];
+
+ fscal dz = hitm1.Z() - T2.z[0];
+
+ if ((staM().timeInfo) && (T2.nTimeMeasurements[0] > 0)) {
+ fscal dt = T2.t[0] + tt * dz - hitm.T();
+ if (dt * dt > 30. * (T2.C55[0] + hitm1.dT2())) { continue; }
+ }
+
+ fscal dx = T2.x[0] + tx * dz - hitm1.X();
+ if (dx * dx > 20. * (T2.C00[0] + hitm1.dX2())) { continue; }
+
+ fscal dy = T2.y[0] + ty * dz - hitm1.Y();
+ if (dy * dy > 30. * (T2.C11[0] + hitm1.dY2())) { continue; }
+
+ if (fAlgo->fParameters.DevIsSuppressOverlapHitsViaMc()) {
+ int indL = fAlgo->fGridHitStartIndex[fIstaL] + fIhitL;
+ int indM = fAlgo->fGridHitStartIndex[fIstaM] + imh;
+ int indM1 = fAlgo->fGridHitStartIndex[fIstaM] + hitsMtmp[j2];
+ int iMC = fAlgo->GetMcTrackIdForGridHit(indL);
+ if ((iMC != fAlgo->GetMcTrackIdForGridHit(indM)) || (iMC != fAlgo->GetMcTrackIdForGridHit(indM1))) {
+ continue;
+ }
+ }
+ hitm1.SetIsSuppresed(1);
+ }
+ }
+
+ fHitsM_2.push_back(imh);
+ fTracks_2.push_back(T2);
+ } // i2
+}
+
+
+void L1TripletConstructor::FindTriplets()
+{
+ if (fIstaR >= fAlgo->fParameters.GetNstationsActive()) { return; }
+ findRightHit();
+ fitTriplets();
+ storeTriplets();
+}
+
+
+void L1TripletConstructor::findRightHit()
+{
+
+ /// Add the middle hits to parameters estimation. Propagate to right station.
+ /// Find the triplets(right hit). Reformat data in the portion of triplets.
+
+ fFit.SetQp0(fvec(0.));
+
+ if (fIstaR >= fAlgo->fParameters.GetNstationsActive()) return;
+
+ {
+ int maxTriplets = fHitsM_2.size() * fAlgo->fParameters.GetMaxTripletPerDoublets();
+ fHitsM_3.clear();
+ fHitsR_3.clear();
+ fHitsM_3.reserve(maxTriplets);
+ fHitsR_3.reserve(maxTriplets);
+ }
+ // ---- Add the middle hits to parameters estimation. Propagate to right station. ----
+
+ for (unsigned int i2 = 0; i2 < fHitsM_2.size(); i2++) {
+
+ fFit.SetTrack(fTracks_2[i2]);
+ L1TrackPar& T2 = fFit.Tr();
+
+ // extrapolate to the right hit station
+
+ fFit.Extrapolate(staR().fZ, fFldL);
+
+ // ---- Find the triplets(right hit). Reformat data in the portion of triplets. ----
+
+ if (fAlgo->kSts == fAlgo->fTrackingMode && (T2.C44[0] < 0)) { continue; }
+ if (T2.C00[0] < 0 || T2.C11[0] < 0 || T2.C22[0] < 0 || T2.C33[0] < 0 || T2.C55[0] < 0) continue;
+
+ if (fabs(T2.tx[0]) > fAlgo->fMaxSlope) continue;
+ if (fabs(T2.ty[0]) > fAlgo->fMaxSlope) continue;
+
+ int iMC = -1;
+ if (fAlgo->fParameters.DevIsMatchTripletsViaMc()) {
+ int indL = fAlgo->fGridHitStartIndex[fIstaL] + fIhitL;
+ int indM = fAlgo->fGridHitStartIndex[fIstaM] + fHitsM_2[i2];
+ iMC = fAlgo->GetMcTrackIdForGridHit(indL);
+ if (iMC < 0 || iMC != fAlgo->GetMcTrackIdForGridHit(indM)) { continue; }
+ }
+
+ L1Vector<L1HitIndex_t> collectedHits;
+ CollectHits(T2, fIstaR, fAlgo->fTripletChi2Cut, iMC, collectedHits, fAlgo->fParameters.GetMaxTripletPerDoublets());
+
+ if (collectedHits.size() >= fAlgo->fParameters.GetMaxTripletPerDoublets()) {
+ LOG(debug) << "L1: GetMaxTripletPerDoublets==" << fAlgo->fParameters.GetMaxTripletPerDoublets()
+ << " reached in findTripletsStep0()";
+ // reject already created triplets for this doublet
+ collectedHits.clear();
+ }
+
+ for (unsigned int ih = 0; ih < collectedHits.size(); ih++) {
+ L1HitIndex_t irh = collectedHits[ih];
+ const L1HitPoint& hitr = fHitsR[irh];
+ if (hitr.IsSuppressed()) continue;
+
+ // pack the triplet
+
+ fHitsM_3.push_back(fHitsM_2[i2]);
+ fHitsR_3.push_back(irh);
+
+ } // search area
+ } // i2
+}
+
+
+void L1TripletConstructor::fitTriplets()
+{
+ constexpr bool dumpTriplets = 0;
+
+ constexpr int nIterations = 2;
+
+ Tindex n3 = fHitsM_3.size();
+ assert(fHitsM_3.size() == fHitsR_3.size());
+
+ fTracks_3.clear();
+ fTracks_3.reset(n3, L1TrackPar());
+
+ /// Refit Triplets
+ if (dumpTriplets) { ca::tools::Debugger::Instance().AddNtuple("triplets", "ev:mc:sta:p:vx:vy:vz:chi2:ndf"); }
+
+ L1Fit fit;
+ fit.SetMask(fmask::One());
+
+ if (fAlgo->fpCurrentIteration->GetElectronFlag()) { fit.SetParticleMass(L1Constants::phys::kElectronMass); }
+ else {
+ fit.SetParticleMass(fAlgo->GetDefaultParticleMass());
+ }
+ const int NHits = 3; // triplets
+
+ // prepare data
+ int ista[NHits] = {fIstaL, fIstaM, fIstaR};
+
+ L1Station sta[3];
+ for (int is = 0; is < NHits; ++is) {
+ sta[is] = fAlgo->fParameters.GetStation(ista[is]);
+ };
+
+ bool isMomentumFitted = ((staL().fieldStatus != 0) || (staM().fieldStatus != 0) || (staR().fieldStatus != 0));
+ bool isTimeFitted = ((staL().timeInfo != 0) || (staM().timeInfo != 0) || (staR().timeInfo != 0));
+ fvec ndf = -4; // straight line
+ ndf += isMomentumFitted ? -1 : 0;
+ ndf += isTimeFitted ? -1 : 0;
+
+ for (int i3 = 0; i3 < n3; ++i3) {
+
+ // prepare data
+ L1HitIndex_t ihit[NHits] = {fAlgo->fGridHitIds[fIhitL + fAlgo->fGridHitStartIndex[ista[0]]],
+ fAlgo->fGridHitIds[fHitsM_3[i3] + fAlgo->fGridHitStartIndex[ista[1]]],
+ fAlgo->fGridHitIds[fHitsR_3[i3] + fAlgo->fGridHitStartIndex[ista[2]]]};
+
+ if (fAlgo->fParameters.DevIsMatchTripletsViaMc()) {
+ int mc1 = fAlgo->GetMcTrackIdForCaHit(ihit[0]);
+ int mc2 = fAlgo->GetMcTrackIdForCaHit(ihit[1]);
+ int mc3 = fAlgo->GetMcTrackIdForCaHit(ihit[2]);
+ if ((mc1 != mc2) || (mc1 != mc3)) { continue; }
+ }
+
+ fscal x[NHits], y[NHits], z[NHits], t[NHits];
+ fscal dt2[NHits];
+ L1XYMeasurementInfo cov[NHits];
+
+ for (int ih = 0; ih < NHits; ++ih) {
+ const L1Hit& hit = fAlgo->fInputData.GetHit(ihit[ih]);
+ x[ih] = hit.x;
+ y[ih] = hit.y;
+ z[ih] = hit.z;
+ t[ih] = hit.t;
+ cov[ih].C00 = hit.dx * hit.dx;
+ cov[ih].C10 = hit.dxy;
+ cov[ih].C11 = hit.dy * hit.dy;
+ dt2[ih] = hit.dt2;
+ };
+
+ // find the field along the track
+
+ L1FieldValue B[3] _fvecalignment;
+ L1FieldRegion fld _fvecalignment;
+ L1FieldRegion fldTarget _fvecalignment;
+
+ fvec tx[3] = {(x[1] - x[0]) / (z[1] - z[0]), (x[2] - x[0]) / (z[2] - z[0]), (x[2] - x[1]) / (z[2] - z[1])};
+ fvec ty[3] = {(y[1] - y[0]) / (z[1] - z[0]), (y[2] - y[0]) / (z[2] - z[0]), (y[2] - y[1]) / (z[2] - z[1])};
+ for (int ih = 0; ih < NHits; ++ih) {
+ fvec dz = (sta[ih].fZ - z[ih]);
+ sta[ih].fieldSlice.GetFieldValue(x[ih] + tx[ih] * dz, y[ih] + ty[ih] * dz, B[ih]);
+ };
+
+ fld.Set(B[0], sta[0].fZ, B[1], sta[1].fZ, B[2], sta[2].fZ);
+ fldTarget.Set(fAlgo->fTargB, fAlgo->fTargZ, B[0], sta[0].fZ, B[1], sta[1].fZ);
+
+ L1TrackPar& T = fit.Tr();
+
+ // initial parameters
+ {
+ fvec dz01 = 1. / (z[1] - z[0]);
+ T.tx = (x[1] - x[0]) * dz01;
+ T.ty = (y[1] - y[0]) * dz01;
+ T.qp = 0.;
+ T.vi = 0.;
+ }
+
+ // repeat several times in order to improve the precision
+ for (int iiter = 0; iiter < nIterations; ++iiter) {
+ // fit forward
+ {
+ fit.SetQp0(T.qp);
+ fit.Qp0()(fit.Qp0() > fAlgo->GetMaxInvMom()) = fAlgo->GetMaxInvMom();
+ fit.Qp0()(fit.Qp0() < -fAlgo->GetMaxInvMom()) = -fAlgo->GetMaxInvMom();
+
+ T.qp = 0.;
+ T.vi = 0.;
+ int ih0 = 0;
+ T.x = x[ih0];
+ T.y = y[ih0];
+ T.z = z[ih0];
+ T.t = t[ih0];
+
+ T.ResetErrors(1., 1., 1., 1., 100., (sta[ih0].timeInfo ? dt2[ih0] : 1.e6), 1.e2);
+ T.C00 = cov[ih0].C00;
+ T.C10 = cov[ih0].C10;
+ T.C11 = cov[ih0].C11;
+
+ T.NDF = ndf;
+
+ // add the target constraint
+ fit.AddTargetToLine(fAlgo->fTargX, fAlgo->fTargY, fAlgo->fTargZ, fAlgo->TargetXYInfo, fldTarget);
+
+ for (int ih = 1; ih < NHits; ++ih) {
+ fit.Extrapolate(z[ih], fld);
+ auto radThick = fAlgo->fParameters.GetMaterialThickness(ista[ih], T.x, T.y);
+ fit.AddMsInMaterial(radThick);
+ fit.EnergyLossCorrection(radThick, fvec(-1.f));
+ fit.FilterXY(cov[ih], x[ih], y[ih]);
+ fit.FilterTime(t[ih], dt2[ih], sta[ih].timeInfo);
+ }
+ }
+
+ if (iiter == nIterations - 1) break;
+
+ // fit backward
+ {
+ fit.SetQp0(T.qp);
+ fit.Qp0()(fit.Qp0() > fAlgo->GetMaxInvMom()) = fAlgo->GetMaxInvMom();
+ fit.Qp0()(fit.Qp0() < -fAlgo->GetMaxInvMom()) = -fAlgo->GetMaxInvMom();
+
+ T.NDF = ndf;
+ T.qp = 0.;
+ T.vi = 0.;
+ int ih0 = NHits - 1;
+ T.x = x[ih0];
+ T.y = y[ih0];
+ T.z = z[ih0];
+ T.t = t[ih0];
+
+ T.ResetErrors(1., 1., 1., 1., 100., (sta[ih0].timeInfo ? dt2[ih0] : 1.e6), 1.e2);
+ T.C00 = cov[ih0].C00;
+ T.C10 = cov[ih0].C10;
+ T.C11 = cov[ih0].C11;
+
+ T.NDF = ndf;
+
+ for (int ih = NHits - 2; ih >= 0; --ih) {
+ fit.Extrapolate(z[ih], fld);
+ auto radThick = fAlgo->fParameters.GetMaterialThickness(ista[ih], T.x, T.y);
+ fit.AddMsInMaterial(radThick);
+ fit.EnergyLossCorrection(radThick, fvec(1.f));
+ fit.FilterXY(cov[ih], x[ih], y[ih]);
+ fit.FilterTime(t[ih], dt2[ih], sta[ih].timeInfo);
+ }
+ }
+ } // for iiter
+
+ fTracks_3[i3] = T;
+
+ if (dumpTriplets) {
+ int mc1 = fAlgo->GetMcTrackIdForCaHit(ihit[0]);
+ int mc2 = fAlgo->GetMcTrackIdForCaHit(ihit[1]);
+ int mc3 = fAlgo->GetMcTrackIdForCaHit(ihit[2]);
+ if ((mc1 >= 0) && (mc1 == mc2) && (mc1 == mc3)) {
+ const CbmL1MCTrack& mctr = CbmL1::Instance()->GetMcTracks()[mc1];
+ float ev = 0;
+ float chi2 = T.chi2[0];
+ float nd = T.NDF[0];
+ ca::tools::Debugger::Instance().FillNtuple("triplets", ev, mc1, fIstaL, mctr.p, mctr.x, mctr.y, mctr.z, chi2,
+ nd);
+ }
+ }
+ } //i3
+
+} // findTripletsStep2
+
+
+void L1TripletConstructor::storeTriplets()
+{
+ /// Selects good triplets and saves them into fTriplets.
+ /// Finds neighbouring triplets at the next station.
+
+ Tindex n3 = fHitsM_3.size();
+
+ bool isMomentumFitted = ((staL().fieldStatus != 0) || (staM().fieldStatus != 0) || (staR().fieldStatus != 0));
+
+ fTriplets.clear();
+ fTriplets.reserve(n3);
+
+ for (Tindex i3 = 0; i3 < n3; ++i3) {
+
+ L1TrackPar& T3 = fTracks_3[i3];
+
+ fscal chi2 = T3.chi2[0]; // / T3.NDF[0];
+
+ const L1HitIndex_t ihitl = fIhitL + fAlgo->fGridHitStartIndex[fIstaL];
+ const L1HitIndex_t ihitm = fHitsM_3[i3] + fAlgo->fGridHitStartIndex[fIstaM];
+ const L1HitIndex_t ihitr = fHitsR_3[i3] + fAlgo->fGridHitStartIndex[fIstaR];
+
+ L1_ASSERT(ihitl < fAlgo->fGridHitStopIndex[fIstaL], ihitl << " < " << fAlgo->fGridHitStopIndex[fIstaL]);
+ L1_ASSERT(ihitm < fAlgo->fGridHitStopIndex[fIstaM], ihitm << " < " << fAlgo->fGridHitStopIndex[fIstaM]);
+ L1_ASSERT(ihitr < fAlgo->fGridHitStopIndex[fIstaR], ihitr << " < " << fAlgo->fGridHitStopIndex[fIstaR]);
+
+ if (!fAlgo->fpCurrentIteration->GetTrackFromTripletsFlag()) {
+ if (chi2 > fAlgo->fTripletFinalChi2Cut) { continue; }
+ }
+
+ // assert(std::isfinite(chi2) && chi2 > 0);
+
+ if (!std::isfinite(chi2) || chi2 < 0) { continue; }
+ //if (!T3.IsEntryConsistent(false, 0)) { continue; }
+
+ fscal qp = T3.qp[0];
+
+ //TODO: why sqrt's? Wouldn't it be faster to skip sqrt() here and
+ //TODO: compare the squared differences dqr*dqp later?
+
+ fscal Cqp = T3.C44[0];
+
+ // TODO: SG: a magic correction that comes from the legacy code
+ // removing it leads to a higher ghost ratio
+ Cqp += 0.001;
+
+ fTriplets.emplace_back(ihitl, ihitm, ihitr, fIstaL, fIstaM, fIstaR, 0, 0, 0, chi2, qp, Cqp, T3.tx[0], T3.C22[0],
+ T3.ty[0], T3.C33[0], isMomentumFitted);
+ }
+}
+
+
+void L1TripletConstructor::CollectHits(const L1TrackPar& Tr, const int iSta, const double chi2Cut, const int iMC,
+ L1Vector<L1HitIndex_t>& collectedHits, int maxNhits)
+{
+ /// Collect hits on a station
+
+ collectedHits.clear();
+ collectedHits.reserve(maxNhits);
+
+ const L1Station& sta = fAlgo->fParameters.GetStation(iSta);
+ const L1HitPoint* hits = &(fAlgo->fGridPoints[0]) + fAlgo->fGridHitStartIndex[iSta];
+ const L1HitIndex_t nHits = fAlgo->fGridHitStopIndex[iSta] - fAlgo->fGridHitStartIndex[iSta];
+
+ fFit.SetTrack(Tr);
+ L1TrackPar& T = fFit.Tr();
+
+ // if make it bigger the found hits will be rejected later because of the chi2 cut.
+ const fvec Pick_m22 = (fvec(chi2Cut) - T.chi2);
+
+ const fscal iz = 1.f / (T.z[0] - fAlgo->fParameters.GetTargetPositionZ()[0]);
+ const fscal timeError2 = T.C55[0];
+ const fscal time = T.t[0];
+
+ L1HitAreaTime areaTime(fAlgo->vGridTime[iSta], T.x[0] * iz, T.y[0] * iz,
+ (sqrt(Pick_m22 * T.C00) + fAlgo->fMaxDx[iSta] + fAlgo->fMaxDZ * abs(T.tx))[0] * iz,
+ (sqrt(Pick_m22 * T.C11) + fAlgo->fMaxDy[iSta] + fAlgo->fMaxDZ * abs(T.ty))[0] * iz, time,
+ sqrt(timeError2) * 5);
+
+ for (L1HitIndex_t ih = 0; (ih < nHits) && ((int) collectedHits.size() < maxNhits);
+ ih++) { // loop over all station hits
+
+ if (!fAlgo->fParameters.DevIsIgnoreHitSearchAreas()) { // only loop over the hits in the area
+ if (!areaTime.GetNext(ih)) { break; }
+ }
+
+ const L1HitPoint& hit = hits[ih];
+ if (hit.IsSuppressed()) continue;
+
+ if (iMC >= 0) { // match via MC
+ if (iMC != fAlgo->GetMcTrackIdForGridHit(fAlgo->fGridHitStartIndex[iSta] + ih)) { continue; }
+ }
+
+ // check y-boundaries
+ //TODO: move hardcoded cuts to parameters
+ if ((sta.timeInfo) && (T.nTimeMeasurements[0] > 0)) {
+ if (fabs(time - hit.T()) > sqrt(timeError2 + hit.dT2()) * 5) continue;
+ if (fabs(time - hit.T()) > 30) continue;
+ }
+
+ // - check whether hit belong to the window ( track position +\- errors ) -
+ const fscal z = hit.Z();
+
+ fvec y, C11;
+ fFit.ExtrapolateYC11Line(z, y, C11);
+
+ fscal dy_est2 = Pick_m22[0] * fabs(C11[0] + hit.dY2());
+
+ fscal xm = hit.X();
+ fscal ym = hit.Y();
+
+ const fscal dY = ym - y[0];
+
+ if (dY * dY > dy_est2) continue;
+
+ // check x-boundaries
+ fvec x, C00;
+ fFit.ExtrapolateXC00Line(z, x, C00);
+
+ fscal dx_est2 = Pick_m22[0] * fabs(C00[0] + hit.dX2());
+
+ const fscal dX = xm - x[0];
+
+ if (dX * dX > dx_est2) continue;
+
+ // check chi2
+
+ fvec C10;
+ fFit.ExtrapolateC10Line(z, C10);
+
+ auto [chi2x, chi2u] = L1Fit::GetChi2XChi2U(hit.X(), hit.Y(), hit.dX2(), hit.dXY(), hit.dY2(), x, y, C00, C10, C11);
+
+ // TODO: adjust the cut, cut on chi2x & chi2u separately
+ if (!fAlgo->fpCurrentIteration->GetTrackFromTripletsFlag()) {
+ if (chi2x[0] > chi2Cut) continue;
+ if (chi2x[0] + chi2u[0] > chi2Cut) continue;
+ }
+
+ collectedHits.push_back(ih);
+
+ } // loop over the hits in the area
+}
diff --git a/reco/L1/L1Algo/L1TripletConstructor.h b/reco/L1/L1Algo/L1TripletConstructor.h
new file mode 100644
index 0000000000..6cf5a6b31d
--- /dev/null
+++ b/reco/L1/L1Algo/L1TripletConstructor.h
@@ -0,0 +1,116 @@
+/* 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 "L1Algo.h"
+#include "L1Field.h"
+#include "L1Fit.h"
+#include "L1HitPoint.h"
+#include "L1Station.h"
+#include "L1TrackPar.h"
+#include "L1Triplet.h"
+#include "L1Vector.h"
+
+
+/// 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, L1HitIndex_t ihl);
+
+ const L1Vector<L1Triplet>& 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 L1TrackPar& Tr, const int iSta, const double chi2Cut, const int iMC,
+ L1Vector<L1HitIndex_t>& collectedHits, int maxNhits);
+
+private:
+ /// left station
+ const L1Station& staL() { return *fStaL; }
+ /// middle station
+ const L1Station& staM() { return *fStaM; }
+ /// right station
+ const L1Station& 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 L1Station* fStaL {nullptr}; ///< left station
+ const L1Station* fStaM {nullptr}; ///< mid station
+ const L1Station* fStaR {nullptr}; ///< right station
+
+ L1HitPoint* fHitsL {nullptr}; ///< hits on the left station
+ L1HitPoint* fHitsM {nullptr}; ///< hits on the middle station
+ L1HitPoint* fHitsR {nullptr}; ///< hits on the right station
+
+ const L1Station* fFld0Sta[2]; // two stations for approximating the field between the target and the left hit
+ const L1Station* fFld1Sta[3]; // three stations for approximating the field between the left and the right hit
+
+ L1HitIndex_t fIhitL;
+ L1TrackPar fTrL;
+ L1FieldRegion fFldL;
+
+ L1Vector<L1HitIndex_t> fHitsM_2 {"L1TripletConstructor::fHitsM_2"};
+ L1Vector<L1TrackPar> fTracks_2 {"L1TripletConstructor::fTracks_2"};
+
+ L1Vector<L1TrackPar> fTracks_3 {"L1TripletConstructor::fTracks_3"};
+ L1Vector<L1HitIndex_t> fHitsM_3 {"L1TripletConstructor::fHitsM_3"};
+ L1Vector<L1HitIndex_t> fHitsR_3 {"L1TripletConstructor::fHitsR_3"};
+
+ L1Vector<L1Triplet> fTriplets {"L1TripletConstructor::fTriplets"};
+
+ L1Fit fFit;
+
+} _fvecalignment;
+
+#endif
--
GitLab