From 9cf4a05a08e48b8020e91619082995fe532cc81d Mon Sep 17 00:00:00 2001
From: Dominik Smith <d.smith@gsi.de>
Date: Thu, 11 Jul 2024 16:38:21 +0200
Subject: [PATCH] Cleaned up the core tracking classes. Summary of changes:
- Removed access to ca::Framework from ca::TrackExtender, ca::TrackFitter
ca::TripletConstructor, ca::CloneMerger, ca::TrackFinder and ca::TrackFinderWindow.
- Simplified the interfaces of these classes.
- General instruction-wise cleanup of these classes (applied std functions in
some places etc.). In particular, streamlining of ca::TrackFinder.
- Made debugging functions in ca::Framework static.
- Moved TrackingMode enum from CaFramework.h to CaParameters.h (needed
to avoid cyclic header file dependency).
- Moved some magnetic-field dependent variables to ca::TripletConstructor (they are
only used there).
- Restricted the scope of many variables, i.e. turned class-fields into
local variables in many instances (e.g. fvMonitorDataThread, fvTrackFinderWindow,
fTrackingMode, fvWData, fvWindowStartThread, fvWindowEndThread,
fvStatNwindows, fvStatNhitsProcessed, fvHitFirstTriplet and fvHitNofTriplets).
- Eliminated the Init() functions of some classes. Moved initialization
to the constructors.
- Improved distinction between class-parameters and working data.
Class parameters are passed to constructor, working data is passed as
arguments to function calls.
---
algo/ca/TrackingChain.cxx | 4 +-
algo/ca/core/pars/CaParameters.h | 11 +-
algo/ca/core/tracking/CaCloneMerger.cxx | 25 +-
algo/ca/core/tracking/CaCloneMerger.h | 18 +-
algo/ca/core/tracking/CaFramework.cxx | 47 +-
algo/ca/core/tracking/CaFramework.h | 121 +----
algo/ca/core/tracking/CaTrackExtender.cxx | 74 +--
algo/ca/core/tracking/CaTrackExtender.h | 13 +-
algo/ca/core/tracking/CaTrackFinder.cxx | 513 ++++++++----------
algo/ca/core/tracking/CaTrackFinder.h | 41 +-
algo/ca/core/tracking/CaTrackFinderWindow.cxx | 286 +++++-----
algo/ca/core/tracking/CaTrackFinderWindow.h | 84 ++-
algo/ca/core/tracking/CaTrackFitter.cxx | 54 +-
algo/ca/core/tracking/CaTrackFitter.h | 17 +-
.../ca/core/tracking/CaTripletConstructor.cxx | 156 +++---
algo/ca/core/tracking/CaTripletConstructor.h | 16 +-
reco/L1/CbmL1.cxx | 20 +-
reco/L1/CbmL1.h | 14 +-
reco/L1/CbmL1Performance.cxx | 4 +-
reco/L1/L1Algo/utils/L1AlgoDraw.cxx | 38 +-
20 files changed, 763 insertions(+), 793 deletions(-)
diff --git a/algo/ca/TrackingChain.cxx b/algo/ca/TrackingChain.cxx
index 39178f8573..6481750481 100644
--- a/algo/ca/TrackingChain.cxx
+++ b/algo/ca/TrackingChain.cxx
@@ -86,7 +86,7 @@ void TrackingChain::Init()
fCaFramework.SetNofThreads(Opts().NumOMPThreads() == std::nullopt ? openmp::GetMaxThreads()
: *(Opts().NumOMPThreads()));
fCaFramework.ReceiveParameters(std::move(parameters));
- fCaFramework.Init(ca::Framework::TrackingMode::kMcbm);
+ fCaFramework.Init(ca::TrackingMode::kMcbm);
// ------ Initialize QA modules
if (fQa.IsSenderDefined()) {
@@ -110,7 +110,7 @@ TrackingChain::Output_t TrackingChain::Run(Input_t recoResults)
// ----- Run reconstruction ------------------------------------------------------------------------------------------
fCaFramework.SetMonitorData(fCaMonitorData);
- fCaFramework.fTrackFinder.FindTracks();
+ fCaFramework.FindTracks();
fCaMonitorData = fCaFramework.GetMonitorData();
// ----- Init output data --------------------------------------------------------------------------------------------
diff --git a/algo/ca/core/pars/CaParameters.h b/algo/ca/core/pars/CaParameters.h
index a7809c1155..30a0496d7c 100644
--- a/algo/ca/core/pars/CaParameters.h
+++ b/algo/ca/core/pars/CaParameters.h
@@ -35,8 +35,15 @@ namespace cbm::algo::ca
/// Type definitions for used containers
using IterationsContainer_t = cbm::algo::ca::Vector<cbm::algo::ca::Iteration>;
template<typename DataT>
- using StationsContainer_t = std::array<ca::Station<DataT>, constants::size::MaxNstations>;
- using MaterialContainer_t = std::array<ca::MaterialMap, constants::size::MaxNstations>;
+ using StationsContainer_t = std::array<ca::Station<DataT>, constants::size::MaxNstations>;
+ using MaterialContainer_t = std::array<ca::MaterialMap, constants::size::MaxNstations>;
+
+ enum TrackingMode
+ {
+ kSts,
+ kGlobal,
+ kMcbm
+ };
/// \class cbm::algo::ca::Parameters
/// \brief A container for all external parameters of the CA tracking algorithm
diff --git a/algo/ca/core/tracking/CaCloneMerger.cxx b/algo/ca/core/tracking/CaCloneMerger.cxx
index 53c1eba4bf..5755539265 100644
--- a/algo/ca/core/tracking/CaCloneMerger.cxx
+++ b/algo/ca/core/tracking/CaCloneMerger.cxx
@@ -22,7 +22,7 @@ using namespace cbm::algo;
// ---------------------------------------------------------------------------------------------------------------------
//
-CloneMerger::CloneMerger(const ca::Framework& algo) : frAlgo(algo) {}
+CloneMerger::CloneMerger(const ca::Parameters<fvec>& pars, const float mass) : fParameters(pars), fDefaultMass(mass) {}
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -30,8 +30,11 @@ CloneMerger::~CloneMerger() {}
// ---------------------------------------------------------------------------------------------------------------------
//
-void CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extRecoHits)
+void CloneMerger::Exec(const ca::InputData& input, WindowData& wData)
{
+ Vector<Track>& extTracks = wData.RecoTracks();
+ Vector<ca::HitIndex_t>& extRecoHits = wData.RecoHitIndices();
+
Vector<unsigned short>& firstStation = fTrackFirstStation;
Vector<unsigned short>& lastStation = fTrackLastStation;
Vector<ca::HitIndex_t>& firstHit = fTrackFirstHit;
@@ -65,10 +68,10 @@ void CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extReco
for (int iTr = 0; iTr < nTracks; iTr++) {
firstHit[iTr] = start_hit;
- firstStation[iTr] = frAlgo.GetInputData().GetHit(extRecoHits[start_hit]).Station();
+ firstStation[iTr] = input.GetHit(extRecoHits[start_hit]).Station();
start_hit += extTracks[iTr].fNofHits - 1;
lastHit[iTr] = start_hit;
- lastStation[iTr] = frAlgo.GetInputData().GetHit(extRecoHits[start_hit]).Station();
+ lastStation[iTr] = input.GetHit(extRecoHits[start_hit]).Station();
start_hit++;
isStored[iTr] = false;
@@ -78,12 +81,12 @@ void CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extReco
}
ca::TrackFit fitB;
- fitB.SetParticleMass(frAlgo.GetDefaultParticleMass());
+ fitB.SetParticleMass(fDefaultMass);
fitB.SetMask(fmask::One());
fitB.SetQp0(fvec(0.));
ca::TrackFit fitF;
- fitF.SetParticleMass(frAlgo.GetDefaultParticleMass());
+ fitF.SetParticleMass(fDefaultMass);
fitF.SetMask(fmask::One());
fitF.SetQp0(fvec(0.));
@@ -93,7 +96,7 @@ void CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extReco
ca::FieldRegion<fvec> fld _fvecalignment;
// Max length for merging
- unsigned char maxLengthForMerge = static_cast<unsigned char>(frAlgo.GetParameters().GetNstationsActive() - 3);
+ unsigned char maxLengthForMerge = static_cast<unsigned char>(fParameters.GetNstationsActive() - 3);
for (int iTr = 0; iTr < nTracks; iTr++) {
if (extTracks[iTr].fNofHits > maxLengthForMerge) continue;
@@ -119,8 +122,8 @@ void CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extReco
unsigned short stam;
- fBf = frAlgo.GetParameters().GetStation(staf).fieldSlice.GetFieldValue(Tf.X(), Tf.Y());
- fBb = frAlgo.GetParameters().GetStation(stab).fieldSlice.GetFieldValue(Tb.X(), Tb.Y());
+ fBf = fParameters.GetStation(staf).fieldSlice.GetFieldValue(Tf.X(), Tf.Y());
+ fBb = fParameters.GetStation(stab).fieldSlice.GetFieldValue(Tb.X(), Tb.Y());
unsigned short dist = firstStation[iTr] - lastStation[jTr];
@@ -129,10 +132,10 @@ void CloneMerger::Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>& extReco
else
stam = staf - 1;
- fvec zm = frAlgo.GetParameters().GetStation(stam).fZ;
+ fvec zm = fParameters.GetStation(stam).fZ;
fvec xm = fvec(0.5) * (Tf.GetX() + Tf.Tx() * (zm - Tf.Z()) + Tb.GetX() + Tb.Tx() * (zm - Tb.Z()));
fvec ym = fvec(0.5) * (Tf.Y() + Tf.Ty() * (zm - Tf.Z()) + Tb.Y() + Tb.Ty() * (zm - Tb.Z()));
- fBm = frAlgo.GetParameters().GetStation(stam).fieldSlice.GetFieldValue(xm, ym);
+ fBm = fParameters.GetStation(stam).fieldSlice.GetFieldValue(xm, ym);
fld.Set(fBb, Tb.Z(), fBm, zm, fBf, Tf.Z());
fvec zMiddle = fvec(0.5) * (Tb.Z() + Tf.Z());
diff --git a/algo/ca/core/tracking/CaCloneMerger.h b/algo/ca/core/tracking/CaCloneMerger.h
index 0cef1c285d..d063c32f08 100644
--- a/algo/ca/core/tracking/CaCloneMerger.h
+++ b/algo/ca/core/tracking/CaCloneMerger.h
@@ -9,11 +9,13 @@
#pragma once // include this header only once per compilation unit
-#include "CaHit.h" // For ca::HitIndex_t
+#include "CaHit.h" // For ca::HitIndex_t
+#include "CaInputData.h"
+#include "CaParameters.h"
#include "CaSimd.h" // TEMPORARY FOR fvec, fscal
#include "CaTrack.h"
#include "CaVector.h"
-
+#include "CaWindowData.h"
namespace cbm::algo::ca
{
@@ -31,7 +33,7 @@ namespace cbm::algo::ca
class CloneMerger {
public:
/// Default constructor
- CloneMerger(const ca::Framework& algo);
+ CloneMerger(const ca::Parameters<fvec>& pars, const float mass);
/// Destructor
~CloneMerger();
@@ -51,9 +53,10 @@ namespace cbm::algo::ca
/// Registers
/// Executes track clones merging algorithm and updates input containers
- /// \param extTracks Reference to the external container of reconstructed tracks
- /// \param extRecoHits Reference to the external container of reconstructed hit indexes
- void Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>&);
+ /// \param input Reference to input data
+ /// \param wData Reference to the external container of reconstructed tracks and hit indices
+ //void Exec(Vector<Track>& extTracks, Vector<ca::HitIndex_t>&, const ca::InputData& input);
+ void Exec(const ca::InputData& input, WindowData& wData);
private:
// ***************
@@ -122,7 +125,8 @@ namespace cbm::algo::ca
Vector<ca::HitIndex_t> fRecoHitsNew{"CaCloneMerger::fRecoHitsNew"}; ///< vector of track hits after the merge
- const ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class
+ const Parameters<fvec>& fParameters; ///< Object of Framework parameters class
+ float fDefaultMass{constants::phys::MuonMass}; ///< mass of the propagated particle [GeV/c2]
};
} // namespace cbm::algo::ca
diff --git a/algo/ca/core/tracking/CaFramework.cxx b/algo/ca/core/tracking/CaFramework.cxx
index 43a7085dbe..ed3d8d2007 100644
--- a/algo/ca/core/tracking/CaFramework.cxx
+++ b/algo/ca/core/tracking/CaFramework.cxx
@@ -1,12 +1,19 @@
/* 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 */
+ Authors: Igor Kulakov [committer], Maksym Zyzak, Valentina Akishina, Ivan Kisel, Sergey Gorbunov, Sergei Zharko, Dominik Smith */
#include "CaFramework.h"
#include "CaGridEntry.h"
+#include "CaTimer.h"
+#include "CaTrack.h"
// #include "CaToolsDebugger.h"
+#include <chrono>
+#include <fstream>
+#include <sstream>
+#include <thread>
+
using namespace cbm::algo::ca;
namespace
@@ -14,31 +21,24 @@ namespace
using namespace cbm::algo; // to keep ca:: prefices in the code
}
-Framework::Framework()
-{
-
-}
-
using cbm::algo::ca::ECounter; // monitor counter key type
using cbm::algo::ca::EDetectorID;
using cbm::algo::ca::ETimer; // monitor timer key type
using cbm::algo::ca::KfFramework_t;
using cbm::algo::ca::KfParameters_t;
+using cbm::algo::ca::Track;
+using constants::phys::ProtonMassD;
+using constants::phys::SpeedOfLightInv;
+using constants::phys::SpeedOfLightInvD;
//using cbm::ca::tools::Debugger;
// ---------------------------------------------------------------------------------------------------------------------
//
void Framework::Init(const TrackingMode mode)
{
- fTrackingMode = mode;
- for (int iThread = 0; iThread < fNofThreads; ++iThread) {
- fvWData.emplace_back();
- fvMonitorDataThread.emplace_back();
- }
- for (int iThread = 0; iThread < fNofThreads; ++iThread) {
- fvTrackFinderWindow.emplace_back(*this, fvWData[iThread], fvMonitorDataThread[iThread]);
- }
- fTrackFinder.Init();
+ fpTrackFinder =
+ std::make_unique<ca::TrackFinder>(fParameters, fDefaultMass, mode, fMonitorData, fNofThreads, fCaRecoTime);
+
fKf.Init();
}
@@ -62,20 +62,9 @@ void Framework::ReceiveInputData(InputData&& inputData)
//
void Framework::ReceiveParameters(Parameters<fvec>&& parameters)
{
- fParameters = std::move(parameters);
-
+ fParameters = std::move(parameters);
fNstationsBeforePipe = fParameters.GetNstationsActive(static_cast<EDetectorID>(0));
- // FIXME: SZh 24.08.2022
- // This approach is suitable only for a case, when all the stations inside a magnetic field come right before
- // all the stations outside of the field!
- fNfieldStations = std::lower_bound(fParameters.GetStations().cbegin(),
- fParameters.GetStations().cbegin() + fParameters.GetNstationsActive(),
- 0, // we are looking for the first zero element
- [](const ca::Station<fvec>& s, int edge) { return bool(s.fieldStatus) > edge; })
- - fParameters.GetStations().cbegin();
-
-
ca::FieldRegion<fvec>::ForceUseOfOriginalField(fParameters.DevIsUseOfOriginalField());
}
@@ -90,7 +79,7 @@ void Framework::ReceiveKfSettings(KfParameters_t&& pars, KfSetup_t&& setup)
// ---------------------------------------------------------------------------------------------------------------------
//
-int Framework::GetMcTrackIdForCaHit(int /*iHit*/) const
+int Framework::GetMcTrackIdForCaHit(int /*iHit*/)
{
return -1;
/*
@@ -101,7 +90,7 @@ int Framework::GetMcTrackIdForCaHit(int /*iHit*/) const
*/
}
-int Framework::GetMcTrackIdForWindowHit(int /*iHit*/) const
+int Framework::GetMcTrackIdForWindowHit(int /*iHit*/)
{
return -1;
/*
diff --git a/algo/ca/core/tracking/CaFramework.h b/algo/ca/core/tracking/CaFramework.h
index 409b165cb9..4b20fae8e2 100644
--- a/algo/ca/core/tracking/CaFramework.h
+++ b/algo/ca/core/tracking/CaFramework.h
@@ -14,14 +14,12 @@
#include "CaInputData.h"
#include "CaMeasurementXy.h"
#include "CaParameters.h"
+#include "CaSimd.h"
#include "CaStation.h"
#include "CaTimer.h"
#include "CaTimesliceHeader.h"
#include "CaTrack.h"
-#include "CaTrackExtender.h"
#include "CaTrackFinder.h"
-#include "CaTrackFinderWindow.h"
-#include "CaTrackFitter.h"
#include "CaTrackParam.h"
#include "CaTrackingMonitor.h"
#include "CaTriplet.h"
@@ -38,7 +36,7 @@
namespace cbm::algo::ca
{
class TripletConstructor;
-
+ class Track;
namespace
{
@@ -60,31 +58,11 @@ namespace cbm::algo::ca
using CaMaterialArray_t = std::array<ca::MaterialMap, constants::size::MaxNstations>;
using Tindex = int; // TODO: Replace with ca::HitIndex_t, if suitable
-
/// Main class of CA track finder algorithm
+ /// Class implements a clones merger algorithm for the CA track finder
///
class Framework {
public:
- // *************************
- // ** Friend classes list **
- // *************************
-
- friend class ca::TripletConstructor;
- friend class ca::TrackExtender;
- friend class ca::TrackFitter;
- friend class ca::TrackFinderWindow;
- friend class ca::TrackFinder;
-
- ///---------------------------
- /// Internal enumerations
-
- enum TrackingMode
- {
- kSts,
- kGlobal,
- kMcbm
- };
-
// **********************************
// ** Member functions declaration **
// **********************************
@@ -92,7 +70,7 @@ namespace cbm::algo::ca
// ** Constructors and destructor
/// Constructor
- Framework();
+ Framework(){};
/// Copy constructor
Framework(const Framework&) = delete;
@@ -109,24 +87,6 @@ namespace cbm::algo::ca
/// Destructor
~Framework() = default;
-
- // ** Functions, which pack and unpack indexes of station and triplet **
-
- // TODO: move to ca::Triplet class (S.Zharko)
- /// Packs station and triplet indices to an unique triplet ID
- /// \param iStation Index of station in the active stations array
- /// \param iTriplet Index of triplet
- static unsigned int PackTripletId(unsigned int iStation, unsigned int iTriplet);
-
- /// Unpacks the triplet ID to its station index
- /// \param id Unique triplet ID
- static unsigned int TripletId2Station(unsigned int id);
-
- /// Unpacks the triplet ID to its triplet index
- /// \param id Unique triplet ID
- static unsigned int TripletId2Triplet(unsigned int id);
-
-
/// Sets Framework parameters object
/// \param other - reference to the Parameters object
void SetParameters(const Parameters<fvec>& other) { fParameters = other; }
@@ -187,18 +147,25 @@ namespace cbm::algo::ca
/// \brief Sets monitor data
void SetMonitorData(const TrackingMonitorData& monitorData) { fMonitorData = monitorData; }
+ TrackingMode GetTrackingMode() { return fpTrackFinder->GetTrackingMode(); }
+
+ const std::vector<ca::WindowData>& GetWData() { return fpTrackFinder->GetWData(); }
+
public:
+ void FindTracks()
+ {
+ std::tie(fRecoTracks, fRecoHits) = fpTrackFinder->FindTracks(fInputData, fTsHeader);
+ return;
+ };
+
/// Gets number of stations before the pipe (MVD stations in CBM)
// int GetNstationsBeforePipe() const { return fNstationsBeforePipe; }
- /// Gets number of stations situated in field region (MVD + STS in CBM)
- int GetNfieldStations() const { return fNfieldStations; }
-
/// Get mc track ID for a hit (debug tool)
- int GetMcTrackIdForCaHit(int iHit) const;
+ static int GetMcTrackIdForCaHit(int iHit);
/// Get mc track ID for a hit (debug tool)
- int GetMcTrackIdForWindowHit(int iGridHit) const;
+ static int GetMcTrackIdForWindowHit(int iGridHit);
// const CbmL1MCTrack* GetMcTrackForWindowHit(int iHit) const;
@@ -223,7 +190,6 @@ namespace cbm::algo::ca
KfFramework_t fKf; ///< KF framework instance
int fNstationsBeforePipe{0}; ///< number of stations before pipe (MVD stations in CBM)
- int fNfieldStations{0}; ///< number of stations in the field region
float fDefaultMass{constants::phys::MuonMass}; ///< mass of the propagated particle [GeV/c2]
TimesliceHeader fTsHeader; ///< current timeslice header
@@ -237,13 +203,16 @@ namespace cbm::algo::ca
Vector<unsigned char> fvHitKeyFlags{
"Framework::fvHitKeyFlags"}; ///< List of key flags: has been this hit or cluster already used
- TrackingMonitorData fMonitorData{}; ///< Tracking monitor data (statistics per call)
- std::vector<TrackingMonitorData> fvMonitorDataThread; ///< Tracking monitor data per thread
+ TrackingMonitorData fMonitorData{}; ///< Tracking monitor data (statistics per call)
int fNofThreads = 1; ///< Number of threads to execute the track-finder
+ ///
+ bool checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2) const;
+
+ std::unique_ptr<ca::TrackFinder> fpTrackFinder; ///< Track finder steer class for the entire time slice
public:
- Vector<CaHitTimeInfo> fHitTimeInfo;
+ // Vector<CaHitTimeInfo> fHitTimeInfo;
double fCaRecoTime{0.}; // time of the track finder + fitter
@@ -254,52 +223,6 @@ namespace cbm::algo::ca
fvec Err{0.f};
public:
- std::vector<ca::WindowData> fvWData; ///< Intrnal data processed in a time-window
- std::vector<ca::TrackFinderWindow> fvTrackFinderWindow; ///< Track finder algorithm for the time window
- ca::TrackFinder fTrackFinder{*this}; ///< Track finder steer class for the entire time slice
- TrackingMode fTrackingMode{kSts};
-
- private:
- /// ================================= DATA PART =================================
-
- /// ----- Different parameters of CATrackFinder -----
-
- bool fIsTargetField{false}; ///< is the magnetic field present at the target
-
-
} _fvecalignment;
-
- // ********************************************
- // ** Inline member functions implementation **
- // ********************************************
-
- // ---------------------------------------------------------------------------------------------------------------------
- //
- [[gnu::always_inline]] inline unsigned int Framework::PackTripletId(unsigned int iStation, unsigned int iTriplet)
- {
-#ifndef FAST_CODE
- assert(iStation < constants::size::MaxNstations);
- assert(iTriplet < constants::size::MaxNtriplets);
-#endif
- constexpr unsigned int kMoveStation = constants::size::TripletBits;
- return (iStation << kMoveStation) + iTriplet;
- }
-
- // ---------------------------------------------------------------------------------------------------------------------
- //
- [[gnu::always_inline]] inline unsigned int Framework::TripletId2Station(unsigned int id)
- {
- constexpr unsigned int kMoveStation = constants::size::TripletBits;
- return id >> kMoveStation;
- }
-
- // ---------------------------------------------------------------------------------------------------------------------
- //
- [[gnu::always_inline]] inline unsigned int Framework::TripletId2Triplet(unsigned int id)
- {
- constexpr unsigned int kTripletMask = (1u << constants::size::TripletBits) - 1u;
- return id & kTripletMask;
- }
-
} // namespace cbm::algo::ca
diff --git a/algo/ca/core/tracking/CaTrackExtender.cxx b/algo/ca/core/tracking/CaTrackExtender.cxx
index 6ece464b2e..a7070603b4 100644
--- a/algo/ca/core/tracking/CaTrackExtender.cxx
+++ b/algo/ca/core/tracking/CaTrackExtender.cxx
@@ -9,6 +9,7 @@
#include "CaDefines.h"
#include "CaFramework.h"
#include "CaGridArea.h"
+#include "CaInputData.h"
#include "CaTrack.h"
#include "CaTrackFit.h"
#include "CaTrackParam.h"
@@ -23,7 +24,10 @@ using namespace cbm::algo::ca;
// ---------------------------------------------------------------------------------------------------------------------
//
-TrackExtender::TrackExtender(const ca::Framework& algo, WindowData& wData) : frAlgo(algo), frWData(wData) {}
+TrackExtender::TrackExtender(const ca::Parameters<fvec>& pars, const float mass) : fParameters(pars), fDefaultMass(mass)
+{
+}
+
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -38,7 +42,7 @@ void TrackExtender::FitBranchFast(const ca::Branch& t, TrackParamV& Tout, const
CBMCA_DEBUG_ASSERT(t.NHits >= 3);
ca::TrackFit fit;
- fit.SetParticleMass(frAlgo.GetDefaultParticleMass());
+ fit.SetParticleMass(fDefaultMass);
fit.SetMask(fmask::One());
fit.SetTrack(Tout);
TrackParamV& T = fit.Tr();
@@ -51,17 +55,17 @@ void TrackExtender::FitBranchFast(const ca::Branch& t, TrackParamV& Tout, const
const int iFirstHit = (upstream) ? nHits - 1 : 0;
const int iLastHit = (upstream) ? 0 : nHits - 1;
- const ca::Hit& hit0 = frAlgo.GetInputData().GetHit(hits[iFirstHit]);
- const ca::Hit& hit1 = frAlgo.GetInputData().GetHit(hits[iFirstHit + step]);
- const ca::Hit& hit2 = frAlgo.GetInputData().GetHit(hits[iFirstHit + 2 * step]);
+ const ca::Hit& hit0 = fInputData->GetHit(hits[iFirstHit]);
+ const ca::Hit& hit1 = fInputData->GetHit(hits[iFirstHit + step]);
+ const ca::Hit& hit2 = fInputData->GetHit(hits[iFirstHit + 2 * step]);
int ista0 = hit0.Station();
int ista1 = hit1.Station();
int ista2 = hit2.Station();
- const ca::Station<fvec>& sta0 = frAlgo.GetParameters().GetStation(ista0);
- const ca::Station<fvec>& sta1 = frAlgo.GetParameters().GetStation(ista1);
- const ca::Station<fvec>& sta2 = frAlgo.GetParameters().GetStation(ista2);
+ const ca::Station<fvec>& sta0 = fParameters.GetStation(ista0);
+ const ca::Station<fvec>& sta1 = fParameters.GetStation(ista1);
+ const ca::Station<fvec>& sta2 = fParameters.GetStation(ista2);
fvec x0 = hit0.X();
fvec y0 = hit0.Y();
@@ -109,13 +113,13 @@ void TrackExtender::FitBranchFast(const ca::Branch& t, TrackParamV& Tout, const
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
for (int i = iFirstHit + step; step * i <= step * iLastHit; i += step) {
- const ca::Hit& hit = frAlgo.GetInputData().GetHit(hits[i]);
- int ista = hit.Station();
- const ca::Station<fvec>& sta = frAlgo.GetParameters().GetStation(ista);
+ const ca::Hit& hit = fInputData->GetHit(hits[i]);
+ int ista = hit.Station();
+ const ca::Station<fvec>& sta = fParameters.GetStation(ista);
fit.Extrapolate(hit.Z(), fld);
fit.FilterHit(sta, hit);
- auto radThick = frAlgo.GetParameters().GetMaterialThickness(ista, T.X(), T.Y());
+ auto radThick = fParameters.GetMaterialThickness(ista, T.X(), T.Y());
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, fvec(upstream ? 1. : -1.));
@@ -146,29 +150,29 @@ void TrackExtender::FitBranch(const ca::Branch& t, TrackParamV& T, const bool up
void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool upstream, const fvec qp0)
{
Vector<ca::HitIndex_t> newHits{"ca::TrackExtender::newHits"};
- newHits.reserve(frAlgo.GetParameters().GetNstationsActive());
+ newHits.reserve(fParameters.GetNstationsActive());
ca::TrackFit fit;
- fit.SetParticleMass(frAlgo.GetDefaultParticleMass());
+ fit.SetParticleMass(fDefaultMass);
fit.SetMask(fmask::One());
fit.SetTrack(Tout);
fit.SetQp0(qp0);
const signed short int step = -2 * static_cast<int>(upstream) + 1; // increment for station index
const int iFirstHit = (upstream) ? 2 : t.NofHits() - 3;
- // int ista = frAlgo.GetInputData().GetHit(t.Hits[iFirstHit]).iSt + 2 * step; // current station. set to the end of track
+ // int ista = fInputData->GetHit(t.Hits[iFirstHit]).iSt + 2 * step; // current station. set to the end of track
- const ca::Hit& hit0 = frAlgo.GetInputData().GetHit(t.Hits()[iFirstHit]); // optimize
- const ca::Hit& hit1 = frAlgo.GetInputData().GetHit(t.Hits()[iFirstHit + step]);
- const ca::Hit& hit2 = frAlgo.GetInputData().GetHit(t.Hits()[iFirstHit + 2 * step]);
+ const ca::Hit& hit0 = fInputData->GetHit(t.Hits()[iFirstHit]); // optimize
+ const ca::Hit& hit1 = fInputData->GetHit(t.Hits()[iFirstHit + step]);
+ const ca::Hit& hit2 = fInputData->GetHit(t.Hits()[iFirstHit + 2 * step]);
const int ista0 = hit0.Station();
const int ista1 = hit1.Station();
const int ista2 = hit2.Station();
- const ca::Station<fvec>& sta0 = frAlgo.GetParameters().GetStation(ista0);
- const ca::Station<fvec>& sta1 = frAlgo.GetParameters().GetStation(ista1);
- const ca::Station<fvec>& sta2 = frAlgo.GetParameters().GetStation(ista2);
+ const ca::Station<fvec>& sta0 = fParameters.GetStation(ista0);
+ const ca::Station<fvec>& sta1 = fParameters.GetStation(ista1);
+ const ca::Station<fvec>& sta2 = fParameters.GetStation(ista2);
fvec x0 = hit0.X();
fvec y0 = hit0.Y();
@@ -192,14 +196,14 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
int ista = ista2 + 2 * step; // skip one station. if there would be hit it has to be found on previous stap
- if (ista2 == frWData.CurrentIteration()->GetFirstStationIndex()) ista = ista2 + step;
+ if (ista2 == frWData->CurrentIteration()->GetFirstStationIndex()) ista = ista2 + step;
- const float pickGather = frWData.CurrentIteration()->GetPickGather();
+ const float pickGather = frWData->CurrentIteration()->GetPickGather();
const fvec pickGather2 = pickGather * pickGather;
- const fvec maxDZ = frWData.CurrentIteration()->GetMaxDZ();
- for (; (ista < frAlgo.GetParameters().GetNstationsActive()) && (ista >= 0); ista += step) { // CHECKME why ista2?
+ const fvec maxDZ = frWData->CurrentIteration()->GetMaxDZ();
+ for (; (ista < fParameters.GetNstationsActive()) && (ista >= 0); ista += step) { // CHECKME why ista2?
- const ca::Station<fvec>& sta = frAlgo.GetParameters().GetStation(ista);
+ const ca::Station<fvec>& sta = fParameters.GetStation(ista);
fit.Extrapolate(sta.fZ, fld);
@@ -208,12 +212,12 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
TrackParamV& tr = fit.Tr();
- const auto& grid = frWData.Grid(ista);
+ const auto& grid = frWData->Grid(ista);
ca::GridArea area(grid, tr.X()[0], tr.Y()[0],
(sqrt(pickGather * tr.C00()) + grid.GetMaxRangeX() + maxDZ * abs(tr.Tx()))[0],
(sqrt(pickGather * tr.C11()) + grid.GetMaxRangeY() + maxDZ * abs(tr.Ty()))[0]);
- if (frAlgo.GetParameters().DevIsIgnoreHitSearchAreas()) {
+ if (fParameters.DevIsIgnoreHitSearchAreas()) {
area.DoLoopOverEntireGrid();
}
@@ -221,10 +225,10 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
while (area.GetNextObjectId(ih)) { // loop over the hits in the area
- if (frWData.IsHitSuppressed(ih)) {
+ if (frWData->IsHitSuppressed(ih)) {
continue;
}
- const ca::Hit& hit = frWData.Hit(ih);
+ const ca::Hit& hit = frWData->Hit(ih);
if (sta.timeInfo && tr.NdfTime()[0] > -2.) {
fscal dt = hit.T() - tr.Time()[0];
@@ -233,7 +237,7 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
//if (GetFUsed((*fStripFlag)[hit.f] | (*fStripFlag)[hit.b])) continue; // if used
- if (frWData.IsHitKeyUsed(hit.FrontKey()) || frWData.IsHitKeyUsed(hit.BackKey())) {
+ if (frWData->IsHitKeyUsed(hit.FrontKey()) || frWData->IsHitKeyUsed(hit.BackKey())) {
continue;
}
@@ -262,12 +266,12 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
if (iHit_best < 0) break;
- const ca::Hit& hit = frWData.Hit(iHit_best);
+ const ca::Hit& hit = frWData->Hit(iHit_best);
newHits.push_back(hit.Id());
fit.Extrapolate(hit.Z(), fld);
fit.FilterHit(sta, hit);
- auto radThick = frAlgo.GetParameters().GetMaterialThickness(ista, tr.X(), tr.Y());
+ auto radThick = fParameters.GetMaterialThickness(ista, tr.X(), tr.Y());
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, upstream ? fvec(1.f) : fvec(-1.f));
@@ -304,8 +308,10 @@ void TrackExtender::FindMoreHits(ca::Branch& t, TrackParamV& Tout, const bool up
} // void ca::Framework::FindMoreHits
/// Try to extrapolate and find additional hits on other stations
-fscal TrackExtender::ExtendBranch(ca::Branch& t)
+fscal TrackExtender::ExtendBranch(ca::Branch& t, const InputData& input, WindowData& wData)
{
+ fInputData = &input;
+ frWData = &wData;
// const unsigned int minNHits = 3;
TrackParamV T;
diff --git a/algo/ca/core/tracking/CaTrackExtender.h b/algo/ca/core/tracking/CaTrackExtender.h
index e5762fb67e..875db7cb35 100644
--- a/algo/ca/core/tracking/CaTrackExtender.h
+++ b/algo/ca/core/tracking/CaTrackExtender.h
@@ -11,16 +11,17 @@
#include "CaBranch.h"
#include "CaHit.h"
+#include "CaParameters.h"
#include "CaSimd.h"
#include "CaTrackParam.h"
#include "CaVector.h"
#include "CaWindowData.h"
-
namespace cbm::algo::ca
{
class Track;
class Framework;
+ class InputData;
namespace
{
@@ -33,7 +34,7 @@ namespace cbm::algo::ca
class TrackExtender {
public:
/// Default constructor
- TrackExtender(const ca::Framework& algo, WindowData& wData);
+ TrackExtender(const ca::Parameters<fvec>& pars, const float mass);
/// Destructor
~TrackExtender();
@@ -53,7 +54,7 @@ namespace cbm::algo::ca
/// Find additional hits for existing track both downstream and upstream
/// \return chi2
- fscal ExtendBranch(ca::Branch& t);
+ fscal ExtendBranch(ca::Branch& t, const ca::InputData& input, WindowData& wData);
private:
///-------------------------------
@@ -88,8 +89,10 @@ namespace cbm::algo::ca
///-------------------------------
/// Data members
- const ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class
- WindowData& frWData;
+ const Parameters<fvec>& fParameters; ///< Object of Framework parameters class
+ const InputData* fInputData; ///< Tracking input data
+ WindowData* frWData;
+ float fDefaultMass{constants::phys::MuonMass}; ///< mass of the propagated particle [GeV/c2]
};
} // namespace cbm::algo::ca
diff --git a/algo/ca/core/tracking/CaTrackFinder.cxx b/algo/ca/core/tracking/CaTrackFinder.cxx
index 1892ea5fbb..6afe5987e1 100644
--- a/algo/ca/core/tracking/CaTrackFinder.cxx
+++ b/algo/ca/core/tracking/CaTrackFinder.cxx
@@ -18,39 +18,66 @@
*
*/
-#include "CaFramework.h"
+#include "CaTrackFinder.h"
+
#include "CaTimer.h"
#include "CaTrack.h"
+#include "CaTriplet.h"
#include <chrono>
#include <fstream>
#include <sstream>
#include <thread>
+
using namespace cbm::algo::ca;
-using cbm::algo::ca::ECounter;
-using cbm::algo::ca::ETimer;
+using cbm::algo::ca::ECounter; // monitor counter key type
+using cbm::algo::ca::ETimer; // monitor timer key type
using cbm::algo::ca::Track;
+using constants::phys::ProtonMassD;
+using constants::phys::SpeedOfLightInv;
+using constants::phys::SpeedOfLightInvD;
// ---------------------------------------------------------------------------------------------------------------------
//
-TrackFinder::TrackFinder(ca::Framework& algo) : frAlgo(algo) {}
-// ---------------------------------------------------------------------------------------------------------------------
-//
-TrackFinder::~TrackFinder() {}
+TrackFinder::TrackFinder(const ca::Parameters<fvec>& pars, const float mass, const ca::TrackingMode& mode,
+ TrackingMonitorData& monitorData, int nThreads, double& recoTime)
+ : fParameters(pars)
+ , fDefaultMass(mass)
+ , fTrackingMode(mode)
+ , fMonitorData(monitorData)
+ , fvMonitorDataThread(nThreads)
+ , fvWData(nThreads)
+ , fNofThreads(nThreads)
+ , fCaRecoTime(recoTime)
+ , fvRecoTracks(nThreads)
+ , fvRecoHitIndices(nThreads)
+ , fWindowLength((ca::TrackingMode::kMcbm == mode) ? 500 : 10000)
+{
+ assert(fNofThreads > 0);
+
+ for (int iThread = 0; iThread < fNofThreads; ++iThread) {
+ fvRecoTracks[iThread].SetName(std::string("TrackFinder::fvRecoTracks_") + std::to_string(iThread));
+ fvRecoHitIndices[iThread].SetName(std::string("TrackFinder::fvRecoHitIndices_") + std::to_string(iThread));
+ }
+}
// ---------------------------------------------------------------------------------------------------------------------
+//CBM Level 1 4D Reconstruction
+//Finds tracks using the Cellular Automaton algorithm
//
-void TrackFinder::FindTracks()
+TrackFinder::Output_t TrackFinder::FindTracks(const InputData& input, TimesliceHeader& tsHeader)
{
- frAlgo.fRecoTracks.clear();
- frAlgo.fRecoHits.clear();
- if (frAlgo.GetInputData().GetNhits() < 1) {
+ Output_t output;
+ Vector<Track>& recoTracks = output.first; //reconstructed tracks
+ Vector<ca::HitIndex_t>& recoHits = output.second; //packed hits of reconstructed tracks
+
+ if (input.GetNhits() < 1) {
LOG(warn) << "No hits were passed to the ca::TrackFinder. Stopping the routine";
- return;
+ return output;
}
//
@@ -58,74 +85,60 @@ void TrackFinder::FindTracks()
// It splits the input data into sub-timeslices
// and runs the track finder over the sub-slices
//
- frAlgo.fMonitorData.StartTimer(ETimer::Tracking);
- frAlgo.fMonitorData.StartTimer(ETimer::PrepareTimeslice);
- frAlgo.fMonitorData.IncrementCounter(ECounter::TrackingCall);
- frAlgo.fMonitorData.IncrementCounter(ECounter::RecoHit, frAlgo.fInputData.GetNhits());
+ fMonitorData.StartTimer(ETimer::Tracking);
+ fMonitorData.StartTimer(ETimer::PrepareTimeslice);
+ fMonitorData.IncrementCounter(ECounter::TrackingCall);
+ fMonitorData.IncrementCounter(ECounter::RecoHit, input.GetNhits());
auto timerStart = std::chrono::high_resolution_clock::now();
- auto& wDataThread0 = frAlgo.fvWData[0]; // NOTE: Thread 0 must be always defined
+ auto& wDataThread0 = fvWData[0]; // NOTE: Thread 0 must be always defined
// ----- Reset data arrays -------------------------------------------------------------------------------------------
- wDataThread0.HitKeyFlags().reset(frAlgo.fInputData.GetNhitKeys(), 0);
-
- frAlgo.fHitTimeInfo.reset(frAlgo.fInputData.GetNhits());
+ wDataThread0.HitKeyFlags().reset(input.GetNhitKeys(), 0);
+ fHitTimeInfo.reset(input.GetNhits());
// TODO: move these values to Parameters namespace (S.Zharko)
// length of sub-TS
- fscal minProtonMomentum = 0.1;
-
- fscal targX = frAlgo.GetParameters().GetTargetPositionX()[0];
- fscal targY = frAlgo.GetParameters().GetTargetPositionY()[0];
- fscal targZ = frAlgo.GetParameters().GetTargetPositionZ()[0];
+ const fscal minProtonMomentum = 0.1;
+ const fscal preFactor = sqrt(1. + ProtonMassD * ProtonMassD / (minProtonMomentum * minProtonMomentum));
+ const fscal targX = fParameters.GetTargetPositionX()[0];
+ const fscal targY = fParameters.GetTargetPositionY()[0];
+ const fscal targZ = fParameters.GetTargetPositionZ()[0];
fStatTsStart = std::numeric_limits<fscal>::max(); // end time of the TS
fStatTsEnd = 0.; // end time of the TS
fStatNhitsTotal = 0;
- // calculate possible event time for the hits (frAlgo.fHitTimeInfo array)
-
- const int nDataStreams = frAlgo.fInputData.GetNdataStreams();
+ // calculate possible event time for the hits (fHitTimeInfo array)
+ const int nDataStreams = input.GetNdataStreams();
for (int iStream = 0; iStream < nDataStreams; ++iStream) {
fscal maxTimeBeforeHit = std::numeric_limits<fscal>::lowest();
-
- int nStreamHits = frAlgo.fInputData.GetStreamNhits(iStream);
+ const int nStreamHits = input.GetStreamNhits(iStream);
fStatNhitsTotal += nStreamHits;
for (int ih = 0; ih < nStreamHits; ++ih) {
- ca::HitIndex_t caHitId = frAlgo.fInputData.GetStreamStartIndex(iStream) + ih;
- const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId);
-
- const ca::Station<fvec>& st = frAlgo.GetParameters().GetStation(h.Station());
-
- fscal dx = h.X() - targX;
- fscal dy = h.Y() - targY;
- fscal dz = h.Z() - targZ;
- fscal l = sqrt(dx * dx + dy * dy + dz * dz);
-
- fscal timeOfFlightMin = l * constants::phys::SpeedOfLightInv;
- fscal timeOfFlightMax =
- 1.5 * l
- * sqrt(1. + constants::phys::ProtonMassD * constants::phys::ProtonMassD / minProtonMomentum / minProtonMomentum)
- * constants::phys::SpeedOfLightInvD;
+ ca::HitIndex_t caHitId = input.GetStreamStartIndex(iStream) + ih;
+ const ca::Hit& h = input.GetHit(caHitId);
+ const ca::Station<fvec>& st = fParameters.GetStation(h.Station());
+ const fscal dx = h.X() - targX;
+ const fscal dy = h.Y() - targY;
+ const fscal dz = h.Z() - targZ;
+ const fscal l = sqrt(dx * dx + dy * dy + dz * dz);
+ const fscal timeOfFlightMin = l * SpeedOfLightInv;
+ const fscal timeOfFlightMax = 1.5 * l * preFactor * SpeedOfLightInvD;
+ const fscal dt = h.RangeT();
// TODO: Is it possible, that the proton mass selection affects the search of heavier particles?
- fscal dt = h.RangeT();
-
- CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId];
- info.fEventTimeMin = h.T() - dt - timeOfFlightMax;
- info.fEventTimeMax = h.T() + dt - timeOfFlightMin;
- if (!st.timeInfo) {
- info.fEventTimeMin = -1.e10;
- info.fEventTimeMax = 1.e10;
- }
+ CaHitTimeInfo& info = fHitTimeInfo[caHitId];
+ info.fEventTimeMin = st.timeInfo ? (h.T() - dt - timeOfFlightMax) : -1.e10;
+ info.fEventTimeMax = st.timeInfo ? (h.T() + dt - timeOfFlightMin) : 1.e10;
// NOTE: if not a MT part, use wDataThread0.IsHitKeyUsed, it will be later copied to other threads
if (info.fEventTimeMin > 500.e6 || info.fEventTimeMax < -500.) { // cut hits with bogus start time > 500 ms
@@ -134,33 +147,23 @@ void TrackFinder::FindTracks()
LOG(error) << "CATrackFinder: skip bogus hit " << h.ToString();
continue;
}
-
- if (maxTimeBeforeHit < info.fEventTimeMax) {
- maxTimeBeforeHit = info.fEventTimeMax;
- }
+ maxTimeBeforeHit = std::max(maxTimeBeforeHit, info.fEventTimeMax);
info.fMaxTimeBeforeHit = maxTimeBeforeHit;
-
- if (fStatTsStart > info.fEventTimeMax - fWindowMargin) {
- fStatTsStart = info.fEventTimeMax - fWindowMargin; // 5 ns margin
- }
- if (fStatTsEnd < info.fEventTimeMin) {
- fStatTsEnd = info.fEventTimeMin;
- }
+ fStatTsStart = std::min(fStatTsStart, info.fEventTimeMax - fWindowMargin); // 5 ns margin
+ fStatTsEnd = std::max(fStatTsEnd, info.fEventTimeMin);
}
fscal minTimeAfterHit = std::numeric_limits<fscal>::max();
// loop in the reverse order to fill CaHitTimeInfo::fMinTimeAfterHit fields
for (int ih = nStreamHits - 1; ih >= 0; --ih) {
- ca::HitIndex_t caHitId = frAlgo.fInputData.GetStreamStartIndex(iStream) + ih;
- const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId);
+ ca::HitIndex_t caHitId = input.GetStreamStartIndex(iStream) + ih;
+ const ca::Hit& h = input.GetHit(caHitId);
if (wDataThread0.IsHitKeyUsed(h.FrontKey()) || wDataThread0.IsHitKeyUsed(h.BackKey())) {
continue;
} // the hit is skipped
- CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId];
- if (minTimeAfterHit > info.fEventTimeMin) {
- minTimeAfterHit = info.fEventTimeMin;
- }
+ CaHitTimeInfo& info = fHitTimeInfo[caHitId];
+ minTimeAfterHit = std::min(minTimeAfterHit, info.fEventTimeMin);
info.fMinTimeAfterHit = minTimeAfterHit;
}
@@ -170,12 +173,12 @@ void TrackFinder::FindTracks()
tmp++;
LOG(warning) << "\n\n stream " << iStream << " hits " << nStreamHits << "\n\n";
for (int ih = 0; (ih < nStreamHits) && (tmp < 10000); ++ih) {
- ca::HitIndex_t caHitId = frAlgo.fInputData.GetStreamStartIndex(iStream) + ih;
- const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId);
+ ca::HitIndex_t caHitId = input.GetStreamStartIndex(iStream) + ih;
+ const ca::Hit& h = input.GetHit(caHitId);
if (wDataThread0.IsHitKeyUsed(h.FrontKey()) || wDataThread0.IsHitKeyUsed(h.BackKey())) {
continue;
} // the hit is skipped
- CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId];
+ CaHitTimeInfo& info = fHitTimeInfo[caHitId];
if (h.Station() < 4) {
tmp++;
LOG(warning) << " hit sta " << h.Station() << " stream " << iStream << " time " << h.T() << " event time "
@@ -187,14 +190,8 @@ void TrackFinder::FindTracks()
}
}
-
- if (fStatTsEnd < fStatTsStart) { // all hits belong to one sub-timeslice
- fStatTsEnd = fStatTsStart;
- }
-
- if (fStatTsEnd > fStatTsStart + 500.e6) { // 500 ms maximal length of the TS
- fStatTsEnd = fStatTsStart + 500.e6;
- }
+ fStatTsEnd = std::max(fStatTsEnd, fStatTsStart); // all hits belong to one sub-timeslice
+ fStatTsEnd = std::min(fStatTsEnd, fStatTsStart + 500.e6f); // 500 ms maximal length of the TS
LOG(debug) << "CA tracker process time slice " << fStatTsStart * 1.e-6 << " -- " << fStatTsEnd * 1.e-6
<< " [ms] with " << fStatNhitsTotal << " hits";
@@ -206,25 +203,26 @@ void TrackFinder::FindTracks()
nWindows = 1;
}
- int nWindowsThread = nWindows / frAlgo.GetNofThreads();
+ int nWindowsThread = nWindows / fNofThreads;
float windowDelta = fWindowLength - fWindowOverlap - fWindowMargin;
//LOG(info) << "CA: estimated number of time windows: " << nWindows;
+ std::vector<float> vWindowStartThread(fNofThreads), vWindowEndThread(fNofThreads);
- fvWindowStartThread[0] = fStatTsStart;
+ vWindowStartThread[0] = fStatTsStart;
{ // Estimation of number of hits in time windows
//Timer time;
//time.Start();
- int nSt = frAlgo.GetParameters().GetNstationsActive();
+ int nSt = fParameters.GetNstationsActive();
// Count number of hits per window and station
- double a = fStatTsStart + fWindowOverlap + fWindowMargin;
- double b = fWindowLength - fWindowOverlap - fWindowMargin;
- HitIndex_t nHitsTot = frAlgo.fInputData.GetNhits();
+ const double a = fStatTsStart + fWindowOverlap + fWindowMargin;
+ const double b = fWindowLength - fWindowOverlap - fWindowMargin;
+ HitIndex_t nHitsTot = input.GetNhits();
std::vector<HitIndex_t> nHitsWindowSta(nWindows * nSt, 0);
for (HitIndex_t iHit = 0; iHit < nHitsTot; ++iHit) {
- const auto& hit = frAlgo.fInputData.GetHit(iHit);
- //const auto& timeInfo = frAlgo.fHitTimeInfo[iHit];
+ const auto& hit = input.GetHit(iHit);
+ //const auto& timeInfo = fHitTimeInfo[iHit];
size_t iWindow = static_cast<size_t>((hit.T() - a) / b);
if (iWindow >= static_cast<size_t>(nWindows)) {
LOG(error) << "ca: Hit out of range: iHit = " << iHit << ", time = " << hit.T() * 1.e-6
@@ -236,7 +234,7 @@ void TrackFinder::FindTracks()
}
// Remove hits from the "monster events"
- if (ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) {
+ if (ca::TrackingMode::kMcbm == fTrackingMode) {
auto maxNofHitsSta = static_cast<HitIndex_t>(50 * fWindowLength / 1.e3);
for (auto& content : nHitsWindowSta) {
if (content > maxNofHitsSta) {
@@ -248,22 +246,20 @@ void TrackFinder::FindTracks()
// Integrate number of hits per window
std::vector<HitIndex_t> nHitsWindow;
nHitsWindow.reserve(nWindows);
- auto windowStaIt = nHitsWindowSta.begin();
HitIndex_t nHitsCollected = 0;
- while (windowStaIt != nHitsWindowSta.end()) {
- nHitsCollected = nHitsWindow.emplace_back(std::accumulate(windowStaIt, windowStaIt + nSt, nHitsCollected));
- windowStaIt = windowStaIt + nSt;
+ for (auto it = nHitsWindowSta.begin(); it != nHitsWindowSta.end(); it += nSt) {
+ nHitsCollected = nHitsWindow.emplace_back(std::accumulate(it, it + nSt, nHitsCollected));
}
// Get time range for threads
- HitIndex_t nHitsPerThread = nHitsCollected / frAlgo.GetNofThreads();
+ HitIndex_t nHitsPerThread = nHitsCollected / fNofThreads;
auto windowIt = nHitsWindow.begin();
size_t iWbegin = 0;
- for (int iTh = 1; iTh < frAlgo.GetNofThreads(); ++iTh) {
- windowIt = std::lower_bound(windowIt, nHitsWindow.end(), iTh * nHitsPerThread);
- iWbegin = std::distance(nHitsWindow.begin(), windowIt) + 1;
- fvWindowStartThread[iTh] = fStatTsStart + iWbegin * windowDelta;
- fvWindowEndThread[iTh - 1] = fvWindowStartThread[iTh];
+ for (int iTh = 1; iTh < fNofThreads; ++iTh) {
+ windowIt = std::lower_bound(windowIt, nHitsWindow.end(), iTh * nHitsPerThread);
+ iWbegin = std::distance(nHitsWindow.begin(), windowIt) + 1;
+ vWindowStartThread[iTh] = fStatTsStart + iWbegin * windowDelta;
+ vWindowEndThread[iTh - 1] = vWindowStartThread[iTh];
}
//time.Stop();
//LOG(info) << "Thread boarders estimation time: " << time.GetTotalMs() << " ms";
@@ -272,120 +268,114 @@ void TrackFinder::FindTracks()
// cut data into sub-timeslices and process them one by one
- //for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) {
- // fvWindowStartThread[iThread] = fStatTsStart + iThread * nWindowsThread * windowDelta;
- // fvWindowEndThread[iThread] =
- // fvWindowStartThread[iThread] + nWindowsThread * windowDelta + fWindowOverlap + fWindowMargin;
+ //for (int iThread = 0; iThread < fNofThreads; ++iThread) {
+ // vWindowStartThread[iThread] = fStatTsStart + iThread * nWindowsThread * windowDelta;
+ // vWindowEndThread[iThread] =
+ // vWindowStartThread[iThread] + nWindowsThread * windowDelta + fWindowOverlap + fWindowMargin;
//}
- fvWindowEndThread[frAlgo.GetNofThreads() - 1] = fStatTsEnd;
+ vWindowEndThread[fNofThreads - 1] = fStatTsEnd;
- for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) {
- double start = fvWindowStartThread[iThread] * 1.e-6;
- double end = fvWindowEndThread[iThread] * 1.e-6;
+ for (int iThread = 0; iThread < fNofThreads; ++iThread) {
+ double start = vWindowStartThread[iThread] * 1.e-6;
+ double end = vWindowEndThread[iThread] * 1.e-6;
LOG(trace) << "Thread: " << iThread << " from " << start << " ms to " << end << " ms (delta = " << end - start
<< " ms)";
}
- frAlgo.fMonitorData.StopTimer(ETimer::PrepareTimeslice);
+ // Statistics for monitoring
+ std::vector<int> vStatNwindows(fNofThreads), vStatNhitsProcessed(fNofThreads);
+
+ fMonitorData.StopTimer(ETimer::PrepareTimeslice);
// Save tracks
- if (frAlgo.GetNofThreads() == 1) {
- this->FindTracksThread(0);
- frAlgo.fMonitorData.StartTimer(ETimer::StoreTracksFinal);
- frAlgo.fRecoTracks = std::move(fvRecoTracks[0]);
- frAlgo.fRecoHits = std::move(fvRecoHitIndices[0]);
- frAlgo.fMonitorData.StopTimer(ETimer::StoreTracksFinal);
+ if (fNofThreads == 1) {
+ this->FindTracksThread(input, 0, std::ref(vWindowStartThread[0]), std::ref(vWindowEndThread[0]),
+ std::ref(vStatNwindows[0]), std::ref(vStatNhitsProcessed[0]));
+ fMonitorData.StartTimer(ETimer::StoreTracksFinal);
+ recoTracks = std::move(fvRecoTracks[0]);
+ recoHits = std::move(fvRecoHitIndices[0]);
+ fMonitorData.StopTimer(ETimer::StoreTracksFinal);
}
else {
std::vector<std::thread> vThreadList;
- vThreadList.reserve(frAlgo.GetNofThreads());
- for (int iTh = 0; iTh < frAlgo.GetNofThreads(); ++iTh) {
- vThreadList.emplace_back(&TrackFinder::FindTracksThread, this, iTh);
+ vThreadList.reserve(fNofThreads);
+ for (int iTh = 0; iTh < fNofThreads; ++iTh) {
+ vThreadList.emplace_back(&TrackFinder::FindTracksThread, this, std::ref(input), iTh,
+ std::ref(vWindowStartThread[iTh]), std::ref(vWindowEndThread[iTh]),
+ std::ref(vStatNwindows[iTh]), std::ref(vStatNhitsProcessed[iTh]));
}
for (auto& th : vThreadList) {
if (th.joinable()) {
th.join();
}
}
- frAlgo.fMonitorData.StartTimer(ETimer::StoreTracksFinal);
+ fMonitorData.StartTimer(ETimer::StoreTracksFinal);
auto Operation = [](size_t acc, const auto& v) { return acc + v.size(); };
int nRecoTracks = std::accumulate(fvRecoTracks.begin(), fvRecoTracks.end(), 0, Operation);
int nRecoHits = std::accumulate(fvRecoHitIndices.begin(), fvRecoHitIndices.end(), 0, Operation);
- frAlgo.fRecoTracks.reserve(nRecoTracks);
- frAlgo.fRecoHits.reserve(nRecoHits);
- for (int iTh = 0; iTh < frAlgo.GetNofThreads(); ++iTh) {
- frAlgo.fRecoTracks.insert(frAlgo.fRecoTracks.end(), fvRecoTracks[iTh].begin(), fvRecoTracks[iTh].end());
- frAlgo.fRecoHits.insert(frAlgo.fRecoHits.end(), fvRecoHitIndices[iTh].begin(), fvRecoHitIndices[iTh].end());
+ recoTracks.reserve(nRecoTracks);
+ recoHits.reserve(nRecoHits);
+ for (int iTh = 0; iTh < fNofThreads; ++iTh) {
+ recoTracks.insert(recoTracks.end(), fvRecoTracks[iTh].begin(), fvRecoTracks[iTh].end());
+ recoHits.insert(recoHits.end(), fvRecoHitIndices[iTh].begin(), fvRecoHitIndices[iTh].end());
}
- frAlgo.fMonitorData.StopTimer(ETimer::StoreTracksFinal);
+ fMonitorData.StopTimer(ETimer::StoreTracksFinal);
}
- frAlgo.fMonitorData.IncrementCounter(ECounter::RecoTrack, frAlgo.fRecoTracks.size());
- frAlgo.fMonitorData.IncrementCounter(ECounter::RecoHitUsed, frAlgo.fRecoHits.size());
+ fMonitorData.IncrementCounter(ECounter::RecoTrack, recoTracks.size());
+ fMonitorData.IncrementCounter(ECounter::RecoHitUsed, recoHits.size());
- auto timerEnd = std::chrono::high_resolution_clock::now();
- frAlgo.fCaRecoTime = (double) (std::chrono::duration<double>(timerEnd - timerStart).count());
+ auto timerEnd = std::chrono::high_resolution_clock::now();
+ fCaRecoTime = (double) (std::chrono::duration<double>(timerEnd - timerStart).count());
// Add thread monitors to the main monitor
- for (auto& monitor : frAlgo.fvMonitorDataThread) {
- frAlgo.fMonitorData.AddMonitorData(monitor, true);
- //frAlgo.fMonitorData.AddMonitorData(monitor);
+ for (auto& monitor : fvMonitorDataThread) {
+ fMonitorData.AddMonitorData(monitor, true);
+ //fMonitorData.AddMonitorData(monitor);
monitor.Reset();
}
- int statNhitsProcessedTotal = 0;
- int statNwindowsTotal = 0;
- for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) {
- statNhitsProcessedTotal += fvStatNhitsProcessed[iThread];
- statNwindowsTotal += fvStatNwindows[iThread];
- }
+ const int statNhitsProcessedTotal = std::accumulate(vStatNhitsProcessed.begin(), vStatNhitsProcessed.end(), 0);
+ const int statNwindowsTotal = std::accumulate(vStatNwindows.begin(), vStatNwindows.end(), 0);
// Filling TS headear
- auto& tsHeader = frAlgo.fTsHeader;
tsHeader.Start() = fStatTsStart;
tsHeader.End() = fStatTsEnd;
- frAlgo.fMonitorData.StopTimer(ETimer::Tracking);
+ fMonitorData.StopTimer(ETimer::Tracking);
- LOG(debug) << "CA tracker: time slice finished. Reconstructed " << frAlgo.fRecoTracks.size() << " tracks with "
- << frAlgo.fRecoHits.size() << " hits. Processed " << statNhitsProcessedTotal << " hits in "
- << statNwindowsTotal << " time windows. Reco time " << frAlgo.fCaRecoTime / 1.e9 << " s";
+ LOG(debug) << "CA tracker: time slice finished. Reconstructed " << recoTracks.size() << " tracks with "
+ << recoHits.size() << " hits. Processed " << statNhitsProcessedTotal << " hits in " << statNwindowsTotal
+ << " time windows. Reco time " << fCaRecoTime / 1.e9 << " s";
+ return output;
}
// ---------------------------------------------------------------------------------------------------------------------
//
-void TrackFinder::FindTracksThread(int iThread)
+void TrackFinder::FindTracksThread(const InputData& input, int iThread, float& windowStart, float& windowEnd,
+ int& statNwindows, int& statNhitsProcessed)
{
//std::stringstream filename;
//filename << "./dbg_caTrackFinder::FindTracksThread_" << iThread << ".txt";
//std::ofstream out(filename.str());
- auto& monitor = frAlgo.fvMonitorDataThread[iThread];
+ auto& monitor = fvMonitorDataThread[iThread];
Timer timer;
timer.Start();
monitor.StartTimer(ETimer::TrackingThread);
monitor.StartTimer(ETimer::PrepareThread);
// Init vectors
- auto& wData = frAlgo.fvWData[iThread];
+ auto& wData = fvWData[iThread];
{
- int nStations = frAlgo.GetParameters().GetNstationsActive();
- size_t nHitsTot = frAlgo.fInputData.GetNhits();
- size_t nHitKeys = frAlgo.fInputData.GetNhitKeys();
- size_t nHitsExpected = 2 * nHitsTot;
- size_t nTracksExpected = 2 * nHitsTot / nStations;
- size_t nThreads = frAlgo.GetNofThreads();
+ const int nStations = fParameters.GetNstationsActive();
+ const size_t nHitsTot = input.GetNhits();
+ const size_t nHitsExpected = 2 * nHitsTot;
+ const size_t nTracksExpected = 2 * nHitsTot / nStations;
fvRecoTracks[iThread].clear();
- fvRecoTracks[iThread].reserve(nTracksExpected / nThreads);
+ fvRecoTracks[iThread].reserve(nTracksExpected / fNofThreads);
fvRecoHitIndices[iThread].clear();
- fvRecoHitIndices[iThread].reserve(nHitsExpected / nThreads);
- frAlgo.fvTrackFinderWindow[iThread].InitTimeslice();
+ fvRecoHitIndices[iThread].reserve(nHitsExpected / fNofThreads);
if (iThread != 0) {
- auto& wData0 = frAlgo.fvWData[0];
- // slow ...
- wData.HitKeyFlags().clear();
- wData.HitKeyFlags().reserve(nHitKeys);
- for (const auto& flag : wData0.HitKeyFlags()) {
- wData.HitKeyFlags().emplace_back(flag);
- }
+ wData.HitKeyFlags() = fvWData[0].HitKeyFlags();
}
for (int iS = 0; iS < nStations; ++iS) {
wData.TsHitIndices(iS).clear();
@@ -393,78 +383,76 @@ void TrackFinder::FindTracksThread(int iThread)
}
}
- const int nDataStreams = frAlgo.fInputData.GetNdataStreams();
- bool areUntouchedDataLeft = true; // is the whole TS processed
+ const int nDataStreams = input.GetNdataStreams();
std::vector<HitIndex_t> sliceFirstHit(nDataStreams, 0);
std::vector<HitIndex_t> sliceLastHit(nDataStreams, 0);
// Define first hit, skip all the hits, which are before the first window
for (int iStream = 0; iStream < nDataStreams; ++iStream) {
- sliceFirstHit[iStream] = frAlgo.fInputData.GetStreamStartIndex(iStream);
- sliceLastHit[iStream] = frAlgo.fInputData.GetStreamStopIndex(iStream);
+ sliceFirstHit[iStream] = input.GetStreamStartIndex(iStream);
+ sliceLastHit[iStream] = input.GetStreamStopIndex(iStream);
for (HitIndex_t caHitId = sliceFirstHit[iStream]; caHitId < sliceLastHit[iStream]; ++caHitId) {
- CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId];
- const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId);
+ const ca::Hit& h = input.GetHit(caHitId);
if (wData.IsHitKeyUsed(h.FrontKey()) || wData.IsHitKeyUsed(h.BackKey())) {
continue;
}
- if (info.fMaxTimeBeforeHit < fvWindowStartThread[iThread]) {
+ const CaHitTimeInfo& info = fHitTimeInfo[caHitId];
+ if (info.fMaxTimeBeforeHit < windowStart) {
sliceFirstHit[iStream] = caHitId + 1;
}
- if (info.fMinTimeAfterHit > fvWindowEndThread[iThread]) {
+ if (info.fMinTimeAfterHit > windowEnd) {
sliceLastHit[iStream] = caHitId;
}
}
}
- fvStatNwindows[iThread] = 0;
- fvStatNhitsProcessed[iThread] = 0;
-
int statLastLogTimeChunk = -1;
+ // Track finder algorithm for the time window
+ ca::TrackFinderWindow trackFinderWindow(fParameters, fDefaultMass, fTrackingMode, fvMonitorDataThread[iThread]);
+ trackFinderWindow.InitTimeslice(input.GetNhitKeys());
+
monitor.StopTimer(ETimer::PrepareThread);
+ bool areUntouchedDataLeft = true; // is the whole TS processed
+
while (areUntouchedDataLeft) {
- frAlgo.fvMonitorDataThread[iThread].IncrementCounter(ECounter::SubTS);
+ fvMonitorDataThread[iThread].IncrementCounter(ECounter::SubTS);
// select the sub-slice hits
- for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); ++iS) {
- frAlgo.fvWData[iThread].TsHitIndices(iS).clear();
+ for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
+ wData.TsHitIndices(iS).clear();
}
-
areUntouchedDataLeft = false;
// TODO: SG: skip empty regions and start the subslice with the earliest hit
- fvStatNwindows[iThread]++;
- //out << fvStatNwindows[iThread] << ' ';
+ statNwindows++;
+ //out << statNwindows << ' ';
int statNwindowHits = 0;
- frAlgo.fvMonitorDataThread[iThread].StartTimer(ETimer::PrepareWindow);
+ fvMonitorDataThread[iThread].StartTimer(ETimer::PrepareWindow);
for (int iStream = 0; iStream < nDataStreams; ++iStream) {
for (ca::HitIndex_t caHitId = sliceFirstHit[iStream]; caHitId < sliceLastHit[iStream]; ++caHitId) {
- const CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId];
- const ca::Hit& h = frAlgo.fInputData.GetHit(caHitId);
+ const CaHitTimeInfo& info = fHitTimeInfo[caHitId];
+ const ca::Hit& h = input.GetHit(caHitId);
if (wData.IsHitKeyUsed(h.FrontKey()) || wData.IsHitKeyUsed(h.BackKey())) {
// the hit is already reconstructed
continue;
}
- if (info.fEventTimeMin
- > fvWindowStartThread[iThread] + fWindowLength) { // the hit is too late for the sub slice
+ if (info.fEventTimeMin > windowStart + fWindowLength) { // the hit is too late for the sub slice
areUntouchedDataLeft = true;
- if (info.fMinTimeAfterHit > fvWindowStartThread[iThread] + fWindowLength) {
+ if (info.fMinTimeAfterHit > windowStart + fWindowLength) {
// this hit and all later hits are out of the sub-slice
break;
}
}
- else {
- if (fvWindowStartThread[iThread] <= info.fEventTimeMax) { // the hit belongs to the sub-slice
- frAlgo.fvWData[iThread].TsHitIndices(h.Station()).push_back(caHitId);
- statNwindowHits++;
- if (info.fMaxTimeBeforeHit < fvWindowStartThread[iThread] + fWindowLength - fWindowOverlap) {
- // this hit and all hits before are before the overlap
- sliceFirstHit[iStream] = caHitId + 1;
- }
+ else if (windowStart <= info.fEventTimeMax) { // the hit belongs to the sub-slice
+ wData.TsHitIndices(h.Station()).push_back(caHitId);
+ statNwindowHits++;
+ if (info.fMaxTimeBeforeHit < windowStart + fWindowLength - fWindowOverlap) {
+ // this hit and all hits before are before the overlap
+ sliceFirstHit[iStream] = caHitId + 1;
}
} // else the hit has been alread processed in previous sub-slices
}
@@ -472,40 +460,38 @@ void TrackFinder::FindTracksThread(int iThread)
//out << statNwindowHits << ' ';
//if (statNwindowHits == 0) { // Empty window
- // frAlgo.fvMonitorDataThread[iThread].StopTimer(ETimer::PrepareWindow);
+ // fvMonitorDataThread[iThread].StopTimer(ETimer::PrepareWindow);
// out << 0 << ' ' << 0 << ' ' << 0 << '\n';
// continue;
//}
- if (ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) {
+ if (ca::TrackingMode::kMcbm == fTrackingMode) {
// cut at 50 hits per station per 1 us.
int maxStationHits = (int) (50 * fWindowLength / 1.e3);
- for (int ista = 0; ista < frAlgo.GetParameters().GetNstationsActive(); ++ista) {
- int nHitsSta = static_cast<int>(frAlgo.fvWData[iThread].TsHitIndices(ista).size());
+ for (int ista = 0; ista < fParameters.GetNstationsActive(); ++ista) {
+ int nHitsSta = static_cast<int>(wData.TsHitIndices(ista).size());
if (nHitsSta > maxStationHits) {
statNwindowHits -= nHitsSta;
- frAlgo.fvWData[iThread].TsHitIndices(ista).clear();
+ wData.TsHitIndices(ista).clear();
}
}
}
- fvStatNhitsProcessed[iThread] += statNwindowHits;
+ statNhitsProcessed += statNwindowHits;
// print the LOG for every 10 ms of data processed
if constexpr (0) {
- int currentChunk = (int) ((fvWindowStartThread[iThread] - fStatTsStart) / 10.e6);
+ int currentChunk = (int) ((windowStart - fStatTsStart) / 10.e6);
if (!areUntouchedDataLeft || currentChunk > statLastLogTimeChunk) {
statLastLogTimeChunk = currentChunk;
- double dataRead =
- 100. * (fvWindowStartThread[iThread] + fWindowLength - fStatTsStart) / (fStatTsEnd - fStatTsStart);
+ double dataRead = 100. * (windowStart + fWindowLength - fStatTsStart) / (fStatTsEnd - fStatTsStart);
if (dataRead > 100.) {
dataRead = 100.;
}
- LOG(info) << "CA tracker process sliding window N " << fvStatNwindows[iThread] << ": time "
- << fvWindowStartThread[iThread] / 1.e6 << " ms + " << fWindowLength / 1.e3 << " us) with "
- << statNwindowHits << " hits. "
+ LOG(info) << "CA tracker process sliding window N " << statNwindows << ": time " << windowStart / 1.e6
+ << " ms + " << fWindowLength / 1.e3 << " us) with " << statNwindowHits << " hits. "
<< " Processing " << dataRead << " % of the TS time and "
- << 100. * fvStatNhitsProcessed[iThread] / fStatNhitsTotal << " % of TS hits."
+ << 100. * statNhitsProcessed / fStatNhitsTotal << " % of TS hits."
<< " Already reconstructed " << fvRecoTracks[iThread].size() << " tracks on thread #" << iThread;
}
}
@@ -513,107 +499,68 @@ void TrackFinder::FindTracksThread(int iThread)
//out << statNwindowHits << ' ';
- frAlgo.fvMonitorDataThread[iThread].StopTimer(ETimer::PrepareWindow);
+ fvMonitorDataThread[iThread].StopTimer(ETimer::PrepareWindow);
//Timer trackingInWindow; //DBG
//trackingInWindow.Start();
- frAlgo.fvMonitorDataThread[iThread].StartTimer(ETimer::TrackingWindow);
- frAlgo.fvTrackFinderWindow[iThread].CaTrackFinderSlice();
- frAlgo.fvMonitorDataThread[iThread].StopTimer(ETimer::TrackingWindow);
+ fvMonitorDataThread[iThread].StartTimer(ETimer::TrackingWindow);
+ trackFinderWindow.CaTrackFinderSlice(input, wData);
+ fvMonitorDataThread[iThread].StopTimer(ETimer::TrackingWindow);
//trackingInWindow.Stop();
//out << trackingInWindow.GetTotalMs() << ' ';
// save reconstructed tracks with no hits in the overlap region
- //if (fvWindowStartThread[iThread] > 13.23e6 && fvWindowStartThread[iThread] < 13.26e6) {
- fvWindowStartThread[iThread] += fWindowLength - fWindowOverlap;
+ //if (windowStart > 13.23e6 && windowStart < 13.26e6) {
+ windowStart += fWindowLength - fWindowOverlap;
// we should add hits from reconstructed but not stored tracks to the new sub-timeslice
// we do it in a simple way by extending the tsStartNew
// TODO: only add those hits from the region before tsStartNew that belong to the not stored tracks
- //out << frAlgo.fvWData[iThread].RecoHitIndices().size() << ' ';
- //out << frAlgo.fvWData[iThread].RecoTracks().size() << '\n';
-
- frAlgo.fvMonitorDataThread[iThread].StartTimer(ETimer::StoreTracksWindow);
- int trackFirstHit = 0;
- for (const auto& track : frAlgo.fvWData[iThread].RecoTracks()) {
- bool isTrackCompletelyInOverlap = true;
- for (int ih = 0; ih < track.fNofHits; ih++) {
- int caHitId = frAlgo.fvWData[iThread].RecoHitIndex(trackFirstHit + ih);
- CaHitTimeInfo& info = frAlgo.fHitTimeInfo[caHitId];
- if (info.fEventTimeMax < fvWindowStartThread[iThread]) { // this hit is before the overlap
- isTrackCompletelyInOverlap = false;
- break;
- }
- }
+ //out << fvWData[iThread].RecoHitIndices().size() << ' ';
+ //out << fvWData[iThread].RecoTracks().size() << '\n';
- if (!areUntouchedDataLeft) { // don't reject tracks in the overlap when no more data are left
- isTrackCompletelyInOverlap = 0;
- }
+ fvMonitorDataThread[iThread].StartTimer(ETimer::StoreTracksWindow);
+ auto trackFirstHit = wData.RecoHitIndices().begin();
+
+ for (const auto& track : wData.RecoTracks()) {
+
+ const bool isTrackCompletelyInOverlap =
+ std::all_of(trackFirstHit, trackFirstHit + track.fNofHits, [&](int caHitId) {
+ CaHitTimeInfo& info = fHitTimeInfo[caHitId];
+ return info.fEventTimeMax >= windowStart;
+ });
- if (isTrackCompletelyInOverlap) {
- //
- // Don't save tracks from the overlap region, since they might have additional hits in the next subslice
- //
-
- // release the track hits
- for (int i = 0; i < track.fNofHits; i++) {
- int caHitId = frAlgo.fvWData[iThread].RecoHitIndex(trackFirstHit + i);
- const auto& h = frAlgo.fInputData.GetHit(caHitId);
- wData.IsHitKeyUsed(h.FrontKey()) = 0;
- wData.IsHitKeyUsed(h.BackKey()) = 0;
+ // Don't save tracks from the overlap region, since they might have additional hits in the next subslice.
+ // Don't reject tracks in the overlap when no more data are left
+ const bool useFlag = !isTrackCompletelyInOverlap || !areUntouchedDataLeft;
+
+ for (int i = 0; i < track.fNofHits; i++) {
+ const int caHitId = *(trackFirstHit + i);
+ const auto& h = input.GetHit(caHitId);
+ wData.IsHitKeyUsed(h.FrontKey()) = static_cast<int>(useFlag);
+ wData.IsHitKeyUsed(h.BackKey()) = static_cast<int>(useFlag);
+
+ if (useFlag) {
+ fvRecoHitIndices[iThread].push_back(caHitId);
}
}
- else { // save the track
+ if (useFlag) {
fvRecoTracks[iThread].push_back(track);
- // mark the track hits as used
- for (int i = 0; i < track.fNofHits; i++) {
- int caHitId = frAlgo.fvWData[iThread].RecoHitIndex(trackFirstHit + i);
- const auto& h = frAlgo.fInputData.GetHit(caHitId);
- wData.IsHitKeyUsed(h.FrontKey()) = 1;
- wData.IsHitKeyUsed(h.BackKey()) = 1;
- fvRecoHitIndices[iThread].push_back(caHitId);
- }
}
+
trackFirstHit += track.fNofHits;
} // sub-timeslice tracks
- frAlgo.fvMonitorDataThread[iThread].StopTimer(ETimer::StoreTracksWindow);
+ fvMonitorDataThread[iThread].StopTimer(ETimer::StoreTracksWindow);
- if (fvWindowStartThread[iThread] > fvWindowEndThread[iThread]) {
+ if (windowStart > windowEnd) {
break;
}
else {
- fvWindowStartThread[iThread] -= fWindowMargin; // do 5 ns margin
+ windowStart -= fWindowMargin; // do 5 ns margin
}
}
- frAlgo.fvMonitorDataThread[iThread].StopTimer(ETimer::TrackingThread);
+ fvMonitorDataThread[iThread].StopTimer(ETimer::TrackingThread);
//timer.Stop();
//LOG(info) << "CA: finishing tracking on thread " << iThread << " (time: " << timer.GetTotalMs() << " ms, "
- // << "hits processed: " << fvStatNhitsProcessed[iThread] << ", "
+ // << "hits processed: " << statNhitsProcessed << ", "
// << "hits used: " << fvRecoHitIndices[iThread].size() << ')';
}
-
-// ---------------------------------------------------------------------------------------------------------------------
-//
-void TrackFinder::Init()
-{
- assert(frAlgo.GetNofThreads() > 0);
-
- fvWindowStartThread.clear();
- fvWindowStartThread.resize(frAlgo.GetNofThreads());
- fvWindowEndThread.clear();
- fvWindowEndThread.resize(frAlgo.GetNofThreads());
- fvStatNwindows.clear();
- fvStatNwindows.resize(frAlgo.GetNofThreads());
- fvStatNhitsProcessed.clear();
- fvStatNhitsProcessed.resize(frAlgo.GetNofThreads());
- fvRecoTracks.clear();
- fvRecoTracks.resize(frAlgo.GetNofThreads());
- fvRecoHitIndices.clear();
- fvRecoHitIndices.resize(frAlgo.GetNofThreads());
-
- for (int iThread = 0; iThread < frAlgo.GetNofThreads(); ++iThread) {
- fvRecoTracks[iThread].SetName(std::string("TrackFinder::fvRecoTracks_") + std::to_string(iThread));
- fvRecoHitIndices[iThread].SetName(std::string("TrackFinder::fvRecoHitIndices_") + std::to_string(iThread));
- }
-
- fWindowLength = (ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) ? 500 : 10000;
-}
diff --git a/algo/ca/core/tracking/CaTrackFinder.h b/algo/ca/core/tracking/CaTrackFinder.h
index ffe5a856a2..7f270d633b 100644
--- a/algo/ca/core/tracking/CaTrackFinder.h
+++ b/algo/ca/core/tracking/CaTrackFinder.h
@@ -9,6 +9,10 @@
#include "CaBranch.h"
#include "CaHit.h"
#include "CaSimd.h"
+#include "CaTimesliceHeader.h"
+#include "CaTrackExtender.h"
+#include "CaTrackFinderWindow.h"
+#include "CaTrackFitter.h"
#include "CaTrackParam.h"
#include "CaVector.h"
@@ -29,11 +33,13 @@ namespace cbm::algo::ca
///
class TrackFinder {
public:
- /// Default constructor
- TrackFinder(ca::Framework& algo);
+ typedef std::pair<Vector<Track>, Vector<ca::HitIndex_t>> Output_t;
+ /// Default constructora
+ TrackFinder(const ca::Parameters<fvec>& pars, const float mass, const ca::TrackingMode& mode,
+ TrackingMonitorData& monitorData, int nThreads, double& recoTime);
/// Destructor
- ~TrackFinder();
+ ~TrackFinder() = default;
/// Copy constructor
TrackFinder(const TrackFinder&) = delete;
@@ -47,31 +53,36 @@ namespace cbm::algo::ca
/// Move assignment operator
TrackFinder& operator=(TrackFinder&&) = delete;
- void FindTracks();
- void FindTracksThread(int iThread);
- void Init();
+ Output_t FindTracks(const InputData& input, TimesliceHeader& tsHeader);
const auto& GetRecoTracksContainer(int iThread) const { return fvRecoTracks[iThread]; }
const auto& GetRecoHitIndicesContainer(int iThread) const { return fvRecoHitIndices[iThread]; }
+ TrackingMode GetTrackingMode() const { return fTrackingMode; }
+ const std::vector<ca::WindowData>& GetWData() const { return fvWData; }
private:
// -------------------------------
// Private methods
+ void FindTracksThread(const InputData& input, int iThread, float& windowStart, float& windowEnd, int& statNwindows,
+ int& statNhitsProcessed);
+ // bool checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2) const;
- bool checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2) const;
-
-
- private:
// -------------------------------
// Data members
- ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class
+ Vector<CaHitTimeInfo> fHitTimeInfo;
+
+ const Parameters<fvec>& fParameters; ///< Object of Framework parameters class
+ float fDefaultMass{constants::phys::MuonMass}; ///< mass of the propagated particle [GeV/c2]
+ ca::TrackingMode fTrackingMode;
+
+ TrackingMonitorData& fMonitorData; ///< Tracking monitor data (statistics per call)
+ std::vector<TrackingMonitorData> fvMonitorDataThread; ///< Tracking monitor data per thread
- std::vector<float> fvWindowStartThread;
- std::vector<float> fvWindowEndThread;
+ std::vector<ca::WindowData> fvWData; ///< Intrnal data processed in a time-window
- std::vector<int> fvStatNwindows;
- std::vector<int> fvStatNhitsProcessed;
+ int fNofThreads; ///< Number of threads to execute the track-finder
+ double& fCaRecoTime; // time of the track finder + fitter
std::vector<Vector<Track>> fvRecoTracks; ///< reconstructed tracks
std::vector<Vector<HitIndex_t>> fvRecoHitIndices; ///< packed hits of reconstructed tracks
diff --git a/algo/ca/core/tracking/CaTrackFinderWindow.cxx b/algo/ca/core/tracking/CaTrackFinderWindow.cxx
index 3ad43cadb0..ff9b0a5b8e 100644
--- a/algo/ca/core/tracking/CaTrackFinderWindow.cxx
+++ b/algo/ca/core/tracking/CaTrackFinderWindow.cxx
@@ -46,13 +46,15 @@ using namespace cbm::algo::ca;
// ---------------------------------------------------------------------------------------------------------------------
//
-TrackFinderWindow::TrackFinderWindow(ca::Framework& algo, ca::WindowData& wData, ca::TrackingMonitorData& monitorData)
- : frAlgo(algo)
- , frWData(wData)
+TrackFinderWindow::TrackFinderWindow(const ca::Parameters<fvec>& pars, const float mass, const ca::TrackingMode& mode,
+ ca::TrackingMonitorData& monitorData)
+ : fParameters(pars)
+ , fDefaultMass(mass)
+ , fTrackingMode(mode)
, frMonitorData(monitorData)
- , fTrackExtender(frAlgo, frWData)
- , fCloneMerger(frAlgo)
- , fTrackFitter(frAlgo, frWData)
+ , fTrackExtender(pars, mass)
+ , fCloneMerger(pars, mass)
+ , fTrackFitter(pars, mass, mode)
{
for (unsigned int i = 0; i < constants::size::MaxNstations; ++i) {
fvTriplets[i].SetName(std::stringstream() << "TrackFinderWindow::fTriplets[" << i << "]");
@@ -68,7 +70,8 @@ TrackFinderWindow::~TrackFinderWindow() {}
//
-bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2) const
+bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2,
+ WindowData& wData) const
{
dchi2 = 1.e20;
@@ -79,7 +82,7 @@ bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triple
if (r.GetLSta() != l.GetMSta()) return false;
- const float tripletLinkChi2 = frWData.CurrentIteration()->GetTripletLinkChi2();
+ const float tripletLinkChi2 = wData.CurrentIteration()->GetTripletLinkChi2();
if (r.IsMomentumFitted()) {
assert(l.IsMomentumFitted());
@@ -120,7 +123,6 @@ bool TrackFinderWindow::checkTripletMatch(const ca::Triplet& l, const ca::Triple
return true;
}
-void TrackFinderWindow::InitTimeslice() { fvHitKeyToTrack.reset(frAlgo.fInputData.GetNhitKeys(), -1); }
// **************************************************************************************************
// * *
@@ -130,65 +132,62 @@ void TrackFinderWindow::InitTimeslice() { fvHitKeyToTrack.reset(frAlgo.fInputDat
// ---------------------------------------------------------------------------------------------------------------------
//
-void TrackFinderWindow::ReadWindowData()
+void TrackFinderWindow::ReadWindowData(const ca::InputData& input, WindowData& wData)
{
int nHits = 0;
- for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); iS++) {
- frWData.HitStartIndexOnStation(iS) = nHits;
- frWData.NofHitsOnStation(iS) = frWData.TsHitIndices(iS).size();
- nHits += frWData.NofHitsOnStation(iS);
+ for (int iS = 0; iS < fParameters.GetNstationsActive(); iS++) {
+ wData.HitStartIndexOnStation(iS) = nHits;
+ wData.NofHitsOnStation(iS) = wData.TsHitIndices(iS).size();
+ nHits += wData.NofHitsOnStation(iS);
}
- frWData.HitStartIndexOnStation(frAlgo.GetParameters().GetNstationsActive()) = nHits;
- frWData.ResetHitData(nHits);
-
- for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); iS++) {
- int iFstHit = frWData.HitStartIndexOnStation(iS);
- for (ca::HitIndex_t ih = 0; ih < frWData.TsHitIndices(iS).size(); ++ih) {
- ca::Hit h = frAlgo.fInputData.GetHit(frWData.TsHitIndex(iS, ih));
- h.SetId(frWData.TsHitIndex(iS, ih));
- frWData.Hit(iFstHit + ih) = h;
+ wData.HitStartIndexOnStation(fParameters.GetNstationsActive()) = nHits;
+ wData.ResetHitData(nHits);
+
+ for (int iS = 0; iS < fParameters.GetNstationsActive(); iS++) {
+ int iFstHit = wData.HitStartIndexOnStation(iS);
+ for (ca::HitIndex_t ih = 0; ih < wData.TsHitIndices(iS).size(); ++ih) {
+ ca::Hit h = input.GetHit(wData.TsHitIndex(iS, ih));
+ h.SetId(wData.TsHitIndex(iS, ih));
+ wData.Hit(iFstHit + ih) = h;
}
}
if constexpr (fDebug) {
LOG(info) << "===== Sliding Window hits: ";
for (int i = 0; i < nHits; ++i) {
- LOG(info) << " " << frWData.Hit(i).ToString();
+ LOG(info) << " " << wData.Hit(i).ToString();
}
LOG(info) << "===== ";
}
- fvHitFirstTriplet.reset(nHits, 0);
- fvHitNofTriplets.reset(nHits, 0);
+ wData.RecoTracks().clear();
+ wData.RecoTracks().reserve(2 * nHits / fParameters.GetNstationsActive());
- frWData.RecoTracks().clear();
- frWData.RecoTracks().reserve(2 * nHits / frAlgo.GetParameters().GetNstationsActive());
-
- frWData.RecoHitIndices().clear();
- frWData.RecoHitIndices().reserve(2 * nHits);
+ wData.RecoHitIndices().clear();
+ wData.RecoHitIndices().reserve(2 * nHits);
fvTrackCandidates.clear();
fvTrackCandidates.reserve(nHits / 10);
for (unsigned int iS = 0; iS < constants::size::MaxNstations; iS++) {
- int nHitsStation = frWData.TsHitIndices(iS).size();
+ int nHitsStation = wData.TsHitIndices(iS).size();
fvTriplets[iS].clear();
fvTriplets[iS].reserve(2 * nHitsStation);
}
}
-void TrackFinderWindow::CaTrackFinderSlice()
+void TrackFinderWindow::CaTrackFinderSlice(const ca::InputData& input, WindowData& wData)
{
/********************************/ /**
* CATrackFinder routine setting
***********************************/
frMonitorData.StartTimer(ETimer::InitWindow);
- ReadWindowData();
+ ReadWindowData(input, wData);
frMonitorData.StopTimer(ETimer::InitWindow);
frMonitorData.StartTimer(ETimer::PrepareGrid);
- for (int iS = 0; iS < frAlgo.GetParameters().GetNstationsActive(); ++iS) {
+ for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
bool timeUninitialised = 1;
fscal lasttime = 0;
@@ -200,8 +199,8 @@ void TrackFinderWindow::CaTrackFinderSlice()
fscal gridMinY = -0.1;
fscal gridMaxY = 0.1;
- for (ca::HitIndex_t ih = 0; ih < frWData.TsHitIndices(iS).size(); ++ih) {
- const ca::Hit& h = frAlgo.fInputData.GetHit(frWData.TsHitIndex(iS, ih));
+ for (ca::HitIndex_t ih = 0; ih < wData.TsHitIndices(iS).size(); ++ih) {
+ const ca::Hit& h = input.GetHit(wData.TsHitIndex(iS, ih));
if (h.X() < gridMinX) {
gridMinX = h.X();
@@ -230,7 +229,7 @@ void TrackFinderWindow::CaTrackFinderSlice()
// TODO: changing the grid also changes the result. Investigate why it happens.
// TODO: find the optimal grid size
- int nSliceHits = frWData.TsHitIndices(iS).size();
+ int nSliceHits = wData.TsHitIndices(iS).size();
fscal sizeY = gridMaxY - gridMinY;
fscal sizeX = gridMaxX - gridMinX;
int nBins2D = 1 + nSliceHits;
@@ -238,20 +237,20 @@ void TrackFinderWindow::CaTrackFinderSlice()
fscal yStep = 0.3 * sizeY / sqrt(nBins2D);
fscal xStep = 0.8 * sizeX / sqrt(nBins2D);
- fscal scale = frAlgo.GetParameters().GetStation(iS).GetZ<float>() - frAlgo.GetParameters().GetTargetPositionZ()[0];
+ fscal scale = fParameters.GetStation(iS).GetZ<float>() - fParameters.GetTargetPositionZ()[0];
if (yStep > 0.3 * scale) yStep = 0.3 * scale;
if (yStep < 0.01 * scale) yStep = 0.01 * scale;
if (xStep > 0.3 * scale) xStep = 0.3 * scale;
if (xStep < 0.01 * scale) xStep = 0.01 * scale;
- auto& grid = frWData.Grid(iS);
+ auto& grid = wData.Grid(iS);
grid.BuildBins(gridMinX, gridMaxX, gridMinY, gridMaxY, xStep, yStep);
/* clang-format off */
- grid.StoreHits(frWData.Hits(),
- frWData.HitStartIndexOnStation(iS),
- frWData.NofHitsOnStation(iS),
- frWData.HitKeyFlags());
+ grid.StoreHits(wData.Hits(),
+ wData.HitStartIndexOnStation(iS),
+ wData.NofHitsOnStation(iS),
+ wData.HitKeyFlags());
/* clang-format on */
}
frMonitorData.StopTimer(ETimer::PrepareGrid);
@@ -260,30 +259,30 @@ void TrackFinderWindow::CaTrackFinderSlice()
//
frMonitorData.StartTimer(ETimer::FindTracks);
// TODO: Get rid of fCurrentTierationIndex. If it is needed of logs, please used caIteration.GetName()
- for (int iIteration = 0; iIteration < frAlgo.GetParameters().GetNcaIterations(); ++iIteration) {
+ for (int iIteration = 0; iIteration < fParameters.GetNcaIterations(); ++iIteration) {
// ----- Prepare iteration
//
frMonitorData.StartTimer(ETimer::PrepareIteration);
- const auto& caIteration = frAlgo.GetParameters().GetCAIterations()[iIteration];
- frWData.SetCurrentIteration(&caIteration);
- frWData.SetCurrentIterationIndex(iIteration);
+ const auto& caIteration = fParameters.GetCAIterations()[iIteration];
+ wData.SetCurrentIteration(&caIteration);
+ wData.SetCurrentIterationIndex(iIteration);
if (iIteration > 0) {
- for (int ista = 0; ista < frAlgo.GetParameters().GetNstationsActive(); ++ista) {
- frWData.Grid(ista).RemoveUsedHits(frWData.Hits(), frWData.HitKeyFlags());
+ for (int ista = 0; ista < fParameters.GetNstationsActive(); ++ista) {
+ wData.Grid(ista).RemoveUsedHits(wData.Hits(), wData.HitKeyFlags());
}
}
- int nHits = frWData.Hits().size();
+ int nHits = wData.Hits().size();
// --> frAlgo.fIsWindowHitSuppressed.reset(frAlgo.fWindowHits.size(), 0);
- frWData.ResetHitSuppressionFlags(); // TODO: ??? No effect?
- for (int j = 0; j < frAlgo.GetParameters().GetNstationsActive(); j++) {
+ wData.ResetHitSuppressionFlags(); // TODO: ??? No effect?
+ for (int j = 0; j < fParameters.GetNstationsActive(); j++) {
fvTriplets[j].clear();
}
- fvHitFirstTriplet.reset(nHits, 0);
- fvHitNofTriplets.reset(nHits, 0);
+ Vector<int> vHitFirstTriplet(nHits, 0); /// link hit -> first triplet { hit, *, *}
+ Vector<int> vHitNofTriplets(nHits, 0); /// link hit ->n triplets { hit, *, *}
{
// #pragma omp task
@@ -295,21 +294,19 @@ void TrackFinderWindow::CaTrackFinderSlice()
// Select magnetic field. For primary tracks - fVtxFieldValue, for secondary tracks - st.fieldSlice
if (caIteration.GetPrimaryFlag()) {
- frWData.TargB() = frAlgo.GetParameters().GetVertexFieldValue();
+ wData.TargB() = fParameters.GetVertexFieldValue();
}
else {
- frWData.TargB() = frAlgo.GetParameters().GetStation(0).fieldSlice.GetFieldValue(0, 0);
+ wData.TargB() = fParameters.GetStation(0).fieldSlice.GetFieldValue(0, 0);
} // NOTE: calculates field frAlgo.fTargB in the center of 0th station
- frAlgo.fIsTargetField = (fabs(frWData.TargB().y[0]) > 0.001);
-
- frWData.TargetMeasurement().SetX(frAlgo.GetParameters().GetTargetPositionX());
- frWData.TargetMeasurement().SetY(frAlgo.GetParameters().GetTargetPositionY());
- frWData.TargetMeasurement().SetDx2(SigmaTargetX * SigmaTargetX);
- frWData.TargetMeasurement().SetDxy(0);
- frWData.TargetMeasurement().SetDy2(SigmaTargetY * SigmaTargetY);
- frWData.TargetMeasurement().SetNdfX(1);
- frWData.TargetMeasurement().SetNdfY(1);
+ wData.TargetMeasurement().SetX(fParameters.GetTargetPositionX());
+ wData.TargetMeasurement().SetY(fParameters.GetTargetPositionY());
+ wData.TargetMeasurement().SetDx2(SigmaTargetX * SigmaTargetX);
+ wData.TargetMeasurement().SetDxy(0);
+ wData.TargetMeasurement().SetDy2(SigmaTargetY * SigmaTargetY);
+ wData.TargetMeasurement().SetNdfX(1);
+ wData.TargetMeasurement().SetNdfY(1);
}
}
frMonitorData.StopTimer(ETimer::PrepareIteration);
@@ -318,20 +315,20 @@ void TrackFinderWindow::CaTrackFinderSlice()
//
frMonitorData.StartTimer(ETimer::ConstructTriplets);
- ca::TripletConstructor constructor(frAlgo, frWData);
+ ca::TripletConstructor constructor(fParameters, input, wData, fDefaultMass, fTrackingMode);
// indices of the two neighbouring station, taking into account allowed gaps
std::vector<std::pair<int, int>> staPattern;
- for (int gap = 0; gap <= frWData.CurrentIteration()->GetMaxStationGap(); gap++) {
+ for (int gap = 0; gap <= wData.CurrentIteration()->GetMaxStationGap(); gap++) {
for (int i = 0; i <= gap; i++) {
staPattern.push_back(std::make_pair(1 + i, 2 + gap));
}
}
- const int iStFirst = frWData.CurrentIteration()->GetFirstStationIndex();
- for (int istal = frAlgo.GetParameters().GetNstationsActive() - 2; istal >= iStFirst; istal--) {
+ const int iStFirst = wData.CurrentIteration()->GetFirstStationIndex();
+ for (int istal = fParameters.GetNstationsActive() - 2; istal >= iStFirst; istal--) {
// start with downstream chambers
- const auto& grid = frWData.Grid(istal);
+ const auto& grid = wData.Grid(istal);
Tindex nGridEntriesL = grid.GetEntries().size();
for (Tindex iel = 0; iel < nGridEntriesL; ++iel) {
ca::HitIndex_t ihitl = grid.GetEntries()[iel].GetObjectId();
@@ -341,8 +338,8 @@ void TrackFinderWindow::CaTrackFinderSlice()
fvTriplets[istal].insert(fvTriplets[istal].end(), constructor.GetTriplets().begin(),
constructor.GetTriplets().end());
}
- fvHitFirstTriplet[ihitl] = frAlgo.PackTripletId(istal, oldSize);
- fvHitNofTriplets[ihitl] = fvTriplets[istal].size() - oldSize;
+ vHitFirstTriplet[ihitl] = PackTripletId(istal, oldSize);
+ vHitNofTriplets[ihitl] = fvTriplets[istal].size() - oldSize;
}
} // istal
frMonitorData.StopTimer(ETimer::ConstructTriplets);
@@ -351,7 +348,7 @@ void TrackFinderWindow::CaTrackFinderSlice()
//
frMonitorData.StartTimer(ETimer::SearchNeighbours);
- for (int istal = frAlgo.GetParameters().GetNstationsActive() - 2; istal >= iStFirst; istal--) {
+ for (int istal = fParameters.GetNstationsActive() - 2; istal >= iStFirst; istal--) {
// start with downstream chambers
for (unsigned int it = 0; it < fvTriplets[istal].size(); ++it) {
@@ -360,15 +357,14 @@ void TrackFinderWindow::CaTrackFinderSlice()
tr.SetFNeighbour(0);
tr.SetNNeighbours(0);
- unsigned int nNeighbours = fvHitNofTriplets[tr.GetMHit()];
-
- unsigned int neighLocation = fvHitFirstTriplet[tr.GetMHit()];
- unsigned int neighStation = frAlgo.TripletId2Station(neighLocation);
- unsigned int neighTriplet = frAlgo.TripletId2Triplet(neighLocation);
+ unsigned int nNeighbours = vHitNofTriplets[tr.GetMHit()];
+ unsigned int neighLocation = vHitFirstTriplet[tr.GetMHit()];
+ unsigned int neighStation = TripletId2Station(neighLocation);
+ unsigned int neighTriplet = TripletId2Triplet(neighLocation);
if (nNeighbours > 0) {
assert((int) neighStation >= istal + 1
- && (int) neighStation <= istal + 1 + frWData.CurrentIteration()->GetMaxStationGap());
+ && (int) neighStation <= istal + 1 + wData.CurrentIteration()->GetMaxStationGap());
}
unsigned char level = 0;
@@ -378,7 +374,7 @@ void TrackFinderWindow::CaTrackFinderSlice()
ca::Triplet& neighbour = fvTriplets[neighStation][neighTriplet];
fscal dchi2 = 0.;
- if (!checkTripletMatch(tr, neighbour, dchi2)) continue;
+ if (!checkTripletMatch(tr, neighbour, dchi2, wData)) continue;
if (tr.GetFNeighbour() == 0) {
tr.SetFNeighbour(neighLocation);
@@ -407,7 +403,7 @@ void TrackFinderWindow::CaTrackFinderSlice()
// min level to start triplet. So min track length = min_level+3.
int min_level = std::min(caIteration.GetMinNhits(), caIteration.GetMinNhitsStation0()) - 3;
- if (frWData.CurrentIteration()->GetTrackFromTripletsFlag()) {
+ if (wData.CurrentIteration()->GetTrackFromTripletsFlag()) {
min_level = 0;
}
@@ -419,48 +415,47 @@ void TrackFinderWindow::CaTrackFinderSlice()
// collect consequtive: the longest tracks, shorter, more shorter and so on
- for (int firstTripletLevel = frAlgo.GetParameters().GetNstationsActive() - 3; firstTripletLevel >= min_level;
+ for (int firstTripletLevel = fParameters.GetNstationsActive() - 3; firstTripletLevel >= min_level;
firstTripletLevel--) {
// choose length in triplets number - firstTripletLevel - the maximum possible triplet level among all triplets in the searched candidate
// how many levels to check
- int nlevel = (frAlgo.GetParameters().GetNstationsActive() - 2) - firstTripletLevel + 1;
+ int nlevel = (fParameters.GetNstationsActive() - 2) - firstTripletLevel + 1;
const unsigned char min_best_l =
(firstTripletLevel > min_level) ? firstTripletLevel + 2 : min_level + 3; // loose maximum
fvTrackCandidates.clear();
- fvTrackCandidates.reserve(frWData.Hits().size() / 10);
+ fvTrackCandidates.reserve(wData.Hits().size() / 10);
- for (const auto& h : frWData.Hits()) {
+ for (const auto& h : wData.Hits()) {
fvHitKeyToTrack[h.FrontKey()] = -1;
fvHitKeyToTrack[h.BackKey()] = -1;
}
//== Loop over triplets with the required level, find and store track candidates
- for (int istaF = iStFirst; istaF <= frAlgo.GetParameters().GetNstationsActive() - 3 - firstTripletLevel;
- ++istaF) {
+ for (int istaF = iStFirst; istaF <= fParameters.GetNstationsActive() - 3 - firstTripletLevel; ++istaF) {
if (--nlevel == 0) break; //TODO: SG: this is not needed
for (Tindex itrip = 0; itrip < (Tindex) fvTriplets[istaF].size(); ++itrip) {
ca::Triplet& first_trip = (fvTriplets[istaF][itrip]);
- const auto& fstTripLHit = frWData.Hit(first_trip.GetLHit());
- if (frWData.IsHitKeyUsed(fstTripLHit.FrontKey()) || frWData.IsHitKeyUsed(fstTripLHit.BackKey())) {
+ const auto& fstTripLHit = wData.Hit(first_trip.GetLHit());
+ if (wData.IsHitKeyUsed(fstTripLHit.FrontKey()) || wData.IsHitKeyUsed(fstTripLHit.BackKey())) {
continue;
}
// skip track candidates that are too short
- int minNhits = frWData.CurrentIteration()->GetMinNhits();
+ int minNhits = wData.CurrentIteration()->GetMinNhits();
if (fstTripLHit.Station() == 0) {
- minNhits = frWData.CurrentIteration()->GetMinNhitsStation0();
+ minNhits = wData.CurrentIteration()->GetMinNhitsStation0();
}
- if (frWData.CurrentIteration()->GetTrackFromTripletsFlag()) {
+ if (wData.CurrentIteration()->GetTrackFromTripletsFlag()) {
minNhits = 0;
}
@@ -482,8 +477,8 @@ void TrackFinderWindow::CaTrackFinderSlice()
best_tr = curr_tr;
- CAFindTrack(istaF, best_tr, &first_trip, curr_tr, min_best_l,
- new_tr); /// reqursive func to build a tree of possible track-candidates and choose the best
+ /// reqursive func to build a tree of possible track-candidates and choose the best
+ CAFindTrack(istaF, best_tr, &first_trip, curr_tr, min_best_l, new_tr, wData);
if (best_tr.NofHits() < firstTripletLevel + 2) continue; // loose maximum one hit
@@ -497,17 +492,16 @@ void TrackFinderWindow::CaTrackFinderSlice()
// TODO: automatize the NDF calculation
- if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode
- || ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) {
+ if (ca::TrackingMode::kGlobal == fTrackingMode || ca::TrackingMode::kMcbm == fTrackingMode) {
ndf = best_tr.NofHits() * 2 - 4;
}
best_tr.SetChi2(best_tr.Chi2() / ndf);
- if (frAlgo.GetParameters().GetGhostSuppression()) {
+ if (fParameters.GetGhostSuppression()) {
if (3 == best_tr.NofHits()) {
- if (!frWData.CurrentIteration()->GetPrimaryFlag() && (istaF != 0))
+ if (!wData.CurrentIteration()->GetPrimaryFlag() && (istaF != 0))
continue; // too /*short*/ non-MAPS track
- if (frWData.CurrentIteration()->GetPrimaryFlag() && (best_tr.Chi2() > 5.0)) continue;
+ if (wData.CurrentIteration()->GetPrimaryFlag() && (best_tr.Chi2() > 5.0)) continue;
}
}
fvTrackCandidates.push_back(best_tr);
@@ -517,10 +511,10 @@ void TrackFinderWindow::CaTrackFinderSlice()
tr.SetAlive(true);
if constexpr (fDebug) {
std::stringstream s;
- s << "iter " << frWData.CurrentIterationIndex() << ", track candidate " << fvTrackCandidates.size() - 1
+ s << "iter " << wData.CurrentIterationIndex() << ", track candidate " << fvTrackCandidates.size() - 1
<< " found, L = " << best_tr.NofHits() << " chi2= " << best_tr.Chi2() << " hits: ";
for (auto hitId : tr.Hits()) {
- s << hitId << " (mc " << frAlgo.GetMcTrackIdForCaHit(hitId) << ") ";
+ s << hitId << " (mc " << ca::Framework::GetMcTrackIdForCaHit(hitId) << ") ";
}
LOG(info) << s.str();
}
@@ -545,7 +539,7 @@ void TrackFinderWindow::CaTrackFinderSlice()
continue;
}
for (auto& hitId : tr.Hits()) {
- const ca::Hit& h = frAlgo.fInputData.GetHit(hitId);
+ const ca::Hit& h = input.GetHit(hitId);
bool isAlive = true;
{ // front strip
auto& stripF = fvHitKeyToTrack[h.FrontKey()];
@@ -596,14 +590,14 @@ void TrackFinderWindow::CaTrackFinderSlice()
tr.SetAlive(true);
for (int iHit = 0; tr.IsAlive() && (iHit < (int) tr.Hits().size()); ++iHit) {
- const ca::Hit& h = frAlgo.fInputData.GetHit(tr.Hits()[iHit]);
+ const ca::Hit& h = input.GetHit(tr.Hits()[iHit]);
tr.SetAlive(tr.IsAlive() && (fvHitKeyToTrack[h.FrontKey()] == tr.Id())
&& (fvHitKeyToTrack[h.BackKey()] == tr.Id()));
}
if (!tr.IsAlive()) { // release strips
for (auto hitId : tr.Hits()) {
- const ca::Hit& h = frAlgo.fInputData.GetHit(hitId);
+ const ca::Hit& h = input.GetHit(hitId);
if (fvHitKeyToTrack[h.FrontKey()] == tr.Id()) {
fvHitKeyToTrack[h.FrontKey()] = -1;
}
@@ -627,34 +621,34 @@ void TrackFinderWindow::CaTrackFinderSlice()
ca::Branch& tr = fvTrackCandidates[iCandidate];
if constexpr (fDebug) {
- LOG(info) << "iter " << frWData.CurrentIterationIndex() << ", track candidate " << iCandidate
+ LOG(info) << "iter " << wData.CurrentIterationIndex() << ", track candidate " << iCandidate
<< ": alive = " << tr.IsAlive();
}
if (!tr.IsAlive()) continue;
- if (frWData.CurrentIteration()->GetExtendTracksFlag()) {
- if (tr.NofHits() < frAlgo.GetParameters().GetNstationsActive()) {
- fTrackExtender.ExtendBranch(tr);
+ if (wData.CurrentIteration()->GetExtendTracksFlag()) {
+ if (tr.NofHits() < fParameters.GetNstationsActive()) {
+ fTrackExtender.ExtendBranch(tr, input, wData);
}
}
for (auto iHit : tr.Hits()) {
- const ca::Hit& hit = frAlgo.fInputData.GetHit(iHit);
+ const ca::Hit& hit = input.GetHit(iHit);
/// used strips are marked
- frWData.IsHitKeyUsed(hit.FrontKey()) = 1;
- frWData.IsHitKeyUsed(hit.BackKey()) = 1;
- frWData.RecoHitIndices().push_back(iHit);
+ wData.IsHitKeyUsed(hit.FrontKey()) = 1;
+ wData.IsHitKeyUsed(hit.BackKey()) = 1;
+ wData.RecoHitIndices().push_back(iHit);
}
Track t;
t.fNofHits = tr.NofHits();
- frWData.RecoTracks().push_back(t);
+ wData.RecoTracks().push_back(t);
if (0) { // SG debug
std::stringstream s;
s << "store track " << iCandidate << " chi2= " << tr.Chi2() << "\n";
s << " hits: ";
for (unsigned int ih = 0; ih < tr.Hits().size(); ih++) {
- s << frAlgo.GetMcTrackIdForCaHit(tr.Hits()[ih]) << " ";
+ s << ca::Framework::GetMcTrackIdForCaHit(tr.Hits()[ih]) << " ";
}
LOG(info) << s.str();
}
@@ -666,11 +660,11 @@ void TrackFinderWindow::CaTrackFinderSlice()
frMonitorData.StartTimer(ETimer::SuppressHitKeys);
// suppress strips of suppressed hits
- for (unsigned int ih = 0; ih < frWData.Hits().size(); ih++) {
- if (frWData.IsHitSuppressed(ih)) {
- const ca::Hit& hit = frWData.Hit(ih);
- frWData.IsHitKeyUsed(hit.FrontKey()) = 1;
- frWData.IsHitKeyUsed(hit.BackKey()) = 1;
+ for (unsigned int ih = 0; ih < wData.Hits().size(); ih++) {
+ if (wData.IsHitSuppressed(ih)) {
+ const ca::Hit& hit = wData.Hit(ih);
+ wData.IsHitKeyUsed(hit.FrontKey()) = 1;
+ wData.IsHitKeyUsed(hit.BackKey()) = 1;
}
}
frMonitorData.StopTimer(ETimer::SuppressHitKeys);
@@ -679,16 +673,16 @@ void TrackFinderWindow::CaTrackFinderSlice()
// Fit tracks
frMonitorData.StartTimer(ETimer::FitTracks);
- fTrackFitter.FitCaTracks();
+ fTrackFitter.FitCaTracks(input, wData);
frMonitorData.StopTimer(ETimer::FitTracks);
// Merge clones
frMonitorData.StartTimer(ETimer::MergeClones);
- fCloneMerger.Exec(frWData.RecoTracks(), frWData.RecoHitIndices());
+ fCloneMerger.Exec(input, wData);
frMonitorData.StopTimer(ETimer::MergeClones);
frMonitorData.StartTimer(ETimer::FitTracks);
- fTrackFitter.FitCaTracks();
+ fTrackFitter.FitCaTracks(input, wData);
frMonitorData.StopTimer(ETimer::FitTracks);
}
@@ -703,7 +697,7 @@ void TrackFinderWindow::CaTrackFinderSlice()
****************************************************************/
void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Triplet* curr_trip, ca::Branch& curr_tr,
- unsigned char min_best_l, ca::Branch* new_tr)
+ unsigned char min_best_l, ca::Branch* new_tr, WindowData& wData)
/// recursive search for tracks
/// input: @ista - station index, @&best_tr - best track for the privious call
/// output: @&NCalls - number of function calls
@@ -713,14 +707,14 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
{
// -- finish with current track
// add rest of hits
- const auto& hitM = frWData.Hit(curr_trip->GetMHit());
- const auto& hitR = frWData.Hit(curr_trip->GetRHit());
+ const auto& hitM = wData.Hit(curr_trip->GetMHit());
+ const auto& hitR = wData.Hit(curr_trip->GetRHit());
- if (!(frWData.IsHitKeyUsed(hitM.FrontKey()) || frWData.IsHitKeyUsed(hitM.BackKey()))) {
+ if (!(wData.IsHitKeyUsed(hitM.FrontKey()) || wData.IsHitKeyUsed(hitM.BackKey()))) {
curr_tr.AddHit(hitM.Id());
}
- if (!(frWData.IsHitKeyUsed(hitR.FrontKey()) || frWData.IsHitKeyUsed(hitR.BackKey()))) {
+ if (!(wData.IsHitKeyUsed(hitR.FrontKey()) || wData.IsHitKeyUsed(hitR.BackKey()))) {
curr_tr.AddHit(hitR.Id());
}
@@ -728,13 +722,12 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
// TODO: SZh 21.08.2023: Replace hard-coded value with a parameter
int ndf = curr_tr.NofHits() * 2 - 5;
- if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode
- || ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) {
+ if (ca::TrackingMode::kGlobal == fTrackingMode || ca::TrackingMode::kMcbm == fTrackingMode) {
ndf = curr_tr.NofHits() * 2 - 4;
}
- if (curr_tr.Chi2() > frWData.CurrentIteration()->GetTrackChi2Cut() * ndf) {
+ if (curr_tr.Chi2() > wData.CurrentIteration()->GetTrackChi2Cut() * ndf) {
return;
}
@@ -754,16 +747,16 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
unsigned int ID = curr_trip->GetFNeighbour() + in;
- unsigned int Station = frAlgo.TripletId2Station(ID);
- unsigned int Triplet = frAlgo.TripletId2Triplet(ID);
+ unsigned int Station = TripletId2Station(ID);
+ unsigned int Triplet = TripletId2Triplet(ID);
const ca::Triplet& new_trip = fvTriplets[Station][Triplet];
fscal dchi2 = 0.;
- if (!checkTripletMatch(*curr_trip, new_trip, dchi2)) continue;
+ if (!checkTripletMatch(*curr_trip, new_trip, dchi2, wData)) continue;
- const auto& hitL = frWData.Hit(new_trip.GetLHit()); // left hit of new triplet
- if (frWData.IsHitKeyUsed(hitL.FrontKey()) || frWData.IsHitKeyUsed(hitL.BackKey())) { //hits are used
+ const auto& hitL = wData.Hit(new_trip.GetLHit()); // left hit of new triplet
+ if (wData.IsHitKeyUsed(hitL.FrontKey()) || wData.IsHitKeyUsed(hitL.BackKey())) { //hits are used
// no used hits allowed -> compare and store track
if ((curr_tr.NofHits() > best_tr.NofHits())
|| ((curr_tr.NofHits() == best_tr.NofHits()) && (curr_tr.Chi2() < best_tr.Chi2()))) {
@@ -777,12 +770,13 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
if constexpr (0) { //SGtrd2d debug!!
- int mc01 = frAlgo.GetMcTrackIdForWindowHit(curr_trip->GetLHit());
- int mc02 = frAlgo.GetMcTrackIdForWindowHit(curr_trip->GetMHit());
- int mc03 = frAlgo.GetMcTrackIdForWindowHit(curr_trip->GetRHit());
- int mc11 = frAlgo.GetMcTrackIdForWindowHit(new_trip.GetLHit());
- int mc12 = frAlgo.GetMcTrackIdForWindowHit(new_trip.GetMHit());
- int mc13 = frAlgo.GetMcTrackIdForWindowHit(new_trip.GetRHit());
+ int mc01 = ca::Framework::GetMcTrackIdForWindowHit(curr_trip->GetLHit());
+ int mc02 = ca::Framework::GetMcTrackIdForWindowHit(curr_trip->GetMHit());
+ int mc03 = ca::Framework::GetMcTrackIdForWindowHit(curr_trip->GetRHit());
+ int mc11 = ca::Framework::GetMcTrackIdForWindowHit(new_trip.GetLHit());
+ int mc12 = ca::Framework::GetMcTrackIdForWindowHit(new_trip.GetMHit());
+ int mc13 = ca::Framework::GetMcTrackIdForWindowHit(new_trip.GetRHit());
+
if ((mc01 == mc02) && (mc02 == mc03)) {
LOG(info) << " sta " << ista << " mc0 " << mc01 << " " << mc02 << " " << mc03 << " mc1 " << mc11 << " "
<< mc12 << " " << mc13 << " chi2 " << curr_tr.Chi2() / (2 * (curr_tr.NofHits() + 2) - 4)
@@ -794,17 +788,17 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
int ndf = 2 * (new_L + 2) - 5;
- if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode) { //SGtrd2d!!!
+ if (ca::TrackingMode::kGlobal == fTrackingMode) { //SGtrd2d!!!
ndf = 2 * (new_L + 2) - 4;
}
- else if (ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) {
+ else if (ca::TrackingMode::kMcbm == fTrackingMode) {
ndf = 2 * (new_L + 2) - 4;
}
else {
ndf = new_L; // TODO: SG: 2 * (new_L + 2) - 5
}
- if (new_chi2 > frWData.CurrentIteration()->GetTrackChi2Cut() * ndf) continue;
+ if (new_chi2 > wData.CurrentIteration()->GetTrackChi2Cut() * ndf) continue;
// add new hit
new_tr[ista] = curr_tr;
@@ -813,7 +807,7 @@ void TrackFinderWindow::CAFindTrack(int ista, ca::Branch& best_tr, const ca::Tri
const int new_ista = ista + new_trip.GetMSta() - new_trip.GetLSta();
new_tr[ista].SetChi2(new_chi2);
- CAFindTrack(new_ista, best_tr, &new_trip, new_tr[ista], min_best_l, new_tr);
+ CAFindTrack(new_ista, best_tr, &new_trip, new_tr[ista], min_best_l, new_tr, wData);
} // add triplet to track
} // for neighbours
} // level = 0
diff --git a/algo/ca/core/tracking/CaTrackFinderWindow.h b/algo/ca/core/tracking/CaTrackFinderWindow.h
index d26cb0000c..5e6aa450a0 100644
--- a/algo/ca/core/tracking/CaTrackFinderWindow.h
+++ b/algo/ca/core/tracking/CaTrackFinderWindow.h
@@ -38,8 +38,8 @@ namespace cbm::algo::ca
class TrackFinderWindow {
public:
/// Default constructor
- TrackFinderWindow(ca::Framework& algo, ca::WindowData& wData, ca::TrackingMonitorData& monitorData);
-
+ TrackFinderWindow(const ca::Parameters<fvec>& pars, const float mass, const ca::TrackingMode& mode,
+ ca::TrackingMonitorData& monitorData);
/// Destructor
~TrackFinderWindow();
@@ -55,20 +55,37 @@ namespace cbm::algo::ca
/// Move assignment operator
TrackFinderWindow& operator=(TrackFinderWindow&&) = delete;
- void ReadWindowData();
- void CaTrackFinderSlice();
-
- void CAFindTrack(int ista, ca::Branch& best_tr, const ca::Triplet* curr_trip, ca::Branch& curr_tr,
- unsigned char min_best_l, ca::Branch* new_tr);
+ void CaTrackFinderSlice(const ca::InputData& input, WindowData& wData);
/// \note The function initializes global arrays for a given thread
- void InitTimeslice();
+ void InitTimeslice(size_t nHitKeys) { fvHitKeyToTrack.reset(nHitKeys, -1); }
private:
///-------------------------------
/// Private methods
- bool checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2) const;
+ void ReadWindowData(const ca::InputData& input, WindowData& wData);
+
+ void CAFindTrack(int ista, ca::Branch& best_tr, const ca::Triplet* curr_trip, ca::Branch& curr_tr,
+ unsigned char min_best_l, ca::Branch* new_tr, WindowData& wData);
+
+ bool checkTripletMatch(const ca::Triplet& l, const ca::Triplet& r, fscal& dchi2, WindowData& wData) const;
+
+ // ** Functions, which pack and unpack indexes of station and triplet **
+
+ // TODO: move to ca::Triplet class (S.Zharko)
+ /// Packs station and triplet indices to an unique triplet ID
+ /// \param iStation Index of station in the active stations array
+ /// \param iTriplet Index of triplet
+ static unsigned int PackTripletId(unsigned int iStation, unsigned int iTriplet);
+
+ /// Unpacks the triplet ID to its station index
+ /// \param id Unique triplet ID
+ static unsigned int TripletId2Station(unsigned int id);
+
+ /// Unpacks the triplet ID to its triplet index
+ /// \param id Unique triplet ID
+ static unsigned int TripletId2Triplet(unsigned int id);
private:
@@ -82,20 +99,53 @@ namespace cbm::algo::ca
/// \note The candidates may share any amount of hits.
Vector<ca::Branch> fvTrackCandidates{"TrackFinderWindow::fTrackCandidates"};
- Vector<int> fvHitFirstTriplet{"TrackFinderWindow::fvHitFirstTriplet"}; /// link hit -> first triplet { hit, *, *}
- Vector<int> fvHitNofTriplets{"TrackFinderWindow::fvHitNofTriplets"}; /// link hit ->n triplets { hit, *, *}
-
/// \note Global array for a given thread
Vector<int> fvHitKeyToTrack{"TrackFinderWindow::fvHitKeyToTrack"};
- ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class
static constexpr bool fDebug = false; // print debug info
- WindowData& frWData;
+ const Parameters<fvec>& fParameters; ///< Object of Framework parameters class
+ float fDefaultMass{constants::phys::MuonMass}; ///< mass of the propagated particle [GeV/c2]
+ ca::TrackingMode fTrackingMode;
+
TrackingMonitorData& frMonitorData; ///< Reference to monitor data
- TrackExtender fTrackExtender; ///< Object of the track extender algorithm
- CloneMerger fCloneMerger; ///< Object of the clone merger algorithm
- TrackFitter fTrackFitter; ///< Object of the track extender algorithm
+ TrackExtender fTrackExtender; ///< Object of the track extender algorithm
+ CloneMerger fCloneMerger; ///< Object of the clone merger algorithm
+ TrackFitter fTrackFitter; ///< Object of the track extender algorithm
};
+ // ********************************************
+ // ** Inline member functions implementation **
+ // ********************************************
+
+ // ---------------------------------------------------------------------------------------------------------------------
+ //
+ [[gnu::always_inline]] inline unsigned int TrackFinderWindow::PackTripletId(unsigned int iStation,
+ unsigned int iTriplet)
+ {
+#ifndef FAST_CODE
+ assert(iStation < constants::size::MaxNstations);
+ assert(iTriplet < constants::size::MaxNtriplets);
+#endif
+ constexpr unsigned int kMoveStation = constants::size::TripletBits;
+ return (iStation << kMoveStation) + iTriplet;
+ }
+
+ // ---------------------------------------------------------------------------------------------------------------------
+ //
+ [[gnu::always_inline]] inline unsigned int TrackFinderWindow::TripletId2Station(unsigned int id)
+ {
+ constexpr unsigned int kMoveStation = constants::size::TripletBits;
+ return id >> kMoveStation;
+ }
+
+ // ---------------------------------------------------------------------------------------------------------------------
+ //
+ [[gnu::always_inline]] inline unsigned int TrackFinderWindow::TripletId2Triplet(unsigned int id)
+ {
+ constexpr unsigned int kTripletMask = (1u << constants::size::TripletBits) - 1u;
+ return id & kTripletMask;
+ }
+
+
} // namespace cbm::algo::ca
diff --git a/algo/ca/core/tracking/CaTrackFitter.cxx b/algo/ca/core/tracking/CaTrackFitter.cxx
index 26a0f22a15..00a599fb12 100644
--- a/algo/ca/core/tracking/CaTrackFitter.cxx
+++ b/algo/ca/core/tracking/CaTrackFitter.cxx
@@ -20,7 +20,12 @@ using namespace cbm::algo;
// ---------------------------------------------------------------------------------------------------------------------
//
-TrackFitter::TrackFitter(ca::Framework& algo, WindowData& wData) : frAlgo(algo), frWData(wData) {}
+TrackFitter::TrackFitter(const ca::Parameters<fvec>& pars, const float mass, const ca::TrackingMode& mode)
+ : fParameters(pars)
+ , fDefaultMass(mass)
+ , fTrackingMode(mode)
+{
+}
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -28,10 +33,10 @@ TrackFitter::~TrackFitter() {}
// ---------------------------------------------------------------------------------------------------------------------
//
-void TrackFitter::FitCaTracks()
+void TrackFitter::FitCaTracks(const ca::InputData& input, WindowData& wData)
{
// LOG(info) << " Start CA Track Fitter ";
- int start_hit = 0; // for interation in frWData.RecoHitIndices()
+ int start_hit = 0; // for interation in wData.RecoHitIndices()
// static ca::FieldValue fldB0, fldB1, fldB2 _fvecalignment;
// static ca::FieldRegion fld _fvecalignment;
@@ -42,17 +47,17 @@ void TrackFitter::FitCaTracks()
ca::FieldValue<fvec> fldB01, fldB11, fldB21 _fvecalignment;
ca::FieldRegion<fvec> fld1 _fvecalignment;
- const int nStations = frAlgo.GetParameters().GetNstationsActive();
+ const int nStations = fParameters.GetNstationsActive();
int nTracks_SIMD = fvec::size();
ca::TrackFit fit; // fit parameters coresponding to the current track
TrackParamV& tr = fit.Tr();
- fit.SetParticleMass(frAlgo.GetDefaultParticleMass());
+ fit.SetParticleMass(fDefaultMass);
fit.SetDoFitVelocity(true);
Track* t[fvec::size()]{nullptr};
- const ca::Station<fvec>* sta = frAlgo.GetParameters().GetStations().begin();
+ const ca::Station<fvec>* sta = fParameters.GetStations().begin();
// Spatial-time position of a hit vs. station and track in the portion
@@ -102,19 +107,19 @@ void TrackFitter::FitCaTracks()
mxy[ista].SetCov(1., 0., 1.);
}
- unsigned short N_vTracks = frWData.RecoTracks().size(); // number of tracks processed per one SSE register
+ unsigned short N_vTracks = wData.RecoTracks().size(); // number of tracks processed per one SSE register
for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvec::size()) {
if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) nTracks_SIMD = N_vTracks - itrack;
for (int i = 0; i < nTracks_SIMD; i++) {
- t[i] = &frWData.RecoTrack(itrack + i);
+ t[i] = &wData.RecoTrack(itrack + i);
}
// fill the rest of the SIMD vectors with something reasonable
for (uint i = nTracks_SIMD; i < fvec::size(); i++) {
- t[i] = &frWData.RecoTrack(itrack);
+ t[i] = &wData.RecoTrack(itrack);
}
// get hits of current track
@@ -133,7 +138,7 @@ void TrackFitter::FitCaTracks()
for (int ih = 0; ih < nHitsTrack; ih++) {
- const ca::Hit& hit = frAlgo.GetInputData().GetHit(frWData.RecoHitIndex(start_hit++));
+ const ca::Hit& hit = input.GetHit(wData.RecoHitIndex(start_hit++));
const int ista = hit.Station();
//if (sta[ista].fieldStatus) { isFieldPresent[iVec] = true; }
@@ -197,16 +202,15 @@ void TrackFitter::FitCaTracks()
}
for (int ih = nHitsTrack - 1; ih >= 0; ih--) {
- const int ista = iSta[ih];
- const ca::Station<fvec>& st = frAlgo.GetParameters().GetStation(ista);
- By[ista][iVec] = st.fieldSlice.cy[0][0];
+ const int ista = iSta[ih];
+ const ca::Station<fvec>& st = fParameters.GetStation(ista);
+ By[ista][iVec] = st.fieldSlice.cy[0][0];
}
}
fit.GuessTrack(z_end, x, y, z, time, By, w, w_time, nStations);
- if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode
- || ca::Framework::TrackingMode::kMcbm == frAlgo.fTrackingMode) {
+ if (ca::TrackingMode::kGlobal == fTrackingMode || ca::TrackingMode::kMcbm == fTrackingMode) {
tr.Qp() = fvec(1. / 1.1);
}
@@ -259,8 +263,8 @@ void TrackFitter::FitCaTracks()
fit.SetMask(initialised);
fit.Extrapolate(z[ista], fld1);
- fit.MultipleScattering(frAlgo.GetParameters().GetMaterialThickness(ista, tr.X(), tr.Y()));
- fit.EnergyLossCorrection(frAlgo.GetParameters().GetMaterialThickness(ista, tr.X(), tr.Y()), fvec(1.f));
+ fit.MultipleScattering(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()));
+ fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()), fvec(1.f));
fit.SetMask(initialised && w[ista]);
fit.FilterXY(mxy[ista]);
@@ -279,12 +283,12 @@ void TrackFitter::FitCaTracks()
{
fitpv.SetMask(fmask::One());
- ca::MeasurementXy<fvec> vtxInfo = frWData.TargetMeasurement();
+ ca::MeasurementXy<fvec> vtxInfo = wData.TargetMeasurement();
vtxInfo.SetDx2(1.e-8);
vtxInfo.SetDxy(0.);
vtxInfo.SetDy2(1.e-8);
- if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode) {
+ if (ca::TrackingMode::kGlobal == fTrackingMode) {
ca::FieldRegion<fvec> fldFull _fvecalignment;
fldFull.SetUseOriginalField(true);
fitpv.SetMaxExtrapolationStep(1.);
@@ -292,13 +296,13 @@ void TrackFitter::FitCaTracks()
fitpv.SetQp0(fitpv.Tr().Qp());
fitpv.Tr() = fit.Tr();
fitpv.Tr().Qp() = fitpv.Qp0();
- fitpv.Extrapolate(frAlgo.GetParameters().GetTargetPositionZ(), fldFull);
+ fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fldFull);
fitpv.FilterXY(vtxInfo);
}
}
else {
fitpv.SetQp0(fitpv.Tr().Qp());
- fitpv.Extrapolate(frAlgo.GetParameters().GetTargetPositionZ(), fld);
+ fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fld);
}
}
@@ -308,7 +312,7 @@ void TrackFitter::FitCaTracks()
//fit.FilterVi(TrackParamV::kClightNsInv);
TrackParamV Tf = fit.Tr();
- if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode) {
+ if (ca::TrackingMode::kGlobal == fTrackingMode) {
Tf.Qp() = fitpv.Tr().Qp();
}
@@ -363,8 +367,8 @@ void TrackFitter::FitCaTracks()
fit.SetMask(initialised);
fit.Extrapolate(z[ista], fld);
- fit.MultipleScattering(frAlgo.GetParameters().GetMaterialThickness(ista, tr.X(), tr.Y()));
- fit.EnergyLossCorrection(frAlgo.GetParameters().GetMaterialThickness(ista, tr.X(), tr.Y()), fvec(-1.f));
+ fit.MultipleScattering(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()));
+ fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.X(), tr.Y()), fvec(-1.f));
fit.SetMask(initialised && w[ista]);
fit.FilterXY(mxy[ista]);
fit.FilterTime(time[ista], dt2[ista], sta[ista].timeInfo);
@@ -381,7 +385,7 @@ void TrackFitter::FitCaTracks()
//fit.FilterVi(TrackParamV::kClightNsInv);
TrackParamV Tl = fit.Tr();
- if (ca::Framework::TrackingMode::kGlobal == frAlgo.fTrackingMode) {
+ if (ca::TrackingMode::kGlobal == fTrackingMode) {
Tl.Qp() = fitpv.Tr().Qp();
}
diff --git a/algo/ca/core/tracking/CaTrackFitter.h b/algo/ca/core/tracking/CaTrackFitter.h
index 2b451763b3..d452716908 100644
--- a/algo/ca/core/tracking/CaTrackFitter.h
+++ b/algo/ca/core/tracking/CaTrackFitter.h
@@ -10,6 +10,8 @@
#include "CaBranch.h"
#include "CaHit.h"
+#include "CaInputData.h"
+#include "CaParameters.h"
#include "CaSimd.h"
#include "CaTrackParam.h"
#include "CaVector.h"
@@ -19,20 +21,18 @@
namespace cbm::algo::ca
{
class Track;
- class Framework;
namespace
{
using namespace cbm::algo; // to keep ca:: prefices in the code
}
-
/// Class implements a track fit the CA track finder
///
class TrackFitter {
public:
/// Default constructor
- TrackFitter(ca::Framework& algo, WindowData& wData);
+ TrackFitter(const ca::Parameters<fvec>& pars, const float mass, const ca::TrackingMode& mode);
/// Destructor
~TrackFitter();
@@ -49,20 +49,17 @@ namespace cbm::algo::ca
/// Move assignment operator
TrackFitter& operator=(TrackFitter&&) = delete;
-
/// Fit tracks, found by the CA tracker
- void FitCaTracks();
+ void FitCaTracks(const ca::InputData& input, WindowData& wData);
- private:
- ///-------------------------------
- /// Private methods
private:
///-------------------------------
/// Data members
- ca::Framework& frAlgo; ///< Reference to the main track finder algorithm class
- ca::WindowData& frWData;
+ const Parameters<fvec>& fParameters; ///< Object of Framework parameters class
+ float fDefaultMass{constants::phys::MuonMass}; ///< mass of the propagated particle [GeV/c2]
+ ca::TrackingMode fTrackingMode;
};
} // namespace cbm::algo::ca
diff --git a/algo/ca/core/tracking/CaTripletConstructor.cxx b/algo/ca/core/tracking/CaTripletConstructor.cxx
index 9476f0f8d6..e10dea41f4 100644
--- a/algo/ca/core/tracking/CaTripletConstructor.cxx
+++ b/algo/ca/core/tracking/CaTripletConstructor.cxx
@@ -19,7 +19,26 @@
using namespace cbm::algo::ca;
-TripletConstructor::TripletConstructor(ca::Framework& algo, WindowData& wData) : frAlgo(algo), frWData(wData) {}
+TripletConstructor::TripletConstructor(const ca::Parameters<fvec>& pars, const ca::InputData& input, WindowData& wData,
+ const float mass, const ca::TrackingMode& mode)
+ : fParameters(pars)
+ , fInputData(input)
+ , frWData(wData)
+ , fDefaultMass(mass)
+ , fTrackingMode(mode)
+{
+ // FIXME: SZh 24.08.2022
+ // This approach is suitable only for a case, when all the stations inside a magnetic field come right before
+ // all the stations outside of the field!
+ fNfieldStations = std::lower_bound(fParameters.GetStations().cbegin(),
+ fParameters.GetStations().cbegin() + fParameters.GetNstationsActive(),
+ 0, // we are looking for the first zero element
+ [](const ca::Station<fvec>& s, int edge) { return bool(s.fieldStatus) > edge; })
+ - fParameters.GetStations().cbegin();
+
+ fIsTargetField = (fabs(frWData.TargB().y[0]) > 0.001);
+}
+
void TripletConstructor::InitStations(int istal, int istam, int istar)
{
@@ -29,29 +48,29 @@ void TripletConstructor::InitStations(int istal, int istam, int istar)
fIstaM = istam;
fIstaR = istar;
- if (fIstaM >= frAlgo.fParameters.GetNstationsActive()) {
+ if (fIstaM >= fParameters.GetNstationsActive()) {
return;
}
- fStaL = &frAlgo.fParameters.GetStation(fIstaL);
- fStaM = &frAlgo.fParameters.GetStation(fIstaM);
- fStaR = &frAlgo.fParameters.GetStation(fIstaR);
+ fStaL = &fParameters.GetStation(fIstaL);
+ fStaM = &fParameters.GetStation(fIstaM);
+ fStaR = &fParameters.GetStation(fIstaR);
{ // two stations for approximating the field between the target and the left hit
int sta1 = fIstaL;
- if (sta1 >= frAlgo.fNfieldStations) {
- sta1 = frAlgo.fNfieldStations - 1;
+ if (sta1 >= fNfieldStations) {
+ sta1 = fNfieldStations - 1;
}
if (sta1 < 1) {
sta1 = 1;
}
int sta0 = sta1 / 2; // station between fIstaL and the target
- assert(0 <= sta0 && sta0 < sta1 && sta1 < frAlgo.fParameters.GetNstationsActive());
+ assert(0 <= sta0 && sta0 < sta1 && sta1 < fParameters.GetNstationsActive());
- fFld0Sta[0] = &frAlgo.fParameters.GetStation(sta0);
- fFld0Sta[1] = &frAlgo.fParameters.GetStation(sta1);
+ fFld0Sta[0] = &fParameters.GetStation(sta0);
+ fFld0Sta[1] = &fParameters.GetStation(sta1);
}
{ // three stations for approximating the field between the left and the right hit
@@ -59,7 +78,7 @@ void TripletConstructor::InitStations(int istal, int istam, int istar)
int sta0 = fIstaL;
int sta1 = fIstaM;
int sta2 = fIstaM + 1;
- if (sta2 >= frAlgo.fParameters.GetNstationsActive()) {
+ if (sta2 >= fParameters.GetNstationsActive()) {
sta2 = fIstaM;
sta1 = fIstaM - 1;
if (sta0 >= sta1) {
@@ -69,16 +88,16 @@ void TripletConstructor::InitStations(int istal, int istam, int istar)
sta0 = 0;
}
}
- if (frAlgo.fParameters.GetNstationsActive() >= 3) {
- assert(0 <= sta0 && sta0 < sta1 && sta1 < sta2 && sta2 < frAlgo.fParameters.GetNstationsActive());
+ if (fParameters.GetNstationsActive() >= 3) {
+ assert(0 <= sta0 && sta0 < sta1 && sta1 < sta2 && sta2 < fParameters.GetNstationsActive());
}
- fFld1Sta[0] = &frAlgo.fParameters.GetStation(sta0);
- fFld1Sta[1] = &frAlgo.fParameters.GetStation(sta1);
- fFld1Sta[2] = &frAlgo.fParameters.GetStation(sta2);
+ fFld1Sta[0] = &fParameters.GetStation(sta0);
+ fFld1Sta[1] = &fParameters.GetStation(sta1);
+ fFld1Sta[2] = &fParameters.GetStation(sta2);
}
- fFit.SetParticleMass(frAlgo.GetDefaultParticleMass());
+ fFit.SetParticleMass(fDefaultMass);
if (frWData.CurrentIteration()->GetElectronFlag()) {
fFit.SetParticleMass(constants::phys::ElectronMass);
}
@@ -95,7 +114,7 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c
if (!fIsInitialized) return;
- fIhitL = iHitL;
+ fIhitL = iHitL;
const auto& hitl = frWData.Hit(fIhitL);
// fit a straight line through the target and the left hit.
@@ -117,10 +136,10 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c
fvec zl = hitl.Z();
fvec time = hitl.T();
fvec timeEr2 = hitl.dT2();
- const fvec dzli = 1.f / (zl - frAlgo.GetParameters().GetTargetPositionZ());
+ const fvec dzli = 1.f / (zl - fParameters.GetTargetPositionZ());
- const fvec tx = (xl - frAlgo.GetParameters().GetTargetPositionX()) * dzli;
- const fvec ty = (yl - frAlgo.GetParameters().GetTargetPositionY()) * dzli;
+ const fvec tx = (xl - fParameters.GetTargetPositionX()) * dzli;
+ const fvec ty = (yl - fParameters.GetTargetPositionY()) * dzli;
TrackParamV& T = fFit.Tr();
@@ -157,7 +176,7 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c
{
ca::FieldValue<fvec> B0 = fFld0Sta[0]->fieldSlice.GetFieldValueForLine(T);
ca::FieldValue<fvec> B1 = fFld0Sta[1]->fieldSlice.GetFieldValueForLine(T);
- fld0.Set(frWData.TargB(), frAlgo.GetParameters().GetTargetPositionZ(), B0, fFld0Sta[0]->fZ, B1, fFld0Sta[1]->fZ);
+ fld0.Set(frWData.TargB(), fParameters.GetTargetPositionZ(), B0, fFld0Sta[0]->fZ, B1, fFld0Sta[1]->fZ);
}
{ // field, made by the left hit, the middle station and the right station
@@ -170,9 +189,9 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c
// add the target constraint
- fFit.FilterWithTargetAtLine(frAlgo.GetParameters().GetTargetPositionZ(), frWData.TargetMeasurement(), fld0);
+ fFit.FilterWithTargetAtLine(fParameters.GetTargetPositionZ(), frWData.TargetMeasurement(), fld0);
- fFit.MultipleScattering(frAlgo.fParameters.GetMaterialThickness(fIstaL, T.GetX(), T.GetY()));
+ fFit.MultipleScattering(fParameters.GetMaterialThickness(fIstaL, T.GetX(), T.GetY()));
// extrapolate to the middle hit
@@ -190,17 +209,16 @@ void TripletConstructor::CreateTripletsForHit(int istal, int istam, int istar, c
void TripletConstructor::FindDoublets()
{
/// Find the doublets. Reformat data into portions of doublets.
-
int iMC = -1;
- if (frAlgo.fParameters.DevIsMatchDoubletsViaMc()) {
- iMC = frAlgo.GetMcTrackIdForWindowHit(fIhitL);
+ if (fParameters.DevIsMatchDoubletsViaMc()) {
+ iMC = ca::Framework::GetMcTrackIdForWindowHit(fIhitL);
if (iMC < 0) {
return;
}
}
CollectHits(fTrL, fIstaM, frWData.CurrentIteration()->GetDoubletChi2Cut(), iMC, fHitsM_2,
- frAlgo.fParameters.GetMaxDoubletsPerSinglet());
+ fParameters.GetMaxDoubletsPerSinglet());
}
@@ -215,7 +233,7 @@ void TripletConstructor::FitDoublets()
fTracks_2.clear();
fTracks_2.reserve(hitsMtmp.size());
- bool isMomentumFitted = (frAlgo.fIsTargetField || (staL().fieldStatus != 0) || (staM().fieldStatus != 0));
+ bool isMomentumFitted = (fIsTargetField || (staL().fieldStatus != 0) || (staM().fieldStatus != 0));
double maxQp = frWData.CurrentIteration()->GetMaxQp();
for (unsigned int i2 = 0; i2 < hitsMtmp.size(); i2++) {
@@ -246,12 +264,12 @@ void TripletConstructor::FitDoublets()
fFit.SetQp0(isMomentumFitted ? fFit.Tr().GetQp() : maxQp);
- fFit.MultipleScattering(frAlgo.fParameters.GetMaterialThickness(fIstaM, T2.GetX(), T2.Y()));
+ fFit.MultipleScattering(fParameters.GetMaterialThickness(fIstaM, T2.GetX(), T2.Y()));
fFit.SetQp0(fFit.Tr().Qp());
// check if there are other hits close to the doublet on the same station
- if (ca::Framework::kMcbm != frAlgo.fTrackingMode) {
+ if (ca::kMcbm != fTrackingMode) {
// TODO: SG: adjust cuts, put them to parameters
fscal tx = T2.Tx()[0];
@@ -282,12 +300,14 @@ void TripletConstructor::FitDoublets()
continue;
}
- if (frAlgo.fParameters.DevIsSuppressOverlapHitsViaMc()) {
- int iMC = frAlgo.GetMcTrackIdForWindowHit(fIhitL);
- if ((iMC != frAlgo.GetMcTrackIdForWindowHit(indM)) || (iMC != frAlgo.GetMcTrackIdForWindowHit(indClone))) {
+ if (fParameters.DevIsSuppressOverlapHitsViaMc()) {
+ int iMC = ca::Framework::GetMcTrackIdForWindowHit(fIhitL);
+ if ((iMC != ca::Framework::GetMcTrackIdForWindowHit(indM))
+ || (iMC != ca::Framework::GetMcTrackIdForWindowHit(indClone))) {
continue;
}
}
+
frWData.SuppressHit(indClone);
}
}
@@ -300,7 +320,7 @@ void TripletConstructor::FitDoublets()
void TripletConstructor::FindTriplets()
{
- if (fIstaR >= frAlgo.fParameters.GetNstationsActive()) {
+ if (fIstaR >= fParameters.GetNstationsActive()) {
return;
}
FindRightHit();
@@ -317,10 +337,10 @@ void TripletConstructor::FindRightHit()
fFit.SetQp0(fvec(0.));
- if (fIstaR >= frAlgo.fParameters.GetNstationsActive()) return;
+ if (fIstaR >= fParameters.GetNstationsActive()) return;
{
- int maxTriplets = fHitsM_2.size() * frAlgo.fParameters.GetMaxTripletPerDoublets();
+ int maxTriplets = fHitsM_2.size() * fParameters.GetMaxTripletPerDoublets();
fHitsM_3.clear();
fHitsR_3.clear();
fHitsM_3.reserve(maxTriplets);
@@ -345,8 +365,8 @@ void TripletConstructor::FindRightHit()
int ih0 = ihit[0];
int ih1 = ihit[1];
- const ca::Hit& h0 = frAlgo.fInputData.GetHit(ih0);
- const ca::Hit& h1 = frAlgo.fInputData.GetHit(ih1);
+ const ca::Hit& h0 = fInputData.GetHit(ih0);
+ const ca::Hit& h1 = fInputData.GetHit(ih1);
LOG(info) << "\n====== Extrapolated Doublet : "
<< " iter " << frWData.CurrentIteration()->GetName() << " hits: {" << fIstaL << "/" << ih0 << " "
@@ -360,7 +380,7 @@ void TripletConstructor::FindRightHit()
bool isDoubletGood = false;
int iMC = -1;
do {
- if (frAlgo.kSts == frAlgo.fTrackingMode && (T2.C44()[0] < 0)) {
+ if (ca::TrackingMode::kSts == fTrackingMode && (T2.C44()[0] < 0)) {
break;
}
if (T2.C00()[0] < 0 || T2.C11()[0] < 0 || T2.C22()[0] < 0 || T2.C33()[0] < 0 || T2.C55()[0] < 0) {
@@ -372,10 +392,10 @@ void TripletConstructor::FindRightHit()
if (fabs(T2.Ty()[0]) > maxSlope) {
break;
}
- if (frAlgo.fParameters.DevIsMatchTripletsViaMc()) {
+ if (fParameters.DevIsMatchTripletsViaMc()) {
int indM = fHitsM_2[i2];
- iMC = frAlgo.GetMcTrackIdForWindowHit(fIhitL);
- if (iMC < 0 || iMC != frAlgo.GetMcTrackIdForWindowHit(indM)) {
+ iMC = ca::Framework::GetMcTrackIdForWindowHit(fIhitL);
+ if (iMC < 0 || iMC != ca::Framework::GetMcTrackIdForWindowHit(indM)) {
break;
}
}
@@ -397,10 +417,10 @@ void TripletConstructor::FindRightHit()
}
Vector<ca::HitIndex_t> collectedHits;
- CollectHits(T2, fIstaR, tripletChi2Cut, iMC, collectedHits, frAlgo.fParameters.GetMaxTripletPerDoublets());
+ CollectHits(T2, fIstaR, tripletChi2Cut, iMC, collectedHits, fParameters.GetMaxTripletPerDoublets());
- if (collectedHits.size() >= frAlgo.fParameters.GetMaxTripletPerDoublets()) {
- LOG(debug) << "Ca: GetMaxTripletPerDoublets==" << frAlgo.fParameters.GetMaxTripletPerDoublets()
+ if (collectedHits.size() >= fParameters.GetMaxTripletPerDoublets()) {
+ LOG(debug) << "Ca: GetMaxTripletPerDoublets==" << fParameters.GetMaxTripletPerDoublets()
<< " reached in findTripletsStep0()";
// reject already created triplets for this doublet
collectedHits.clear();
@@ -412,7 +432,7 @@ void TripletConstructor::FindRightHit()
ca::HitIndex_t irh = collectedHits[ih];
if constexpr (fDebugDublets) {
ca::HitIndex_t ihit = frWData.Hit(irh).Id();
- const ca::Hit& h = frAlgo.fInputData.GetHit(ihit);
+ const ca::Hit& h = fInputData.GetHit(ihit);
LOG(info) << " hit " << ihit << " " << h.ToString();
}
if (frWData.IsHitSuppressed(irh)) {
@@ -455,7 +475,7 @@ void TripletConstructor::FitTriplets()
fit.SetParticleMass(constants::phys::ElectronMass);
}
else {
- fit.SetParticleMass(frAlgo.GetDefaultParticleMass());
+ fit.SetParticleMass(fDefaultMass);
}
const int NHits = 3; // triplets
@@ -464,7 +484,7 @@ void TripletConstructor::FitTriplets()
ca::Station<fvec> sta[3];
for (int is = 0; is < NHits; ++is) {
- sta[is] = frAlgo.fParameters.GetStation(ista[is]);
+ sta[is] = fParameters.GetStation(ista[is]);
};
bool isMomentumFitted = ((staL().fieldStatus != 0) || (staM().fieldStatus != 0) || (staR().fieldStatus != 0));
@@ -480,10 +500,10 @@ void TripletConstructor::FitTriplets()
ca::HitIndex_t ihit[NHits] = {frWData.Hit(iwhit[0]).Id(), frWData.Hit(iwhit[1]).Id(), frWData.Hit(iwhit[2]).Id()};
- if (frAlgo.fParameters.DevIsMatchTripletsViaMc()) {
- int mc1 = frAlgo.GetMcTrackIdForCaHit(ihit[0]);
- int mc2 = frAlgo.GetMcTrackIdForCaHit(ihit[1]);
- int mc3 = frAlgo.GetMcTrackIdForCaHit(ihit[2]);
+ if (fParameters.DevIsMatchTripletsViaMc()) {
+ int mc1 = ca::Framework::GetMcTrackIdForCaHit(ihit[0]);
+ int mc2 = ca::Framework::GetMcTrackIdForCaHit(ihit[1]);
+ int mc3 = ca::Framework::GetMcTrackIdForCaHit(ihit[2]);
if ((mc1 != mc2) || (mc1 != mc3)) {
continue;
}
@@ -494,7 +514,7 @@ void TripletConstructor::FitTriplets()
ca::MeasurementXy<fvec> mxy[NHits];
for (int ih = 0; ih < NHits; ++ih) {
- const ca::Hit& hit = frAlgo.fInputData.GetHit(ihit[ih]);
+ const ca::Hit& hit = fInputData.GetHit(ihit[ih]);
mxy[ih] = ca::MeasurementXy<fvec>(hit.X(), hit.Y(), hit.dX2(), hit.dY2(), hit.dXY(), fvec::One(), fvec::One());
x[ih] = hit.X();
@@ -518,7 +538,7 @@ void TripletConstructor::FitTriplets()
};
fld.Set(B[0], sta[0].fZ, B[1], sta[1].fZ, B[2], sta[2].fZ);
- fldTarget.Set(frWData.TargB(), frAlgo.GetParameters().GetTargetPositionZ(), B[0], sta[0].fZ, B[1], sta[1].fZ);
+ fldTarget.Set(frWData.TargB(), fParameters.GetTargetPositionZ(), B[0], sta[0].fZ, B[1], sta[1].fZ);
TrackParamV& T = fit.Tr();
@@ -556,11 +576,11 @@ void TripletConstructor::FitTriplets()
T.NdfTime() = sta[ih0].timeInfo ? 0 : -1;
// add the target constraint
- fit.FilterWithTargetAtLine(frAlgo.GetParameters().GetTargetPositionZ(), frWData.TargetMeasurement(), fldTarget);
+ fit.FilterWithTargetAtLine(fParameters.GetTargetPositionZ(), frWData.TargetMeasurement(), fldTarget);
for (int ih = 1; ih < NHits; ++ih) {
fit.Extrapolate(z[ih], fld);
- auto radThick = frAlgo.fParameters.GetMaterialThickness(ista[ih], T.X(), T.Y());
+ auto radThick = fParameters.GetMaterialThickness(ista[ih], T.X(), T.Y());
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, fvec(-1.f));
fit.FilterXY(mxy[ih]);
@@ -594,7 +614,7 @@ void TripletConstructor::FitTriplets()
for (int ih = NHits - 2; ih >= 0; --ih) {
fit.Extrapolate(z[ih], fld);
- auto radThick = frAlgo.fParameters.GetMaterialThickness(ista[ih], T.X(), T.Y());
+ auto radThick = fParameters.GetMaterialThickness(ista[ih], T.X(), T.Y());
fit.MultipleScattering(radThick);
fit.EnergyLossCorrection(radThick, fvec(1.f));
fit.FilterXY(mxy[ih]);
@@ -609,13 +629,14 @@ void TripletConstructor::FitTriplets()
int ih0 = ihit[0];
int ih1 = ihit[1];
int ih2 = ihit[2];
- int mc1 = frAlgo.GetMcTrackIdForCaHit(ih0);
- int mc2 = frAlgo.GetMcTrackIdForCaHit(ih1);
- int mc3 = frAlgo.GetMcTrackIdForCaHit(ih2);
+ int mc1 = ca::Framework::GetMcTrackIdForCaHit(ih0);
+ int mc2 = ca::Framework::GetMcTrackIdForCaHit(ih1);
+ int mc3 = ca::Framework::GetMcTrackIdForCaHit(ih2);
+
if (1 || (mc1 >= 0) && (mc1 == mc2) && (mc1 == mc3)) {
- const ca::Hit& h0 = frAlgo.fInputData.GetHit(ih0);
- const ca::Hit& h1 = frAlgo.fInputData.GetHit(ih1);
- const ca::Hit& h2 = frAlgo.fInputData.GetHit(ih2);
+ const ca::Hit& h0 = fInputData.GetHit(ih0);
+ const ca::Hit& h1 = fInputData.GetHit(ih1);
+ const ca::Hit& h2 = fInputData.GetHit(ih2);
//const CbmL1MCTrack& mctr = CbmL1::Instance()->GetMcTracks()[mc1];
LOG(info) << "== fitted triplet: "
<< " iter " << frWData.CurrentIterationIndex() << " hits: {" << fIstaL << "/" << ih0 << " " << fIstaM
@@ -700,7 +721,7 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
collectedHits.clear();
collectedHits.reserve(maxNhits);
- const ca::Station<fvec>& sta = frAlgo.fParameters.GetStation(iSta);
+ const ca::Station<fvec>& sta = fParameters.GetStation(iSta);
fFit.SetTrack(Tr);
TrackParamV& T = fFit.Tr();
@@ -721,7 +742,7 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
<< (sqrt(Pick_m22 * T.C00()) + grid.GetMaxRangeX() + maxDZ * abs(T.Tx()))[0] << " "
<< (sqrt(Pick_m22 * T.C11()) + grid.GetMaxRangeY() + maxDZ * abs(T.Ty()))[0];
}
- if (frAlgo.fParameters.DevIsIgnoreHitSearchAreas()) {
+ if (fParameters.DevIsIgnoreHitSearchAreas()) {
area.DoLoopOverEntireGrid();
}
@@ -740,7 +761,8 @@ void TripletConstructor::CollectHits(const TrackParamV& Tr, const int iSta, cons
LOG(info) << " check the hit.. ";
}
if (iMC >= 0) { // match via MC
- if (iMC != frAlgo.GetMcTrackIdForWindowHit(ih)) {
+
+ if (iMC != ca::Framework::GetMcTrackIdForWindowHit(ih)) {
if constexpr (fDebugCollectHits) {
LOG(info) << " hit mc does not match";
}
diff --git a/algo/ca/core/tracking/CaTripletConstructor.h b/algo/ca/core/tracking/CaTripletConstructor.h
index de18a4d718..162e0c5279 100644
--- a/algo/ca/core/tracking/CaTripletConstructor.h
+++ b/algo/ca/core/tracking/CaTripletConstructor.h
@@ -34,7 +34,8 @@ namespace cbm::algo::ca
/// Constructor
/// \param nThreads Number of threads for multi-threaded mode
- TripletConstructor(ca::Framework& algo, ca::WindowData& wData);
+ TripletConstructor(const ca::Parameters<fvec>& pars, const ca::InputData& input, WindowData& wData,
+ const float mass, const ca::TrackingMode& mode);
/// Copy constructor
TripletConstructor(const TripletConstructor&) = delete;
@@ -89,8 +90,18 @@ namespace cbm::algo::ca
const ca::Station<fvec>& staR() { return *fStaR; }
private:
+ const Parameters<fvec>& fParameters; ///< Object of Framework parameters class
+ const InputData& fInputData; ///< Tracking input data
+ WindowData& frWData;
+ float fDefaultMass{constants::phys::MuonMass}; ///< mass of the propagated particle [GeV/c2]
+ ca::TrackingMode fTrackingMode;
+
+ bool fIsTargetField{false}; ///< is the magnetic field present at the target
+
+ // Number of stations situated in field region (MVD + STS in CBM)
+ int fNfieldStations{0};
+
bool fIsInitialized{false}; ///< is initialized;
- ca::Framework& frAlgo; ///< reference to ca::Framework object
int fIstaL{-1}; ///< left station index
int fIstaM{-1}; ///< middle station index
@@ -119,7 +130,6 @@ namespace cbm::algo::ca
Vector<ca::Triplet> fTriplets{"TripletConstructor::fTriplets"};
ca::TrackFit fFit;
- WindowData& frWData;
private:
static constexpr bool fDebugDublets = false; // print debug info for dublets
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index a2d72bd140..59e4629e60 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -197,7 +197,7 @@ try {
fUseTRD = false;
fUseTOF = false;
- if (ca::Framework::TrackingMode::kSts == fTrackingMode) {
+ if (ca::TrackingMode::kSts == fTrackingMode) {
fUseMVD = true;
fUseSTS = true;
fUseMUCH = false;
@@ -208,7 +208,7 @@ try {
if (findTask) fUseMVD = findTask->MvdUsage();
}
- if (ca::Framework::TrackingMode::kMcbm == fTrackingMode) {
+ if (ca::TrackingMode::kMcbm == fTrackingMode) {
fUseMVD = false;
fUseSTS = true;
fUseMUCH = true;
@@ -217,7 +217,7 @@ try {
// fInitManager.DevSetIgnoreHitSearchAreas(true); // uncomment for debug
}
- if (ca::Framework::TrackingMode::kGlobal == fTrackingMode) {
+ if (ca::TrackingMode::kGlobal == fTrackingMode) {
//at the moment trd2d tracking only
fUseMVD = false;
fUseSTS = false;
@@ -264,9 +264,9 @@ try {
std::string mainConfig = std::string(gSystem->Getenv("VMCWORKDIR")) + "/macro/L1/configs/";
switch (fTrackingMode) {
- case ca::Framework::TrackingMode::kSts: mainConfig += "ca_params_sts.yaml"; break;
- case ca::Framework::TrackingMode::kMcbm: mainConfig += "ca_params_mcbm.yaml"; break;
- case ca::Framework::TrackingMode::kGlobal: mainConfig += "ca_params_global.yaml"; break;
+ case ca::TrackingMode::kSts: mainConfig += "ca_params_sts.yaml"; break;
+ case ca::TrackingMode::kMcbm: mainConfig += "ca_params_mcbm.yaml"; break;
+ case ca::TrackingMode::kGlobal: mainConfig += "ca_params_global.yaml"; break;
}
fInitManager.SetConfigMain(mainConfig);
}
@@ -372,7 +372,7 @@ try {
// TODO: SZh 15.08.2022: Replace station type with ca::EDetectorID
stationInfo.SetStationType(1); // MVD
stationInfo.SetTimeInfo(!bDisableTime && mvdInterface->IsTimeInfoProvided(iSt));
- stationInfo.SetFieldStatus(fTrackingMode == ca::Framework::TrackingMode::kMcbm ? 0 : 1);
+ stationInfo.SetFieldStatus(fTrackingMode == ca::TrackingMode::kMcbm ? 0 : 1);
stationInfo.SetZref(mvdInterface->GetZref(iSt));
stationInfo.SetZmin(mvdInterface->GetZmin(iSt));
stationInfo.SetZmax(mvdInterface->GetZmax(iSt));
@@ -395,7 +395,7 @@ try {
// TODO: SZh 15.08.2022: Replace station type with ca::EDetectorID
stationInfo.SetStationType(0); // STS
stationInfo.SetTimeInfo(!bDisableTime && stsInterface->IsTimeInfoProvided(iSt));
- stationInfo.SetFieldStatus(ca::Framework::TrackingMode::kMcbm == fTrackingMode ? 0 : 1);
+ stationInfo.SetFieldStatus(ca::TrackingMode::kMcbm == fTrackingMode ? 0 : 1);
stationInfo.SetZref(stsInterface->GetZref(iSt));
stationInfo.SetZmin(stsInterface->GetZmin(iSt));
stationInfo.SetZmax(stsInterface->GetZmax(iSt));
@@ -450,7 +450,7 @@ try {
stationInfo.SetXmax(std::max(fabs(trdInterface->GetXmax(iSt)), fabs(trdInterface->GetXmin(iSt))));
stationInfo.SetYmax(std::max(fabs(trdInterface->GetYmax(iSt)), fabs(trdInterface->GetYmin(iSt))));
- if (ca::Framework::TrackingMode::kGlobal == fTrackingMode) {
+ if (ca::TrackingMode::kGlobal == fTrackingMode) {
stationInfo.SetTimeInfo(false);
}
stationInfo.SetTrackingStatus(true);
@@ -645,7 +645,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
LOG_IF(info, fVerbose > 0) << "\n======= Ca Track finder: process timeslice ...";
}
- fpAlgo->fTrackFinder.FindTracks();
+ fpAlgo->FindTracks();
// IdealTrackFinder();
fTrackingTime = fpAlgo->fCaRecoTime; // TODO: remove (not used)
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 76bde7f973..2d6c6d1d57 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -258,9 +258,9 @@ class CbmL1 : public FairTask {
}
void SetExtrapolateToTheEndOfSTS(bool b) { fExtrapolateToTheEndOfSTS = b; }
- void SetStsOnlyMode() { fTrackingMode = ca::Framework::TrackingMode::kSts; }
- void SetMcbmMode() { fTrackingMode = ca::Framework::TrackingMode::kMcbm; }
- void SetGlobalMode() { fTrackingMode = ca::Framework::TrackingMode::kGlobal; }
+ void SetStsOnlyMode() { fTrackingMode = ca::TrackingMode::kSts; }
+ void SetMcbmMode() { fTrackingMode = ca::TrackingMode::kMcbm; }
+ void SetGlobalMode() { fTrackingMode = ca::TrackingMode::kGlobal; }
ca::TrackingMonitor fMonitor{}; ///< Tracking monitor
@@ -440,7 +440,7 @@ class CbmL1 : public FairTask {
ca::Framework* fpAlgo = nullptr; ///< Pointer to the L1 track finder algorithm
- ca::Framework::TrackingMode fTrackingMode = ca::Framework::TrackingMode::kSts; ///< Tracking mode
+ ca::TrackingMode fTrackingMode = ca::TrackingMode::kSts; ///< Tracking mode
ca::Vector<CbmL1Track> fvRecoTracks = {"CbmL1::fvRecoTracks"}; ///< Reconstructed tracks container
@@ -468,7 +468,7 @@ class CbmL1 : public FairTask {
///
int fSTAPDataMode = 4; ///< Option to work with files for the standalone mode
- TString fSTAPDataDir = "."; ///< Name of input/output directory for running in a STAP mode
+ TString fSTAPDataDir = "."; ///< Name of input/output directory for running in a STAP mode
TString fSTAPDataPrefix = "test"; ///< Name of input/output file prefix. The prefix is defined by output TTree file
TString fSTAPParamFile = ""; ///< Name of the parameter file (generated automatically, if not provided manually)
@@ -487,8 +487,8 @@ class CbmL1 : public FairTask {
static constexpr std::string_view kSTAPAlgoIDataDir = "input_hits";
- Int_t fTrackingLevel = 2; // currently not used
- Double_t fMomentumCutOff = 0.1; // currently not used
+ Int_t fTrackingLevel = 2; // currently not used
+ Double_t fMomentumCutOff = 0.1; // currently not used
bool fUseMVD = false; ///< if Mvd data should be processed
bool fUseSTS = false; ///< if Mvd data should be processed
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 2ee6f2bd06..66be5f0a40 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -1133,7 +1133,7 @@ void CbmL1::TrackFitPerformance()
{"pvi", "Pull Vi [residual/estimated_error]", 100, -10., 10.},
{"distrVi", "Vi distribution [1/c]", 100, 0., 4.}};
- if (ca::Framework::kGlobal == fpAlgo->fTrackingMode) {
+ if (ca::TrackingMode::kGlobal == fpAlgo->GetTrackingMode()) {
Table[4].l = -100.;
Table[4].r = 100.;
}
@@ -1166,7 +1166,7 @@ void CbmL1::TrackFitPerformance()
{"pvi", "Pull Vi [residual/estimated_error]", 100, -10., 10.},
{"distrVi", "Vi distribution [1/c]", 100, -10., 10.}};
- if (ca::Framework::kGlobal == fpAlgo->fTrackingMode) {
+ if (ca::TrackingMode::kGlobal == fpAlgo->GetTrackingMode()) {
TableVertex[4].l = -1.;
TableVertex[4].r = 1.;
}
diff --git a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
index a203bc8bfc..56e657724c 100644
--- a/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
+++ b/reco/L1/L1Algo/utils/L1AlgoDraw.cxx
@@ -79,18 +79,18 @@ void L1AlgoDraw::InitL1Draw(ca::Framework* algo_)
vHits.clear();
vHitsQa.clear();
- auto nHits = algo->fvWData[0].Hits().size();
+ auto nHits = algo->GetWData()[0].Hits().size();
vHits.reserve(nHits);
vHitsQa.reserve(nHits);
- for (const auto& hit : algo->fvWData[0].Hits()) {
+ for (const auto& hit : algo->GetWData()[0].Hits()) {
vHits.push_back(hit);
int iQaHit = hit.Id();
vHitsQa.push_back(CbmL1::Instance()->GetQaHits()[iQaHit]);
}
NStations = algo->GetParameters().GetNstationsActive();
for (int i = 0; i < NStations; i++) {
- HitsStartIndex[i] = algo->fvWData[0].HitStartIndexOnStation(i);
- HitsStopIndex[i] = HitsStartIndex[i] + algo->fvWData[0].NofHitsOnStation(i);
+ HitsStartIndex[i] = algo->GetWData()[0].HitStartIndexOnStation(i);
+ HitsStopIndex[i] = HitsStartIndex[i] + algo->GetWData()[0].NofHitsOnStation(i);
vStations[i] = algo->GetParameters().GetStation(i);
}
}
@@ -244,12 +244,12 @@ void L1AlgoDraw::DrawRecoTracks()
int NRecTracks = 0;
// CbmL1 &L1 = *CbmL1::Instance();
- int curRecoHit = 0;
- cbm::algo::ca::Vector<ca::HitIndex_t>& recoHits = algo->fvWData[0].RecoHitIndices();
- for (vector<Track>::iterator it = algo->fvWData[0].RecoTracks().begin(); it != algo->fvWData[0].RecoTracks().end();
- ++it) {
- Track& T = *it;
- int nHits = T.fNofHits;
+ int curRecoHit = 0;
+ const cbm::algo::ca::Vector<ca::HitIndex_t>& recoHits = algo->GetWData()[0].RecoHitIndices();
+ for (vector<Track>::const_iterator it = algo->GetWData()[0].RecoTracks().begin();
+ it != algo->GetWData()[0].RecoTracks().end(); ++it) {
+ const Track& T = *it;
+ int nHits = T.fNofHits;
// if (nHits > 5) continue; // draw clones
// YZ->cd(); YZ->Update();
// XZ->cd(); XZ->Update();
@@ -485,8 +485,8 @@ void L1AlgoDraw::DrawDoublet(int il, int ir)
void L1AlgoDraw::DrawInfo()
{
cout << " vHits.size = " << algo->GetInputData().GetNhits() << endl;
- cout << " vRecoHits.size = " << algo->fvWData[0].RecoHitIndices().size() << endl;
- cout << " vTracks.size = " << algo->fvWData[0].RecoTracks().size() << endl;
+ cout << " vRecoHits.size = " << algo->GetWData()[0].RecoHitIndices().size() << endl;
+ cout << " vTracks.size = " << algo->GetWData()[0].RecoTracks().size() << endl;
}
void L1AlgoDraw::DrawTarget()
@@ -571,10 +571,10 @@ void L1AlgoDraw::DrawInputHits()
for (int ista = NStations - 1; ista >= 0; ista--) { // //start downstream chambers
ca::Station<ca::fvec>& st = vStations[ista];
- Int_t n_poly = 0;
- Int_t n_poly_fake = 0;
+ Int_t n_poly = 0;
+ Int_t n_poly_fake = 0;
for (int ih = HitsStartIndex[ista]; ih < HitsStopIndex[ista]; ih++) {
- ca::Hit& h = vHits[ih];
+ ca::Hit& h = vHits[ih];
const auto& hQa = vHitsQa[ih];
int iMC = hQa.GetBestMcPointId();
//if( (vSFlag[h.f] | vSFlagB[h.b] )&0x02 ) continue; // if used
@@ -704,11 +704,11 @@ void L1AlgoDraw::DrawRestHits(ca::HitIndex_t* StsRestHitsStartIndex, ca::HitInde
for (int ista = NStations - 1; ista >= 0; ista--) { // //start downstream chambers
ca::Station<ca::fvec>& st = vStations[ista];
- Int_t n_poly = 0;
- Int_t n_poly_fake = 0;
+ Int_t n_poly = 0;
+ Int_t n_poly_fake = 0;
for (ca::HitIndex_t iRestHit = StsRestHitsStartIndex[ista]; iRestHit < StsRestHitsStopIndex[ista]; iRestHit++) {
- int ih = realIHit[iRestHit];
- ca::Hit& h = vHits[ih];
+ int ih = realIHit[iRestHit];
+ ca::Hit& h = vHits[ih];
const auto& qaHit = vHitsQa[ih];
int iMC = qaHit.GetBestMcPointId();
//if( (vSFlag[h.f] | vSFlagB[h.b] )&0x02 ) continue; // if used
--
GitLab